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

Research paper

# Multiphysics modelling of 3D concrete printing: From material model to process simulation and optimisation

Maxime Pierre<sup>a, ID,*</sup>, Siavash Ghabezloo<sup>a, ID</sup>, Patrick Dangla<sup>a, ID</sup>, Romain Mesnil<sup>a,b, ID</sup>, Matthieu Vandamme<sup>a, ID</sup>, Jean-François Caron<sup>a, ID</sup>

## ARTICLE INFO

### Keywords:

- 3D concrete printing
- Numerical simulation
- Multiphysics couplings
- Process optimisation

## ABSTRACT

Predictive simulation of 3D concrete printing is important to warrant printability and durability of print pieces and to optimise printing parameters, yet tedious due to the complexity of the material behaviour and printing process. From a constitutive model allowing a continuous description of the coupled thermo-poro-mechanical behaviour of cement-based materials from the early-age to the hardened state, a comprehensive finite element simulation framework is designed. It aims at modelling extrusion-based 3D printing processes, taking into account the sequential deposition of material. Study of the onset of plastic collapse on specific geometries at different printing speeds show the complexity of collapse prediction as well as the importance of process-related effects. An optimisation scheme is proposed to determine optimal printing speed modulations from numerical simulations with the perspective of increasing productivity in 3D concrete printing. The model shows good predicting capabilities when compared with experimental printing failures, and is able to extrapolate to other accelerator dosages without model re-calibration.

## 1. Introduction

The construction industry is currently faced with several challenges. It has been shown that its productivity has stagnated, if not decreased, in the past 50 years [1], due to building techniques that have rarely evolved for the majority of a century. Growing concerns about environmental issues and energy consumption also put pressure on the buildings and construction industry responsible for 36% of the final use of energy and 39% of carbon dioxide emissions related to energy and processes [2]. Cement production alone accounts for roughly 8% of global CO₂ emissions (as of 2012 [3]), hence the necessity to reduce our cement consumption to meet emission goals. To face these challenges, industrialisation of the construction sector is imperative, and 3D printing with cement-based materials, especially through extrusion technologies, has emerged as one of the promising off-site processes. Now a fast-growing field in both the industry and academia, it yields hope for increasing productivity, along with a reduced impact from optimised structures, as opposed to traditional, “monolithic” casting methods.

There are still several issues that need to be addressed with regard to extrusion-based 3D printing of cementitious materials before this technology becomes widespread in the industry. The first crucial challenge is to ensure the printability, or buildability, of a given print piece.

This requires a good understanding of how the strength and stiffness of the material develop throughout the hydration process during printing, hence the need for a material model based on hydration and its induced evolution of material properties. Indeed, many external factors can influence the hydration kinetics, in particular temperature and relative humidity, rendering a model relying on time-dependent rather than hydration-dependent properties of limited generality as it would be tied to specific printing conditions.

The perspective of printing in uncontrolled environments also raises the question of the multiphysics aspect of material behaviour. Cementitious materials are porous materials, and as such, the effect of the pore fluid cannot be overlooked: this is even more true in that the water inside the pore is a critical component of the hydration reaction. Understanding the mechanisms behind shrinkage, for instance, be it from self-desiccation or external drying, requires the modelling of poromechanical couplings as well as fluid transport in the medium. Shrinkage can significantly alter the geometry of a 3D printed object, making it deviate from the intended geometry, which can be detrimental in many applications, especially in prefabrication where different elements may need to interlock. Thorough modelling of shrinkage phenomena could help counterbalance this effect by printing a slightly different geometry

<footnote><sup>a</sup> Navier, ENPC, Institut Polytechnique de Paris, Univ Gustave Eiffel, CNRS, Marne-la-Vallée, France</footnote>

<footnote><sup>b</sup> Build'in, Ecole des Ponts, Marne-la-Vallée, France</footnote>

<footnote>* Corresponding author.
E-mail address: maxime.pierre@enpc.fr (M. Pierre).</footnote>
from the intended final geometry. Drying can, in turn, have detrimental consequences on durability, through cracking or *shark skin* effect at the early age, or by reducing the extent of hydration in dry areas. The latter may be a factor in the deterioration of the interlayer bond strength known as *cold joint* [4,5], leading to preferential mechanical failure at the interface between layers [6]. It also implies that material properties are likely to be heterogeneous in the lace section, as the air-exposed periphery will tend to be dryer and thus potentially less hydrated than the centre of the section, which should be taken into account when estimating the mechanical properties of a printed object [7]. Furthermore, the long-term durability of 3D printed cement-based objects with respect to chloride ingress or freeze-thaw damage, for instance, is also affected by heterogeneities in the material's structure which stem from the printing process [8–10]. Comprehensive modelling of cement-based 3D-printed structures can thus only be achieved using a multiphysics material model taking into account the different couplings between mechanics, pore fluid, thermal and chemical processes.

Finally, 3D printing is meant to produce complex geometries otherwise costly if not impossible to realise by conventional casting methods. This also entails that analytical methods are not suitable in a majority of real-life cases, and that the development of simulation tools, based, for instance, on the finite element method (FEM), is required.

In the present paper, we propose a simulation framework for extrusion-based 3D printing processes with cementitious materials, backed by a multiphysics material model, which describes cementitious materials as porous media with hydration-dependent properties. The model formulation is based on previous work performed in our laboratory on the behaviour of oil-well cements. Quite similar to 3D concrete printing, the simulation of the response of a cement sheath in a deep well needs a constitutive law for the cement paste from the early age to the hardened state, taking into account the couplings between hydraulic, thermal, and mechanical phenomena. The thermoporomechanical behaviour of the hardened and hydrating cement paste was studied in a vast experimental programme [11–16]. The results were then associated with detailed microstructure characterisation [17,18] and micromechanical modelling [19–22] to develop a modelling framework for cement paste from early age to hardened state [23–26]. This model is enriched to take into account thermal effects as well as the possibility of expansion during hydration, and adapted to 3D concrete printing materials. The material model is calibrated using different experiments, and its ability to predict the evolution of materials with different accelerator contents using a single set of parameters is shown. It is then implemented in a finite element code and used to perform 2D axisymmetric and 3D simulations of the concrete printing process. After analysing the mechanisms of plastic collapse on specific geometries, we propose to leverage such a simulation framework towards process optimisation, in particular with the modulation of printing speed during printing to achieve maximum productivity on a given print-piece. Finally, a quantitative comparison with 3 different experimental cases is conducted as validation for the complete simulation framework.

List of symbols

<table><tr><td>σ</td><td>Stress tensor</td><td>Pa</td></tr><tr><td>ε</td><td>Strain tensor</td><td></td></tr><tr><td>T</td><td>Temperature</td><td>K</td></tr><tr><td>S</td><td>Lagrangian density of entropy</td><td>J K<sup>-1</sup> m<sup>-3</sup></td></tr><tr><td>M<sub>lw</sub></td><td>Lagrangian fluid mass flow vector</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>ρ<sub>lw</sub></td><td>Water density</td><td>kg m<sup>3</sup></td></tr><tr><td>p<sub>l</sub></td><td>Liquid pressure</td><td>Pa</td></tr><tr><td>S<sub>s</sub></td><td>Lagrangian density of solid entropy</td><td>J K<sup>-1</sup> m<sup>3</sup></td></tr><tr><td>ν<sub>lw</sub></td><td>Water stoichiometric coefficient</td><td></td></tr></table>

<table><tr><td>M<sub>lw</sub></td><td>Water molar mass</td><td>kg mol<sup>-1</sup></td></tr><tr><td>ξ</td><td>Lagrangian density of extent of reaction</td><td>mol m<sup>3</sup></td></tr><tr><td>φ</td><td>Porosity</td><td></td></tr><tr><td>φ<sub>0</sub></td><td>Initial porosity</td><td></td></tr><tr><td>φ</td><td>Mechanical contribution to porosity variation (φ - φ<sub>0</sub>)</td><td></td></tr><tr><td>S<sub>l</sub></td><td>Lagrangian liquid saturation degree</td><td></td></tr><tr><td>ν<sub>s</sub></td><td>Solid volume fraction</td><td></td></tr><tr><td>p<sub>c</sub></td><td>Capillary pressure</td><td>Pa</td></tr><tr><td>p<sup>*</sup></td><td>Equivalent pore pressure</td><td>Pa</td></tr><tr><td>Θ</td><td>Stiffness tensor</td><td>Pa</td></tr><tr><td>b</td><td>Biot's effective stress coefficient</td><td></td></tr><tr><td>α<sub>s</sub></td><td>Solid thermal expansion</td><td>K<sup>-1</sup></td></tr><tr><td>L</td><td>Latent heat of hydration</td><td>J mol<sup>-1</sup></td></tr><tr><td>C<sub>αs</sub></td><td>Volumetric heat capacity of the solid at constant stress</td><td>J K<sup>-1</sup> m<sup>-3</sup></td></tr><tr><td>A</td><td>Chemical affinity</td><td>J mol<sup>-1</sup></td></tr><tr><td>ε<sub>∞</sub></td><td>Theoretical maximum lagrangian density of hydration extent</td><td>mol m<sup>-3</sup></td></tr><tr><td>α</td><td>Hydration degree (ε/ε<sub>∞</sub>)</td><td></td></tr><tr><td>R<sub>H</sub></td><td>Relative humidity</td><td></td></tr><tr><td>σ'</td><td>Mean effective stress</td><td>Pa</td></tr><tr><td>τ</td><td>Deviatoric stress</td><td>Pa</td></tr><tr><td>σ<sub>c</sub></td><td>Yield stress in compression</td><td>Pa</td></tr><tr><td>σ<sub>f</sub></td><td>Yield stress in traction</td><td>Pa</td></tr><tr><td>k<sub>l</sub></td><td>Liquid permeability</td><td>m<sup>2</sup></td></tr><tr><td>η<sub>lw</sub></td><td>Water dynamic viscosity</td><td>Pa s</td></tr><tr><td>k<sub>T</sub></td><td>Thermal conductivity</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>E</td><td>Young's modulus of cement paste</td><td>Pa</td></tr><tr><td>K</td><td>Bulk modulus of cement paste</td><td>Pa</td></tr><tr><td>E<sub>ε</sub></td><td>Chemical stiffening modulus</td><td>Pa</td></tr><tr><td>h<sub>ε</sub></td><td>Chemical hardening modulus</td><td>Pa</td></tr><tr><td>h<sub>c</sub></td><td>Isotropic hardening coefficient</td><td>Pa</td></tr></table>

## 2. Thermoporomechanical material model

The present model constitutes an extension of the work of Coussy on porous media [27,28] to cement-based materials, taking into account the hydration process and the associated evolution of material properties, developed by Samudio [24]. A detailed derivation of the following constitutive equations can be found in [26]. This model takes classical equations of partially saturated thermoporomechanics and enriches them by adding evolution laws with regard to the progressive advancement of the chemical reaction of cement hydration. It can therefore be used to model the behaviour of cement-based materials from their very early-age, right after mixing, all the way to their hardened state, while taking into account heat and fluid transport in the medium and full thermo-hydro-mechanical couplings. It features elasto-plasticity with a chemically hardening yield surface, an hydration-dependent retention-curve, and also some newer elements specific to accelerated mortars for 2K printing systems, such as early age expansion due to the growth of ettringite crystals. Compared with material models traditionally used for modelling 3D concrete printing, which are mostly limited to mechanics and assume a time-based evolution of the material properties, this model provides more generalisability by considering the underlying chemical reactions responsible for cement hardening, and also incorporates couplings with other physics which are of high importance at the early age, such as hydromechanical couplings which enable to understand self-desiccation and drying shrinkage of cement-based materials.

### 2.1. Main constitutive equations

The following equations relate the evolution of the stress $\sigma$, porosity variation $\phi = \phi - \phi_0$ and solid entropy $S_s$ with the variations of the
variables: elastic strain $\xi^{el}$, liquid pressure $p_l$, temperature $T$ and extent of the hydration reaction $\xi$:

$$ 
\begin{aligned}
d\sigma|_{S_l} &= \mathcal{C}(\xi) : d\xi^{el} - S_l \cdot b(\xi) \frac{1}{\xi} d\rho_l - 3\kappa(\xi) \frac{1}{\xi} dT, \\
d\varphi|_{S_l} &= b(\xi) \frac{1}{\xi} : d\xi^{el} + \frac{b(\xi) - \phi_0(\xi)}{K_s} S_l d\rho_l - 3\alpha_s(b(\xi) - \phi_0(\xi))dT. \\
dS_s|_{S_l} &= 3\kappa(\xi) \frac{1}{\xi} : d\xi^{el} - 3\alpha_s(b(\xi) - \phi_0(\xi))S_l d\rho_l + \frac{(1 - \phi_0(\xi))C_{\sigma_s}}{T} dT - \frac{L}{T} d\xi.
\end{aligned}
\quad (1) $$

