Extrusion process simulation and layer shape prediction during 3D-concrete-printing using the Particle Finite Element Method

To enable purposeful design and implementation of automated concrete technologies, precise assessment and prediction of the complex material flow at various stages of the process chain are necessary. This paper investigates the intermediate stage of the extrusion and deposition processes in extrusion-based 3D-concreteprinting, using a numerical model based on the Particle Finite Element Method (PFEM). In PFEM, due to the Lagrangian description of motion, remeshing algorithms and the alpha shape method are used to track the free surface during large deformation scenarios. The Bingham constitutive model was used for describing the rheological behaviour of fresh concrete. This model is validated by comparing the numerically predicted layer geometries with those obtained from laboratory 3D printing experiments. Extensive parametric studies were then conducted using the numerical simulation, delineating the influence of process and material parameters on the layer geometries, the dynamic surface forces generated under the extrusion nozzle and the inter-layer interactions.


Introduction
Additive manufacturing techniques have found rich applications in various industrial manufacturing processes, in particular of high-end components as required in the biomedical and aeronautical industry [1]. Possibilities of these techniques for manufacturing complex shapes and structures with apparently no limitation in freedom of architectural design have also attracted the interest of the construction industry in 5 the last decade [2,3,4,5,6,7,8]. Developments in this field have the potential of initiating a paradigm change in the practice of construction by enabling novel, highly optimized multifunctional designs of structural components and by improving productivity, efficiency and safety of construction processes. Different approaches have been proposed in the past decade, including contour crafting [2], concrete printing [9], automated slip forming [6], CONPrint3D R [10] and selective activation [11]. For comprehensive classifications 10 and comparisons of additive construction approaches refer to [12,13,14,15]. These approaches vary in scale and comprise applications of both on site manufacturing and prefabrication in plants that operate under controlled environment. Currently, extrusion-based approaches are in majority among additive construction techniques; they are suitable for wide-ranging implementation and for manufacturing of large-scale structures [10]. In extrusion-based 3D-concrete-printing, the material is deposited layer-wise using an automated 15 control of the extrusion nozzle by a robot arm. Extrusion-based 3D-concrete-printing process chains can be divided into production, transportation/pumping, extrusion and deposition stages [16]. Each of these stages requires a specific range of optimal rheological and mechanical properties of the print material, which in time to enable faster deposition rates. The complex material requirements on printable concretes are elaborated in [6,16,17,18,19,20].
In the context of this paper, we look into the intermediate stage at the end of extrusion stage and the beginning of deposition stage. Since the rheological properties required for efficient extrusion and deposition stages vary drastically, the critical properties in this intermediate stage and their implications on both 25 extrusion and deposition are of high significance.
There has been significant progress made in understanding the deposition stage during 3D-concreteprinting. Two mechanisms have been recognized as causes of structure collapse in extrusion-based, layered 3D-concrete-printing during manufacturing: Material failure (plastic failure) and loss of structural stability (buckling) [17,19,21]. In consideration of material failure and static yield stress as the key material 30 performance parameter, Perrot et al. [21] proposed a prediction of critical failure timee after deposition at which the concrete element ruptures if vertical load is further increased due to continued material deposition. Wolfs et al. [19] developed and validated a numerical model for predicting the failure of 3D-printed concrete elements due to loss of structural stability. There are numerous other experimental or numerical works focusing on the printable concrete behavior after the deposition. In contrast to the deposition stage, the 35 research on extrusion stage is scarce. Blockage free, precise extrusion of concrete through varying geometry nozzles at different printing rates is a very complex process. While these challenges are well recognized [22,23], the are dealt with so far primarily using empirical approaches. Considering the fact, that printable concretes have a vast range of time-dependent material properties, and are printed with varying machine configurations, trial-and-error approaches for studying the extrusion process are inadequate. 40 The concrete properties during extrusion process, the machine geometry, machine-material interactions and process parameters (e.g. printing rate, nozzle shape) define the geometry of the printed layers. The printed material undergoes large deformations resulting in shapes, which highly depend on the (timedependent) material behavior, the printing process parameters (e.g. material flow rate or printing speed) and geometrical configuration (e.g. nozzle geometry), with exception of the cases where material character- 45 ized by a very high yield stress is extruded through a horizontal-deposition nozzle, see [10,15,20]. In all cases, any change in the rheological properties during the pumping or the extrusion processes will result in a variation of the shape and size of the extruded layers. Such variations shall be prevented or well-controlled to achieve a consistent layer thickness and desired geometry of the printed structure. It is noteworthy that any changes in layer geometry could also influence the stability of printed structure. The effective geometry 50 of printed layers is not considered in existing models for 3D printed concrete as yet [17,19,21]. The authors believe that it is of outmost importance to delineate the links between concrete properties, nozzle geometry, process parameters and the resulting shape and size of the printed layers. Such an insight into extrusion processes will allow i) to determine the resulting layer geometry and therefore to predict the structural stability of the printed structure, and ii) to dynamically adjust the machine and process parameters to 55 compensate the changes in material properties so that a consistent layer geometry with minimal variations can be ensured.
To allow for a precise assessment of the fabrication process and an optimized machine control, the existing complex material-machine interactions need to be predicted and interpreted correctly. To this end, a computational model is presented for the simulation of extrusion processes in 3D-concrete-printing which 60 allows to analyze the resulting geometries, the time-variant forces and the structural response. Furthermore, extensive numerical analysis was carried out to obtain insight into the forces evolving between the substrate layer and the material being deposited right under the extrusion nozzle during the extrusion. The significance of predicting such forces during the extrusion process is corroborated by preliminary findings in [24,25], leading to the conclusion, that these forces can be few times larger than the actual self-weight of a deposited The presented constitutive law for print material is incorporated into a framework based on the Particle Finite Element Method (PFEM) [28,29]. PFEM utilizes standard finite elements in an updated Lagrangian description of motion to avoid dealing with convective terms in the balance equation. In conjunction with frequent re-meshing, to circumvent element distortions in large deformation scenarios, PFEM is a reliable numerical tool, that showed its robustness in many different applications covering fluid and solid 75 mechanics. In the context of this work, constrained Delaunay triangulations during re-meshing is found to be advantageous for simulating extrusion processes of 3D-concrete-printing. The flow of the freshly printed concrete is modelled with a regularized Bingham model [30]. The advantage of this model lies in the simplicity by approximating the unyielded (solid) regime with a large viscosity. Therefore, the theoretical framework is based on a fluid dynamic formulation. After the introduction of the important theoretical and 80 numerical features of the model, this work focuses on the validation and demonstration of the capabilities of the presented model with respect to different scenarios of extrusion processes of 3D-concrete-printing.
The paper is structured as follows: Section 2 gives a brief overview of the theoretical and numerical framework of the model, covering the constitutive law, the strong form of the governing equations and the basic PFEM algorithm. In Section 3 a one layer 3D-concrete-printing experiment is conducted and used as 85 a validation example for the numerical model. In Section 4 the results from two different numerical studies are discussed to analyze the resulting extruded geometry and the extrusion forces of one printed concrete layer. In Section 5 a third numerical analysis is conducted for studying the dynamic interactions during printing of a second layer on top of an existing layer. Section 6 gives a general discussion on the practical consequences of the findings of this work. Finally, Section 7 provides a summary of the paper and the most 90 important results.

