Article

# Additive Manufacturing of Cementitious Materials by Selective Paste Intrusion: Numerical Modeling of the Flow Using a 2D Axisymmetric Phase Field Method

Alexandre Pierre<fnref n="1" />, Daniel Weger<fnref n="2" />,*, Arnaud Perrot<fnref n="3" />, and Dirk Lowke<fnref n="4" />

Received: 15 October 2020; Accepted: 5 November 2020; Published: 7 November 2020

<figure><img src="image_1.png" /></figure>

check for
updates

**Abstract:** The 3D printing of concrete has now entered a new era and a transformation of the construction sector is expected to reshape fabrication with concrete. This work focuses on the selective paste intrusion method, which consists of bonding dry particles of aggregate with a cement paste. This innovative technique could lead to the production of very precise component for specific applications. The main obstacle to tackle in order to reach a high shape accuracy of high mechanical performances of 3D printing elements by selectively activating the material is the control of the distribution of the cement paste through the particle bed. With the aim to better understand the path followed by the solution as it penetrates a cut-section of the granular packing, two-dimensional numerical modeling is carried out using Comsol software. A phase-field method combined with a continuous visco-plastic model has been used to study the influence of the average grain diameter, the contact angle, and the rheological properties of cement pastes on the penetration depth. We compare the numerical modeling results to existing experimental results from 3D experiments and a one-dimensional analytical model. We then highlight that the proposed numerical approach is reliable to predict the final penetration of the cement pastes.

**Keywords:** 3D printing; selective paste intrusion; yield stress; cement; numerical simulation

## 1. Introduction

Digital construction using the 3D printing of concrete has been democratized in many countries around the world as large scale projects have emerged and scientific knowledge have grown rapidly in the last few years [1,2]. The free-form possibility allowed by 3D-printers is one of the opportunities among other expected advantages, but several challenges regarding both materials and processes are remain in the digital concrete industry [3,4].

The most common method is the deposit of the cementitious material using digitally controlled extrusion process [4–10]. The basic principle is to combine a well-defined and precise process using robots with a material for which the rheological behaviour and structural build-up are well characterized. For extrusion-deposition techniques, the material must be fluid enough to be extruded from the nozzle, but the success of the printing process mainly depends on the fine-tuning of the construction rate according to the material structural build-up [9]. The rheological requirements for printable concrete deal with several physical properties, such as yield stress, thixotropic effect, critical strain, or elastic modulus,

<footnote n="1">¹ L2MGC EA4114, CY Cergy Paris University, 5 mail Gay-Lussac—Neuville-sur-Oise, 95031 Cergy-Pontoise, France; alexandre.pierre@cyu.fr</footnote>

<footnote n="2">² Centre for Building Materials (CBM), Technical University of Munich, Baumbachstraße 7, 81245 Munich, Germany; weger@ib-schiessl.de</footnote>

<footnote n="3">³ IRDL, UMR CNRS 6027, Université Bretagne Sud, F-56100 Lorient, France</footnote>

<footnote n="4">⁴ Institute of Building Materials, Concrete Construction and Fire Safety (iBMB), University of Braunschweig, Beethovenstr. 52, 38106 Braunschweig, Germany; D.Lowke@ibmb.tu-bs.de</footnote>

<footnote>* Correspondence: arnaud.perrot@univ-ubs.fr; Tel.: +33-(0)1-34-25-69-36</footnote>
as reported in [11]. In addition to the social, economic, and durability requirements, the scientific
challenge is huge and recent work shows that the rheological requirements of the cement based materials
used are dependent on both the process and object shape [11].

An alternative is to use an additive layer-wise process based on the a selectively activated powder-bed by injecting a fluid or an activator [2,4]. After successive depositions of thin layers of powder, complex shapes and overhang structures can be produced by selectively depositing the fluid activator which ensures the particles binding at the hardened state. Unbounded particles are then removed and could be used for other applications or even 3D printing [2]. Finally, a post-print treatment could be applied to improve the final properties of the structure regarding to thermal or mechanical performances [2].

In principle, particle-bed 3D printing methods can be divided in selective binder activation, selective paste intrusion and binder jetting [2]. Depending on the application, selective binder activation can be used: a mixture of aggregate and binder is poured and the hardening will be activated by spraying or jetting of a solution of water and admixture [2,12–15]. Recent works [2,16–18] used the selective paste intrusion method to build cement based material components: the binder composed of water, cement, and admixtures is applied to a bed of sand particles. Considering the binder jetting, the binder can be a resin which reacts with a constituent of the particle bed, but this last technique is less used for civil engineering applications nowadays.