In Eq. (1), $\mathcal{C}$ designates the 4th-order elasticity tensor, $S_l$ the liquid saturation degree, b Biot's effective stress coefficient, $\phi_0$ the initial porosity, $K_s$ the bulk modulus of the solid, $\alpha_s$ the thermal expansion of the solid, $\kappa = \alpha_s K$ where K is the solid skeleton's bulk modulus, $C_{\sigma_s}$ the heat capacity of the solid under constant stress, and L the latent heat of hydration. The tangent thermoporomechanical coefficients all depend on the extent of the hydration reaction $\xi$, with a dependency that is supposed to be linear. This assumption is supported by the results of micromechanical modelling on hydrating cement paste, where the thermoporomechanical coefficients follow a near-linear trend with hydration after percolation [23,26]. For the materials of interest in 3D printing, we make the hypothesis that percolation has already occurred when the material exits the nozzle. Indeed, as the material can sustain its own weight, a solid skeleton is therefore present. Introducing the change in solid molar volume $\Delta\overline{V}_s$ in the hydration reaction, the increase in solid volume fraction $d\nu_s$ associated to an increment of extent of reaction $d\xi$ is given by:

$$ d\nu_s = \Delta\overline{V}_s d\xi, \quad (2) $$

such that the evolution of the initial porosity $\phi_0$ is:

$$ d\phi_0 = -d\nu_s = -\Delta\overline{V}_s d\xi. \quad (3) $$

Defining a chemical stiffening modulus $E_\xi$, the evolution of stiffness is as follows:

$$ dE = E_\xi \Delta\overline{V}_s d\xi. \quad (4) $$

Biot's effective stress coefficient $b = 1 - K/K_s$ will hereafter be considered constant equal to 1, as this study aims at very early-age materials for which $K \ll K_s$.

## 2.2. Inelastic strains

In addition to elastic strain, we consider plastic strain as well as a so-called expansive strain:

$$ \underline{\underline{\epsilon}} = \underline{\underline{\epsilon}}^{el} + \underline{\underline{\epsilon}}^{p} + \underline{\underline{\epsilon}}^{exp}. \quad (5) $$

Indeed, accelerated mortars using shotcrete accelerators, which promote the precipitation of ettringite at the very early-age, can undergo significant expansion [29]. We will consider it isotropic and linear with respect to the hydration extent $\xi$, so that if $\epsilon_V^\xi$ is the volumetric expansion associated with an hydration extent increment of 1 mol:

$$ \underline{\underline{d\epsilon}}^{exp} = \frac{1}{3} \underline{\underline{\epsilon}}_V^\xi \frac{1}{\xi} d\xi. \quad (6) $$

Given the high stiffness of the solid compared to that of the porous skeleton, this expansive deformation is due in its integrality to a deformation of the pore space:

$$ d\phi^{exp} = \text{tr}(d\underline{\underline{\epsilon}}^{exp}) = \underline{\underline{\epsilon}}_V^\xi d\xi. \quad (7) $$

The plastic behaviour of cement-based materials, at any age, is known to be of importance and can in some cases cause object collapse during printing [30]. The Asymmetric Cam-Clay (ACC) model [24] adopted here proposes a continuous description of the plastic behaviour

<figure><img src="image_2.png" /><figcaption>Fig. 1. Evolution of elastic domain with hydration.</figcaption></figure>

of cement-based materials through chemical hardening laws. The yield function is given by:

$$ f(\sigma', \tau) = \tau^2 \exp \left( k \frac{2\sigma' + \sigma_t - \sigma_c}{\sigma_t + \sigma_c} \right) + M^2(\sigma' + \sigma_t)(\sigma' - \sigma_c), \quad (8) $$

where $\sigma' = \text{tr}(\underline{\underline{\sigma}} + \beta p^*)$ is an effective stress measure with $p^* = S_l p_l$ as the equivalent pore pressure, $\tau$ the Von Mises norm of the stress's deviator, $k$ is a material parameter and $\sigma_c$, $\sigma_t$ and $M$ are variables evolving with the hydration extent $\xi$ and the plastic strain $\underline{\underline{\epsilon}}^p$. Their dependencies are as follows:

$$ \frac{\partial \sigma_c}{\partial \xi} = -h_{\underline{\underline{\xi}}} \frac{d\phi_0}{d\underline{\underline{\xi}}} = h_{\underline{\underline{\xi}}} \Delta \overline{V}_s, \quad (9) $$

$$ dM = M_{\underline{\underline{\xi}}} \Delta \overline{V}_s d\underline{\underline{\xi}}, \quad (10) $$

with $h_{\underline{\underline{\xi}}}$ the chemical hardening modulus and $M_{\underline{\underline{\xi}}}$ the chemical shear hardening coefficient. The corresponding evolution of the yield surface with hydration is depicted in Fig. 1. The evolution of $\sigma_c$ and $\sigma_t$ is complete with the isotropic hardening rule:

$$ \frac{\partial \sigma_c}{\partial \underline{\underline{\epsilon}}^p} = \sigma_c h_{\underline{\underline{\epsilon}}} (1 - \phi_0(\underline{\underline{\epsilon}})) \underline{\underline{\underline{1}}}, \quad (11) $$

$$ \sigma_t = f_t \sigma_c, \quad (12) $$

where $h_{\underline{\underline{\epsilon}}}$ is the isotropic hardening coefficient and $f_t$ the tensile strength ratio.

The plastic variation of porosity is given by:

$$ d\phi^p = \beta \text{tr}(d\underline{\underline{\epsilon}}^p), \quad (13) $$

which results in $\beta = 1$ under the assumption of a plastically incompressible matrix [27].

## 2.3. Complementary evolution laws

### 2.3.1. Hydration kinetics

The hydration kinetics is driven by the constitutive model as many parameters depend on it. We assume that the hydration rate is proportional to the chemical affinity and that the reaction is thermo-activated through an Arrhenius-style law:

$$ \frac{d\underline{\underline{\xi}}}{dt} = \frac{\mathcal{A}}{\eta_{\underline{\underline{\xi}}}} \exp \left( -\frac{E_{\mathcal{A}}}{RT} \right), \quad (14) $$

with $\mathcal{A}$ the chemical affinity, $E_{\mathcal{A}}$ the activation energy of the chemical reaction forming hydrates, and $\eta_{\underline{\underline{\xi}}}$ the viscosity due to the diffusion of water through the layer of already formed hydrates. The chemical affinity of the reaction also depends on the relative humidity in the pore space, as it has been reported both from a theoretical and experimental standpoint that hydration reactions stop under a certain relative humidity (typically around 80% for alite) [31–33]. Here, we will use the empirical formula proposed by Wyrzykowski et al. [33]. Combining it with the forms of affinity and viscosity from Tomas [34] yields:

$$ \mathcal{A}(\underline{\underline{\xi}}, R_H) = \min \left( 1, 3.7 (R_H)^{2.5} \right) k_{\underline{\underline{\xi}}} \left( \frac{\mathcal{A}_0}{k_{\underline{\underline{\xi}}} \underline{\underline{\xi}}_{\infty}} + \underline{\underline{\xi}} \right) (\underline{\underline{\xi}}_{\infty} - \underline{\underline{\xi}}), \quad (15) $$
$$ \eta_{\xi}(\xi) = \eta_{\xi,0} \exp\left(\bar{\eta} \frac{\xi}{\xi_{\infty}}\right), \qquad (16) $$

where $R_H$ is the relative humidity, $A_0$, $k_\xi$, $\eta_{\xi,0}$ and $\bar{\eta}$ are material parameters and $\xi_\infty$ is the hydration extent corresponding to a complete hydration.

### 2.3.2. Water retention curve

The retention curve relates the saturation degree $S_l$ to the capillary pressure $p_c = p_g - p_l \approx -p_l$ as we consider the gas pressure to be negligible. The Van Genuchten model is adopted [35], in a modified form [24] which depends on the hydration extent $\xi$:

$$ S_l = \left(1 + \left(\frac{p_c}{p_r(\xi)}\right)^n\right)^{\frac{1-n}{n}}, \qquad (17) $$

where $n > 1$ is a real exponent and $p_r$ a reference pressure dependent on $\xi$ through:

$$ p_r(\xi) = p_{r,0} \exp\left(k_{sl} \frac{\xi}{\xi_{\infty}}\right), \qquad (18) $$

where $k_{sl}$ is a real material parameter, describing the shift of the retention curve with hydration due to the change of the pore structure.

## 2.4. Field equations and transport properties

In the absence of body forces other than self-weight induced by gravity's acceleration $\gamma$, the mechanical equilibrium reads:

$$ \underline{\nabla} \cdot \underline{\sigma} = \rho \underline{\gamma}. \qquad (19) $$

The field equation for the fluid is obtained by writing the conservation of the fluid mass, and hence:

$$ \frac{dm_w}{dt} = -\underline{\nabla} \cdot \underline{M}_w + \nu_w M_w \frac{d\xi}{dt}, \qquad (20) $$

where $M_w$ is the water mass flow vector, $\nu_w$ is the stoichiometric coefficient of water in the hydration reaction and $M_w$ the molar mass of water. The water mass flow can be obtained from Darcy's law:

$$ \underline{M}_w = -\rho_w \frac{k_l(S_l, \phi)}{\eta_w} \underline{\nabla} \underline{P}_l, \qquad (21) $$

with $\eta_w$ the dynamic viscosity of water (≈ 1 mPa·s at 20 °C). The liquid permeability $k_l$ is not considered constant, but is a function of the saturation degree $S_l$ and the porosity $\phi$. Permeability decreases with decreasing saturation degree, and here we use the widely adopted model from Mualem [36]:

$$ \frac{k_l(S_l)}{k_l(S_l = 1)} = \sqrt{S_l} \left( 1 - \left( 1 - S_l^{\frac{n}{n-1}} \right)^{\frac{n-1}{n}} \right)^2, \qquad (22) $$

where exponent $n$ refers to the same exponent as the one from the water retention curve (17). The decrease of the permeability with decreasing porosity has also been the focus of several models, the most used being that of Kozeny and Carman [37], and that of Taylor [38]. Following the work of Picandet et al. [39] on cement pastes, we will adopt Taylor's model, linking the logarithm of the permeability to the porosity:

$$ \frac{\partial}{\partial \phi} (\ln k_l) = \frac{C_{k_l}}{(1 - \phi)^2} \qquad (23) $$

The thermal field equation is obtained by writing the heat conservation:

$$ T \left( \frac{\partial S}{\partial t} - \underline{\nabla} \cdot \left( s_w \underline{M}_w \right) \right) = - \underline{\nabla} \cdot \underline{q}. \qquad (24) $$

Combining this equation with Fourier's law:

$$ \underline{q} = -k_T \underline{\nabla} T, \qquad (25) $$

<table><caption>**Table 1**
Mix 2 proportions for 1 kg before addition of the accelerator.</caption><thead><tr><td>Material</td><td>Quantity</td></tr></thead><tbody><tr><td>Premix</td><td>908 g</td></tr><tr><td>Water</td><td>88 g</td></tr><tr><td>Superplasticizer</td><td>4 g</td></tr><tr><td>Accelerator</td><td>6 g</td></tr></tbody></table>

<table><caption>**Table 2**
Mix 2 proportions for 1 kg before addition of the accelerator.</caption><thead><tr><td>Material</td><td>Quantity</td></tr></thead><tbody><tr><td>Premix</td><td>896 g</td></tr><tr><td>Water</td><td>102 g</td></tr><tr><td>Superplasticizer</td><td>2 g</td></tr><tr><td>Accelerator</td><td>4.7 g or 3 g</td></tr></tbody></table>

gives the following form for the thermal field equation:

$$ T \left( \frac{\partial S}{\partial t} - \underline{\nabla} \cdot \left( s_w \underline{M}_w \right) \right) = k_T \Delta T. \qquad (26) $$

# 3. Model parameters for two component 3D-printing

Although the framework presented above is meant to be as general as possible and could be applicable to a wide range of cement-based materials, from here on we will focus on more specific mix designs suitable for 3D printing with a two-component strategy. This strategy relies on using a flowable mortar with relatively low yield stress (typically around 200 Pa [40]) which is easy to pump and mixed with an accelerator in the printer head, right before the nozzle. We focus thereafter on the fine determination of the most significant parameters when considering plastic collapse of 3D-printed objects, namely hydration kinetics, evolution of mechanical properties and partially saturated behaviour with the evolution of the retention curve.

## 3.1. Material and methods