Bingham Model
As was shown in many studies [31,32,33,34,35], the flow of homogenized fresh concrete can be suitably modelled with viscoplastic fluid models on the macro scale. The material response of these models 95 is controlled by a certain yield stress that must be exceeded to initiate a viscous flow of the material. Typically, the regime before yielding is neglected and approximated by a rigid behavior. The Bingham model is a viscoplastic fluid, where the viscous flow is characterized by a constant plastic viscosity. The constitutive relation of a Bingham model is defined in one dimensional shear as τ xy = τ 0 + µγ when |τ xy | ≥ τ 0 , 100 with the deviatoric stress component τ xy , the yield stress τ 0 , the plastic viscosity µ and the shear strain ratė γ. Due to a rigid material behavior in the unyielded regime, difficulties of an unbounded apparent viscosity for strain rates being zero arise in numerical implementations. This issue is generally addressed by using regularized models, where the rigid material response is approximated by a high viscosity. In this paper the regularized Bingham model of [30] is applied, which, formulated in 3D as where τ is the deviatoric Cauchy stress tensor, |γ| = √ 2D : D is the equivalent strain rate, D is the deformation rate tensor and m is a regularization parameter. The deformation rate tensor is defined by the symmetric part of the velocity gradient tensor as 110 The regularization parameter m must be calibrated cautiously. It should be chosen large enough, so that the rigid response is captured accurately, but small enough, so that the problem remains numerically solvable.
Typically, a value between 1000 and 3000 has been found to be a reasonable approximation. Furthermore, the equivalent stress is defined as the square root of the second invariant of the deviatoric stresses

115
The equivalent stress is used for the 3D extension of the yield criterion in (1) as a measure to distinguish between yielded (τ eq ≥ τ 0 ) and unyielded (τ eq < τ 0 ) regimes of the material in the studied examples.

Governing equations
Using the regularized Bingham model (2) as the constitutive equation, the flow of fresh concrete is described by the balance of momentum and mass, with all essential boundary conditions. The local form of 120 the balance of momentum for a continuum description of the material flow in a time interval [0, T ], covering the domain Ω, is given as with the density ρ and the body forces ρb i acting within the domain. Furthermore, displacements v * i are prescribed on the Dirichlet boundary Γ v as and tractions t * i are acting on the Neumann boundary Γ t as with n j as the normal vector. The balance of mass is given, assuming a slight compressibility of the material, via [36] 130ṗ with the volumetric part of the deformation rate tensor D v = D ii , the pressure p and the compressibility K. The Cauchy stress tensor is split into a deviatoric and a volumetric part as with the second order identity tensor I ij = δ ij . 135