Advantages of the particle-bed 3D printing techniques are that the final shape of the element or the structure only depends on the injection and the penetration of the binder as the sand particle bed is mechanically stable at rest. Furthermore, cantilevers or other overhanging elements can be built without the need for additional support structures due to the presence of unbound particles. This technique aims to be resource efficient compared to the extrusion method as fewer polymer additives are needed to produce complex structural elements. The selective paste intrusion method has already allowed to produce high strength concrete parts [2,16–18]. Nevertheless, the selective paste intrusion has not yet been fully transferred to the concrete construction industry. However, good material properties depend on a sufficient penetration of the layers. As the future strength and final shape accuracy is related to the degree of paste penetration, it is a major concern for the additive manufacturing of cementitious materials. Therefore, the penetration of the binder should be predicted depending on several factors, such as the rheological properties of the liquid binder, the particle-bed morphology, and the interface properties between fluid (cement suspensions) and solid (particles). A way to optimize the process could be the application of numerical modeling to predict quantitative results of the final penetration depth.

The numerical modeling of fresh concrete and cementitious material flows has gained interest in the last decade and a special technical letter concerning the new insight and challenges about this subject has been recently published [19]. For 3D printing with particle-based methods, the challenge is to capture the evolution of free surface with time as the moving interfaces of two-phase flows to capture the fluid flowing through air content of the porosity of the particle-bed. Three different approaches can be used: interface tracking where materials flows through a fixed mesh (Eulerian), interface capturing where the motion of materials particle is tracked (Lagrangian), or ALE (arbitrary Lagragian–Eulerian), which combines both but where free surface movement is limited [20]. In Eulerian tracking methods, an additional field variable is used to describe the interface between two fluids or a solid and a fluid. As reported by Mokbel et al. [21], several interface capturing approach can be used: level-set [22], volume of fluid [23], and phase field method [24]. Volume-of-fluid (VOF) could have been used, but it can suffer from inaccurate computations of the interface properties as reported in [25]. The VOF method is then less accurate for 3D printing with selective paste intrusion where surface tension could be the dominant effect. Also, the phase-field model is a front capturing method where a field variable is used to indicate each phases of both fluids while dealing with diphasic 3D printing process. This method includes the mass conservation and transport stabilization, which is not the case while dealing with
level-set methods or other Eulerian methods. Consequently, phase field methods allow for coupling
interface advection and flow equations in order to prevent time step restrictions of the used solver.

It is commonly admitted that cement pastes are viscoplastic fluids which behave as Bingham
fluid. Its flow through sand particle layer depends on its rheological properties: yield stress, viscosity,
the morphological properties of the sand particle-bed (layer height, medium grain diameter, shape,
roughness, void volume fraction), the sand wettability, and absorption properties. Theoretical work
validated by experimental results established the general relation between pressure drop as a function
of the penetration rate for a yield stress fluid flowing through a porous medium [26]. This physical
approach has recently been adapted to develop and validate an analytical model to evaluate the
penetration of cement paste through sand particle-beds assuming no-slip between aggregates and
cement paste for selective cement paste activation [27]. In addition, a recent analytical and digital
modeling study [28] presents four scenarios that may arise when a drop of Newtonian fluid is deposited
into a particle bed: (1) the drop does not bond the grains due to a lack of penetration, (2) the drop
remains trapped between the grains, (3) the drop descends to the bottom of the layer, and finally, (4) the
drop splits and spreads around the particles. These four scenarios can occur for non-Newtonian fluid
like cement paste until the pressure applied to the fluid overcomes the magnitude of its yield stress.
Within the above frame, the influence of the rheological properties that affect the penetration of the
binder through a sand layer bed will be discussed in this study. In addition, special attention will be
paid to the influence of the relationship between surface tension and the forces of gravity.

Then, the aim is to provide new information about the influence of surface tension effect and
contact angle on the path followed by the solution as it penetrates the granular matrix. It is necessary
to ask questions about the parameters to be controlled to produce homogeneous structures. In this
study, we first present the Fourier–Voronoi-based generation from the work of Mollon et al. [29] to
generate three sand bed packings which differ in their the average particle diameter and which have
been used previously in the experimental measurements of Pierre et al. [27]. Then, the phase field
method and the modeling of the rheological behaviour following a continuous visco-plastic model are
presented to follow the path of the suspensions through the powder-bed. It is shown that the contact
angle at the fluid–particle interface, the average grain diameter, and the rheological properties of the
suspensions govern the kinetic of the flow and the final penetration depth. We then highlight that
the numerical modeling is able to provide a quantitative estimation of the penetration final depth.
A comparison of the final penetration depth with existing 3D experimental results and a 1D analytical
approach is also carried out.

**2. Materials and Methods**

2.1. *Sand Packing Generation*