We make use of different experiments in order to determine some of our model parameters. In particular, measuring the p-wave velocities during hydration gives us a proxy measurement of the hydration kinetics. Performing mercury intrusion porosimetry on samples that have been freeze-dried at different times following the introduction of the accelerating agent allows to follow the evolution of the total porosity and of the pore size distribution, and therefore the water retention curve, with hydration. For monitoring the evolution of the yield stress of the material, we use a pocket shear vane, as suggested by Demont et al. [41].

### 3.1.1. Mortar mixes

We use two mortar mixes which have been employed at the printing facility at Ecole nationale des ponts et chaussées (ENPC). Mix 1 is comprised of a premix (Ductal NaG3 from Lafarge-Holcim) mixed with water and a superplasticizer (CHRYSO®Fluid Optima 100), to which an alkali-free accelerator (Floquat ASL) is added. Mix proportions are given in Table 1. The mortar was mixed in a oil-well cement blender for 5 min prior to the addition of the accelerator, after which it was mixed for another 10 s.

Mix 2 is composed of a premix developed in-house (portland cement, sand and limestone filler), a superplasticizer mix and the same alkali-free accelerator as Mix 1. Table 2 gives the mass proportions.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Fitting of the hydration kinetics model on P-wave velocity data for Mix 1.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Fitting of the hydration kinetics model from P-wave velocity data for Mix 2 with 7 mL·L⁻¹ and 11 mL·L⁻¹ accelerator dosage.</figcaption></figure>

### 3.1.2. P-wave velocity measurements

Measurements of P-wave velocities were performed in the early age on a sample of 2 cm thickness and 8 cm diameter, up to 3 h after adding the accelerator to the fresh mortar. A paper mould was used to prevent the P-waves from travelling through the mould, which could occur in the case of a more rigid mould and alter the results. The velocity of the P-waves is obtained by measuring the travel time through the sample and two 2 mm thick PVC plates that isolate the wave sensors from the fresh mortar. The travel time of the P-waves through these PVC plates was calibrated prior to the experiment. The sensors operated at a frequency of 54 kHz and a tension of 500 V, and the acquisition was carried out with a precision of 0.5 µs, which translates into data quantification as seen in Fig. 2.

### 3.1.3. Mercury intrusion porosimetry (MIP)

Fresh mortar was deposited in hermetic containers immediately after mixing to avoid evaporation to occur and influence the hydration process. The containers were opened at different times to measure time evolution of the porosity, and several cubes of approximately 1 cm were cut from the mortar. Then they were freeze dried in liquid nitrogen under vacuum. Dried samples were used to perform the MIP

<figure><img src="image_5.png" /><figcaption>Fig. 4. Evolution of total porosity in the first 90 min following accelerator introduction (MIP measurements).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Water retention curves from MIP measurements at different mortar ages and dual Van Genuchten fits.</figcaption></figure>

measurements. The total porosity measured on two distinct mortar batches is shown in Fig. 4.

### 3.2. Hydration reaction

The accelerator is an alkali-free accelerator composed of aluminium sulphate, which is known to greatly improve cement setting time by stimulating the production of ettringites (AFt) (28) from dissolved alite (C₃S) (27) [42].

$$ \mathrm{Ca_3SiO_5} + 5\mathrm{H_2O} \rightarrow 3\mathrm{Ca^{2+}} + 6\mathrm{OH^{-}} + \mathrm{H_4SiO_4} \quad (27) \\ \mathrm{C_3S} \text{ (Alite)} $$

$$ 2[\mathrm{Al(OH)_4}^{-}] + 6\mathrm{Ca^{2+}} + 4\mathrm{OH^{-}} + 3\mathrm{SO_4^{2-}} + 26\mathrm{H_2O} \rightarrow 3\mathrm{CaO} \cdot \mathrm{Al_2O_3} \cdot 3\mathrm{CaSO_4} \cdot 32\mathrm{H_2O} \quad (28) \\ \text{ÅFt (Ettringite)} $$

