Penetration of Cement Pastes into Particle-Beds: A Comparison of Penetration Models

For the selective paste intrusion (SPI) method, thin layers of aggregate are locally bound by cement paste where the structure shall arise. After completion of the printing process, the structure is excavated from the particle-bed and the unbound particles are removed. However, for a sufficient layer bonding and shape accuracy, the rheology of the cement paste must be adapted to the flow resistance of the particle-bed. For practical application, that means mostly time and material consuming “trial and error” tests. To prevent that, analytical models can help to predict the penetration of the cement paste. This paper presents four analytical models to calculate the penetration depth of a cement paste into a particle packing. Based on Darcy’s law, an already existing model is slightly modified (model A+) and a generalized (model C), an advanced generalized (model D) as well as a simplified model (model B/B+) are developed. Compared to conducted tests on the penetration depth, model B showed good accuracy (deviation <1.5 mm) for pastes with a yield stress ≥8.2 Pa, model A+/B+/C for ≥ 5.4 Pa and model D even for <5.4 Pa. Finally, an application guide for each model for practical use will be given.


Introduction
Due to the growing scarcity of resources, especially of aggregates for concrete production, and the high energy requirements as well as CO 2 emissions of cement production, lightweight construction and material saving are becoming increasingly important [1][2][3].
Consequently, the design principles in civil engineering must move away from a simple, compact construction to a more complex design that needs less material displaying multifunctionalities. First projects point out the enormous potential of saving resources by using additive manufacturing processes [4][5][6][7][8]. In addition, there are further ecological and economic advantages, such as faster production, elimination of formwork (and thus savings in waste and costs) [9][10][11][12][13][14][15] and the possibility of realizing new ways in the design of buildings (e.g., the principle that form follows force) [16][17][18][19] or new reinforcement concepts [20][21][22][23].
Many of the digital construction projects using additive manufacturing deal with depositing processes, also called extrusion. These enable one to produce large concrete components within a short time. In addition, there is usually no limitation by a construction space. A drawback of the depositing process is the lack of freedom of form, since strong overhangs can only be realized with the help of supporting structures [24][25][26]. An exception is shotcrete 3D printing technology, which enables also strong overhanging structures [27]. However, these techniques could require the use of a combination of additives and admixtures which increase the cost and the complexity of the mix-design of the material.
Overhangs of almost any complexity combined with much higher surface resolutions can be realized with selective binding processes [28], which were first used in 1995 by Pegna [29,30] for the production of concrete elements. A variant (which is applied in this paper) is the Selective (Cement) Paste Intrusion-SPI method [31,32]. Here, a particle-bed of aggregate is bound with cement paste where the component is to be produced; see Figure 1a-d. Due to the supporting function of the unbound particles, almost any overhanging construction with high surface resolution, high strength (up to 78 MPa) and high durability as well as density comparable to ordinary concrete (by filling up a present particle network) [26,[33][34][35] can be created, Figure 2. In contrast, the selective binding method Selective Cement Activation (SCA) uses an aggregate-binder (cement) mixture which is locally hardened by a water-based activator [28,[36][37][38][39]. The SCA method has already been in focus for building lunar outposts [40] and can be used for the production of structures made of lightweight concrete [41].
SCA offers a high surface resolution [42,43]. The surface resolution of SPI is in between the extrusion depositing processes and SCA. However, compared to SCA, SPI can achieve a faster construction speed due its coarser particles and its associated higher layer thickness.
The ability to control the penetration of the cement paste into the particle-bed is crucial for the sufficient layer bonding and the successful application of SPI (see Figure 1). This requires a profound knowledge of the interplay between the rheology of the cement paste and the flow resistance of the particle-bed.
In case of a material change or a batch changeover, it is currently still necessary to carry out material-and time-consuming "trial and error" tests to prove the success of the cement paste (particle-bed) layer thickness combination selected for production. This can be very uneconomical, especially when using large printers. In order to avoid this procedure in future, prediction of the penetration depth of the cement paste by means of simple material tests and analytical models is required. A first approach of an analytical model to calculate the penetration of cement pastes into particle-beds by Pierre et al. [44] showed good results to predict a penetration of the cement paste into the particle-bed under the assumption that the penetration ends at the border to the next layer of the particle-bed. Moreover, numerical simulation appears to be an alternative for the prediction of the penetration depth [45,46]. This paper presents analytical models which are able to calculate penetration depths which are even deeper than a particle-bed layer. Furthermore, the models in this paper try to take the flow behavior during the application of the cement paste on the particlebed into account, before the penetration starts, which enables a slight modification of a model already existing by Pierre et al. [44] (model A) to be used for penetration depths greater than the layer height (model A+). Additionally, a new approach for an analytical implementation of the pore structure of the particle-bed and of the process-parameters of the printer (applied volume as well as nozzle diameter) are considered. Therefore, a generalized and an advanced generalized model as well as a simplified model (with easier determination of the input parameters) will be developed based on Darcy s law [44,47,48]. The results of the analytical models will be compared to experimentally determined penetration depths which were varied by the rheological properties of the cement paste, the size of the aggregates in the particle-bed with constant packing density (flow resistance) and the humidity of the aggregates (dry and wet).

Selective Paste Intrusion (SPI) Method
The Selective Paste Intrusion (SPI) is able to create a link between the depositing processes with their high strengths and the particle-bed based processes with their high surface resolution combined with high durability. In SPI, only the aggregate is applied as dry material which is bound by a local intrusion of cement paste.
The experimental setup of the printer used for the investigations in this paper consists of a conical, round nozzle with an inner diameter b nz of 2.0 mm, which moves with a height h nz of 15.0 mm over the particle-bed with dimensions of 0.305 m to 0.375 m by means of an x-y gantry system. The height of the building space is 0.250 m, see Figure 3. The cement paste is pumped from a reservoir to the nozzle by a peristaltic pump controlled by a stepping motor. After the local application of the cement paste on the particle-bed, the building platform, on which the particle-bed is applied, is lowered by the set layer thickness and a new particle layer is applied (here 3.0 mm). After the printing is finished, the component is excavated from the building space.
A major advantage is the almost dust-free unpacking process and the good flowability of the unbound particles, which results in easy unpacking. The unconsolidated aggregate can be reused for a new print.

Cement Paste
An Ordinary Portland Cement (OPC, cement type: CEM I 42.5 R) was used as cement for all investigations. The grain size distribution of the cement can be taken from Figure 4. Demineralized water was utilized as mixing water. Furthermore, a polycarboxylate etherbased superplasticizer (PCE) with a solid content of 35.1 wt.%. was applied. The water content of the PCE was charged in the mixing water. The PCE was mixed into the mixing water. The cement was placed in the mixer and the solution of water with the PCE was added over the first 30 s of the mixing process. The first mixing section lasted 90 s, followed by a 120 s pause in order to return cement adhering to the bottom of the mixing container and to the mixing tool into the mixing process. Then, the mixing process was continued for another 90 s. The device used was an intensive mixer with a star type rotor (Eirich, type R 02) at fastest speed (stage 2 of 2 for mixing tool and container). The mixing water was pre-cooled to 1.5 • C using a cryostat in order to achieve a cement paste temperature of 20 ± 1 • C.
The investigations were performed using cement pastes with different water to cement ratios (w/c-ratio 0.30, 0.35 and 0.40) as well as mini slump flows (250 mm, 300 mm, 350 mm and 400 mm). An overview of the cement pastes with the superplasticizer contents used which allows to reach the targeted mini slump flows can be found in Table 1.  [50,51]. In addition, the viscosity (as well as the thixotropy) of the cement paste change with the w/c-ratio for the same mini slump flow. The density and the results of the rheological measurements of the cement paste can be found in Appendix A, Table A1.

Aggregates
Three fractions of a sieved and fire-dried quartz sand from the same deposit with a grain size of 0.7-1.2 mm (medium grain size d 50 = 1 mm), 1.0-2.2 mm (d 50 = 1.6 mm) and 2.0-3.2 mm (d 50 = 2.6 mm) were used. All fractions exhibit comparable densities as well as bulk densities and thus similar porosity contents (ratio between cavity volume to total volume), see Table 2. This causes a comparable requirement of cement paste in order to fill the cavities. Assuming a same surface roughness and moisture content, the flow resistance, is thus only determined by the d 50 or by the size of the gaps between the particles, respectively. The yield stress and the viscosity were determined in a rotational double-plate measuring system (diameter 50 mm, gap distance 1 mm, rheometer: Anton Paar MCR 502) 900 s after addition of the mixing water. The plates had a surface roughness with a depth of 0.5 mm to prevent wall slip. The profile started with an average shear rate of 40 s −1 for 10 s to achieve a complete structural break up (for the given shear history). Afterwards, the measurement was carried out with 19 descending steps of 80 s −1 to 0.02 s −1 with a step duration of 6 s. Preliminary tests showed that with this profile an equilibrium (an approximately constant value of the shear stress in every step) could be achieved in the steps between 60 s −1 and 2.5 s −1 for all measured pastes. Below 2.5 s −1 , the effect of thixotropy was already recognizable by a rising value of the shear stress.
The yield stress and the viscosity were calculated using the model of Herschel-Bulkley, which is a common model to describe the flow behavior of cement pastes [52][53][54][55] (see Equations (1) and (2)).
τ 0,HB (Pa) is the (Herschel-Bulkley) yield stress and k (Pa·s n ) is the consistency factor. In addition, the flow index n [-] describes a shear-thinning behavior in the range 0 < n < 1 and a shear-thickening behavior in the range 1 < n < ∞. If n = 1, it becomes a Bingham fluid.
The viscosity of Herschel-Bulkley is always given as a function of shear rate. Consequently, the shear rate-dependent, (Herschel-Bulkley) viscosity η( . γ) (Pa·s) is calculated according to Equation (2).

Penetration Tests
For the measurement of the penetration depth noted e, two strands of cement paste with a length of each 0.25 m per particle fraction were applied (see Figures 5 and 6). Thus, the cement paste could penetrate the particle-bed without reaching the bottom of the building space.
In order to quantify the effect of the process parameters, the penetration depth was investigated for constant cement paste properties (w/c-ratio 0.3 and mini slump flow of 300 mm) by changing (a) the velocity of the nozzle (by gantry movement speed) at the same volume output per length and by changing (b) the volume output per length at constant gantry movement speed (see Appendix B, Table A2).
All tests were carried out on both dry and pre-wetted (wet) aggregates. Pierre et al. have recently shown that aggregates water content affects the penetration depth [45]. For pre-wetting, aggregates were stored under water for at least 12 h and removed from the water before printing. The wet aggregate surface was dabbed with paper towels to quantify the dehydration effect on the cement paste by a dry aggregate. water before printing. The wet aggregate surface was dabbed with paper towels to quantify the dehydration effect on the cement paste by a dry aggregate.
The hardened strands were poured into epoxy resin, sawed apart lengthwise in the middle and sprayed with phenolphthalein in order to clearly distinguish the cement matrix from the aggregates. Then, the penetration depth e of the strands was measured every 0.01 m which results in 50 values for the penetration depth per series (see

Concept of Model Development
A first approach of an analytical model to calculate the penetration of cement pastes into particle-beds was proposed by Pierre et al. [44] (model A). Based on Darcy´s law and the Green and Ampt equation [44,47,48], model A showed good results to predict a penetration of the cement paste into the particle-bed under the assumption that the penetration ends at the border to the next layer of the particle-bed, which makes it impossible to generally use it in the SPI-printing process.
Therefore, this paper presents analytical models which are able to calculate penetration depths which are even deeper than a particle-bed layer. Furthermore, the models in this paper take the flow behavior during the application of the cement paste on the particle-bed into account, before the penetration starts, which enables model A also to be used for penetration depths greater than the layer height (model A+). Additionally, a new approach for an analytical implementation of the pore structure of the particle-bed and of the process-parameters (applied volume as well as nozzle diameter) are inserted directly in the model or indirectly by the calculation of the spread characteristic of the strands.
Thus, a generalized model (model C), an advanced generalized model (model D) as well as a simplified model (model B/B+), with easier determination of the input parameters, are developed based on Darcy´s law [44,47,48]. The results of the analytical models will be compared to experimentally determined penetration depths which were varied by The hardened strands were poured into epoxy resin, sawed apart lengthwise in the middle and sprayed with phenolphthalein in order to clearly distinguish the cement matrix from the aggregates. Then, the penetration depth e of the strands was measured every 0.01 m which results in 50 values for the penetration depth per series (see Figures 5 and 6).

Concept of Model Development
A first approach of an analytical model to calculate the penetration of cement pastes into particle-beds was proposed by Pierre et al. [44] (model A). Based on Darcy s law and the Green and Ampt equation [44,47,48], model A showed good results to predict a penetration of the cement paste into the particle-bed under the assumption that the penetration ends at the border to the next layer of the particle-bed, which makes it impossible to generally use it in the SPI-printing process.
Therefore, this paper presents analytical models which are able to calculate penetration depths which are even deeper than a particle-bed layer. Furthermore, the models in this paper take the flow behavior during the application of the cement paste on the particle-bed into account, before the penetration starts, which enables model A also to be used for penetration depths greater than the layer height (model A+). Additionally, a new approach for an analytical implementation of the pore structure of the particle-bed and of the processparameters (applied volume as well as nozzle diameter) are inserted directly in the model or indirectly by the calculation of the spread characteristic of the strands.
Thus, a generalized model (model C), an advanced generalized model (model D) as well as a simplified model (model B/B+), with easier determination of the input parameters, are developed based on Darcy s law [44,47,48]. The results of the analytical models will be compared to experimentally determined penetration depths which were varied by the rheological properties of the cement paste, the size of the aggregates in the particle-bed (flow resistance) and the humidity of the aggregates (dry and wet).  Models C and D use a novel approach to calculate the factors α and β, which are usually values determined by experiments to describe the pore structure [47,48]. This is done by an assumption of an effective capillary pore system. Model C describes a simpler variant, since only the yield stress is used as a rheological input parameter. In model D, the flow behavior in the porous medium is extended by the shear rate occurring in the pores . γ and the Herschel-Bulkley parameters n and k. This can improve the results especially for very flowable pastes.
All presented models assume a pressureless deposition of the cement paste on the particle-bed (laminar flow) and thus only the dead weight of the paste (hydrostatic pressure) as driving force of the penetration against the flow resistance of the particle-bed.
If additional pressure is applied via the nozzles, it has to be considered in addition to the hydrostatic pressure. The applied pressure can furthermore modify the flow regime itself and thus affect the reliability of the whole penetration process. For very high pressures, the type of flow could change from laminar to turbulent and the resulting loss of pressure could make it necessary to apply, e.g., the Forchheimer equation [56][57][58]. However, this is not part of the presented models.
The presented models assume that the penetration of the cement paste is considerably shorter than the setting and hardening time. Therefore, the penetration must be at least 10 times faster than the setting/hardening kinetics. This is also valid for the thixotropy of the cement pastes used in this investigation, which is neglected in this approach.
The change of process-technological aspects like (a) the velocity of the nozzle (by gantry movement speed) at the same volume output changing and (b) the volume output at constant gantry movement speed (see Section 2.3.2 and Appendix B, Table A2) showed no significant effect on the penetration depths in the selected test setup; see also [34]. Therefore, the velocity of the nozzle as an influencing factor is neglected for the models presented in this paper. However, the deposited volume of material could have an effect and is nevertheless implemented.

Theoretical Background of the Models
Many studies of a flow through porous media of non-Newtonian fluids concentrate mainly on the description of viscous power laws fluids without having a yield stress [59][60][61][62][63][64]. However, especially when considering the flow behavior of concretes or cement pastes, e.g., in the field of oil well cementing, soil injection or the simulation of the flow of Self Compacting Concrete (SCC) through rebars [65][66][67], the consideration of the yield stress is decisive.
Chevalier et al. [47,48] combined Darcy's law with the Herschel-Bulkley model assuming spherical particles in the particle-bed in order to model the flow of a fluid with a yield stress through a porous medium, see Equation (3).
In which D (m) is an undefined unit of length that characterizes the porous medium and ∇p (Pa/m) is the pressure loss over this unit of length. α (here = 5.5) and β (here = 85) are two unknown parameters which must be determined experimentally by flow measurements.
α depends on the size ratio of the particles in the fluid and the widest path between the particles of the particle-bed. β describes the flow resistance of the entire particle-bed depending on the pore size distribution and structure. β also considers the permeability K (m 2 ) of the particle-bed. Assuming a Newtonian fluid (τ 0 = 0 and n = 1), the following relationship can be obtained using Darcy's law and the Kozeny-Carman equation [44,68], see Equation (4).
The ratio of the flow velocity of the fluid v (m/s)/D (m) is equivalent to the shear rate . γ (s −1 ) occurring between the particles in the porous medium. In Equation (3), the first term including the yield stress describes the mobilization of the fluid in the widest path through the porous medium. As the shear rate increases, the fluid also flows in finer pores (second term). Thus, the ratio of the pressure drop increases with the activation of smaller pores [47].
However, when using this model, it should be noted that the particles in the fluid should be much smaller than the particles of the particle-bed.
However, the literature [47,56,63,69] presents very large ranges of α (0.98-5.5) and β (1.23-102.7) and mutual dependencies. Therefore, the coefficients should always be determined experimentally. However, in the case of cement-based suspensions, this is very difficult due to the time-dependent flow properties. Consequently, this paper shows a new assumption for the calculation of these coefficients.

Spread Characteristic of the Cement Paste During Printing
The models of this paper assume that during the printing process the cement paste does not only flow vertically from the nozzle into the particle-bed, but spreads also horizontally on the particle-bed before it penetrates (see Figures 8 and 9). Application of the cement paste during printing process).
The applied cement paste strand with an initial height H 0 (m) and width b 0 (m) will spread aside to a width b eff (m) and a related effective height H 0,eff (m) which determines the effective hydrostatic pressure as driving force for the penetration.
By using fine nozzles and thus thin strands of cement paste comparing to the applied length l (m), we can assume a rectangular cross section of the strand. Following [50,70,71] and under the condition that the cement paste is not penetrating the particle-bed before the strand has spread; we can further assume that the final spread b eff is depending on the yield stress of the cement paste. Thus, the yield stress determines the final height H 0 . eff (and therefore the effective hydrostatic pressure) of the strand. The applied cement paste strand with an initial height H0 (m) and width b0 (m) will spread aside to a width beff (m) and a related effective height H0,eff (m) which determines the effective hydrostatic pressure as driving force for the penetration. By using fine nozzles and thus thin strands of cement paste comparing to the applied length l (m), we can assume a rectangular cross section of the strand. Following [50,70,71] and under the condition that the cement paste is not penetrating the particle-bed before the strand has spread; we can further assume that the final spread beff is depending on the yield stress of the cement paste. Thus, the yield stress determines the final height H0.eff (and therefore the effective hydrostatic pressure) of the strand. Therefore, Equation (5) describes the final condition (h = Heff,0) of the spread before penetrating the particle-bed.
Here, τ0 (Pa) is the yield stress, ρp (kg/m 3 ) is the density of the cement paste, g (m/s 2 ) is the gravity and bnz (m) the diameter of the nozzle.
However, the effective width beff of the strand is still not known. Again, under assumption of a rectangular cross section and already achieved final spread beff (bnz = 0), the  The applied cement paste strand with an initial height H0 (m) and width b0 (m) will spread aside to a width beff (m) and a related effective height H0,eff (m) which determines the effective hydrostatic pressure as driving force for the penetration. By using fine nozzles and thus thin strands of cement paste comparing to the applied length l (m), we can assume a rectangular cross section of the strand. Following [50,70,71] and under the condition that the cement paste is not penetrating the particle-bed before the strand has spread; we can further assume that the final spread beff is depending on the yield stress of the cement paste. Thus, the yield stress determines the final height H0.eff (and therefore the effective hydrostatic pressure) of the strand. Therefore, Equation (5) describes the final condition (h = Heff,0) of the spread before penetrating the particle-bed.
Here, τ0 (Pa) is the yield stress, ρp (kg/m 3 ) is the density of the cement paste, g (m/s 2 ) is the gravity and bnz (m) the diameter of the nozzle.
However, the effective width beff of the strand is still not known. Again, under assumption of a rectangular cross section and already achieved final spread beff (bnz = 0), the  Therefore, Equation (5) describes the final condition (h = H eff,0 ) of the spread before penetrating the particle-bed.
Here, τ 0 (Pa) is the yield stress, ρ p (kg/m 3 ) is the density of the cement paste, g (m/s 2 ) is the gravity and b nz (m) the diameter of the nozzle.
However, the effective width b eff of the strand is still not known. Again, under assumption of a rectangular cross section and already achieved final spread b eff (b nz = 0), the applied volume V l (m 3 ) can be calculated with the length l (m) of the strand and H 0,eff using Equation (7).
This leads to Equation (8) for calculation of b eff .
Consequently, the yield stress τ 0 can be calculated following equation (9), when b eff is known.
However, the final height H 0,eff can also simplified be approximated by the final height h(R) of the cement paste after the mini slump flow test following [50,51], which was used in [34] (see Appendix C).

Modification of Model A to Model A+
A detailed derivation of model A (see Equation (10)), can be found in [44].
Here, ϕ is the solid fraction of the porous medium, h lay (m) the layer height of the applied particles and κ the Janssen parameter [72].
In order to make the calculation also valid for penetration depths greater than the layer height h lay , model A has to become independent to this parameter (see [34]). Therefore, all parts of Equation (10) which describe the height of the cement paste depending on h lay have to be replaced by Equation (6) which leads to Equation (11).

Theoretical Background of Model B/B+
For a first and simple approach in model B/B+ the β-term associated to the viscosity effect will be neglected. Thus, Equation (3) simplifies to Equation (10).
Consequently, as limiting boundary condition for a penetration of the cement paste, the yield stress must firstly be overcome [44] (see Equation (12)).
For model B/B+, it is assumed that the medium particle diameter d 50 (m) can be used for D as the determining factor for the flow resistance of the particle-bed. According to Chevalier et al. [47], coefficient α is set to 5.5. Furthermore, the effective hydrostatic pressure decreases over time. This is taken into account by a linear regression [73], which can be expressed by dividing the effective hydrostatic pressure H 0,eff from Equation (6) by 2.
Following this assumption and that the hydrostatic pressure turns to zero after reaching the final penetration depth e (∇p = ρ p · g H 0,eff/ e), Equation (12) turns into model B/B+ in Equation (14).
The difference between models B and B+ is only the measurement of the yield stress. Model B does not need a rheological measurement by a rheometer but uses the assumption of Roussel [50] to determine the yield stress only by the slump flow value. In contrast, model B+ utilizes the results of the rheometric measurements following the model of Herschel-Bulkley, which probably enables more accurate results.

Theoretical Background of Model C and D
Since model B/B+ use for α values from literature or difficult to determine experimental values, models C and D calculate the values for α and β according to the pore structure.
Therefore, we assume that the porous medium of the particle-bed can be described by a capillary pore structure between the particles. Consequently, the determinant factor for the flow resistance turns to the effective pore diameter d eff (m) (=D).
Chevalier et al. [47] describe that α defines the penetration ability of a fluid which depends on the widest path between the particles in the particle-bed or on the porosity of the particle-bed ε, respectively. Consequently, α can be assumed according to Equation (15).
In contrast, β describes the flow resistance of the entire particle-bed as a function of pore size distribution and structure. Thus, β designates the mobilization of the fluid in the pores of the particle-bed (porosity ε) depending on the limiting yield stress τ 0,lim (Pa) (depending on the pore size distribution and structure) as well as the actual yield stress of the fluid τ 0 . This results in the following assumption for factor β, see Equation (17).
Implementing Equations (15) and (17) in Equation (3) under the further assumption of the boundary conditions of model C gives model D in Equation (18).

Determination of the Effective Pore Structure and Shear Rate in the Particle-Bed
Assuming that the particle-bed consists of a dense, hexagonally arranged packing (hexagonally densest packing or cubically densest packing), tetrahedral or octahedral gaps occur between the particles [63,74]. If the tetrahedral gaps are approximated as an equilateral triangle and the octahedral gap as a square or two isosceles triangles, the area and circumference can be determined (see Figure 10).
If the area of the octahedral gap A octa (m 2 ) is derived from two isosceles triangles (see Equation (19)), the circumference of the octahedral gap U octa (m) can be determined in Equation (20).
The area of the equilateral triangle of the tetrahedral gap A tetra (m 2 ) can be determined according to Equation (21). The circumference of the tetrahedral gap U tetra (m) is calculated according to Equation (22).
If the area and circumference of the pore is known, an equivalent effective pore radius r eff [m] of a round capillary can be calculated according to Lopez et al. [63] in Equation (23).
G can be calculated using Poiseuille's law according to Equation (24) and form factor χ [75], see Table 3.  Table 3. Values of the form factor χ for different cross sections of the pore according to [75]. With the effective equivalent radius r eff an equivalent capillary pore system can now be approximated. In a pressureless capillary pore, the shear rate . γ can be given as a function of the effective pore radius r eff and the shear rate-dependent viscosity η [73] (see Equation (25)).
There, K (m 2 /Darcy) is the absolute permeability coefficient. This coefficient can be calculated by analogies for the simulation of flow movements in porous media. Tamayol et al. [69,81] describe the pore structure of a three-dimensional porous medium simplified with a three-dimensional cubic configuration of a lattice with a flow resistance that can be expressed in a lattice permeability coefficient K G [m 2 ] (see Equation (27) and Figure 11a. The constant 0.08 describes the arrangement of the 3D lattice for a wide range of porosities and was found by Tamayol and Bahrami [69] by comparison with numerical data from Higdon and Ford [82]. S G [m] is the centric distance of the non-permeable part of the lattice and d G [m] the diameter of the lattice bars. The parameters S G and d G in three-dimensional structures can be expressed by the relationship in Equation (28) [69,82].
In this case, S G (m) and d G (m) are unknown parameters and ϕ is the solid ratio. However, if we now assume that the lattice system itself can be flowed through, and that the sections between the lattice no longer have porosity, the three-dimensional capillary system of the particle-bed results (see Figures 11b and 12). Now the porosity of the particle bed K 3D (m 2 ) for octahedral d eff , octa (m) or tetrahedral gaps d eff,tetra (m) can be represented according to Equations (29) and (30).
According to Hagen-Poiseuille, a limiting boundary condition for the penetration of the cement paste can now be given in this capillary system analogous to Equation (13) in combination with the largest effective radius of the pore system r eff,octa (see Equation (31)).
Since the dynamic viscosity η is dependent on the shear rate . γ, the shear rate according to Equation (25) still cannot be calculated despite L eff being known. However, the limiting yield stress τ 0,lim,C,D according to Equation (31) can be used, which characterizes the transition of the fluid into the flow state. By combining Equations (2), (25) and (31), Equation (32) can be derived for the calculation of the shear rate . γ in the pores of the particle-bed.  The measured penetration depths (white and blue circles) increase with increasing slump flow (decreasing yield stress) and increase slightly with increasing w/c-ratio (decreasing viscosity). Thus, the yield stress exhibits a dominant role in comparison to the viscosity for the penetration of the cement paste. Furthermore, the penetration depth e increases with an increasing d 50 (decreasing flow resistance of the particle-bed) and with increasing humidity of the particle-bed (circles in white = dry and blue = wet). These observations correspond to the results in [25,26,28,31,32,44].

Results and Discussion
Model A+ exhibits a good correlation to the penetration depths up to a mini slump flow of 350 mm with a deviation of ≤0.6 mm for dry aggregates and ≤1.3 mm for wet aggregates. However, the deviation for a slump flow of 400 mm increases to ≤3.9 mm for dry and ≤3.3 mm for wet aggregates.
Model B show only small deviations from the measured penetration depths when considering the dry aggregates in the range of a mini slump flow from 250 mm to 300 mm. For a low yield stress (mini slump flow 350 mm and 400 mm), the deviation of model B strongly increases with increasing particle size d 50 (decreasing flow resistance of the particle-bed). Thus, model B overestimates the penetration depths with up to 10.2 mm (approx. +270%). The results show furthermore, that a determination of the rheolocal properties (model B+) leads to better results than using an approximation of the yield stress by the mini slump flow due to the shear thickening flow behavior of the cement paste (see results in Appendix A, Table A1). The tendency to exhibit an increased deviation with decreasing yield stress and increasing d 50 was also shown in Pierre et al. [44] for model A and model A+. This effect could be caused by neglecting the β-term of Equation (3). Therefore, the effect of a shearthinning or shear-thickening flow behavior (see results in Appendix A, Table A1) cannot be reproduced in Model A and Model B/B+, although its effect increases in a porous medium with decreasing flow resistance.
Model C shows only small absolute deviations of ≤1.3 mm to the measured penetration depths for the wet aggregate up to a mini slump flow of 350 mm. However, model C shows a better agreement with the penetration depths using dry aggregate with absolute deviations of ≤0.7 mm for the same mini slump flow range. Furthermore, considering the effective pore structure in the α-term leads to better results compared to model B/B+. Comparing the results for model A+ and C, the calculated penetrations depths are almost congruent, though the approach of both models is different. However, both models seem to be able to describe the structure of the porous medium quite accurately. However, neglecting the β-term also leads to relatively high deviations for cement pastes with low yield stress (mini slump flow 400 mm) for model A+ and model C.
Model D exhibits only slight deviations of ≤1.4 mm for dry and ≤1.3 mm for wet aggregates for all penetration measurements. By applying the β-term from Equation (3), model D succeeds in eliminating the major deviations even for the flowable cement pastes with 400 mm mini slump flow.
Although the penetration depths continue to increase due to the increasing pressure gradient (τ 0,lim /τ 0 ), they are probably simultaneously reduced by the internal friction of the cement paste (respectively the dynamic viscosity, expressed by k· γ .. n ) and the friction on the particle surfaces. In addition, the flow through the particle-bed strongly depends on the prevailing pressure gradient and the shear rate in the pores. Due to the friction loss, the pressure gradient decreases rapidly and the flow stops more quickly as the pore size decreases. Thus, the penetration depths are additionally reduced. This effect increases with increasing ratio τ 0,lim /τ 0 .
The results for model C and D in Figure 13 were calculated using tetrahedral gaps. The assumption of octahedral gaps showed less accurate results for the penetration depth which can be explained by the shape and bulk density of the used aggregates (see Appendix D).
The velocity of the nozzle as well as the applied volume of the cement paste (see Section 2.3.2 and Appendix B, Table A2) showed no significant influence on the penetration depth for the used printer set up (see Appendix B). However, especially higher amounts of cement paste per length could have an effect on the penetration depth.
The investigations show that model D, in contrast to the other presented models, is capable of validly calculating the penetration of the cement paste across all investigated properties of the cement pastes and particle-bed by using a pore system with tetrahedral gaps. An additional discussion of the theoretical background of the models can be found in [34,[44][45][46].
Finally, the following procedure is recommended as a guideline for a successful application of SPI in practice (see Figure 14).

Conclusions and Outlook
In order to prevent effortful "trial and error" tests to validate new material combinations in the SPI printing process, we presented four models to calculate the penetration depth of a cement paste in a particle-bed based on material and process technological parameters. All models showed good results for cement pastes up to medium yield stress (mini slump flow of 300 mm to 350 mm). However, model C and D use a new approach to calculate the pore structure of the particle-bed. Model D achieved the best correlation to the real penetration measurements, even for very fluid cement paste with low yield stress (mini slump flow of 400 mm) for which viscous effects cannot be neglected.
However, the models should also be validated for other printer setups (e.g., other nozzle diameters or distances to the particle-bed). Furthermore, an expansion of the models for turbulent flow conditions could be interesting as soon as the cement paste is jetted out of the nozzle and is not deposited pressureless like for the investigations in this paper.  Data Availability Statement: Data sharing is not applicable to this article.

Acknowledgments:
The authors gratefully acknowledge the funding by the Bavarian-French University Cooperation Centre (BFHZ-CCUFB) as well as the German Research Foundation (DFG). D.W. thanks the Ingenieurbüro (engineering office) Schiessl Gehlen Sodeikat GmbH in Munich, Germany, for the opportunity to write this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The density of the fresh cement paste was determined volumetrically according to [83]. The rheological properties were calculated according to Heschel-Bulkley. The results can be found in Table A1.

Appendix B
In order to investigate the effect of process-related parameters on the penetration behavior of the cement pastes, a) the speed of the gantry was varied (1000 mm/min, 2000 mm/min, 3000 mm/min) for a constant application volume of 5.54·10 −5 m 3 , as well as b) the application volume (2.92·10 −5 m 3 , 5.54·10 5 m 3 , 8.54·10 −5 m 3 ) for a constant gantry speed of 2000 mm/min (see Table A2). For these tests, a cement paste with a w/c ratio of 0.3 and a mini slump flow of 300 mm was used.  Figure A1a shows the results of the variation of the gantry speed. It can be seen that the gantry speed seems to have no influence on the penetration depth. Only the penetration depth in dry aggregate with a d 50 = 0.0016 m shows a gradient with increasing speed. However, the slope of the penetration depth lies within the scatter of the measurement results and can therefore be neglected. Figure A1b exhibits the penetration depth depending on the applied volume of cement paste per length. Here, an increase of volume from 5.54·10 −5 m 3 to 8.32·10 −5 m 3 does not increase the penetration depth. However, halving the applied volume to 2.92·10 −5 m 3 appears to have an influence on the penetration depth for all aggregate fractions, which increases with increasing medium grain diameter d 50 . However, for the smallest cement paste volume per length no continuous paste application was possible with the printer set up.
This leads to an inconsistent penetration front and a probably changed propagation behavior. Therefore, this value is not considered to be usable, as a continuous application is aimed.
Thus, the gantry speed does not represent an input parameter for the analytical calculation model. However, especially for larger applied amounts of cement paste the used volume per length should nevertheless be taken into account. assumption of too large pore diameters which leads to an overestimation of the shear rates occurring in the pores and therefore an underestimation of the pressure gradient (ratio τ 0,lim /τ 0 ) as well as the flow resistance of the particle-bed. Figure A4 shows the longitudinal cross section through of a penetration measurement. The samples exhibit that the pore system between the particles is dominated by tetrahedral gaps. Consequently, when considering the flow resistance under free flow in the β-term, the smaller and finest pores play a decisive role [47]. In contrast, the largest path between the pores determines whether or not flow begins in the pores (α-term).