In this work, a Fourier–Voronoi-based generation is used to generate the sand packing in 2D
following the open-source program of Mollon and Zhao [29]. Using a Matlab code, three sand packings
layers of 12.5 mm width and 20 mm height were generated to model the experimental conditions of
the measurements performed in Pierre et al. [27]. The void fraction of sand packings, ie. the ratio of
the volume of the air phase occupied in the aggregate layer, was experimentally determined to a value
of 0.46 by Pierre et al. [27]. The results of the three sand packing layers obtained are shown in Figure 1
for the average particle diameter, noted $d_{agg}$, of 1.0 mm, 1.6 mm, and 2.6 mm, respectively.
<figure><img src="image_2.png" /><figcaption>Figure 1. Sand packing ($\Phi_p$ = 0.46) generated with a Fourier–Voronoi-based generation from the program of Mollon et al. [29].</figcaption></figure>

At first sight, the 2D formulation of the proposed method may seemingly limit its practical use for real 3D particle modeling [29]. In order to qualitatively compare theses modeling results to real case, the numerical assembly in two-dimensions is considered as a cut-section of a 3D granular packing, as illustrated in Figure 2. It can be observed that, both in real (Figure 2) and in the modeling cases (Figure 1), most of the sand grains do not have an interconnection with each other. Also, the void fraction is constant for all the sand packings generated as it is computed by the code of Mollon and Zhao [29]. Then, the schematic sand packing used and schematized in Figure 1 will be used to model the flow path of the cement paste.

<figure><img src="image_3.png" /><figcaption>Figure 2. (a) Picture of the sand packing in 3D; (b) Picture of the 2D cut-section of the sand packing assembly.</figcaption></figure>

## 2.2. Modeling of the Rheological Behaviour

A Papanastasiou model [30] is used to describe the Bingham rheological behaviour of the cement pastes as already proposed in a previous study from the authors using a level-set method [31]. The review of Mitsoulis [32] showed that the regularized continuous visco-plastic model of Papanastasiou [30] is reliable to approximate the Bingham model. In this model, the shear rate is linearly linked to the shear stress through the dynamic viscosity of the penetrating fluid that is expressed following Equation (1):

$$ \eta_{paste} = \eta_{exp} |\dot{\gamma}|^{n_{exp}-1} + \frac{\tau_{0,exp}}{|\dot{\gamma}|} [1 - e^{(-m|\dot{\gamma}|)}] \quad (1) $$

The value of the exponent $n_{exp}$ is 1 as the suspensions of cementitious materials behave as a Bingham fluid [11,33,34] ($\tau = \tau_{0,exp} + \eta_{exp}\dot{\gamma}$) and also showed a reliable fitting with the rheological
measurements used in this study and carried out by Pierre et al. [27]. Values of experimental yield stress $\tau_{0,exp}$ and viscosity $\eta_{exp}$ are evaluated from the fitting of the Bingham model on the measured data from rheological measurements on two different cement pastes carried out in a previous work of Pierre et al. [27] are used in this study and can be seen in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3. Rheological behaviour of the cement pastes modeled by the Bingham model from experimental measurements of Pierre et al. [27] approached with the Papanastasiou model.</figcaption></figure>

The shear stress-shear rates flow curves for the Bingham model obtained from experimental measurements of Pierre et al. [27] are plotted with the best fitting curve of the Papanastasiou model. It can be noted that the Papanastasiou model perfectly fits the Bingham rheological behaviour for shear rates higher than 0.01 s⁻¹ with a value of m parameter of 500 s. The value of m parameter is a fitting parameter chosen to adjust the reliable fitting domain [30]. This condition was always fulfilled to evaluate reliable penetration depth for all the numerical simulations as shown by Pierre et al. [31].

## 2.3. Numerical Simulation Procedure

The numerical simulations were carried out using Comsol® software and the phase field method in a 2D axi-symmetric condition. The phase-field model in the conservative form used in this study breaks the Cahn-Hilliard equation into two second-order partial differential equations [35,36]:

$$ \frac{\partial \varphi}{\partial t} + u \nabla \cdot \varphi = \nabla \cdot \frac{\gamma \lambda}{\varepsilon^2} \nabla \psi \qquad (2) $$

$$ \psi = -\nabla \cdot \varepsilon^2 \nabla \varphi + (\varphi^2 - 1)\varphi \qquad (3) $$

where $\varphi$ denotes the phase field variable, $\psi$ an auxiliary variable to separate fourth order equations into two second order equations, $t$ the time, and $\varepsilon$ the parameter determining the thickness of the interface; $\gamma$ is the mobility (m³·s/kg) which determines the relaxation time of interface and timescale of diffusion. This parameter $\gamma$ is determined by the relation $\gamma = \chi \varepsilon^2$ where $\chi$ is a mobility tuning
parameter fixed to a default value of 1 m·s/kg; $\lambda$ is the mixing energy quantity (N), and $u$ is the fluid velocity (m/s).

The maximum mesh element will be taken as the incertitude while evaluating the penetration depth of yield stress fluid through sand packing. The mobility and interface thickness parameters are related to the surface tension coefficient $\sigma$ (N/m) and defined through Equation (4):