The Particle Finite Element Method
The Particle Finite Element Method (PFEM) is a numerical method for solving fluid dynamical problems in an updated Lagrangian framework [28, 29,37]. PFEM has proven its robustness in many different applications in fluid dynamics, for example fluid structure interaction problems [36,38], Bingham fluid problems [39,40] or multi-phase flow [41]. The reliability of PFEM was also demonstrated in applications which is treated be deletion or insertion of particles. Solution variables are then linearly mapped to newly inserted particles.
In the context of this work it was found, that modifying the original re-meshing procedure with a constrained Delaunay triangulation is advantageous with respect to free surface recognition and mass con-155 servation. The advantage of a constrained Delaunay triangulation lies in the fact, that the free surface of the domain is constrained throughout the re-meshing procedure. The free surface must be tracked throughout the simulation process and the discretization must be updated, when element edges are becoming too long or too short by insertion or deletion of particles. This treatment is straightforward in 2D, when self contact is ignored. However, in 3D, to guarantee a well-shaped surface triangulation can be difficult, due 160 to higher geometrical complexity. Therefore, when necessary, free surface recognition is supported by an additional alpha shape triangulation in 3D. Contact between bodies is modelled by activation of contact elements between bodies, that are formed from a second constrained triangulation, only consisting of the free surface and boundary particles. Contact elements are activated, when the distance between two bodies becomes smaller than a given threshold. The material parameters of these contact elements are the same as 165 for the main body, see [28,37] for more details about contact treatment in PFEM. Contact with the ground is modelled by a standard penalty method. The basic steps of the PFEM algorithms are summarized as 1. Define a cloud of pseudo-particles at the current time step t n and define the free surface. 2. Generate the finite element discretization by the constrained Delaunay triangulation of the domain.
3. Discretize and solve the system of governing equations in an updated Lagrangian framework at t n+1 = and discretized by linear finite elements. Consequently, the Ladyzhenskaya-Babuska-Brezzi (LBB) condition is violated with equal order finite elements for the velocity and pressure field. To this end, a direct pressure stabilization technique is adopted to stabilize the mixed formulation. The stabilization technique is based on a polynomial projection of the C 0 approximated pressure field onto a pressure field, which is described by a one order lower approximation [48]. Finally, the system of equations is solved monolithically.

One layer 3D-concrete-printing experiment and numerical validation
In this section the results of the numerical simulation are compared with those obtained in experiments for one 3D-printed concrete layer. A vertical, circular extrusion nozzle is chosen, since it produces more complex shapes of the printed cross-section in comparison to a rectangular extrusion nozzle, thus, making the numerical validation more meaningful. The experimental background covering all relevant information 185 of the concrete mixture and the conducted printing experiment are presented first. Then, the numerical study is explained and the results from the simulations and the laboratory tests are compared.

Laboratory test
A fine grained printable concrete, very similar to Mixture A in [22] is extruded with a gantry printer. Portland cement CEM I 52.5 R according to EN 197-1 [49] was chosen as the main binder. An aqueous highrange, water-reducing admixture (HRWRA) on a polycarboxylate-basis (PCE) was used as superplasticiser. This PCE was developed for ready-mix concrete applications and thus ensures an extended consistency retention, i.e., over 90 min after water addition. Tap water at sub-room temperature was used for all mixtures under investigation. Supplementary cementitious materials (SCMs) were added to reduce the total cement content and, consequently, the CO 2 footprint. At the same time, the SCMs employed were expected 195 to act as rheology modifiers, which increase the structural build-up rate. These SCMs were i) Class F Fly ash (denoted as FA) according to ASTM C618-12a [50] and ii) microsilica suspension (denoted as MSS) with a specific surface of 15-30 m 2 /g and a mean particle size of 0.15 µm; 80 wt% of primary particles had diameters less than 5 µm. The MSS contained 50 wt% microsilica and 50 wt% of water. Furthermore, very fine quartz sand (fraction 1: 0.06 to 0.2 mm) and two fine natural river sands (fraction 2: 0 to 1 mm and 200 fraction 3: 0 to 2 mm) were used as aggregates. The particle size distribution of the sands may be found in [16]. The effective water-to-binder ratio of the mixture was 0.42. The mixture had a binder volume of 52.5 %.
A gantry printer with a custom made printhead was used. Its printhead can extrude concrete with maximum aggregate size of up to 10 mm, see Fig. 2. Further details can be found in [10]. For the validation 205 of numerical model, three seperate layers were produced with a circular nozzle with a diameter of 56 mm. The length of each layer was 50 cm. The deposition drop height was kept constant at 24 mm. The printhead velocity was 50 mm/s, the time interval between the layers was below 2 min.
The rheological properties of the printable concrete were measured using a Viskomat XL rotational rheometer with a vane rotor. Dynamic and static yield stresses were measured at a concrete age of 20 min, 210 which coincides with the concrete age at which the extrusion occurred. Torque and rotational velocities were transformed into shear stress and shear rate values using the Reiner-Riwlin equations. The dynamic yield stress was calculated from the downward curve of a hysteresis test in which the shear rate was raised to 13.04 s −1 (at 80 rpm), with constant shearing applied for 25 s, before dropping back to 0 s −1 . The static yield stress was determined from constant shear rate tests with an applied shear rate of 0.05 s −1 , making 215 sure that the effective shear rates are equal to the applied shear rates and flow onset takes place in every test as described in [51]. Cross-section measurements were taken in both fresh and hardened state of the concrete, by cutting the printed layers laterally. Layer geometries were measured carefully with a Vernier calipers, at the cut cross sections immediately after the extrusion, as well as, after one-day age.