Combining the two reactions, 2 moles of C₃S can produce, in the presence of aluminium sulphate, one mole of ettringite by consuming 36 moles of water in the process. The molar volumes are given in Table 3. The molar variation of the solid volume amounts to ΔV_s = V̄_AFt - 2V̄_C₃S = 574.2 cm³·mol⁻¹, while the net change in molar volume
<table><caption>**Table 3**
Molar volumes in ettringitic reaction.</caption><thead><tr><th>Species</th><th>Molar volume [cm<sup>3</sup>·mol<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Alite</td><td>75.4</td></tr><tr><td>Water</td><td>18.07</td></tr><tr><td>Ettringite</td><td>725</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Yield stress measurements using a pocket shear vane on Mix 2 with 7 mL·L⁻¹ and 11 mL·L⁻¹ dosages of aluminium-sulphate accelerator, compared with model fitting.</figcaption></figure>

in the reaction is $\Delta\overline{V} = \Delta\overline{V}_s - 36\overline{V}_H = -76.3 \text{ cm}^3\cdot\text{mol}^{-1}$, indicating a chemical shrinkage as usually expected in hydration reactions.

### 3.3. Hydration kinetics

In the absence of calorimetric measurements, we used the evolution of the P-wave velocity from ultrasonic measurements to calibrate the hydration kinetics. This is based on various experimental results showing a linear correlation between the P-wave velocity and hydration degree obtained from calorimetric measurements [43–45]. Experimental data was fitted with the hydration model described in Section 2.3.1. The ultimate value of the hydration extent $\xi_{\infty}$ was estimated for each mix based on the quantity of accelerator, which is supposed to be the limiting reactant of ettringite production. From [42], the estimated ultimate hydration extents are $20 \text{ mol} \cdot \text{m}^{-3}$ for Mix 1 and $16 \text{ mol} \cdot \text{m}^{-3}$ and $10 \text{ mol} \cdot \text{m}^{-3}$ for Mix 2 with $7 \text{ mL} \cdot \text{L}^{-1}$ and $11 \text{ mL} \cdot \text{L}^{-1}$ accelerator, respectively.

The results of the fit for Mix 1 can be seen in Fig. 2, and the corresponding values of the parameters are given in the recapitulating Table 4. For Mix 2, where we consider two different accelerator dosages, the model was fitted to the case with accelerator dosage $11 \text{ mL} \cdot \text{L}^{-1}$. The same set of parameters was retained for the lower accelerator dosage at $7 \text{ mL} \cdot \text{L}^{-1}$, and shows good agreement with the kinetics deduced experimentally from the P-waves, as seen in Fig. 3. This model can therefore be used to extrapolate hydration kinetics data to other accelerator dosages by adjusting only the ultimate extent of hydration $\xi_{\infty}$. Model parameters can be found in Table 5.

### 3.4. Porosity variations and chemical stresses

Porosity of cement-based materials is known in general to decrease with hydration. However, at the very early age, this is not what is observed for mortar mixed with aluminium-based accelerators, for

which a significant increase in porosity occurs as can be seen from the results of mercury intrusion porosimetry (MIP) show (Fig. 4).

Indeed, it is known that the formation of ettringite can lead to a significant expansion [29]. Porosity measurements can be used to estimate the volumetric chemical expansion $\epsilon_V^\xi$ under the hypothesis of a solely chemical action on porosity:

$$d\phi = -\Delta\overline{V}_s d\xi + \epsilon_V^\xi d\xi = (\epsilon_V^\xi - \Delta\overline{V}_s)d\xi. \qquad (29)$$

From the porosity measurements (Fig. 4), we consider an increase in porosity of 0.008 between 15 min and 45 min. During this time interval, the hydration extent goes from 10.4 to 16.8, so that the volumetric chemical expansion is:

$$\epsilon_V^\xi = \Delta\overline{V}_s + \frac{0.008}{16.8 - 10.4} \approx 1.8 \cdot 10^{-3} \text{ m}^3\cdot\text{mol}^{-1}. \qquad (30)$$

### 3.5. Water retention curve from mercury intrusion porosimetry (MIP)

The water retention curves of porous materials can be approximated using the results of mercury intrusion porosimetry [53]. Indeed, the progressive intrusion of the medium by a nonwetting fluid (mercury) is a process analogous to that of the drying a medium initially filled with water. We leverage Washburn’s equation linking the capillary radius $r$ to the pressure $p$ and the surface tension $\gamma$ of the fluid, and to the contact angle $\theta$ between the fluid and the solid:

$$r = - \frac{2\gamma \cos(\theta)}{p}, \qquad (31)$$

which we can write for both water and mercury and combine to express the capillary pressure $p_c$ in terms of the mercury intrusion pressure $p_{Hg}$:

$$p_c = - \frac{\gamma_w \cos(\theta_w)}{\gamma_{Hg} \cos(\theta_{Hg})} p_{Hg}, \qquad (32)$$

where subscripts $w$ and $Hg$ refer to water and mercury, respectively. Fig. 5 shows the water retention curves extracted from MIP data after freeze-drying of mortar samples at different times (15, 30, 45 min and at the hardened state). The surface tensions considered were $\gamma_w = 72.8 \text{ mN} \cdot \text{m}^{-1}$ and $\gamma_{Hg} = 436 \text{ mN} \cdot \text{m}^{-1}$, with contact angles $\theta_w = 40^\circ$ and $\theta_{Hg} = 130^\circ$. No significant evolution of the water retention curve can be observed in the first 45 min, a typical time frame for a print. For studies limited to the printing process, it is reasonable to adopt a constant water retention curve ($k_{sl} = 0$). However, studies aimed at longer-term simulations should take into account the dependence of the water retention curve on hydration, as it is apparent in Fig. 5 that there is an evolution between fresh and hardened states [24].

At the fresh state, there are low capillary pressures associated with states close to saturation, which cannot be rendered through a classical Van Genuchten model. Inaccurate modelling of this part of the retention curve can have a significant impact on the model at a very early age: the material is expected to be slightly desaturated and an overestimation of capillary pressure would lead to unrealistic drying shrinkage magnitude. In fact, there is a fraction of macropores that can be rendered using a linear combination of two retention curves, one for macropores and another for micropores, as has already been done for soils with multimodal pore size distributions [54,55]. Introducing a macroporosity $\phi_M$, and a couple of retention curves $S_{l,M}$ and $S_{l,m}$ characteristic of macropores and micropores, respectively, the saturation degree $S_l$ can thus be expressed as:

$$S_l = \frac{\phi_M}{\phi} S_{l,M} + \left(1 - \frac{\phi_M}{\phi}\right) S_{l,m}, \qquad (33)$$

with $S_{l,M}$ and $S_{l,m}$ each given by a Van Genuchten model:

$$S_{l,M} = \left(1 + \left(\frac{p_c}{p_{r,M}}\right)^{n_M}\right)^{\frac{1-n_M}{n_M}}, \qquad (34)$$
<table><caption>**Table 4**
Summary of model parameters for Mix 1.</caption><thead><tr><th>Symbol</th><th>Value</th><th>Unit</th><th>Description</th><th>Determination</th></tr></thead><tbody><tr><td colspan="5"><strong>Thermoporoelastic parameters</strong></td></tr><tr><td>E<sub>0</sub></td><td>1·10<sup>5</sup></td><td>Pa</td><td>Initial Young's modulus (ξ = 0)</td><td>Ref. [46]</td></tr><tr><td>E<sub>z</sub></td><td>4·10<sup>9</sup></td><td>Pa</td><td>Chemical stiffening modulus</td><td>From uniaxial test</td></tr><tr><td>ν</td><td>0.2</td><td>[-]</td><td>Drained Poisson's ratio</td><td>Assumption of a constant ratio typical for mortar</td></tr><tr><td>b</td><td>1</td><td>[-]</td><td>Biot's effective stress coefficient</td><td>Assumption based on the early age (low skeleton stiffness)</td></tr><tr><td>α<sub>s</sub></td><td>1·10<sup>-5</sup></td><td>K<sup>-1</sup></td><td>Solid thermal expansion</td><td>From clinker thermal expansion [47]</td></tr><tr><td>K<sub>s</sub></td><td>1.35·10<sup>11</sup></td><td>Pa</td><td>Solid stiffness</td><td>Stiffness of clinker constituents [48]</td></tr><tr><td>C<sub>σs</sub></td><td>2.4·10<sup>6</sup></td><td>J·K<sup>-1</sup>·m<sup>-3</sup></td><td>Solid heat capacity at constant stress</td><td>For clinker with ρ<sub>s</sub> = 3200 kg·m<sup>-3</sup> [49]</td></tr><tr><td>L</td><td>4.52·10<sup>5</sup></td><td>J·mol<sup>-1</sup></td><td>Latent heat of hydration</td><td>Ref. [50]</td></tr><tr><td>φ<sub>0</sub>(0)</td><td>0.21</td><td>[-]</td><td>Initial porosity</td><td>MIP results</td></tr><tr><td>ΔV&#773;<sub>s</sub></td><td>5.74·10<sup>-4</sup></td><td>m<sup>3</sup>·mol<sup>-1</sup></td><td>Solid molar volume variation</td><td>Ref. [50]</td></tr><tr><td>ε<sup>z</sup><sub>V</sub></td><td>1.8·10<sup>-3</sup></td><td>m<sup>3</sup>·mol<sup>-1</sup></td><td>Expansion from ettringite formation</td><td>From porosity variation (Eq. (30))</td></tr><tr><td colspan="5"><strong>Hydration kinetics</strong></td></tr><tr><td>ε<sup>z</sup><sub>∞</sub></td><td>20</td><td>mol·m<sup>-3</sup></td><td>Ultimate hydration extent</td><td>Based on accelerator dosage from [42]</td></tr><tr><td>E<sub>A</sub></td><td>4.15·10<sup>4</sup></td><td>J·mol<sup>-1</sup></td><td>Hydration activation energy</td><td>Ref. [51] (C<sub>s</sub>S → AFt)</td></tr><tr><td>A<sub>0</sub>/k<sub>z</sub></td><td>18.46</td><td>mol<sup>2</sup>·m<sup>-6</sup></td><td></td><td>Fitting on ultrasonic wave data</td></tr><tr><td>k<sub>z</sub>/η<sub>z,0</sub></td><td>8609</td><td>m<sup>3</sup>·mol<sup>-1</sup>·s<sup>-1</sup></td><td></td><td>Fitting on ultrasonic wave data</td></tr><tr><td>η&#773;</td><td>2.96</td><td>[-]</td><td>Viscosity due to water diffusion through hydrates</td><td>Fitting on ultrasonic wave data</td></tr><tr><td colspan="5"><strong>Retention curve</strong></td></tr><tr><td>p<sub>r,m</sub></td><td>6·10<sup>5</sup></td><td>Pa</td><td>Micropores reference pressure</td><td>Fitting on MIP data</td></tr><tr><td>n<sub>m</sub></td><td>2.7</td><td>[-]</td><td>Exponent of VG model for micropores</td><td>Fitting on MIP data</td></tr><tr><td>k<sub>s,l</sub></td><td>0</td><td>[-]</td><td>WRC logarithmic shift with hydration</td><td>Fitting on MIP data</td></tr><tr><td>p<sub>r,M</sub></td><td>850</td><td>Pa</td><td>Macropores reference pressure</td><td>Fitting on MIP data</td></tr><tr><td>n<sub>M</sub></td><td>1.7</td><td>[-]</td><td>Exponent of VG model for macropores</td><td>Fitting on MIP data</td></tr><tr><td>f<sub>M</sub></td><td>0.25</td><td>[-]</td><td>Fraction of macroporosity from ettringite expansion</td><td>Fitting on MIP data</td></tr><tr><td colspan="5"><strong>Plasticity</strong></td></tr><tr><td>k</td><td>-0.7</td><td>[-]</td><td>ACC asymmetry coefficient</td><td>Ref. [24]</td></tr><tr><td>h<sub>z</sub></td><td>8·10<sup>7</sup></td><td>Pa</td><td>Chemical hardening modulus</td><td>From uniaxial test</td></tr><tr><td>h<sub>z</sub></td><td>100</td><td>[-]</td><td>Isotropic hardening coefficient</td><td>Ref. [24]</td></tr><tr><td>M<sub>z</sub></td><td>35</td><td>[-]</td><td>Chemical friction hardening coefficient</td><td>Assumption based on low shear strength at early age</td></tr><tr><td>M<sub>0</sub></td><td>0.4</td><td>[-]</td><td>Initial friction ratio</td><td>Assumption based on low shear strength at early age</td></tr><tr><td>f<sub>i</sub></td><td>0.1</td><td>[-]</td><td>Isotropic tensile to compressive strength ratio</td><td>Classical mortar assumption for traction-compression ratio</td></tr><tr><td colspan="5"><strong>Fluid and heat transport</strong></td></tr><tr><td>k<sub>T</sub></td><td>1.5</td><td>W·m<sup>-1</sup>·K<sup>-1</sup></td><td>Thermal conductivity</td><td>Ref. [52]</td></tr><tr><td>k<sub>i</sub>(e<sub>0</sub>)</td><td>10<sup>-16</sup></td><td>m<sup>2</sup></td><td>Initial permeability</td><td>Permeability cell measurement</td></tr><tr><td>C<sub>k,i</sub></td><td>0.77</td><td>[-]</td><td>Taylor permeability model coefficient</td><td>Ref. [39]</td></tr><tr><td>ν<sub>w</sub></td><td>36</td><td>[-]</td><td>Water consumption stoichiometric coefficient</td><td>Ref. [50]</td></tr></tbody></table>

$$S_{l,m} = \left( 1 + \left( \frac{p_c}{p_{r,m}} \right)^{\frac{n_m}{n_m}} \right)^{\frac{1-n_m}{n_m}} . \quad (35)$$

MIP results clearly indicate an increase in the macroporosity $\phi_M$ with hydration in the early age, linked to the expansive behaviour due to the formation of ettringite. Under the hypothesis that the expansion mechanism induces the creation of macropores, denoting $f_M$ the fraction of porosity induced by the expansion which is macroporosity, its evolution reads:

$$d\phi_M = f_M d\phi^{exp} = f_M \epsilon_z^z d\epsilon_z. \quad (36)$$

Fitting of the model on experimental data gives the following parameters: $p_{r,m} = 650$ kPa, $n_m = 2.7$, $p_{r,M} = 730$ Pa, $n_M = 1.54$, as well as a fraction $f_M = 0.25$. We want to stress here the fact that $f_M = 0.25$ does not mean that the macroporosity represents 25% of the porosity, which would be very significant. Rather, this coefficient only concerns the additional porosity generated from the expansive ettringite hydration mechanics, so that the macroporosity only represents a few percents of the total porosity.

### 3.6. Shear vane test, slugs measurement and initial yield stress

On Mix 2, we used a pocket shear vane to determine the yield stress evolution of the material [41]. We use this test to determine the chemical hardening parameters $h_z$ and $M_z$ from the plasticity model. We consider that the shear vane test corresponds to a pure shear test with a shear stress $\tau_s$, so that $\sigma' = 0$ and $\tau = \sqrt{3}\tau_s$. This measurement gives one point on the yield surface corresponding to the intersection

with the deviatoric axis $\tau$. Fig. 6 shows the measured yield stress for Mix 2 at 7 mL·L<sup>-1</sup> and 11 mL·L<sup>-1</sup> of accelerating admixture.

The hydration extent at deposition is necessarily non-zero, and is considered as 2% of the ultimate value by analogy with the hydration degree of ordinary portland cement at the end of the induction period, which is generally deemed to be of the order of 2 to 4% [56]. Upon this hypothesis, fitting of the parameters was done on the experimental data for 11 mL·L<sup>-1</sup> accelerator dosage, and gives the following parameters: $h_z = 287$ MPa, $M_0 = 0.058$ and $M_z = 54$. It can be seen from Fig. 6 that the model provides a very good fit in the case of 11 mL·L<sup>-1</sup> of accelerator. In addition, the same parameter set is also able to well reproduce the evolution in yield stress in the case of 7 mL·L<sup>-1</sup> accelerator dosage. This underlines the generalisation capability of the model, which is able to extrapolate the behaviour to different admixture contents.

### 4. 3D-printing simulation framework

Simulating the process of 3D-printing with cementitious materials requires not only a material model as presented through previous sections, but also a modelling framework that defines the way in which the geometry, process, and environmental interactions are handled. The classical finite element method (FEM) has been used in the past in conjunction with time-dependent mechanical constitutive models in 2D-axisymmetric simulations [46] as well as with more complex geometries [57]. Other types of modelling, such as computational fluid dynamics (CFD) based or particle finite element methods have also been used, especially focussing on the deposition method and precise lace geometry [58–60]. Herein, we propose an approach based on FEM,
<table><caption>**Table 5**
Summary of model parameters for Mix 2.</caption><thead><tr><th>Symbol</th><th>Value</th><th>Unit</th><th>Description</th><th>Determination</th></tr></thead><tbody><tr><td colspan="5"><strong>Thermoporoelastic parameters</strong></td></tr><tr><td>E<sub>0</sub></td><td>1·10<sup>5</sup></td><td>Pa</td><td>Initial Young's modulus (ε = 0)</td><td></td></tr><tr><td>E<sub>z</sub></td><td>4·10<sup>9</sup></td><td>Pa</td><td>Chemical stiffening modulus</td><td>From Mix 1</td></tr><tr><td>ν</td><td>0.2</td><td>[-]</td><td>Drained Poisson's ratio</td><td>Assumption of a constant ratio typical for mortar</td></tr><tr><td>b</td><td>1</td><td>[-]</td><td>Biot's effective stress coefficient</td><td>Assumption based on the early age (low skeleton stiffness)</td></tr><tr><td>α<sub>s</sub></td><td>1·10<sup>-5</sup></td><td>K<sup>-1</sup></td><td>Solid thermal expansion</td><td>From clinker thermal expansion [47]</td></tr><tr><td>K<sub>s</sub></td><td>1.35·10<sup>11</sup></td><td>Pa</td><td>Solid stiffness</td><td>Stiffness of clinker constituents [48]</td></tr><tr><td>C<sub>σs</sub></td><td>2.4·10<sup>6</sup></td><td>J·K<sup>-1</sup>·m<sup>-3</sup></td><td>Solid heat capacity at constant stress</td><td>For clinker with ρ<sub>s</sub> = 3200 kg·m<sup>-3</sup> [49]</td></tr><tr><td>L</td><td>4.52·10<sup>5</sup></td><td>J·mol<sup>-1</sup></td><td>Latent heat of hydration</td><td>Ref. [50]</td></tr><tr><td>φ<sub>0</sub>(0)</td><td>0.24</td><td>[-]</td><td>Initial porosity</td><td>From mortar composition</td></tr><tr><td>ΔV&#773;<sub>s</sub></td><td>5.74·10<sup>-4</sup></td><td>m<sup>3</sup>·mol<sup>-1</sup></td><td>Solid molar volume variation</td><td>Ref. [50]</td></tr><tr><td>ε<sup>s</sup><sub>V&#773;</sub></td><td>1.8·10<sup>-3</sup></td><td>m<sup>3</sup>·mol<sup>-1</sup></td><td>Expansion from ettringite formation</td><td>From Mix 1</td></tr><tr><td colspan="5"><strong>Hydration kinetics</strong></td></tr><tr><td>ε<sup>s</sup><sub>∞</sub></td><td>16 or 10</td><td>mol·m<sup>-3</sup></td><td>Ultimate hydration extent</td><td>Based on accelerator dosage from [42]</td></tr><tr><td>E<sub>A</sub></td><td>4.15·10<sup>4</sup></td><td>J·mol<sup>-1</sup></td><td>Hydration activation energy</td><td>Ref. [51] (C<sub>s</sub>S → AFt)</td></tr><tr><td>A<sub>0</sub>/k<sub>z</sub></td><td>4.5</td><td>mol<sup>2</sup>·m<sup>-6</sup></td><td></td><td>Fitting on ultrasonic wave data</td></tr><tr><td>k<sub>z</sub>/η<sub>z,0</sub></td><td>7000</td><td>m<sup>3</sup>·mol<sup>-1</sup>·s<sup>-1</sup></td><td></td><td>Fitting on ultrasonic wave data</td></tr><tr><td>η&#773;</td><td>3</td><td>[-]</td><td>Viscosity due to water diffusion through hydrates</td><td>Fitting on ultrasonic wave data</td></tr><tr><td colspan="5"><strong>Retention curve</strong></td></tr><tr><td>p<sub>r,m</sub></td><td>6·10<sup>5</sup></td><td>Pa</td><td>Micropores reference pressure</td><td>From Mix 1</td></tr><tr><td>n<sub>m</sub></td><td>2.7</td><td>[-]</td><td>Exponent of VG model for micropores</td><td>From Mix 1</td></tr><tr><td>k<sub>s,l</sub></td><td>0</td><td>[-]</td><td>WRC logarithmic shift with hydration</td><td>From Mix 1</td></tr><tr><td>p<sub>r,M</sub></td><td>850</td><td>Pa</td><td>Macropores reference pressure</td><td>From Mix 1</td></tr><tr><td>n<sub>M</sub></td><td>1.7</td><td>[-]</td><td>Exponent of VG model for macropores</td><td>From Mix 1</td></tr><tr><td>f<sub>M</sub></td><td>0.25</td><td>[-]</td><td>Fraction of macroporosity from ettringite expansion</td><td>From Mix 1</td></tr><tr><td colspan="5"><strong>Plasticity</strong></td></tr><tr><td>k</td><td>-0.7</td><td>[-]</td><td>ACC asymmetry coefficient</td><td>Ref. [24]</td></tr><tr><td>h<sub>z</sub></td><td>8·10<sup>7</sup></td><td>Pa</td><td>Chemical hardening modulus</td><td>From shear vane tests</td></tr><tr><td>h<sub>s</sub></td><td>100</td><td>[-]</td><td>Isotropic hardening coefficient</td><td>Ref. [24]</td></tr><tr><td>M<sub>z</sub></td><td>35</td><td>[-]</td><td>Chemical friction hardening coefficient</td><td>From shear vane tests</td></tr><tr><td>M<sub>0</sub></td><td>0.4</td><td>[-]</td><td>Initial friction ratio</td><td>From shear vane tests</td></tr><tr><td>f<sub>i</sub></td><td>0.1</td><td>[-]</td><td>Isotropic tensile to compressive strength ratio</td><td>Classical mortar assumption for traction-compression ratio</td></tr><tr><td colspan="5"><strong>Fluid and heat transport</strong></td></tr><tr><td>k<sub>T</sub></td><td>1.5</td><td>W·m<sup>-1</sup>·K<sup>-1</sup></td><td>Thermal conductivity</td><td>Ref. [52]</td></tr><tr><td>k<sub>i</sub>(ε<sub>0</sub>)</td><td>10<sup>-16</sup></td><td>m<sup>2</sup></td><td>Initial permeability</td><td>From Mix 1</td></tr><tr><td>C<sub>k<sub>i</sub></sub></td><td>0.77</td><td>[-]</td><td>Taylor permeability model coefficient</td><td>Ref. [39]</td></tr><tr><td>ν<sub>w</sub></td><td>36</td><td>[-]</td><td>Water consumption stoichiometric coefficient</td><td>Ref. [50]</td></tr></tbody></table>

with a mesh evolving through time as the material is deposited during printing. First, a suitable mesh must be generated from the printer's toolpath. Then a simulation framework based on FEM is presented that is capable of handling the sequential addition of elements to the current mesh. These new elements represent the newly added material in the printing process. Finally, interactions of the printpiece with the printer nozzle as well as with the environment are addressed.

## 4.1. Mesh definition from a print path

Based on the printing solution offered by the technological platform Build'in from École nationale des ponts et chaussées [61], which allows for 5-axis printing (x, y, z positioning along with x and y rotations, the z rotation serving no purpose since a circular nozzle is used), the toolpath given to print a workpiece is composed of a sequence of planes. Equivalently, the toolpath can be represented through a sequence of positions with an associated normal n to the printing plane. We will use this representation to systematically build a suitable mesh that is well segmented to represent the sequential manner in which material is deposited by the printer. Meshes are generated with the GMSH software [62], used through its API for python.

### 4.1.1. 2D-axisymmetric approach

For simple analysis of pieces obtained by revolution around the z-axis, a 2D-axisymmetric model represents a faster solution for simulation, although it does not give the ability to study dissymmetries in loading such as those induced by punctual loads from the printer's head.

In this case, only one point per layer is needed, with its coordinates (r, z) as well as the normal n to the printing plane. If (e<sub>r</sub>, e<sub>θ</sub>, e<sub>z</sub>) is a cylindrical base and ||n|| = 1, then we can define s = e<sub>θ</sub> × n so that the (e<sub>θ</sub>, s) plane corresponds to the current printed surface. In most cases, the nozzle is vertical (no rotation around e<sub>θ</sub>), in which case s = r. If O<sub>n</sub> is the position of the nozzle at the n<sup>th</sup> layer, using s we can build the edges of the top surface A<sub>n</sub> and B<sub>n</sub> so that O<sub>n</sub>A<sub>n</sub> = -w/2 s and O<sub>n</sub>B<sub>n</sub> = w/2 s with w the width of the lace. The quadrilateral A<sub>n-1</sub>B<sub>n-1</sub>B<sub>n</sub>B<sub>n</sub> then represents the section of the n<sup>th</sup> layer and will be made as a physical surface entity, which means that the mesh will have to conform to the boundaries of this quadrilateral. Fig. 7 illustrates that this meshing method allows for systematic handling of all cases, including objects with radii differing from layer to layer and non-horizontal printing.

### 4.1.2. 3D approach

In the 3D approach, the key step is the reconstruction of the section of the lace, orthogonal to the printing direction, at each point of the printing path. The volume of added material during the nozzle's travelling from one point to the next is the volume linking two consecutive lace sections.

As illustrated in Fig. 8, to construct the section at point O<sub>n</sub>, we first define the print direction that is tangent to the printing path as d = O<sub>n-1</sub>O<sub>n+1</sub>/||O<sub>n-1</sub>O<sub>n+1</sub>||. Defining t = d × n/||d × n||, the planes (O<sub>n</sub>, d, t) and (O<sub>n</sub>, t, n) then contain the top surface of the print and the print section, respectively. Similarly to the 2D case, t allows us to build the two top edges of the lace section, with the two bottom edges being the top edges of the previous layer. However, this supposes that each
<figure><img src="image_8.png" /><figcaption>Fig. 7. Different print path cases and their handling in 2D-axisymmetric conditions.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Mesh region construction principle for printing simulations.</figcaption></figure>

layer has the same number of discretisation points, which is in practice often the case due to the slicing process of the 3D model resulting in the printing path, but is not a prerequisite for printing. Should the projection on the (x, y) plane of the points of two consecutive layers be too far apart, this would also result in ill-shaped volumes, which can cause difficulties in finite element method (FEM) calculations. Nevertheless, it is always possible to re-slice the original 3D model to meet these requirements, so that it does not result in any loss of generality.

The main limitation of this method, which does not occur in Vox-
elPrint for example [63], is the case of overlapping printing paths.
Although dealing with such print paths is not impossible, it would re-
quire substantial additional work to the meshing framework presented
and is out of the scope of this study. We will hence limit our print paths
to single-lace, non-overlapping ones.

4.2. Finite element solver for sequential material deposition

The FEM resolution is done through the FEniCS [64] partial
derivative equations solver, while constitutive law integration is done
using the MFront code generator [65]. Combining these two tools
through the MFront Generic Interface Support (MGIS), fully
coupled thermohydromechanical simulations can be performed. Used
in conjunction with the meshing tool described above, it forms a
framework for the simulation of the printing process of cementitious
materials.

In contrast to usual FEM simulations, the sequential deposition of
material must be handled, which means that the geometry is different

at every time step in the calculation. The approach of this framework
is to extract a submesh from the mesh of the entire geometry and
solve the FEM problem on this submesh. From one time step to the
next, the nodal values and internal variables at the elements have
to be transported from one submesh to the next. However, due to
changes in the connectivity of the nodes from one submesh to another,
the numbering of nodes and elements may be different. In order to
circumvent this problem, we store the values at the end of the time
step in arrays referring to the mesh of the entire geometry, through
maps between nodes and elements of the submesh and entire mesh.
These arrays are initialised with the initial conditions of the material.
At the beginning of a new time step, new maps are set between the new
submesh and the entire mesh, allowing for the extraction of the initial
conditions: in elements and nodes that were already present in the
previous submesh, this pulls the results of the previous time step, and
in elements and nodes that were not present yet, this pulls the initial
conditions of the material with which the arrays have been initiated.
It allows for the continuity of the calculations, despite the usage of a
different mesh at every time step.

4.3. Boundary conditions and loads

4.3.1. Substrate interaction

In reality, the boundary condition at the substrate is not trivial.
Vertical displacements should indeed be zero (setting aside the possi-
bility of the bottom layer partially lifting off of the substrate, which
should only be the case in very specific geometries in which this should
be addressed specifically and will not be considered here in general),
but the horizontal boundary condition is unclear. Indeed, one can reasonably think of it as some kind of frictional contact, with more or less influence of friction depending, namely, on the roughness of the substrate. Frictional contact as a boundary condition is tricky to implement: herein, we will consider a perfect slip boundary condition at the base for the sake of simplicity.

### 4.3.2. Layer pressing by the nozzle

As reported by Carneau et al. [66], in extrusion-based 3D concrete printing using the layer pressing method to shape the layer from a circular nozzle, pressures equivalent to several layers' self-weight can be applied upon passage of the nozzle. It depends on the printing parameters and, more specifically, on the velocity of the printer head relative to that of the mortar flow. This could be critical in studying the stability during printing, especially considering that this pressure is applied very locally and can bring asymmetry in the stress state even in otherwise symmetric structures.

Every section of the printed object will be loaded, upon deposition, with this pressure from the nozzle head. Strictly speaking, at a given time, only the area directly beneath the nozzle head should be loaded. However, including it in this way in the simulation would require a continuous deposition of the material, in contrast to the sequential deposition methods described above. Such a simulation with continuous deposition seems hard to achieve using FEM, and may be better suited for discrete element models (DEM) or computational fluid dynamics (CFD) methods, for instance, which have been used successfully by other authors while investigating at the lace scale [58]. However, the authors' focus here lies on the structural scale, hence the adoption of a coarser approach regarding direct actions on the lace. The nozzle head pressure is then applied to the last current section of the mesh, on its top surface, so that every section is loaded at some stage. As no viscous effects are considered in the material model and considering that the length of the print sections in the mesh is typically of the same order of magnitude as the nozzle's diameter, the inadequacy between the duration of loading in the simulation compared to reality should not influence the results of the simulation.

### 4.3.3. Environmental boundary conditions

For the sake of simplicity, we will consider simplified interactions with the environment. For the temperature boundary conditions, we will impose the external temperature (20 °C), although in practice the mortar surface is at a higher temperature than room temperature and is bound to radiation and convective heat fluxes. As for the fluid, evaporation fluxes will not be taken into account here and are left for further studies. Therefore, a sealed condition ($M_{-w} \cdot \bar{n} = 0$) is enforced at the interface.

## 5. Simulation examples and discussion

In this section, we discuss the results of simulations on two object geometries based on the material properties from Mix 1, a simple cylinder and a slightly more complex vase-like shape thereafter called “corolla”. The aim is to study the case of plastic collapse failure during printing, as well as factors which may influence it, such as layer pressing [66]. With the cylinder, we show how a simple analytical criterion can be safely substituted to simulation in the case of simple geometries. The corolla reveals a greater complexity of failure mechanisms due to subtle interplay between the evolution of stress and strength in the structure. We show how different printing speeds can lead to radically different failure modes. We also compare the results of axisymmetric 2D and 3D simulations incorporating the effect of layer pressure.

### 5.1. Cylinder geometry and analytical solution

The considered geometry is a circular cylinder of 5 cm radius. The lace section is 25 mm in width and 7 mm in height, and we consider

<figure><img src="image_10.png" /><figcaption>Fig. 9. Stability of cylinder print for different printing speeds.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Temperature elevation map relative to room temperature at 40 layers ($v_{-} = 6 \text{ cm} \cdot \text{s}^{-1}$).</figcaption></figure>

the full object to be 60 layers high. The horizontal nozzle speeds $v_{-}$ considered range between 6 and $18 \text{ cm} \cdot \text{s}^{-1}$. Taking into account the geometry of the lace, this corresponds to mortar flows ranging from 0.63 to $1.89 \text{ L} \cdot \text{min}^{-1}$, which are usual flow values for printing. We show here how a simple analytical estimation using an abacus provides similar predictive capacity as the complex simulation.

For a cylinder of current height $H$ with perfect slip conditions at the base, and neglecting the effects of internal pore pressure and thermal loading, the stress state is simply uniaxial:

$$ \bar{\sigma} = \rho g(H - z) \bar{e}_z \otimes \bar{e}_z, \quad z \in [0, H]. \qquad (37) $$

In the plane of volumetric and deviatoric stresses ($\sigma, \tau$), this corresponds to a loading slope of 3. Plasticity will thus occur at the intersection point between the line $\tau = 3\sigma'$ and the yield surface. Let ($\sigma_{crit}, \tau_{crit}$) be this intersection point. Then $\tau_{crit} = 3\sigma_{crit}$ and $\sigma_{crit}$ is so that $f(\sigma_{crit}, 3\sigma_{crit}) = 0$. This stress $\sigma_{crit}$ can also be converted in terms of number of layers through the self-weight of the latter. The volumetric stress from a layer of height $h_0$ is $\rho g h_0 / 3$ ($\sigma = \sigma_{zz} / 3$), so that we can extract a critical number of layers $N_{crit} = 3\sigma_{crit} / \rho g h_0$.

Knowing the hydration kinetics, added to the fact that given the geometry of the cylinder in a plasticity-based analysis, the yielded layer will always be the bottom one, allowing for the simple analysis of
<figure><img src="image_12.png" /><figcaption>Fig. 11. Hydration profiles at different printing stages.</figcaption></figure>

failure. Fig. 9 presents this failure criterion in the form of an abacus: to a given printing speed corresponds a line whose slope is proportional to the printing speed. The intersection of this line with the failure zone represents the failure point (in time and number of layers). In blue is the line determined by $N_{crit}$, compared to the number of layers at the current time for different printing speeds. This analytical method covers only the mechanical aspect of the problem. The temperature, which is expected to be higher at the centre of the lace's section, as well as pore pressure will influence the failure. Given the short time scale of printings presented here, temperature will have a small effect on the hydration degree. However, pore pressure can influence the onset of plasticity, as the yield criterion is evaluated with respect to the effective stress $\sigma + \beta p^* \delta$ (Eq. (8)). A lower bound for the critical number of layers could then be obtained by considering a pore pressure equal to the mean stress, i.e. a null mean effective stress. $\sigma_{crit}$ in this case would be the result of $f(0, 3\sigma_{crit}) = 0$. This lower bound would constitute a safer estimate of the actual critical number of layers, which should lie between the two bounds ($p_l = 0$ and $p_l = \text{tr}(\sigma)/3$). The lower bound is represented in orange in Fig. 9. Regardless of the bound considered (upper bound corresponding to $p_l = 0$ or lower bound corresponding to $\sigma' = 0$), the analytical analysis predicts successful printings for printing speeds $v_\to$ of 6 and $10 \text{ cm} \cdot \text{s}^{-1}$, while printing at 14 and $18 \text{ cm} \cdot \text{s}^{-1}$ should result in failures in layers 20–30 and 16–21, respectively. Printability at a speed of $12 \text{ cm} \cdot \text{s}^{-1}$ is unclear as only the lower bound is crossed at 25 layers.