$$ \sigma = \frac{\sqrt{8} \lambda}{3 \varepsilon} \qquad (4) $$

The phase field method finally solves: the Cahn–Hiliard Equation (2), a continuity equation Equation (5), and the Navier–Stokes Equation (6).

$$ \nabla \cdot u = 0 \qquad (5) $$

$$ \rho \left( \frac{\partial u}{\partial t} + (u \cdot \nabla) u \right) = -\nabla P + \nabla \cdot \left( \eta \left( \nabla u + (\nabla u)^T \right) \right) + \rho g \qquad (6) $$

From Equation (6), it is shown that the unique driving force for the penetration is the gravity, because no initial velocity, i.e., $u = 0$, is considered in this study. The mass conservation is checked for each simulation in this paper as it could be a problem to ensure a careful implementation of the numerical modeling. The density and the dynamic viscosity were computed using the following equations:

$$ \rho = \rho_{air} + (\rho_{paste} - \rho_{air})\varphi \qquad (7) $$

$$ \eta = \eta_{air} + (\eta_{paste} - \eta_{air})\varphi \qquad (8) $$

Figure 4 illustrates the different domains (air in blue and cement paste yield stress fluid in red) and the meshing used for the simulations at the initial time of the numerical simulations.

<figure><img src="image_5.png" /><figcaption>Figure 4. Domains where cement paste is in red and air in blue and meshing used for the numerical simulations using CFD module from Comsol® software.</figcaption></figure>

## 2.4. Contact Angle and Boundary Conditions

Pierre et al. [27] showed from experimental results that the water content of the sand affects the penetration depth during 3D printing process like the selective paste intrusion method. Sand particles
absorb some of the water from the cement paste leading to an evolution of the yield stress and the viscosity during its penetration through the particle-bed. The effect of the surface wetting and absorption of water by dry particles of sand on fluid penetration could be modelled with different boundary conditions at the sand-percolating fluid interface. Considering the numerical simulations, a way to study the fluid-solid contact influence on the penetration is to analyze the influence of the contact angle at the interface on the penetration. The contact angle values studied ranged from 3π/8 to 5π/8 (see Figure 5). At the domain boundaries, the contact angle was kept constant throughout all simulations at a value of 90° to avoid the fluid being attracted to the walls as proposed by Boyce et al. [28]. Also, the no-slip condition is defined as a fluid velocity which tends to 0 at the fluid-solid interface.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Contact angle at the interface cement-air-sand: (a) $\theta = \pi/2$; (b) $\theta = 3\pi/8$; (c) $\theta = 5\pi/8$.</figcaption></figure>

The slip length is the height of the element size taken between 0.025 mm and 0.563 mm depending of the medium diameter of the sand bed particles. The discretization scheme, the solver settings and the time step width are standard setting of Comsol® software. Concerning the boundary conditions, no mass flows across boundaries is allowed.

### 3. Results

#### 3.1. Influence of the Contact Angle

Once the yield stress is overcome, the motion of the cement paste through the granular matrix is governed by two properties: a term due to gravity and a term due to the surface tension as well as the contact angle. This part aims to define the influence of the contact angle on the penetration of the cement paste through the sand particle beds of different average diameter. Snapshots of the penetration of the fluid having a yield stress of 2 Pa illustrate, in Figure 6, the temporal development of the cement paste distribution through the granular beds depending of the value of the contact angle for a value of $d_{agg} = 1.6$ mm. These snapshots are from a revolution of the 2D axisymmetric numerical modeling.

The influence of the values of the contact angle is shown in more details as a function of time in Figures 7 and 8 plotting the penetration depth of respectively pastes of a value of yield stresses of 2 Pa and 4 Pa through a granular packing of $d_{agg} = 1.6$ mm. A value of $d_{agg} = 1.6$ mm of particles was chosen to study the influence of the contact angle as the ratio $\sqrt{\sigma/\rho g}$ is approximatively equal to 1.6 mm.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Snapshots of the penetration of a cement paste through a sand bed with particles of $d_{agg} = 1.6$ mm: influence of the value of the contact angle.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Influence of the contact angle at the interface for the penetration of a yield stress fluid ($\tau_{0,exp} = 2$ Pa; $\eta_{exp} = 0.1$ Pa·s) in a sand bed particles of 20 mm height and $d_{agg} = 1.6$ mm.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Influence of the contact angle at the interface for the penetration of a yield stress fluid ($\tau_{0,exp} = 4$ Pa; $\eta_{exp} = 0.1$ Pa·s) in a sand bed particles of 20 mm height and $d_{agg} = 1.6$ mm.</figcaption></figure>