Computational simulations of a concrete printing experiment 220
In this subsection the results from numerical simulations of the printing process are validated by comparing the computed cross-section and the real cross-section of one layer of the printed material. The flow rate of the material is deduced from the cross-sectional area and the printing speed in the experiment. In the model, the extrusion flow is approximated by prescribing a constant inlet velocity directly at the outlet of the extrusion nozzle in the numerical simulation. The concrete processing inside the 3D-printer is not considered. The inlet flow is modelled by continuously pushing new elements into the domain. Material and process parameters are given in Table 1. Two numerical studies with different yield stress values are performed. The lower value is the dynamic yield stress of the material, measured from the down-curve of a hysteresis loop rheology test. During these measurements the concrete is strongly sheared and it is questionable, if the printed concrete is similarly 230 sheared during the extrusion process, considering that the printed concrete behaves very thixotropic by developing a static yield stress exceeding 1 kPa quickly after extrusion. Therefore, the numerical experiment is also conducted with the larger static yield stress of 1.2 kPa to evaluate the sensitivity of the printing process with respect to the chosen yield stress. Assuming symmetry of the problem along the printing direction, only one half of the cross-section is simulated to reduce the computational time. No-slip boundary conditions are 235 applied along the interface between the printed concrete and the surface.
With a printing speed of v print = 0.05 m/s and a layer length of 0.5 m the printing process takes t = 10 s. The numerically "printed" concrete with the yield stress τ 0 = 0.306 kPa is shown in Fig. 3 at three stages (t = 3 s, t = 6 s and t = 9 s). As can be observed, after a few seconds of printing time, a steady state printing process is reached and a constant cross-section of the layer is obtained. In Fig. 4, the velocity norm �t���s �t���s �t���s Figure 3: Numerical simulation of one layer extrusion: Shape of the printed concrete layer for τ 0 = 0.306 kPa after t = 3 s, t = 6 s and t = 9 s. 240 around the extrusion nozzle after t = 10 s is depicted. Even though a constant inlet flow velocity is applied in this study, the material quickly develops a certain steady state velocity profile after extrusion.
In Fig. 5, the geometry of the numerically "printed" concrete layer, obtained from simulations using the dynamic yield stress τ 0 of 0.306 kPa and the static yield stress τ 0 of 1.2 kPa, is compared to that measured during the laboratory tests. Considering the simplifications made in this study with respect to the inlet flow the height of the layer. The computed and measured areas of the cross-sections differ less than 1 %. This 250 difference is attributed to a slight compressibility considered in the numerical analysis, see (8). It should be noted that the inlet flow velocity or extrusion rate in the simulation is deduced from the experiments. Thus, it is to be expected, that the amount of printed material is the same in the simulations and the experiments. Generally, the geometry of the numerically "printed" layer using the dynamic yield stress τ 0 of 0.306 kPa matches the experimental results better as compared to the results obtained with the static yield stress 255 τ 0 of 1.200 kPa. However, the analysis using the larger yield stress leads to a more rounded shape at the lower lateral parts of the cross-section, which agrees better with visual assessment in the experiment. As stated earlier, the yield stress of the concrete at an intermediate stage during the printing process could be somewhere in between the rheometrically measured dynamic and static yield stresses. The shear history of such concrete is neither as high as that during the downwards curve of a hysteresis test, nor it is as low as 260 during constant shear rate tests with an applied shear rate of 0.05 s −1 . Hence, simulations using a Bingham model with a level of the yield stress between the static and dynamic yield stress is expected to yield an even better correlation of the predicted and the measured geometries. It shall also be noted that the regularized Bingham model does not model neither shear thinning behaviour nor thixotropy of concrete. Consequently, shear thinning fluid models such as the Herschel-Bulkley model may be used to further improve the obtained 265 results. However, this study was not considered in the study at hand, since additional parameter calibration is needed for such models, for which no directly measured experimental data is available, and since the Bingham model yielded satisfactory results.

Numerical study on parameters affecting the printing process of one layer of concrete
The following numerical investigations intend to provide insight into the evolving loadings and topological 270 changes during the concrete printing process. They are divided into three analyses. In the first study (Subsection 4.1), the influence of different process parameters on the printed cross-sections is analyzed. This analysis is designed as a parametric study of a one-layer printing test with a round printing nozzle. The round printing nozzle is again in focus here, as the resulting shapes of the cross-section are more complex and less predictable than, for example, in case of a rectangular cross-section. In the second study 275 (Subsection 4.2), the extrusion forces generated below the printing nozzle are analyzed. These forces are several times higher than that induced by the dead weight of a layer and may have an influence on the deformations of lower layers or interlayer bond strength. In this analysis, the influences of the material and process parameters are examined. As this is an extensive study with an associated vast calculation effort, the study is performed in 2D, adopting a rectangular printing nozzle. In Section 5 the dynamic interactions 280 when printing a second layer on an existing layer will be analyzed. An overview of the numerical studies conducted in this section and in Section 5 is provided in Tab. 2. Finding the optimal process parameters in a 3D-concrete-printing process which result in the desired width and height of the printed concrete is a complex matter due to interactions of many different factors 285 during manufacturing. For example, by varying the printing speed, the flow behavior of the material may change, which vice versa results in a different cross-sectional shape. To this end, numerical models can help to understand these interactions and support the design of the process. Therefore, in the following computational analysis the influence of different parameters controlling a 3D-concrete-printing process is analyzed by evaluating the resulting cross-sectional shapes. The layout of the printing process adopted in 290 the numerical analysis is illustrated in Fig. 6. This analysis constitutes a parameter study of the experiment in the previous section 3. One layer of 50 cm length is printed with a circular nozzle. The material and process parameters are summarized in Tab