The validity of these two estimates is confirmed by performing 2D axisymmetric simulations for the same geometry. The results compiled in Table 6 confirm good agreement between the simulation and the expected result of the analytical results.

## 5.2. Corolla geometry

This second object is defined by the following relationship between radius $r$ and height $z$:

$$r(z) = 0.05 \left( 1 + \left( \frac{z}{0.14} \right)^2 \right) \quad (z \in [0, 0.28]). \quad (38)$$

The geometry of the lace is the same as the cylinder and the full object is comprised of 40 layers.

<table><caption>Table 6
Maximum number of layers for the cylinder, analytical vs simulation.</caption><thead><tr><th rowspan="2">v<sub>−</sub> [cm·s<sup>-1</sup>]</th><th colspan="2">Maximum nb. of layers</th></tr><tr><th>Analytical</th><th>Simulation</th></tr></thead><tbody><tr><td>6</td><td>60</td><td>60</td></tr><tr><td>10</td><td>60</td><td>60</td></tr><tr><td>12</td><td>25–60</td><td>29</td></tr><tr><td>14</td><td>19–29</td><td>19</td></tr><tr><td>18</td><td>15–20</td><td>15</td></tr></tbody></table>

### 5.2.1. Hydration profiles

In an axisymmetric analysis, properties are unchanged by any rotation around the symmetry axis, allowing for a profile study of variables along the $z$ axis. However, in the $(r, z)$ section, properties are also subject to being heterogeneous along the $r$ axis. Indeed, the latent heat of hydration will cause the mortar to heat during hydration, whereas heat dissipation at the mortar surface will lead to a lower temperature at the surface relative to the centre of the lace. In Fig. 10, it can be seen that in the case of $v_\to = 6 \text{ cm} \cdot \text{s}^{-1}$ at the 40th layer, the centre of the lace is up to 0.3 K hotter than the periphery which is kept at room temperature. This is also a lower-bound estimation of the temperature inside the structure, as a simplified constant temperature boundary condition is considered. Heating during pumping and extrusion of the material is also not taken into account in this work, and would tend to further increase the temperature of the material above that of the room. In reality, the surface of the mortar should be at a temperature higher than that of the room, and the cooling is driven mainly by convection and radiation on the surface. As hydration is thermoactivated through the activation energy $E_A$, the hydration rate $\dot{\xi}$ will be higher in the centre of the lace than at its periphery. This implies that from the effect of temperature alone, material properties are expected to be heterogeneous in the section, due to a gradient in the extent of hydration reaction.