It is observed from Figures 7 and 8 that the penetration depth kinetic depends on the contact angle. A small contact angle tends to slow down the penetration depth as function of the time for both values of yield stress. On the contrary, the liquid falls to the bottom of the granular bed in a shorter time for a value of contact angle of $5\pi/8$. In addition to the value of yield stress, the value of the contact angle at the interface between fluid and aggregates is a key parameter for the modeling of the 3D printing selective binding technique.

These results can be linked as the capillary force $F = \sigma \cos(\theta)$ vanishes if the contact angle is equal to $\pi/2$; is negative if the contact angle is $3\pi/8$ and is positive for a contact angle value of $5\pi/8$. Then, the contact angle with a value of $3\pi/8$ tends to pull upward the paste and on the contrary a value of contact angle of $5\pi/8$ tends to affect the motion by pulling downward the paste. If the value of contact angle is approximatively equal to $\pi/2$, the surface tension less affect motion as the force $du$ to surface tension and gravity are balanced and the penetration is only dependent of the rheological parameters of the fluid (yield stress and viscosity). Nevertheless, as the porous media does present surfaces which are not parallel, the resulting interface could be curved and hence result in a net pressure drop. Nevertheless, as we want to compare the numerical modeling results with previous existing results taking into account the rheological properties as well as the morphology and granular distribution of the granular bed, we choose a value of the contact angle of $\pi/2$. Within this angle value, we avoid the fluid being attracted to the grains as proposed by Boyce et al. [28].

## 3.2. Penetration Depth

Figures 9 and 10 show the penetration depth as a function of time for medium diameters of the sand bed particles and the cement pastes presenting yield stresses of 2 Pa and 4 Pa respective and viscosities values of 0.1 Pa s Herein, no scaling is applied for the model as these values refer to previous experimental measurements of Pierre et al. [27]. Results from Figure 9 show that the penetration of the fluid of $\tau_{0,exp} = 2$ Pa is complete for medium diameters of sand bed particles of 1.6 mm and 2.6 mm. Here, the numerical simulation for the average grain diameter of 2.6 mm did not last 20 s as we introduce a stop criterion when the mass conservation loss becomes higher than 5%. The kinetic of the penetration is slower for the fluid penetrating into the sand bed of $d_{agg} = 1.6$ mm than the sand bed of $d_{agg} = 2.6$ mm. The entire penetration through the sand bed of average diameter
of 1 mm is not achieved as it tends to a plateau value of 6 mm. As already shown using a level-set method [31], these results confirm that, if the contact angle is neutral, a balance between the average grain diameter of the sand particles and the rheological properties of the cement paste is crucial to achieve a complete penetration using a 3D printing particle based method process, as previously shown with the establishment of modified Darcy's law for yield stress [36].

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Influence of the average diameter of the sand bed particles on the penetration depth of a cement paste ($\tau_{0,exp} = 2$ Pa; $\eta_{exp} = 0.1$ Pa·s) of 20 mm height as a function of time.</figcaption></figure>

Figure 10 illustrates the penetration depth of the suspension having a yield stress of 4 Pa through the three different sand particles bed. We observe that plateau values occur for the sand bed particles of average grain diameters of 1 and 1.6 mm. The penetration is entirely achieved through a sand bed particle of average grain diameter of 2.6 mm.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Influence of the average diameter of the sand bed particles on the penetration depth of a cement paste ($\tau_{0,exp} = 4$ Pa; $\eta_{exp} = 0.1$ Pa·s) of 20 mm height as a function of time.</figcaption></figure>

### 3.3. *Comparison of Experimental Results and Numerical Simulations*

Figures 11 and 12 compare the penetration depths through particle-beds of different medium diameters obtained by numerical simulations with the analytical model and experimental measurements from Pierre et al. [27]. We plot Figures 11 and 12 the results of these experiments which were carried out
on both dry and wetted sand with a water content of 0.6% in order to assess a qualitative comparison of the 2D fluid invasion with two different cases of real 3D experiments. We should recall here that dry sand can lead to shear-thickening of the penetrating fluid and then lower penetration values than saturated wetted sand because of water absorption and surface wetting.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison of the numerical penetration depth of cement paste $\tau_{0,exp} = 2$ Pa; $\eta_{exp} = 0.1$ Pa·s) and the experimental results and analytical model from Pierre et al. [27].</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Comparison of the numerical penetration depth of cement paste ($\tau_{0,exp} = 4$ Pa; $\eta_{exp} = 0.1$ Pa·s) and the experimental results and analytical model from Pierre et al. [27].</figcaption></figure>
The analytical model of Pierre et al. [27] allows for evaluating the height of penetration in function of paste rheological parameters and density, particles bed average grain diameter and void ratio, solid volume fraction, following Equation (9)

$$h_{pen} = \frac{\rho g d_{agg} (1 - \Phi_s) h_{lay}}{\alpha \tau_{0,exp} - \rho g d_{agg} \Phi_s + \frac{\Phi_s}{1 - \Phi_s} 6 \tau_{0,exp} \kappa} \quad (9)$$