Parameter
Value The shapes of concrete cross-sections for varying printing heights h are shown in Fig. 7a. With increasing printing height also the height of the printed cross-section increases, which is plausible. The width of the cross-section decreases with increasing printing heights, due to a constant flow rate of the material. In 300 this example, the change of the cross-sectional width is approximately two times the change of the height. Hence, if the height of the printed layer needs to be adapted without changing the width of the layer, the flow rate needs to be scaled with the printing height h by the amount of material change in width direction. Fig. 7b shows the results of the numerical simulations for different radii r of the extrusion nozzle. In this example, the extrusion flow is kept constant by scaling the inlet velocity with the extrusion nozzle radius 305 r. Hence, the same amount of material is printed for all extrusion nozzle radii r. As can be observed, by using a larger extrusion nozzle radius, the upper layer surface is flattened and the overall shape becomes more rectangular-like. When a smaller extrusion nozzle radius is used, the shape becomes less rectangular, resulting in a slightly larger layer width and smaller layer height.   Note that the extrusion flow is constant and due to varying printing speeds, different cross-sections of the printed material are obtained. Changing the printing speed affects both, the width and the height. Hence, finding a general scaling rule between printing speed and printed shapes is difficult. For this reason and also for the reason, that printing speeds may change the timing of the whole printing process (e.g. time intervals between layers), tuning the printing 325 speed for adapting the layer shapes is not recommended.
In addition, another study with varying printing speed v print is carried out, where the inlet velocity v inlet is linearly scaled with the printing speed v print to yield the same amount of printed material, see Fig.  9b. As can be observed, larger printing speeds lead to slightly smaller layer heights and widths. However, these differences are small and in general it can be stated, that when the flow rate is linearly scaled with the 330 printing speed, identical cross-sections are produced. A reason for the observed small differences could be viscosity and inertia effects at different flow velocities. This finding supports the fact, that the yield stress is actually the relevant material parameter for controlling the shape of printed layers and that viscosity and inertia effects are almost negligible. The findings from the presented numerical investigation can be summarized up as follows: For fine tuning 335 the width of printed layers, varying the inlet velocity (or flow rate) is the most predictable strategy. If the layer height needs to be adapted, changing the printing height is the only measure. While doing so, the flow rate must be scaled, to keep the layer width constant. Furthermore, by changing the radius of the extrusion nozzle, more drastic changes of the cross-section can be made in terms of layer width and height.
Increasing the extrusion nozzle radius may also flatten the upper layer surface, which is usually necessary, 340 when stacking multiple layers during manufacturing. Changing the printing speed is not recommended to control the shape of the printed layer. However, when the inlet velocity (or flow rate) is scaled with the printing speed, almost identical cross-sections can be produced for different printing speeds.

Forces under the extrusion nozzle
The following numerical study focuses on the analysis of dynamic surface forces induced below the nozzle 345 during the extrusion process of the printed material. As will be seen, these forces are several times higher than the force induced by self-weight of one layer and may cause deformations of substrate layers and/ or influence the interlayer bond strength. Hence, it is important to understand the general characteristics and dependencies of such forces.
The numerical experiments presented in this subsection are characterized by a systematic investigation of 350 material and process parameters to gain deeper insights into the complex interactions during the extrusion process of 3D-concrete-printing. Because of the large computational effort of 3D simulations and since interpretation of 2D results is easier, this study is conducted in 2D. It is based on a one layer test with a layer length of 40 cm, which is printed with different material parameters, geometries and process parameters, see Fig. 10. The material and process parameters and the intervals of certain parameters are summarized 355 in Tab. 4. The process parameters refer to a printing process with a rectangular cross-section, since only this nozzle shape can be realized in 2D. The printing speed v print is set to be equal to the inlet velocity v inlet and the distance between the extrusion nozzle and the ground surface is set equal to the width of the extrusion nozzle plus a certain separation distance d sep to tune the printing process, see The numerically "printed" layer and the distribution of the pressure at three different stages obtained for the reference parameters are shown in Fig. 11. The pressure distribution in the vicinity of the extrusion nozzle is zoomed in. It is observed, that the pressure is substantially increased below the extrusion nozzle. In real printing processes this pressure must be carried by the lower layers. Consequently, it must be guaranteed, that the lower layers are already strong enough to carry such loads, so that an accurate printing    At first, the influence of the separation distance d sep is investigated. In Fig. 12a the computed shape of the printed material around the extrusion nozzle is shown for different separation distances at a position 370 of 30 cm from the beginning of printing, adopting the reference parameters in the numerical model. As can be observed, a "belly" is formed in front of the extrusion nozzle, which becomes more prominent with decreasing separation distance. For small separation distances this "belly" accumulates over the printing time and more material is deposited in the front region. However, if the separation distance becomes large enough, a more or less vertical front is formed below the nozzle. It is expected, that this effect becomes less 375 important in 3D analysis, when the material can simply flow laterally to the sides. A separation distance of d sep = w/20 was found to produce satisfactory results with respect to the formation of this "belly". Hence, this value is used as a fixed value in the remaining studies. The resulting vertical surface tractions at the ground below the extrusion nozzle are depicted in Fig. 12b for different separation distances. For small separation distances, a maximum of vertical surface traction is observed near the rear corner of the nozzle, the level of which decreases with increasing separation. For larger separation distances, the loading distribution changes and the location of the maximum traction moves towards the front corner of the nozzle. The peak values of the surface forces are roughly two times higher in comparison to those induced by the self-weight of one layer. Consequently, in a real printing process lower layer would need to carry twice the amount of the self-weight temporarily and locally under the extrusion nozzle. The next study to be presented focuses on the influence of the structuration of the freshly printed concrete after being extruded. In Fig. 13a, the vertical surface tractions under the nozzle are plotted along the length of the printed concrete layer for different yield stresses. The surface tractions increase with an increasing dynamic yield stress, as more energy and consequently a higher injection pressure is needed to print these materials. A similar trend can be observed for an increasing plastic viscosity, see Fig. 13b, which 390 shows, that the maximum tractions also increases with increasing plastic viscosity. While changing plastic viscosities qualitatively does not change significantly the force distribution, the yield stress has a strong influence on the force distribution. For larger yield stresses, the contact region decreases and the pressure tends to localize near the front end of the nozzle.  In the last part of this analysis, the influence of selected process parameters (printing velocity, width of 395 the nozzle) are investigated. Fig. 14a clearly indicates increasing tractions with increasing printing speed, due to larger inlet velocities v inlet , which are scaled by the printing speed v print . Lower printing speeds are therefore more preferable, when the forces acting on substrate layers shall be reduced. Fig. 14b compares the influence of different extrusion nozzle widths on the vertical surface tractions below the extrusion nozzle. Note that the height is also scaled with the nozzle width via h = w + d sep /20. Printing speed and inlet 400 velocity remain constant. Hence, more material is printed with an increasing nozzle width. A first guess would be, that the tractions simply scale with the nozzle width. However, Fig. 14b elucidates, that the peak tractions below the extrusion nozzle are almost the same, independent from the extrusion nozzle width (and layer heights). This is attributed to the fact, that a more localized distribution of the stresses in the material under the extrusion nozzle is observed for smaller nozzle widths, see Fig. 15, as less material for 405 stress transfer is available. In a final numerical experiment, the analysis is repeated for different nozzle widths, however, now the flow rate is kept constant and the inlet velocity v inlet is scaled by the flow rate of the reference state parameter. For this example the distance between the ground and the nozzle is h = 2.1 cm. Fig. 16 shows the vertical surface tractions for different extrusion nozzle widths. As can be observed, the peak level of the tractions 410 is almost identical for widths between 1 and 2 cm. It is interesting to note, that for a smaller width the forces are distributed over a larger area, while for a larger width the traction forces acting on the supporting ground are more concentrated. For larger widths (w = 3 cm), the forces below the nozzle decrease and the "S-like" shape of the tractions turns into a "L-like" shape, indicating that another flow regime is reached. For this case, the inlet velocity v inlet is sufficiently small, and a certain "pulling effect" of the material starts These "S-like" and "L-like" flow regimes can also be observed in Figs. 13a and 13b for changing material parameters, indicating that the flow regime depends on both material and process parameters.