Fig. 11 shows the hydration profiles of the corolla in the 10th, 20th, 30th and 40th layers for the different printing speeds when plasticity has not yet been reached. Printing time given the size of the object and the considered printing speeds is limited to the first 20 min of the mortar's
<figure><img src="image_13.png" /><figcaption>Fig. 12. Vertical stress, orthoradial stress and pore pressure maps at 40 layers ($v_{\infty} = 6 \text{ cm} \cdot \text{s}^{-1}$).</figcaption></figure>

age. As such, hydration kinetics is mostly linear (Fig. 2). Furthermore, the curvature of the vertical profile of the print is limited in the first layers, as the tangent to the profile at $z = 0$ is vertical ($dr_c/dz(0) = 0$). This results in a mostly linear hydration profile $\xi(z)/\xi_{\infty}$ in most of the simulated printings (Figs. 11(a), 11(b), 11(c)). Obviously, the higher the printing speed, the lower the hydration degree of a given layer. At 40 layers (Fig. 11(d)), time between each layer increases along with the radius of the object. As a consequence, the change in hydration degree from layer to layer also increases, resulting in a profile akin to that of the print section itself.

### 5.2.2. Stress states and failure mechanisms

Figs. 12(a) and 12(b) show the vertical and orthoradial stress states in the corolla at 40 layers, for printing at $v_{\infty} = 6 \text{ cm} \cdot \text{s}^{-1}$. Regarding vertical stress, the lower part of the corolla's section lies mainly in compression, whereas the upper, inner part lies in tension from the corolla's tendency to fall over by opening outward. Similarly, the lower part is on average in orthoradial compression, while the upper part is in tension, again from the geometry's tendency to open and fall over. This shift in stress state between the lower and upper parts of the section will result in more subtle failure mechanisms compared to the straight cylinder, for which plastic failure always occurred at the base. The “jagged” look of the orthoradial stress map originates in the differential swelling of the mortar in neighbouring layers. Indeed, each layer has a different magnitude of swelling because its extent of hydration is different. Their stiffness is also different, leading to the observable discontinuities in the stress at the interface. A layer towards the bottom is more hydrated, and thus is stiffer and has swollen more than a layer situated towards the top, which is consistent with each layer generally lying in traction at its bottom (at the interface with a more hydrated layer) and in compression at the top of the layer (at the interface with a less hydrated layer). Fig. 12(c) shows the repartition of the pore pressure. It is also, as expected, heterogeneous throughout the structure. The hydration reaction consumes water over time, which induces a desaturation of the pores: therefore, lower pore pressures are observed in bottom layers where the mortar is older. This suction is the driver of self-desiccation shrinkage, and as such is of paramount interest for the quality control of the geometry of print-pieces. Moreover, since poromechanical theory leads to consider plasticity as being dependent on an effective stress, the onset of plastic deformation is also a function of the evolving pore pressure.

Finding the onset point of plastic failure amounts to understanding the competition between the development of material's strength versus stresses. Fig. 13 illustrates the plastic failure – or lack thereof – of the different prints with speeds 6 to 14 cm·s⁻¹. In all cases except 6 cm·s⁻¹,

<figure><img src="image_14.png" /><figcaption>Fig. 13. Plasticity onset points and maximum height for different printing speeds.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Evaluation of stress to yield stress ratio.</figcaption></figure>

plastic strains develop before the 40th layer is reached, but at different locations: at the bottommost layer for the highest speeds at 12 and 14 cm·s⁻¹, while it occurs at roughly 8 layers from the top in the case of 8 and 10 cm·s⁻¹. Depending on the geometry, plastic collapse can then occur outside of the bottommost layer, even though it supports the most load.
<figure><img src="image_16.png" /><figcaption>Fig. 15. Profiles of stress to yield stress ratios.</figcaption></figure>

To illustrate this change in failure mechanism, it can be interesting to look at the ratio between the current stress and the yield stress. However, this can only be meaningful in one dimension: to do so, a yield stress is evaluated by extending the current stress path, so that distances are calculated along this direction, as illustrated in Fig. 14. The ratio of the equivalent stress $\tilde{\sigma}_n$ to the equivalent yield stress $\tilde{\sigma}_y$, later referred to as the “stress to yield stress ratio”, is then given by:

$$ \frac{\tilde{\sigma}_n}{\tilde{\sigma}_y} = \frac{\tilde{\sigma}_n}{\tilde{\sigma}_n + \tilde{\sigma}_r} \qquad (39) $$

This ratio varies between 0 and 1, a higher value meaning a stress state closer to the yield surface, with a value of 1 corresponding to the yield of the material ($\tilde{\sigma}_r = 0$). It allows for a simple evaluation of how much of the mortar’s resistance is attained at each point in the structure, regardless of the heterogeneity in yield stresses (Fig. 15). This shows, at each stage before failure occurs, which layer is the most critical from the perspective of plasticity. At the early stages of printing (Fig. 15(a)), the most critical layer is the bottom layer. At 18 and 22 layers, the shape of the curve has clearly shifted, and the layer in which the highest percentage of the yield stress is attained is at midheight of the current object. It is consistent with the points of initiation of plastic deformations illustrated previously in Fig. 13. This behaviour is consistent with that observed experimentally when printing similar geometries, as can be seen in Fig. 16(b): the collapsed corolla at the bottom, for which we can see the remains of the upper layers on the printing surface, has kept its first few layers intact, showing that failure, in this case, was indeed initiated from an intermediate layer rather than from the bottom.

### 5.2.3. 3D modelling and effect of layer pressing

Axisymmetric simulations, although significantly more computationally efficient than fully 3D analysis, cannot render the true complexity of the printing process itself. Even on printed objects displaying an axis of symmetry, the sequential manner in which material is deposited is so that the printpiece is never truly symmetric, be it in its loads or in its material properties. In that sense, the results of the 3D simulation are expected to differ from those of their 2D counterparts.

In fact, 3D simulation of a printing of the same geometry reveals a significant drop in the predicted buildability for printing speeds of 5 cm·s⁻¹ and higher (Fig. 18): for instance, in the case of 6 cm·s⁻¹, the number of printable layers is reduced from 16 to only 10 layers. Plasticity appears locally in the bottom layer as the new layer is starting to be printed, vertically from the freshly deposited sections, as can be seen in Fig. 17. This drop in the number of printable layers can be explained by three factors. Firstly, the partial deposition of the layer means that it does not contribute to the structural rigidity of the print piece in the way a full layer would. Secondly, the mean hydration extent of each layer is lower than in the 2D case where the underlying assumption is that the whole layer is being deposited instantly at the time when the printer starts to print this layer. Thirdly, the self-weight load of the partially deposited layer will induce a nonsymmetric deformation of the printpiece. The accumulation of these three factors thus leads to worse buildability in reality than the 2D analysis could predict.

However, it is worth noting that the importance of partial layer deposition can be expected to be highly dependent on the print-piece’s geometry. Similar analysis on the cylinder piece with a printing speed $v_↓ = 14 \text{ cm} \cdot \text{s}^{-1}$, where the number of buildable layers decreases from 19 (Table 6) to 18 layers, almost no decrease compared to the case of the corolla. The overhang in the corolla is probably more detrimental as it constitutes an eccentric load from a structural point of view.

In addition to the partial deposition of layers, the printer nozzle can, in some cases, apply pressure on the object that is being printed. This is true for instance in processes which leverage layer-pressing in order to control the geometry of the lace, which is common with 2 component systems. As the pressure applied from layer-pressing can be of the same order of magnitude as the self-weight of the lace [66], we propose here to study its effect by taking a pressure value corresponding to the self-weight of one layer. The considered layer-pressing pressure $p_↓$ is then:

$$ p_↓ = \rho g h_0 \approx 160 \text{ Pa.} \qquad (40) $$

Results compiled in Fig. 18 for the corolla geometry show a further reduction in the number of printable layers, although it should be noted that the difference between the predictions of the 2D and
<figure><img src="image_17.png" /><figcaption>Fig. 16. Successful printing with carbon fibre reinforcement (a) and failed printing without reinforcement (b) of a geometry similar to the presented simulations (courtesy of V. Debono, XtreeE). Simulations do not take into account reinforcement.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Local onset of plasticity due to partially deposited top layer, $v_-$ = 10 cm·s⁻¹.</figcaption></figure>

<caption>Fig. 18. Buildability assessment from 2D axisymmetric and 3D (with and without layer pressing) simulations of the corolla.</caption>

3D simulations without considering the layer pressing phenomenon is significantly more pronounced, although the layer-pressing pressure $p_-$ can be higher than 160 Pa in practice, which would enlarge the gap [66].

Overall, whether a particular printing process utilises layer-pressing or not, it appears that 3D simulations have their interest even in the simplest geometries in order to reach a certain level of accuracy. 2D simulations should thus be analysed carefully as they constitute a priori an overestimation of the number of printable layers.