where $h_{lay}$ is the height of the layer of the sand packing bed, $\Phi_s$ is the sand volume fraction, $\tau_{0,exp}$ is the yield stress, $\kappa$ is a parameter that represents the fraction of the sand particles surface area where the fluid is sheared. The coefficient $\alpha$ is a fitting parameter computed from experiments on Carbopol gels and is independent on the fluid properties. Chevalier et al. [36] propose $\alpha = 5.5$ for a spherical particle assembly. The reader can refer to the work of Pierre et al. [27] for more details on Equation (9).

Figures 11 and 12 aim to compare the results from three different approaches: the 2D axi-symmetric numerical modeling presented in this study, the 1D analytical approach of Equation (9) and the 3D experimental results presented in the work of Pierre et al. [27]. Nevertheless, it is observed from both Figures 11 and 12 that the numerical modeling results show a quantitative evolution which follows the experimental measurements carried out on wet sand.

The penetration evaluated from Equation (9) shows higher penetration depth for a fluid with a value of yield stress of 2 Pa whereas a better agreement with the analytical model is found with a value of yield stress of 4 Pa. As a reminder, the experimental results with wet sand showed higher penetration than the dry sand as the water saturation of wet sand the particles prevent the absorption of the cement pastes. At a layer particle bed having an average grain diameter of 1 mm, a gap is observed between all the approaches. This comparison clearly shows that some interactions at the fluid-solid interface are not captured by both analytical 1D approach and developed 2D axi-symmetric numerical modeling for fine grains where interface between fluid and solid may govern the flow path. The gap between the experimental results and numerical simulations always ranges from one size of grains to two sizes of grains. For the sand packing with an average grain diameter of 1.6 mm, the difference of penetration depth values is about 3 mm for yield stress suspension of 4 Pa, which corresponds to two grains. Of course, the value of the penetration can be tailored with the values of the contact angles which govern the flow behaviour of the suspension. It then allows improving the penetration during the printing process through the use of surface tension modifiers that will improve the wettability. It can be noted that numerical modeling can capture the effect on heterogeneous granular packing that may create preferential flow path that cannot be predicted using the one-dimensional analytical modeling of Pierre et al. [27].

## 4. Conclusions

In this paper, we have discussed the influence of the average diameter of the sand particles, the rheological properties of fluids, and the contact angle on the influence of the penetration of yield stress fluids through a sand particles layer using numerical simulations. We emphasized that the contact angle can govern the kinetic of the penetration and the final depth penetration as reported from previous experimental works. Regarding the numerical modeling, we have shown that using a Fourier-Voronoi based generation to model a sand particle-bed and a phase-field method allows following the penetration depth of yield stress fluid through a sand particle bed. We highlighted that the contact angle is predominant in terms of the fluid behaviour traveling down the packing bed. Different scenarios can occur regarding the effect of surface tension, which is dependent on the contact angle: the penetration is accelerated, trapped, or does not infiltrate the sand bed. Also, the numerical simulations results were compared with analytical modeling, taking into account the sand particle bed properties and the rheology of the fluid and experimental results. The numerical results have provided a reliable quantitative prediction of the final depth penetration and more information about the flow path in comparison to previous analytical unidirectional modeling of the flow.
Advances in calculation clusters are expected to improve three-dimensional digital modeling by taking into account the surface tension effects and the modification of the contact angle by additives. Analytical modeling must also make advances in predicting the penetration of fluids into particles, whether through the activation of cement or the intrusion of paste.

**Author Contributions:** Conceptualization, A.P. (Alexandre Pierre), A.P. (Arnaud Perrot); methodology, A.P. (Alexandre Pierre); validation, D.W., A.P. (Alexandre Pierre), A.P. (Arnaud Perrot), and D.L.; formal analysis, A.P. (Alexandre Pierre); investigation, A.P. (Alexandre Pierre); resources, A.P. (Alexandre Pierre), A.P. (Arnaud Perrot); writing—original draft preparation, A.P. (Alexandre Pierre), A.P. (Arnaud Perrot); writing—review and editing, D.W., A.P. (Alexandre Pierre), A.P. (Arnaud Perrot) and D.L.; visualization, A.P. (Alexandre Pierre); supervision, A.P. (Alexandre Pierre), A.P. (Arnaud Perrot) and D.L.; project administration, A.P. (Alexandre Pierre), A.P. (Arnaud Perrot) and D.L.; funding acquisition D.W., A.P. (Alexandre Pierre), A.P. (Arnaud Perrot) and D.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors gratefully acknowledge the funding by the Bavarian-French University Cooperation Centre (BFHZ-CCUFB).

**Acknowledgments:** The authors gratefully acknowledge the funding by the Bavarian-French University Cooperation Centre (BFHZ-CCUFB) as well as the German Research Foundation (DFG).

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