Numerical studies on parameters affecting the interactions between two layers of printed concrete 420
In this section the dynamic interactions during printing of a second layer on top of an existing layer are analyzed. The study focuses on the question, whether and when the stresses in a substrate layer exceed the yield stress. When this is the case, lower layers may get affected due to plastic deformations. The analysis is first carried out in 2D and the obtained results are then compared to results from 3D analyses.

425
The numerical investigation of a printing process on top of an existing layer of printed concrete is performed as follows: The initial state is characterized by an already printed layer of a certain length (taken as 40 times the width of the nozzle). This layer has a substrate yield stress τ 0,s . Subsequently, a second layer is printed on top of the existing layer. This induces forces during the printing process in the substrate layer. The yield stress of the new printed layer will be denoted as the overlay yield stress τ 0,o . Fig. 17 430 illustrates the layout of the numerical experiment and the analyzed printing scenario. The range of material and process parameters are the same as in the previous section, with the exception, that the extrusion nozzle width w only varies between 1 and 2 cm, see Tab. 4. The reference material parameters also remain the same. It was found, that a substrate yield stress lower than 1 kPa did not provide enough capacity for preventing plastic collapse over the entire studied parameter range. Hence, the substrate yield stress is 435 varied between 1 and 2 kPa. In a first analysis, the case when the substrate and overlay yield stresses are both equal to 1 kPa is investigated for two widths of the extrusion nozzle (1 cm and 2 cm). Fig. 18a shows the equivalent stresses induced in the two concrete layers. A similar trend as in the previous analysis can be observed: Larger stress concentrations appear for smaller extrusion nozzle widths. When the equivalent stress exceeds 440 the yield stress, plastic deformations occur and the material yields. Hence, the influence of forces generated during the printing process and acting on the substrate layer can be examined by checking the yield criterion to distinguish between yielded (red) and unyielded (blue) material, see Fig. 18b. As can be seen, lower layers may not be strong enough to fully bear the dynamic loading induced by the extrusion process, due to yielding of the material of the substrate layer under the extrusion nozzle. This effect is more pronounced 445 for the smaller extrusion nozzle width, which is in agreement with the findings from the previous analysis; see Section 4.2. Hence, one may conclude that printing layers with larger extrusion nozzles is preferable. Next, the influence of the yield stress in the substrate layer is evaluated. When the substrate yield stress increases, i.e. when the time interval between layers or the microstructure build-up increases, the substrate layer becomes sufficiently strong to fully carry the loads from the extrusion process. This is visualized in Fig.   450 19a, which shows yielded and unyielded material for different substrate yield stresses in the vicinity of the extrusion nozzle. The overlay yield stress is adopted as τ 0,o = 1 kPa. For τ 0,s > 1.4 kPa the substrate layer does not yield anymore and is strong enough to carry all loads from the extrusion nozzle. Evaluating the total volume of the yielded material is a useful indicator for the influence of forces acting on the substrate layer during the printing process. The volume of the yielded material is evaluated for both layers for varying 455 overlay yield stresses τ 0,o and plotted as a function of the substrate yield stress τ 0,s in Fig. 19b. For a sufficiently large substrate yield stress, the volume of the yielded material approaches a constant value. At this point only the currently printed layer yields, while the substrate layer is strong enough to carry all loads during printing. When the substrate reaches a yield stress of 2 kPa, the lower layer is strong enough for the complete range of overlay yield stresses investigated. For overlay yield stresses lower than 0.650 kPa no 460 influence on lower layers is observed at all.
A similar investigation was performed for varying plastic viscosity. In Fig. 20a the total yielded volume is plotted as a function of the substrate yield stress for different plastic viscosities. For the range of values considered here, the influence the viscosity is smaller than that of the yield stress; cf. Fig. 19b. This observation is in agreement with the results presented in Section 4.2, where it was found, that with increasing 465 plastic viscosity the tractions are more uniformly distributed below the extrusion nozzle in comparison to the tractions for increasing yield stresses. Hence, for an increasing overlay yield stress the stresses in the material are more localized towards the front edge under the nozzle, which leads to local stress states exceeding the yield stress, see Fig. 18b.
Next, the effect of an increased printing speed v print = 0.2 m/s is examined. The simulation results are 470 depicted in Fig. 20b, which is to be compared with Fig. 19b, based on v print = 0.1 m/s. As expected, the total yielded volume becomes larger for a large printing speed. Already with an overlay yield stress of 0.475 kPa, the stresses may exceed the yield stress in the substrate layer. Finally, the yielded volume is analyzed for a smaller extrusion nozzle width of w = 1 cm; see Fig. 21a. Due to less printed material, the yielded volume is smaller in comparison to the results obtained from the 475 larger extrusion nozzle width w = 2 cm in Fig. 19b. For a more meaningful comparison, the obtained yielded volumes for both extrusion nozzle widths must be normalized individually by an equivalent volume. This volume is herein chosen as the yielded volume obtained with a substrate yield stress of τ 0,s = 2 kPa. The normalized yielded volume for the case of an overlay yield stress τ 0,o = 1 kPa is depicted in Fig. 21b. It can be observed, that the yielded volume is more sensitive to the process with the smaller extrusion nozzle width, 480 as the yielded volume of the smaller extrusion nozzle width increases more rapidly for higher substrate yield stresses. The reason for this are the more concentrated stress states for smaller printed layers, which was already found in Section 4.2 and visually assessed in Fig. 18a.