## 6. Applications: predictive capacity and printing optimisation pathways

In this section, we propose a quantitative comparison between numerical model and experiment based on printing with varying speed and accelerator dosage, on the basis of Mix 2, followed by an example of a possible optimisation scheme on printing speed of the robot.

### 6.1. Change of printing speed

Two identical objects are printed with accelerator dosage 11 mL·L⁻¹ and nozzle speeds of 13 cm·s⁻¹ and 20 cm·s⁻¹, respectively. Fig. 19 shows the failure at 13 cm·s⁻¹, where it can be seen that in the experiment, the object failed during the printing of layer 40 with the top 7 layers collectively failing while the bottom 32 still stand after failure. The simulation predicts a failure when printing layer 44. As can be seen in Fig. 19b, around 10 layers are fully plastified in the whole layer width. This indicates that the bottom 33 layers would not fail, which is consistent with the experimental finding.

Fig. 20 shows an identical analysis for printing at 20 cm·s⁻¹. Experimentally, the object failed while printing layer 39, close to the result obtained by printing at 13 cm·s⁻¹, with a similar failure pattern: the top 7 layers fail collectively while the bottom 31 remain intact. The simulation predicts failure while printing layer 42, and shows that plasticity occurs both in the bottom layer and in the upper part of the object, as seen in Fig. 20b. However, the bottom layer is only plastified on the outer side of the object, while some layers in the upper part are
<figure><img src="image_19.png" /><figcaption>**Fig. 19.** Failure of the printed object with Mix 2, accelerator dosage 11 mL·L⁻¹ and printing speed 13 cm·s⁻¹: (a) Failure during experiment at layer 40 with superimposed simulation to scale; (b) Plastified zone in the simulation when predicting failure at layer 44.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 20.** Failure of the printed object with Mix 2, accelerator dosage 11 mL·L⁻¹ and printing speed 20 cm·s⁻¹: (a) Failure during experiment at layer 39 with superimposed simulation to scale; (b) Plastified zone in the simulation when predicting failure at layer 42.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. 21.** Failure of the printed object with Mix 2, accelerator dosage 7 mL·L⁻¹ and printing speed 20 cm·s⁻¹: (a) Failure during experiment at layer 27 with superimposed simulation to scale; (b) Plastified zone in the simulation when predicting failure at layer 27.</figcaption></figure>

fully plastified, which points to a failure occurring in the upper layers similarly to the previous case. Again, fully plastified layers are located above layer 32, which is in good agreement with the experimental observation.

## 6.2. Change of accelerator dosage

An identical object is then printed while maintaining a printing speed of 20 cm·s⁻¹ and reducing the accelerator dosage to 7 mL·L⁻¹. Fig. 21 shows the comparison between the experiment and simulation in this case. Failure occurred in the experiment while printing layer 27, with failure originating from the first printed layer which yields under the weight of the object. The simulation is done with the same set of parameters, except for the ultimate hydration extent $\xi_{\infty}$, which is proportional to accelerator dosage, as explained previously. The simulation predicts a failure at layer 27, which originates from the first layer, as can be seen in Fig. 21b, in perfect agreement with the experimental observation. Overall, on the 3 presented test cases, the model error in terms of predicted number of printable layers is at most 10%.

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

<figure><img src="image_23.png" /><figcaption>**Fig. 22.** Optimised printing compared with printing at a constant $v_{\infty} = 6$ cm·s⁻¹.</figcaption></figure>

## 6.3. Towards print planning: optimising process parameters for each print-piece

Printpath planning has been widely discussed for other technologies, such as plastic extrusion, and aims at finding optimal nozzle trajectories to print objects [67–69]. For concrete 3D printing, printpath planning is also important to broaden the range of printable geometries, for example, by using non-planar nozzle paths [70]. However, many current applications of 3D printing with cementitious materials involve thin-walled, vertical structures for which the range of possible printpaths is scarce. Then of greater interest is the optimisation of the process parameters, such as printing speed, material mix design, or other parameters related to specific printing technologies. Beyond simple buildability assessment, the perspective of developing simulation tools for 3D printing of cement-based materials is to derive optimisation schemes that will ensure that printing is done under optimal conditions.

From the presented simulation framework, we derive a simple optimisation scheme on the printing speed $v_{\infty}$: for each layer, using a dichotomy algorithm, the highest speed such that no plastic deformations occur is selected, and the process is then repeated for each subsequent layer. By doing so, the algorithm ensures both the buildability of the print-piece and an optimal total printing time.

Applied to the object geometry considered in this section, the optimised nozzle speed profile by layer is shown in Fig. 22. A maximum printing speed of 30 cm·s⁻¹ was enforced to reflect the limitations of the printer. From the 887 s necessary to print the full object at a constant speed of 6 cm·s⁻¹, the optimisation reduces this time to 604 s, thus reducing the printing time by 32%.

## 7. Conclusions

Extrusion-based 3D printing with cementitious materials is a promising technology to reduce the carbon footprint of the construction industry. However, difficulties in assessing the buildability of complex geometries as well as predicting long-term behaviours of print pieces remains are some of the challenges that must be overcome in order for this technology to reach its potential. Numerical modelling that takes into account both the printed material and the printing process could help answer these unknowns and make 3D printing more predictable and reproducible, as well as open a pathway for process optimisation.

Herein, we presented a simulation framework aimed at automatised simulations of extrusion-based 3D printing with cement-based materials. Based on a coupled chemo-thermoporomechanical material model, it encompasses all couplings necessary to investigate buildability as
well as drying and temperature related issues, paramount in 3D printing and which will be addressed in future studies. We demonstrated the strength of this approach compared with phenomenological ones as it was able to extrapolate accurately to different accelerator contents with a single set of parameters. The processes side is also well taken into account, through a meshing tool from robot toolpath and a finite element method based solver mimicking the deposition process. This framework allows for the creation of a “digital twin” of the printing process itself, a powerful tool for buildability assessment of complex 3D-printed geometries.

Using this framework on different geometries and with different printing parameters reveals the difficulty of predicting the onset of plastic failure. In addition to coupled phenomena in the material itself, the printing process introduces another layer of complexity with the sequential deposition of the material and potential layer pressing, supporting the need for fully 3D simulations even in the case of simple geometries with a symmetry axis. Consequently, in addition to testing the raw material, measurements of process parameters such as the force applied by the printer head on the print piece seem of great interest in predicting buildability.

Beyond the assessment of buildability, it has been shown how such a framework can be used to optimise process parameters such as printing speed and maximise productivity. This opens the way for new workflows in 3D concrete printing: from the design of a print-piece, using simulations tools for printpath and process parameter planning, in order to operate the printer under safe and optimised conditions.

### CRediT authorship contribution statement

**Maxime Pierre:** Writing – review & editing, Writing – original draft, Software, Methodology, Conceptualization. **Siavash Ghabezloo:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Patrick Dangla:** Writing – review & editing, Methodology, Conceptualization. **Romain Mesnil:** Writing – review & editing, Methodology, Conceptualization. **Matthieu Vandamme:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Jean-François Caron:** Writing – review & editing, Project administration, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### References

[1] C. Carson, M. Abbott, A review of productivity analysis of the New Zealand construction industry, Australas. J. Constr. Econ. Build. 12 (3) (2012) 1–15, http://dx.doi.org/10.5130/AJCEB.v12i3.2584.

[2] IEA, 2019 global status report for buildings and construction, 2019.

[3] S.A. Miller, A. Horvath, P.J. Monteiro, Readily implementable techniques can cut annual CO2 emissions from the production of concrete by over 20%, Environ. Res. Lett. 11 (7) (2016) 074029, http://dx.doi.org/10.1088/1748-9326/11/7/074029.

[4] G.M. Moelich, J. Kruger, R. Combrinck, Modelling the interlayer bond strength of 3D printed concrete with surface moisture, Cem. Concr. Res. 150 (2021) 106559, http://dx.doi.org/10.1016/j.cemconres.2021.106559.

[5] E. Keita, H. Bessaies-Bey, W. Zuo, P. Belin, N. Roussel, Weak bond strength between successive layers in extrusion-based additive manufacturing: measurement and physical origin, Cem. Concr. Res. 123 (2019) 105787, http://dx.doi.org/10.1016/j.cemconres.2019.105787.

[6] J. Xiao, H. Liu, T. Ding, Finite element analysis on the anisotropic behavior of 3D printed concrete under compression and flexure, Addit. Manuf. 39 (2021) 101712, http://dx.doi.org/10.1016/j.addma.2020.101712.

[7] R.M. Ghantous, Y. Valadez-Carranza, S.R. Reese, W.J. Weiss, Impact of drying of 3D printed cementitious pastes on their degree of hydration, in: Third RILEM International Conference on Concrete and Digital Fabrication: Digital Concrete 2022, Springer, 2022, pp. 216–221, http://dx.doi.org/10.1007/978-3-031-06116-5_32.

[8] J. Van Der Putten, M. De Volder, P. Van den Heede, M. Deprez, V. Cnudde, G. De Schutter, K. Van Tittelboom, Transport properties of 3D printed cementitious materials with prolonged time gap between successive layers, Cem. Concr. Res. 155 (2022) 106777, http://dx.doi.org/10.1016/j.cemconres.2022.106777.

[9] Y. Zhang, H. Qiao, R. Qian, C. Xue, Q. Feng, L. Su, Y. Zhang, G. Liu, H. Du, Relationship between water transport behaviour and interlayer voids of 3D printed concrete, Constr. Build. Mater. 326 (2022) 126731, http://dx.doi.org/10.1016/j.conbuildmat.2022.126731.

[10] A. Das, Y. Song, S. Mantellato, T. Wangler, D.A. Lange, R.J. Flatt, Effect of processing on the air void system of 3D printed concrete, Cem. Concr. Res. 156 (2022) 106789, http://dx.doi.org/10.1016/j.cemconres.2022.106789.

[11] S. Ghabezloo, J. Sulem, S. Guédon, F. Martineau, J. Saint-Marc, Poromechanical behaviour of hardened cement paste under isotropic loading, Cem. Concr. Res. 38 (12) (2008) 1424–1437, http://dx.doi.org/10.1016/j.cemconres.2008.06.007.

[12] S. Ghabezloo, J. Sulem, J. Saint-Marc, The effect of undrained heating on a fluid-saturated hardened cement paste, Cem. Concr. Res. 39 (1) (2009) 54–64, http://dx.doi.org/10.1016/j.cemconres.2008.09.004.

[13] S. Ghabezloo, J. Sulem, J. Saint-Marc, Evaluation of a permeability–porosity relationship in a low-permeability creeping material using a single transient test, Int. J. Rock Mech. Min. Sci. 46 (4) (2009) 761–768, http://dx.doi.org/10.1016/j.ijrmms.2008.10.003.

[14] M.-H. Vu, J. Sulem, S. Ghabezloo, J.-B. Laudet, A. Garnier, S. Guédon, Time-dependent behaviour of hardened cement paste under isotropic loading, Cem. Concr. Res. 42 (6) (2012) 789–797, http://dx.doi.org/10.1016/j.cemconres.2012.03.002.

[15] N. Agofack, S. Ghabezloo, J. Sulem, A. Garnier, C. Urbanczyk, Experimental investigation of the early-age mechanical behaviour of oil-well cement paste, Cem. Concr. Res. 117 (2019) 91–102, http://dx.doi.org/10.1016/j.cemconres.2019.01.001.

[16] N. Agofack, S. Ghabezloo, J. Sulem, A. Garnier, C. Urbanczyk, Onset of creation of residual strain during the hydration of oil-well cement paste, Cem. Concr. Res. 116 (2019) 27–37, http://dx.doi.org/10.1016/j.cemconres.2018.10.022.

[17] S. Bahafid, S. Ghabezloo, M. Duc, P. Faure, J. Sulem, Effect of the hydration temperature on the microstructure of class G cement: C-S-H composition and density, Cem. Concr. Res. 95 (2017) 270–281, http://dx.doi.org/10.1016/j.cemconres.2017.02.008.

[18] S. Bahafid, S. Ghabezloo, P. Faure, M. Duc, J. Sulem, Effect of the hydration temperature on the pore structure of cement paste: Experimental investigation and micromechanical modelling, Cem. Concr. Res. 111 (September 2017) (2018) 1–14, http://dx.doi.org/10.1016/j.cemconres.2018.06.014.

[19] S. Ghabezloo, Association of macroscopic laboratory testing and micromechanics modelling for the evaluation of the poroelastic parameters of a hardened cement paste, Cem. Concr. Res. 40 (8) (2010) 1197–1210, http://dx.doi.org/10.1016/j.cemconres.2010.03.016.

[20] S. Ghabezloo, Multiscale modeling of the poroelastic properties of various oil-well cement pastes, J. Multiscale Model. 2 (3–4) (2010) 199–215, http://dx.doi.org/10.1142/S1756973710000412.