## References

1. De Schutter, G.; Lesage, K.; Mechtcherine, V.; Nerella, V.N.; Habert, G.; Agusti-Juan, I. Technical potential of 3D printing, economic and environmental issues. *Cem. Concr. Res.* **2018**, *112*, 25–36. [CrossRef]

2. Lowke, D.; Dini, E.; Perrot, A.; Weger, D.; Gehlen, C.; Dillenburger, B. Particle-bed 3D printing in concrete construction—Possibilities and challenges. *Cem. Concr. Res.* **2018**, *112*, 50–65, ISSN 0008-8846. [CrossRef]

3. Wangler, T.; Lloret, E.; Reiter, L.; Hack, N.; Gramazio, F.; Kohler, M.; Bernhard, M.; Dillenburger, B.; Buchli, J.; Roussel, N.; et al. Digital Concrete: Opportunities and Challenges. *RILEM Tech. Lett.* **2016**, *1*, 67–75. [CrossRef]

4. Perrot, A. 3D *Printing of Concrete: State of the Art and Challenge of the Digital Construction Revolution*; Wiley: Hoboken, NJ, USA, 2019.

5. Khoshnevis, B.; Hwang, D.; Yao, K.T.; Yeh, Z. Mega-scale fabrication by Contour Crafting. *Int. J. Ind. Syst. Eng.* **2006**, *1*, 301–320. [CrossRef]

6. Zhang, J.; Khoshnevis, B. Optimal machine operation planning for construction by Contour Crafting. *Autom. Constr.* **2013**, *29*, 50–67. [CrossRef]

7. Le, T.T.; Austin, S.A.; Lim, S.; Buswell, R.; Gibb, A.G.F.; Thorpe, T. Mix design and fresh properties for high-performance printing concrete. *Mater. Struct.* **2012**, *45*, 1221–1232. [CrossRef]

8. Lim, S.; Buswell, R.A.; Le, T.T.; Austin, S.A.; Gibb, A.G.; Thorpe, T. Developments in construction-scale additive manufacturing processes. *Autom. Constr.* **2012**, *21*, 262–268. [CrossRef]

9. Perrot, A.; Rangeard, D.; Pierre, A. Structural built-up of cement-based materials used for 3D-printing extrusion techniques. *Mater. Struct.* **2015**, *49*, 1213–1220. [CrossRef]

10. Lloret, E.; Shahab, A.R.; Linus, M.; Flatt, R.J.; Gramazio, F.; Kohler, M.; Langenberg, S. Complex concrete structures. *Comput. Des.* **2015**, *60*, 40–49. [CrossRef]

11. Roussel, N. Rheological requirements for printable concretes. *Cem. Concr. Res.* **2018**, *112*, 76–85. [CrossRef]

12. Cesaretti, G.; Dini, E.; De Kestelier, X.; Colla, V.; Pambaguian, L. Building components for an outpost on the Lunar soil by means of a novel 3D printing technology. *Acta Astronaut.* **2014**, *93*, 430–450. [CrossRef]

13. Fromm, A. 3-D-*Printing Zementgebundender Formteile: Grundlagen, Entwicklung und Verwendung*; Kassel University Press GmbH: Kassel, Germany, 2014.

14. Lowke, D.; Weger, D.; Henke, K.; Talke, D.; Winter, S.; Gehlen, C. 3D-Drucken von Betonbauteilen durch selektives Binden mit calciumsilikatbasierten Zementen: Erste Ergebnisse zu betontechnologischen und verfahrenstechnischen Einflüssen. In *Tagungsbericht '19. Internationale Baustofftagung*, 700th ed.; Ludwig, P.D.-I.H.-M., Ed.; Weimar Pub.: Weimar, Germany, 2015.

15. Shakor, P.; Sanjayan, J.; Nazari, A.; Nejadi, S. Modified 3D printed powder to cement-based material and mechanical properties of cement scaffold used in 3D printing. *Constr. Build. Mater.* **2017**, *138*, 398–409. [CrossRef]
16. Weger, D.; Lowke, D.; Gehlen, C. 3D printing of concrete structures using the selective binding method—Effect of concrete technology on contour precision and compressive strength. In Proceedings of the 11th Fib International PhD Symposium in Civil Engineering, Tokyo, Japan, 29–31 August 2016; Maekawa, K., Kasuga, A., Yamazaki, J., Eds.; The University of Tokyo: Tokyo, Japan, 2016; pp. 403–410, ISBN 978-4-9909148-0-6.

17. Weger, D.; Lowke, D.; Gehlen, C. 3D Printing of Concrete Structures with Calcium Silicate based Cements using the Selective Binding Method—Effects of Concrete Technology on Penetration Depth of Cement Paste. In Proceedings of the Hipermat 2016—4th International Symposium on Ultra-High Performance Concrete and High Performance Construction Materials Kassel, Kassel, Germany, 9–11 March 2016; Kassel University Press: Kassel, Germany, 2016.