3D simulations of a 2-layer printing process
The 3D simulations in this subsection serve as verification of the 2D studies and gives further insights � eq ��� � � � eq ��� � w=1 cm w=2 cm Nevertheless, the yielded volume is smaller in comparison to the results obtained in the 2D analysis. It is expected, that increasing the width of the extrusion nozzle in direction of the layer width, the 3D 500 state will even more resemble the 2D state. Even though, the material is not constrained laterally in the 3D simulations, horizontal deformations of the substrate layer are around 0.5 %, which is negligible from the practical perspective. It must be noted that these deformations are caused by plastic flow only. More advanced elasto-viscoplastic models are necessary to capture elastic deformations in the substrate layer during printing. According to [15,26,27], the yielding characteristics of the printed material has consequences for the interlayer bond strength of the material. Layer bonds may get improved, when the substrate layer yields.

510
For the material parameters, however, the analyses presented in this work elucidate a paradox situation: The yielded volume of the substrate layer is the larger the lower its yield stress is. At the same time, the yield stress of the currently printed layer has to be higher in order to increase the yielded volume. However, concrete with a lower yield stress would fill the contours on the substrate layer better and lead to an enhanced bond. Hence, a generalized conclusion on the optimal tuning of the material properties is not yet possible 515 based on the numerical results only.
With respect to the process design, the situation is similarly complex. If only the results presented in this paper are taken into account, increasing the printing speed would increase the yielded volume, and therefore could result in a better bond. However, in [52] it was shown in a micromechanical investigation, that lower printing speeds are preferable for higher surface roughness, which in its turn would increase the interlayer 520 bond strength. Hence, with regard to process parameters, different influences are counteracting, which again makes a generalized statement difficult. Furthermore, increasing the yielded volume for enhanced interlayer bonds is in conflict to the general purpose of reducing the deformations in the substrate layer.
As found in Section 5, plastic deformations can occur in lower layers during printing. In the 3D examples investigated in this paper, such deformations are moderate and negligibly small for real applications.

525
However, it is known, that 3D printed concrete shows elastic behavior at very early age after printing [17,19]. The Young's modulus increases with time due to concrete's thixotropic nature [33]. Hence, during the printing process, also elastic deformations may arise, which are not captured in this work. This is due to the limitations of the adopted viscoplastic constitutive law, where the material response before yielding is approximated by a high viscosity value. In order to capture the elastic effects before yielding, an 530 elasto-viscoplastic model is necessary.
In the analyses it was observed, that the yielded volume increases nonlinearly with decreasing yield stresses; see Fig. 21a. For lower values of the yield stress the yielded volume and, consequently, plastic failure, become highly sensitive to the given yield stress, and even small changes in the values might lead to plastic failure. In the numerical examples the substrate yield stress was chosen high enough, so that plastic 535 collapse during printing was prevented.
The dynamic surface forces generated during extrusion of the material may influence the stability of the printed structure. While stability and plastic failure of 3D-concrete-printed structures due to self-weight has been addressed in the literature (e.g. [17,19,21,53,54]) the dynamic forces induced by extrusion were not considered in these approaches. Since such forces act as an additional load at the top of the currently 540 printed wall segment, this may result in a smaller safety margin w.r.t. collapse and a smaller critical height of the structure, see Fig. 24.