[21] S. Ghabezloo, Effect of the variations of clinker composition on the poroelastic properties of hardened class G cement paste, Cem. Concr. Res. 41 (8) (2011) 920–922, http://dx.doi.org/10.1016/j.cemconres.2011.03.022.

[22] S. Ghabezloo, Micromechanics analysis of thermal expansion and thermal pressurization of a hardened cement paste, Cem. Concr. Res. 41 (5) (2011) 520–532, http://dx.doi.org/10.1016/j.cemconres.2011.01.023.

[23] N. Agofack, S. Ghabezloo, J. Sulem, Chemo-poro-elastoplastic modelling of an oilwell cement paste: Macroscopic shrinkage and stress-strain behaviour, Cem. Concr. Res. 132 (2020) 106046, http://dx.doi.org/10.1016/j.cemconres.2020.106046.

[24] M. Samudio, Modelling of an oil well cement paste from early age to hardened state: hydration kinetics and poromechanical behaviour (Ph.D. thesis), Université Paris-Est, 2017.

[25] M. Pierre, S. Ghabezloo, P. Dangla, R. Mesnil, M. Vandamme, J.-F. Caron, Multi-physics modelling for extrusion-based 3D-printing: Material, process and applications, in: RILEM International Conference on Concrete and Digital Fabrication, Springer, 2024, pp. 449–456, http://dx.doi.org/10.1007/978-3-031-70031-6_52.

[26] M. Pierre, M. Samudio, S. Ghabezloo, P. Dangla, Modelling the poromechanical behaviour of class G cement paste: A multiphysics approach from early age to hardened state, Cem. Concr. Res. 193 (2025) 107852, http://dx.doi.org/10.1016/j.cemconres.2025.107852.

[27] O. Coussy, Poromechanics, John Wiley & Sons, 2004.

[28] O. Coussy, Mechanics and physics of porous solids, John Wiley & Sons, 2011.

[29] D. Min, T. Mingshu, Formation and expansion of ettringite crystals, Cem. Concr. Res. 24 (1) (1994) 119–126, http://dx.doi.org/10.1016/0008-8846(94)90092-2.
[30] N. Roussel, Rheological requirements for printable concretes, Cem. Concr. Res. 112 (2018) 76-85, http://dx.doi.org/10.1016/j.cemconres.2018.04.005.

[31] B. Albert, B. Guy, D. Damidot, Water chemical potential: A key parameter to determine the thermodynamic stability of some hydrated cement phases in concrete? Cem. Concr. Res. 36 (5) (2006) 783-790, http://dx.doi.org/10.1016/j.cemconres.2005.12.016.

[32] R.J. Flatt, G.W. Scherer, J.W. Bullard, Why alite stops hydrating below 80% relative humidity, Cem. Concr. Res. 41 (9) (2011) 987-992, http://dx.doi.org/10.1016/j.cemconres.2011.06.001.

[33] M. Wyrzykowski, P. Lura, Effect of relative humidity decrease due to self-desiccation on the hydration kinetics of cement, Cem. Concr. Res. 85 (2016) 75-81, http://dx.doi.org/10.1016/j.cemconres.2016.04.003.

[34] M. Cervera, J. Oliver, T. Prato, Thermo-chemo-mechanical model for concrete. I: Hydration and aging, J. Eng. Mech. 125 (1999) 1018-1027, http://dx.doi.org/10.1061/(ASCE)0733-9399(1999)125:9(1018).

[35] M.T. Van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Am. J. 44 (5) (1980) 892-898, http://dx.doi.org/10.2136/ssaj1980.0361599500440050002x.

[36] Y. Mualem, A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res. 12 (3) (1976) 513-522, http://dx.doi.org/10.1029/WR012i003p00513.

[37] P.C. Carman, Fluid flow through a granular bed, Trans. Inst. Chem. Eng. Lond. 15 (1937) 150-156.

[38] D.W. Taylor, Fundamentals of soil mechanics, vol. 66, (2) LWW, 1948.

[39] V. Picandet, D. Rangeard, A. Perrot, T. Lecompte, Permeability measurement of fresh cement paste, Cem. Concr. Res. 41 (3) (2011) 330-338, http://dx.doi.org/10.1016/j.cemconres.2010.11.019.

[40] N. Ducoulombier, R. Mesnil, P. Carneau, L. Demont, H. Bessaies-Bey, J.-F. Caron, N. Roussel, The “slugs-test” for extrusion-based additive manufacturing: Protocol, analysis and practical limits, Cem. Concr. Compos. 121 (2021) 104074, http://dx.doi.org/10.1016/j.cemconcomp.2021.104074.

[41] L. Demont, R. Mesnil, N. Ducoulombier, J.-F. Caron, Affordable inline struc- turation measurements of printable mortar with a pocket shear vane, Constr. Build. Mater. 408 (2023) 133602, http://dx.doi.org/10.1016/j.conbuildmat.2023.133602.

[42] X. Liu, B. Ma, H. Tan, B. Gu, T. Zhang, P. Chen, H. Li, J. Mei, Effect of aluminum sulfate on the hydration of Portland cement, tricalcium silicate and tricalcium aluminate, Constr. Build. Mater. 232 (2020) 117179, http://dx.doi.org/10.1016/j.conbuildmat.2019.117179.

[43] A. Boumiz, C. Vernet, F.C. Tenoudji, Mechanical properties of cement pastes and mortars at early ages: Evolution with time and degree of hydration, Adv. Cem. Based Mater. 3 (3-4) (1996) 94-106, http://dx.doi.org/10.1016/S1065-7355(96)00402-5.

[44] M.-H. Vu, Effet des contraintes et de la température sur l'intégrité des ciments des puits pétroliers (Ph.D. thesis), Université Paris-Est, 2012, p. 240.

[45] L. Sun, X. Pang, S. Ghabezloo, H. Wang, J. Sun, Hydration kinetics and strength retrogression mechanism of silica-cement systems in the temperature range of 110 °C-200 °C, Cem. Concr. Res. 167 (2023) 107120, http://dx.doi.org/10.1016/j.cemconres.2023.107120.

[46] R. Wolfs, F. Bos, T. Salet, Early age mechanical behaviour of 3D printed concrete: Numerical modelling and experimental testing, Cem. Concr. Res. 106 (2018) 103-116, http://dx.doi.org/10.1016/j.cemconres.2018.02.001.

[47] H. Kada, M. Lachemi, N. Petrov, O. Bonneau, P.C. Aïtin, Determination of the coefficient of thermal expansion of high performance concrete from initial setting, Mater. Struct. 35 (2002) 35-41, http://dx.doi.org/10.1007/BF02482088.

[48] K. Velez, S. Maximilien, D. Damidot, G. Fantozzi, F. Sorrentino, Determination by nanoindentation of elastic modulus and hardness of pure constituents of Portland cement clinker, Cem. Concr. Res. 31 (4) (2001) 555-561, http://dx.doi.org/10.1016/S0008-8846(00)00505-6.

[49] T. Honorio, B. Bary, F. Benboudjema, Thermal properties of cement-based materials: Multiscale estimations at early-age, Cem. Concr. Compos. 87 (2018) 205-219, http://dx.doi.org/10.1016/j.cemconcomp.2018.01.003.

[50] D.P. Bentz, P.E. Stutzman, F. Zunino, Low-temperature curing strength enhance- ment in cement-based materials containing limestone powder, Mater. Struct. 50 (2017) 1-14, http://dx.doi.org/10.1617/s11527-017-1042-6.

[51] L. D'aloia, G. Chanvillard, Determining the "apparent" activation energy of concrete: Ea—numerical simulations of the heat of hydration of cement, Cem. Concr. Res. 32 (8) (2002) 1277-1289, http://dx.doi.org/10.1016/S0008-8846(02)00791-3.

[52] P. Mounanga, A. Khelidj, G. Bastian, Experimental study and modelling ap- proaches for the thermal conductivity evolution of hydrating cement paste, Adv. Cem. Res. 16 (3) (2004) 95-103, http://dx.doi.org/10.1680/adrc.2004.16.3.95.

[53] S. Prapaharan, A. Altschaeffl, B. Dempsey, Moisture curve of compacted clay: mercury intrusion method, J. Geotech. Eng. 111 (9) (1985) 1139-1143, http://dx.doi.org/10.1061/(ASCE)0733-9410(1985)111:9(1139).

[54] P.J. Ross, K.R. Smettem, Describing soil hydraulic properties with sums of simple functions, Soil Sci. Am. J. 57 (1) (1993) 26-29, http://dx.doi.org/10.2136/ssaj1993.03615995005700010006x.

[55] W. Durner, Hydraulic conductivity estimation for soils with heterogeneous pore structure, Water Resour. Res. 30 (2) (1994) 211-223, http://dx.doi.org/10.1029/93WR02676.

[56] K.L. Scrivener, P. Juilland, P.J. Monteiro, Advances in understanding hydration of portland cement, Cem. Concr. Res. 78 (2015) 38-56, http://dx.doi.org/10.1016/j.cemconres.2015.05.025.

[57] T. Ooms, G. Vantyghem, R. Van Coile, W. De Corte, A parametric modelling strategy for the numerical simulation of 3D concrete printing with complex ge- ometries, Addit. Manuf. 38 (2021) 101743, http://dx.doi.org/10.1016/j.addma.2020.101743.

[58] R.J. Wolfs, T.A. Salet, N. Roussel, Filament geometry control in extrusion-based additive manufacturing of concrete: The good, the bad and the ugly, Cem. Concr. Res. 150 (2021) 106615, http://dx.doi.org/10.1016/j.cemconres.2021.106615.

[59] J. Reinold, V.N. Nerella, V. Mechtcherine, G. Meschke, Extrusion process simu- lation and layer shape prediction during 3D-concrete-printing using the particle finite element method, Autom. Constr. 136 (2022) 104173, http://dx.doi.org/10.1016/j.autcon.2022.104173.

[60] G. Rizzieri, L. Ferrara, M. Cremonesi, Numerical simulation of the extrusion and layer deposition processes in 3D concrete printing with the particle finite element method, Comput. Mech. (2023) 1-19, http://dx.doi.org/10.1007/s00466-023-02367-y.

[61] Build'in, a regional technological platform for construction, 2016, (Accessed 12 July 2016). https://www.buildin-enpc.fr/.

[62] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities, Internat. J. Numer. Methods Engrg. 79 (11) (2009) 1309-1331, http://dx.doi.org/10.1002/nme.2579.

[63] G. Vantyghem, T. Ooms, W. De Corte, VoxelPrint: A grasshopper plug-in for voxel-based numerical simulation of concrete printing, Autom. Constr. 122 (2021) 103469, http://dx.doi.org/10.1016/j.autcon.2020.103469.

[64] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M.E. Rognes, G.N. Wells, The FEniCS project version 1.5, Arch. Numer. Softw. 3 (100) (2015).

[65] T. Helfer, B. Michel, J.-M. Proix, M. Salvo, J. Sercombe, M. Casella, Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element mod-elling platform, Comput. Math. Appl. 70 (5) (2015) 994-1023, http://dx.doi.org/10.1016/j.camwa.2015.06.027.

[66] P. Carneau, R. Mesnil, O. Baverel, N. Roussel, Layer pressing in concrete extrusion-based 3D-printing: Experiments and analysis, Cem. Concr. Res. 155 (2022) 106741, http://dx.doi.org/10.1016/j.cemconres.2022.106741.

[67] N. Ganganath, C.-T. Cheng, K.-Y. Fok, K.T. Chi, Trajectory planning for 3D printing: A revisit to traveling salesman problem, in: 2016 2nd International Conference on Control, Automation and Robotics, ICCAR, IEEE, 2016, pp. 287-290, http://dx.doi.org/10.1109/ICCAR.2016.7486742.

[68] A.V. Shembekar, Y.J. Yoon, A. Kanyuck, S.K. Gupta, Trajectory planning for conformal 3D printing using non-planar layers, in: International Design Engi- neering Technical Conferences and Computers and Information in Engineering Conference, 51722, American Society of Mechanical Engineers, 2018, http://dx.doi.org/10.1115/DETC2018-85975, V01AT02A026.

[69] R.D. Weeks, R.L. Truby, S.G. Uzel, J.A. Lewis, Embedded 3D printing of multimaterial polymer lattices via graph-based print path planning, Adv. Mater. 35 (5) (2023) 2206958, http://dx.doi.org/10.1002/adma.202206958.

[70] P. Carneau, R. Mesnil, N. Roussel, O. Baverel, Additive manufacturing of cantilever-from masonry to concrete 3D printing, Autom. Constr. 116 (2020) 103184, http://dx.doi.org/10.1016/j.autcon.2020.103184.