18. Weger, D.; Lowke, D.; Gehlen, C. Additive manufacturing of concrete elements using selective cement paste intrusion—Effect of layer orientation on strength and durability. In Proceedings of the RILEM 1st International Conference on Concrete and Digital Fabrication, Zurich, Switzerland, 15 January 2018.

19. Vasilic, K.; Gram, A.; Wallevik, J.E. Numerical simulation of fresh concrete flow: Insight and challenges. *RILEM Tech. Lett.* **2019**, *4*, 57–66. [CrossRef]

20. Cremonesi, M.; Ferrara, L.; Frangi, A.; Perego, U. Simulation of the flow of fresh cement suspensions by a Lagrangian finite element approach. *J. Non-Newtonian Fluid Mech.* **2010**, *165*, 1555–1563. [CrossRef]

21. Mokbel, D.; Abels, H.; Aland, S. A phase-field model for fluid-structure interaction. *J. Comput. Phys.* **2018**, *372*, 823–840. [CrossRef]

22. Osher, S.; Fedkiw, R. *Level Set Methods and Dynamic Implicit Surfaces*; Springer Science and Business Media LLC: Berlin/Heidelberg, Germany, 2003.

23. Hirt, C.; Nichols, B. Volume of fluid (VOF) method for the dynamics of free boundaries. *J. Comput. Phys.* **1981**, *39*, 201–225. [CrossRef]

24. Aland, S.; Voigt, A. Benchmark computations of diffuse interface models for two-dimensional bubble dynamics. *Int. J. Numer. Methods Fluids* **2012**, *69*, 747–761. [CrossRef]

25. Ge, Z.; Loiseau, J.-C.; Tammisola, O.; Brandt, L. An efficient mass-preserving interface-correction levelset/ghost fluid method for droplet suspensions under depletion forces. *J. Comput. Phys.* **2018**, *353*, 435–459. [CrossRef]

26. Thibaud Chevalier. *Ecoulements de Fluides à Seuil en Milieux Confinés. Mécanique des Fluides* [*Physics.class-ph*]; Université Paris-Est: Champs-sur-Marne, France, 2013.

27. Pierre, A.; Weger, D.; Perrot, A.; Lowke, D. Penetration of cement pastes into sand packings during 3D printing: Analytical and experimental study. *Mater. Struct.* **2018**, *51*, 22. [CrossRef]

28. Boyce, C.M.; Ozel, A.; Sundaresan, S. Intrusion of a Liquid Droplet into a Powder under Gravity. *Langmuir* **2016**, *32*, 8631–8640. [CrossRef]

29. Mollon, G.; Zhao, J. Fourier–Voronoi-based generation of realistic samples for discrete modelling of granular materials. *Granul. Matter* **2012**, *14*, 621–638. [CrossRef]

30. Papanastasiou, T.C. Flows of Materials with Yield. *J. Rheol.* **1987**, *31*, 385–404. [CrossRef]

31. Pierre, A.; Weger, D.; Perrot, A.; Lowke, D. 2D Numerical Modelling of Particle-Bed 3D Printing by Selective Paste Intrusion. In *Proceedings of the RILEM Bookseries*; Springer Science and Business Media LLC: Berlin/Heidelberg, Germany, 2020; Volume 28, pp. 342–352.

32. Mitsoulis, E. Rheology reviews. British Society of Rheology: Aberystwyth, UK, 2007; pp. 135–178.

33. Feys, D.; Wallevik, J.E.; Yahia, A.; Khayat, K.H.; Wallevik, O.H. Extension of the Reiner–Riwlin equation to determine modified Bingham parameters measured in coaxial cylinders rheometers. *Mater. Struct.* **2012**, *46*, 289–311. [CrossRef]

34. Ferraris, C.F.; Brower, L.E.; Banfill, P.; Beaupré, D.; Chapdelaine, F.; De Larrard, F.; Domone, P.; Nachbaur, L.; Sedran, T.; Wallevik, O.; et al. *Comparison of Concrete Rheometers*; National Institute of Standards and Technology (NIST): Gaithersburg, MA, USA, 2001.
35. Bai, F.; He, X.; Yang, X.; Zhou, R.; Wang, C. Three dimensional phase-field investigation of droplet formation in microfluidic flow focusing devices with experimental validation. *Int. J. Multiph. Flow* **2017**, *93*, 130–141.
[CrossRef]

36. Chevalier, T.; Clain, X.; Dupla, J.; Canou, J.; Rodts, S.; Coussot, P. Darcy's law for yield stress fluid flowing through a porous medium. *J. Non-Newtonian Fluid Mech.* **2013**, *195*, 57–66. [CrossRef]

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_14.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).