Conclusions
In this paper, results from computational simulations of the extrusion process during 3D-concrete-printing have been presented. The Particle Finite Element Model (PFEM), based on a Lagrangian formulation and 545 continuous re-meshing strategy, was employed in conjunction with a fluid dynamic approach with a regularized Bingham constitutive law to simulate the flow of fresh concrete. The numerical investigations provided insights into the shape of the printed cross-section, the extrusion forces, yielding characteristics of substrate layers and layer interactions during the intermediate stage at the end of extrusion and the beginning of deposition. Due to the very different requirements for efficient extrusion and deposition, extensive numerical 550 experiments were performed to understand the key tendencies with respect to a purposeful choice of process and material parameters.
For the validation of the numerical model a one-layer 3D-concrete-printing experiment with a circular, vertical extrusion nozzle was simulated. Despite the simplifications connected with use of a viscoplastic material formulation with a constant viscosity, the geometry of the numerically "printed" cross-section 555 showed a very good agreement with the measurements' results obtained from the laboratory tests.
To investigate the influence of the process and material parameters on various features of the printing process, computational analyses of the printing of one and of two layers of concrete were performed.
In a first numerical analysis, the influence of various process parameters on the cross-section shape of a one layer 3D-concrete-printing test was investigated. The study revealed, that printing height and flow rate 560 are most effective process parameters to tune the height and width of the resulting concrete layer. Though it is also possible to achieve wider layers with smaller nozzle sizes, such setup increases the unevenness of printed layers. Larger nozzle diameters generally lead to printed layers with a more regular shape of cross-sections, when the space under the nozzle is completely filled and the flow rate is constant. It was also found, that identical cross-sections can be obtained for varying printing speeds, if the flow rate is scaled accordingly. Hence, inertia and viscosity effects play a minor role in affecting the layer shape, leaving yield stress as the major material parameter influencing the shape of printed cross-sections. The presented computational model is helpful to predict the resulting shapes of the deposited layers, without the need for an a priori time consuming trial-and-error testing. This might be relevant in particular, when the printed cross-section is not supposed to be constant, but designed to change during the printing 570 process to enable a flexible control of the layer geometry. Conversely, the developed numerical model provides guidance how to maintain a constant layer cross-section despite variations in material properties, by adapting the process parameters accordingly.
In a second numerical study, the influences of process and material parameter on the dynamic surface forces generated below the extrusion nozzle during 3D-concrete-printing were investigated. It was found 575 that such forces are several times larger than those induced by the self-weight of one layer. Hence, lower layers would need to bear this load temporarily under the extrusion nozzle. Increasing the yield stress or the viscosity will increase the extrusion forces, due to higher forces that are necessary to print these materials. Smaller printing heights and higher printing speeds lead to larger forces. This study also showed, that stress concentrations under the extrusion nozzle are higher for smaller printed layers.

580
Interestingly, varying extrusion nozzle widths with a constant flow rate result in a similar peak level of extrusion forces. In this study two flow regimes were deduced from considering the forces. An "S-like" regime is obtained when more material is deposited in the front region of the nozzle, leading to a "belly"-shaped printed material. In the "S-like" regime the peak level of extrusion forces is located near the rear corner of the extrusion nozzle. The shape of the printed material approaches an "L-like" regime when less material 585 is deposited in the front region of the nozzle and a more or less vertical front is formed below the nozzle. In the "L-like", regime the peak of the extrusion forces is located under the front corner of the extrusion nozzle.
In the third numerical investigation, the dynamic interactions between the printed layer and the substrate layer were analyzed. This study revealed, that for certain combinations of process and material parameters 590 the lower layer might start yielding under the nozzle due to the additional forces from the extrusion process. As an indicator for measuring the influence of the extrusion process on the substrate layer, the volume of the yielded material was calculated. Relating the substrate yield stress to the volume yielded during the printing process, the minimum age-related yield stress necessary to carry the loads from the extrusion process was determined. A substrate yield stress of roughly 2 kPa was found to be sufficient in all printing scenarios 595 under investigation. Similar relationships were found between material and process parameters on one hand and the resulting yielded volume on the other hand as in the analyses of a single layer printing process. Only small plastic deformations of the printed substrate layer were observed. However, depending on the specific combination of process and material parameters, a substrate's yield stress lower than approximately 1 kPa would lead to a collapse during the printing process.

600
The computational experiments presented in this paper provide a deeper insight into the printed geometries, dynamic interactions between the nozzle and deposited material as well as into the material and structural response during 3D-concrete-printing. These findings may help in the design of printing materials and processes and possibly, when transferred into a real-time steering software, may support adapting process parameters during printing.

605
Reflecting the findings from the computational study in view of the experiments with real 3D-printing, it was noted, that microstructural features of 3D concrete printing, such as surface roughness, also affect the printed structure and its properties, and specifically the bond between the layers. Due to limitations of the adopted rigid-viscoplastic model, elastic deformations in the unyielded regime of the printed concrete are not taken into account yet. When considered, such deformations and related stresses during printing 610 may be captured. Finally, the dynamic forces induced during the extrusion process under the nozzle act as a load at the top of the currently printed wall segment, which may influence the early age stability of the printed structure.