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

# Thickness optimisation in 3D printed concrete structures

Romain Mesnil ID*, Pedro Sarkis Rosa, Léo Demont ID

École Nationale des Ponts et Chaussées, Institut Polytechnique de Paris, Marne-La-Vallée, France

## ARTICLE INFO

**Keywords:**

Concrete printing

3DCP

Shell structures

Optimisation

Thickness optimisation

Buckling

## ABSTRACT

Layer pressing in 3D concrete printing (3DCP) allows to continuously modify the thickness of printed laces by changing adequately the robot speed. However, most applications consider a constant thickness throughout the printing and do not leverage all the possibilities from robotic technologies. The aim of this paper is to demonstrate the potential offered by thickness variation to achieve higher structural efficiency and to lower the material usage. To do so, analytical solutions for stress and buckling of tapered heavy column are recalled and highlight a potential of reduction of 25% of material for simple geometries with materials with low structuration rate. Numerical optimisation based on a penalty method and on the finite element simulation with shell elements is then implemented to minimise the volume of printed components for more complex geometries. Promising results are observed and should encourage the 3DCP community to further study this previously unexplored dimension of the process.

## 1. Introduction

### 1.1. Context

3D Concrete Printing (3DCP) is one of the focus points of research on construction automation and is considered as one of the main enablers of construction 4.0 [1]. The possibility to create more complex geometries is often seen as a path towards reduction of environmental footprint without compromising performance [2]. Increasing efficiency of concrete structures is an important topic, since concrete is the most widely produced synthetic material worldwide and contributes to a consequent share of anthropic greenhouse gases emissions [3].

Extrusion is by far the most popular process in 3D Concrete Printing (3DCP). This method, while mature, limits design freedom compared to powder bed printing, where the artifact is supported by material during printing. A key challenge is ensuring structural stability during printing due to stresses from self-weight [4]. This challenge makes simulation of 3DCP a growing research area, essential for advancing the technology [5]. The construction industry's interest in 3DCP is rising, but current scientific progress in modelling and simulation does not meet industry demands [6]. As numerical models become increasingly predominant over analytical models [5], these tools need optimisation for computational efficiency and ease of use to be widely adopted [7]. Real-time feedback and control are crucial for process optimisation [7] but are not yet achievable with existing models.

The 3DCP process faces challenges at various scales [6]: mixing and pumping, extrusion and filament, and the structure or printed element.

Industry demands advances in structural-scale models able to simulate structures during printing, ensuring printability and form finding. Such model must consider time-evolution of the printed material as well as build rate, and layer interaction [5]. The model must predict buckling and plastic collapse, as well as deformation. This would enable the process to be optimised, for example, by defining the optimum printing speed or the minimum structuration rate without the occurrence of failure [7]. In addition, shape optimisation - and hence material savings - would also be possible.

Analytical methods, which provide quick insights, are gaining interest, but are often limited to simple shapes [7]. For example, [8] proposed rheological requirements for printable concrete in terms of layer stability, as well as geometry control, collective buckling failure, and plastic cracking. Furthermore, [9] has developed a parametric model to assess the failure of 3D-printed wall structures in terms of plastic collapse and elastic buckling considering the evolution of material properties. In addition, a rheological model was developed by [10] to quantify the constructability of a printable mortar, considering material strength, construction rate, and failure modes based on a Mohr-Coulomb criterion.

Numerical modelling is required for the evaluation of more complex structures, especially for non-standard geometries. In the reviews carried out by [5,6], two numerical models stand out for their ability to simulate complex structures: [11,12] (the latter was further developed in two subsequent publications [13,14]). Both studies employ volumetric finite element modelling (FEM) using Abaqus software,

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: romain.mesnil@enpc.fr (R. Mesnil).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. 3DCP finite element modelling presented in the literature.</figcaption></figure>

employing eight-node solid elements. The temporal nature of the model is implemented by dividing the calculation into a sequence of steps, with progressive activation of elements at each step (“model change” option in Abaqus), made possible by using a layered mesh. Once activated, time-dependent material variables are assigned to the elements to account for the evolution of the property gradient during printing, using Mohr-Coulomb criteria. In [12], two different models are proposed, which are called CobraPrint and VoxelPrint. Both are based on the model previously proposed by [15], who first demonstrated the application of numerical methods to control structural behaviour during printing [7]. The primary distinction is in how the 3D model is discretised; whereas VoxelPrint uses voxelisation to approximate the 3D shape of the object, CobraPrint builds its 3D models from a prior definition of the tool path, similar to the model in [11]. Fig. 1 illustrates these three models.

This approach may not be scalable for models of components with several hundred layers. A convergence study performed in [16] shows that such models requires hundreds of thousands of degrees of freedom even for simple structures in order to capture the influence of lace geometry. 3D FEM tends to have high computational costs, especially if geometric non-linearity for buckling and material non-linearity for plasticity have to be considered [7]. There is thus room for alternative, lighter, numerical methods for the structural modelling of the 3DCP process. One such alternative is the construction of an equivalent plate model under the assumption of plane stresses ($\sigma_{13} = \sigma_{23} = \sigma_{33} = 0$). A significant advantage of plate model is the notion of generalised forces, which provide a more intuitive tool for visualising and understanding the real phenomena observed.

Modelling 3DCP has already been carried out using plate elements [16]. By calculating the equivalent Kirchhoff-Love plate elastic stiffness tensors as a function of the 3D elastic properties of the unit cell (in this example, the lace), the authors were able to quantify the effect of the local geometry of the lace on the overall anisotropy of the printed object. The results were conclusive, with a ratio of up to 150% between the bending stiffness of the different plate directions, demonstrating the ability of FEM modelling using plate elements to capture phenomena of a three-dimensional nature.

Another critical requirement for fully developing 3DCP is to integrate the numerical modelling process into a comprehensive design protocol. Five factors make up the 3DCP framework [5]: print path, fresh-state properties, mix design, cured-state properties, and design optimisation. These factors are interdependent, i.e., modifying any of them can lead to modification of the others. Although no one has yet succeeded in condensing all these factors into a single and comprehensive numerical model, we must consider optimising communication and feedback between the individual models. The models mentioned above capable of simulating complex geometries also offer their proposals for the 3DCP design flow [11–14,17]. There is scope for further development, given that, in all these cases, the proposal involves two or more independent software. This makes establishing the necessary communication flow for the process optimisation routines more complex.

Despite the exponential growth of industrial developments, structural optimisation considering process constraints, including stability during printing seems to be out of the scope of research at the moment, as most analysis tools are leaning towards computationally expensive approaches. Substantial material savings is however a necessity in order

<table><caption>**Table 1**

Printing parameters. Bold letters highlight parameters that can be independently chosen. Other parameters are considered as process outputs.</caption><thead><tr><th>Category</th><th>Parameter</th><th>Notation</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="2">Kinematics</td><td>Nozzle velocity</td><td>V<sub>r</sub></td><td>mm/s</td></tr><tr><td>Extrusion velocity</td><td>V<sub>e</sub></td><td>mm/s</td></tr><tr><td rowspan="3">Geometric</td><td>Nozzle diameter</td><td>D</td><td>mm</td></tr><tr><td>Layer height</td><td>h</td><td>mm</td></tr><tr><td>Layer width</td><td>B</td><td>mm</td></tr><tr><td rowspan="5">Process</td><td>Initial yield stress</td><td>τ<sub>0</sub></td><td>Pa</td></tr><tr><td>Flow rate</td><td>Q</td><td>L/s</td></tr><tr><td>Initial shear modulus</td><td>G<sub>0</sub></td><td>Pa</td></tr><tr><td>Extrusion force</td><td>F</td><td>N</td></tr><tr><td>Inter-layer time</td><td>Δt</td><td>s</td></tr><tr><td rowspan="3">Material</td><td>Structuration rate</td><td>A<sub>thix</sub></td><td>Pa/s</td></tr><tr><td>Shear stiffening rate</td><td>G<sub>thix</sub></td><td>Pa/s<sup>2</sup></td></tr><tr><td>Density</td><td>ρ<sub>g</sub></td><td>N/m<sup>3</sup></td></tr></tbody></table>

to mitigate the environmental impact of the 3D printing process itself and make it competitive with conventional construction process [18]. Finding the limits of economy of material in the 3DCP by extrusion process is thus an important question that underlies all efforts to reduce its environmental footprint. This article makes a modest contribution in this direction by exploring the potential of varying thickness throughout printing.

## 1.2. Problem statement

The knowledge on 3DCP process increases steadily, and more and more research is aiming at the simulation of the process. Yet, the topic of process optimisation is rarely addressed. The primary objective of this study is to explore the potential of process optimisation to achieve economy of material in 3DCP. Three questions are addressed:

* Is it possible to link process parameters to fast finite element simulation to allow for optimisation of process parameters?

* Is there potential to vary thickness during printing to achieve material savings without loss of performance?

* Which properties should be optimised?

The second section recalls fundamental results that allow to link geometrical parameters of the lace to process parameters, such as extrusion and nozzle velocity. The third section presents a simplified framework for evaluation of stability of 3D-printed components, opening the way towards optimisation. The optimisation procedure is presented in the fourth section. The results of the optimisation are presented in the fifth Section: they show that the manufacturer has to make a compromise between economical performance and environmental impact, and that the proposed optimisation framework is adapted to help him or her to make this decision.

# 2. Identification of process constraints and parameters

## 2.1. Process parametrisation

3D printing by extrusion is a complex task that involves material, kinematics and process parameters. By varying those parameters, different ways of printing the same object can be devised, but with different levels of efficiency. In this study, the parametrisation proposed by [19] is employed, summarised in Table 1. It is essential to note that the reference addresses the layer pressing 3DCP strategy, which is also adopted here. To optimise the process, one can employ the integration of a parametric model in tandem with optimisation algorithms, as will be illustrated.
In reality, it is not possible to choose all parameters independently, it is thus important to differentiate input parameters from output parameters. This work has already been carried out in [19]. This article builds upon this work and adds several assumptions.

It is assumed that structuration rate and shear stiffening rate are constant throughout the process. Although this assumption is legitimate for 1K printing strategy, it has been shown in some studies on 2K strategies that the rate of adjuvantation significantly modifies the structuration rate [20,21]. The assumption made in this paper is thus equivalent to an assumption of constant rate of adjuvantation.

We assume an oblong shape for the lace [19,22], it is possible to derive the relation from Eq. (1) between nozzle velocity, lace width B and height h.

$$V_r = \frac{Q}{h \left( B + \left( \frac{\pi}{4} - 1 \right) h \right)} \quad (1)$$

It is assumed that the nozzle velocity remains inferior to extrusion speed, which is the oriented lace pressing regime: $V_r < V_e$. It is also assumed that the lace height h remains inferior to the nozzle diameter D

$$\left\{ \begin{array}{l} V_r < \frac{4Q}{\pi D^2} \\ \frac{h}{D} < 1 \end{array} \right. \quad (2)$$

The robot velocity is assumed to vary during printing, that is to say, it is a function of time. The robot toolpath consists of a series of positions and velocities of the tool center point: it is possible to vary this velocity during the process, for example by using dedicated software [23]. All other input process parameters (Q, D, h) are supposed to be constant with time. This assumption will be discussed in the last section of the article. The result of these assumptions and of Eq. (1), is that the lace thickness B can be changed by changing robot velocity, according to Eq. (3).

$$B(t) = \frac{Q}{hV_r} + \left(1 - \frac{\pi}{4}\right) h \quad (3)$$

Additionally, the combination of (2) and (3) shows that the lace thickness B remains superior to D. The set of constraints of Eq. (2) can thus be rewritten as in Eq. (4).

$$\left\{ \begin{array}{l} B > D \\ \frac{h}{D} < 1 \end{array} \right. \quad (4)$$

## 2.2. Practical bounds on process parameters

In practice, the overpressing of a layer may result in a local failure, as shown in [19]. Carneau et al. have proposed an analytical formula recalled in Eq. (5) for the pressing force at the nozzle. It has been backed by a series of experiments.

$$F = \frac{\tau_0 D B^2}{2h} \quad (5)$$

A necessary condition for stability is that the stress resultant from the pressing and the self-weight of the last and penultimate layer does not exceed the yield stress of the material of the penultimate layer:

$$\frac{4F}{\pi D^2} + 2\rho gh \le \tau_0 + A_{thix} \cdot \Delta t \quad (6)$$

Using Eq. (5), an upper for the ratio B/D can be found with (7), after introduction of the dimensionless structuration rate $\beta = \frac{A_{thix} \Delta t}{\tau_0}$ and of the dimensionless lace weight $\gamma = \frac{\rho g h}{\tau_0}$.

$$1 < \frac{B}{D} < \sqrt{\frac{\pi h}{D}} \cdot \sqrt{1 + \beta - 2\gamma} \quad (7)$$

This shows that materials with higher structuration rate $\beta$ allow for higher contrast between layer thickness. Therefore, the thickness variation strategy is only compatible with materials with relatively high structuration, typically found in 2K systems.

**Numerical application:** Demont et al. measured $A_{thix}$ in the range of 90Pa/s for a bicomponent strategy [21]. Considering an inter-layer time of 30 s, an initial yield stress of 1000 Pa, a layer height of 8 mm and a nozzle diameter of 20 mm, the application of Eq. (7) results in practical bounds for the lace thickness B.

$$1 < \frac{B}{D} < 2 \quad (8)$$

## 2.3. Numerical modelling

As discussed in the introduction, analytical models, although necessary to understand the scaling of a given physical problem, are usually constrained to simple solutions and are not sufficient in practice. Therefore, it is necessary to implement a numerical model that can be used in a broader context. A large part of this article is devoted to the introduction of a simplified modelling strategy that can be used together with classical optimisation procedures.

# 3. Numerical modelling of 3D printing

## 3.1. Model assumptions

The model developed in this study is based on the following simplifying assumptions:

**Plate model.** the printed wall thickness is assumed to be much smaller than the other characteristic dimensions. The plane stress model is adopted. This assumption is not verified in the immediate vicinity of the nozzle, since the layer pressing results in transverse strains.

**Elastic linear simulation.** if the critical stress is exceeded, global instability is assumed, as there are few means of redistributing the stresses.

**Elastic linear buckling.** it can be dominant on slender structures [24]. A safety factor is used to compensate for second-order effects.

**Material time evolution model.** The shear threshold $\tau_c(t)$ is defined as a linear function of the material's structuration rate $A_{thix}$ [21,25]. The time $t$ since deposition of the material (so that $t = 0$ occurs at the nozzle level) is also called age of the material. The evolution of yield stress is written as follows:

$$\tau_c(t) = \tau_{c0} + A_{thix} t$$

The shear modulus $G(t)$ is defined analogously to $\tau_c$ but employs a quadratic function:

$$G(t) = G_0 + G_{thix} t^2$$

This behavioural law was obtained from preliminary rheometer tests and requires further investigation. However, some theoretical results on thixotropy allow to capture the phenomenology of this evolution, so that quadratic evolution of storage modulus can be viewed as a consequence of the linearity assumption for yield stress evolution with time. Indeed, Chateau et al. derived a relationship between yield stress and shear modulus for suspensions with monodisperse rigid and noncolloidal particles [26]. This relation is shown in Eq. (9):

$$\frac{\tau_c(\phi)}{\tau_c(0)} = \sqrt{\frac{(1-\phi)G(\phi)}{G(0)}} \quad (9)$$

The use of the accelerator leads to the nucleation of ettringite needles, which create a crystalline structure [27]. Even a slight change in solid fraction may lead to drastic change in mechanical properties. Therefore, Eq. (9) can be rewritten by considering that the solid volume
fraction remains close to the solid volume fraction $\phi_0$ before addition of accelerator.

$$ \frac{\tau_c(\phi)}{\tau_c(\phi_0)} = \sqrt{\frac{1-\phi}{1-\phi_0} \cdot \frac{G(\phi)}{G(\phi_0)}} \sim \sqrt{\frac{G(\phi)}{G(\phi_0)}} \quad (10) $$

Thus, the linearity assumption for yield stress results in a quadratic evolution of shear modulus with time.

$$ \left(1 + \frac{A_{thix}}{\tau_{c,0}} \cdot t\right)^2 G_0 = G(t) \quad (11) $$

It should be remembered that this estimation is valid under simplifying assumptions. Eq. (11) allows to understand the fundamental scaling of time evolution of stiffness. Other models can be used to describe the evolution of shear modulus with time, depending on the physical or chemical mechanism at stake in the appearant structuration. For example, Wolfs and Salet have observed a linear evolution of Young's modulus through time for 1K material [15]. Motamedi et al. have postulated exponential increase in stiffness for 3D printed earth structures based on soil mechanics models [28]. The choice of evolution law for stiffness and strength should thus be backed by experiments.

**Von Mises failure criterion:** the yield strength is considered to be $\tau_c$, according to the von Mises criterion. This failure criterion is often used for fresh mortar [29]. Other criteria considering effects of consolidation and difference between tension and compression could be used.

**Layer-pressing strategy:** although some of the model's features can be applied to other printing strategies, such as “infinite brick extrusion”, this model is based on the layer-pressing strategy, as presented in [19].

## 3.2. Model construction and computation

The principle for evaluating load-bearing capacity of a 3D printed-model during printing is illustrated in Fig. 2, shown below. It consists of 7 steps: the first two require input from the user, while the other are fully automated.

1. A surface representing the object to be printed is supplied by the user.

2. The user then defines the number N of subdivisions of the surface, corresponding to the number of horizontal layers into which the object will be divided.

3. A layer-wise mesh is created according to the N subdivisions.

4. N independent models are created, each representing a printing step (a snapshot during the printing process).

5. For each model (step), the average time of deposition of each mesh face is calculated as a function of geometry and printing parameters.

6. Material properties $\tau_c(t)$, $G(t)$, and $E(t)$ are assigned to each mesh slice based on the times of deposition.

7. Finite element calculations of N independent models are performed under self-weight load.

This algorithm will hereafter be referred to as “step-by-step” and will serve as the foundational algorithm upon which the optimisation routine will be built. In practical terms, users have the flexibility to specify which of the N steps they wish to compute. Furthermore, the number of subdivisions N is not necessarily equal to the actual number of layers in the physical item. For each selected step, two key outcomes are obtained: the buckling factor $\lambda_{cr}$ and the Von Mises stress utilisation rate $\eta$, derived through independent calculations (second-order equilibrium and linear buckling analysis). Both results obtained can be interpreted as the safety factors for each step, in terms of plastic collapse and elastic buckling. If one wishes to assess the safety of the print during the entire 3DCP process, one can instruct the program to calculate all the steps.

The time of deposition of the nth layer is determined by computing the volume $V$ of material deposited. This volume is a function of the

layer heights $H_n$, layer length $\ell$ and vertical robot velocity $V_z$, shown in Fig. 3.

$$ t_{deposition}^{n} = \frac{V}{Q} = \frac{1}{Q} \sum_{i=1}^{n} V_z(z_i) H_i \quad (12) $$

The vertical velocity is related to nozzle velocity by the equation $V_z = h/\ell V_r$, where $\ell$ is the average length of the layer. The average length is computed by averaging the total boundary length for each lace (blue curves in Fig. 3). The nozzle velocity is computed as a function of layer height h and target thickness B with Eq. (1). For a given step n, the material age $t_k^n$ of the kth layer is found by subtracting the time of depositions.

$$ t_k^n = t_{deposition}^n - t_{deposition}^k \quad (13) $$

The finite element software used is Karamba, a FEA library interoperable with parametric design environment [30]. The triangular TRIC shell elements are used [31]. Although second-order equilibrium and elastic buckling calculations are performed for each of the N independent models, the suggested modelling is still appealing from the perspective of computational cost. Indeed, compared to a volumetric element calculation with geometric and material non-linearity, the linear elastic plate model is considerably lighter. Whereas nonlinear volumetric modelling can take several minutes (or hours) to converge, each step in the proposed model is solved within seconds. Let us use Fig. 2's shape as an example, when modelled with about 4,000 elements.

## 3.3. Convergence analysis: analytical considerations

When the number of strips N equals the actual number of layers in the printed structures, the Finite Element Model contains all the information about gradient of material properties. However, a large-scale structure may consist of hundreds, if not thousands of layers. As a result, the exact model becomes very computationally intensive. Selecting the appropriate subdivision that remains accurate while keeping computation time as low as possible is thus an important preliminary step. The model's convergence with regard to mesh fineness is explored from two distinct perspectives: (i) the gradient of concrete properties, and (ii) finite element calculations.

Fig. 4 illustrates the modelling assumption and introduces some notations required for an analytical study. The mesh is superimposed with the actual structure, and average material properties are selected based on the estimated age of the element. The variation of age across an element of height $H_{element}$ is defined by Eq. (14).

$$ \nabla t = \frac{H_{element} \cdot \Delta t}{h} \quad (14) $$

For the sake of simplicity, ratio $h/\Delta t$ is replaced by $V_z$ the vertical velocity of the robot. The variations of yield stress across the element are thus written in Eq. (15).

$$ \nabla \tau_c = \frac{H_{element} \cdot A_{thix}}{V_z} \quad (15) $$

The relative variation of yield stress is thus

$$ \epsilon = \frac{\nabla \tau_c}{\tau_c(H)} = \frac{1}{\frac{\tau_0 V_z}{H_{el} A_{thix}} + \frac{(H-z)}{H_{el}}} \quad (16) $$

A necessary condition for the accuracy of the model is that the relative variation remains small: $\epsilon \ll 1$. Two asymptotic cases can be identified:

* $\frac{\tau_0 V_z}{H_{el} A_{thix}} \gg 1$: the material has a build-up across the element that remains small compared to the initial yield stress; In that case, the equation $\epsilon \ll 1$ is always verified

* $\frac{\tau_0 V_z}{H_{el} A_{thix}} < 1$: the equality $\epsilon \ll 1$ is verified only when $(H-z)/H_{el}$, i.e. when the element is sufficiently far away from the last deposited lace.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Construction of *N* independent finite element models representing the *N* calculation steps.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Notations for determining the age of layers in the FEM. The colour gradient corresponds to the age of the material: the lighter, the more recent the deposition.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Notations for the convergence analysis. Darker shades of grey correspond to older material.</figcaption></figure>

This study suggests that a non-uniform meshing could be adopted to limit the meshing error. However, the contribution of the last deposited layers to the overall strain energy is negligible once enough layers have been deposited.

### 3.4. Convergence analysis: numerical validation

The analytical analysis can be validated numerically on a case-study close to the practical considerations (see Fig. 5). In order to study the convergence of finite element calculations, two global factors were selected: the buckling factor $\lambda$ and the potential energy $E_p$. The results are then calculated for the entire printing process, i.e. all model

<figure><img src="image_6.png" /><figcaption>Fig. 5. Different mesh fineness levels.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Buckling factor as a function of mesh density (hourglass).</figcaption></figure>

calculation steps, and plotted against the z-coordinate of the print at each step (Figs. 6 and 7). The results indicate that a mesh density of 1,100 elements per square meter is adequate for modelling the two objects under study. Beyond this threshold, the computational cost escalates disproportionately to the gains in accuracy.

Notably, an exception to the general convergence observed occurs within the first 10 cm of the print, where substantial variability among the different meshes is found in both cases. This variability could stem from the previously mentioned properties gradient error, or from the validity of plate model hypothesis in this region, where the structure is not yet slender. This confirms the analytical considerations that show that the finite element should be refined near the newly deposited layers, in order to accommodate for the gradient in material properties.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Strain energy as a function of mesh density (hourglass).</figcaption></figure>

### 3.5. Conclusion on numerical modelling

This section proposed a modelling strategy for the printing process. Unlike most existing approaches adopted in the literature, this modelling approach focuses on speed, with the use of shell elements coupled with linear elasticity. From a practical perspective, we advocate for this simplified approach rather than plasticity, since plastic failure introduces major defects in practice, and that the post-failure behaviour of structures is not relevant for manufacturers. Additionally, flow rule for plasticity has yet to be measured on 3D printed mortars, and the accuracy of more complex models has yet to be demonstrated with physical experiments. This section has also shown that the mesh size can be decoupled from layer size without significant loss of accuracy, except at the nozzle level. This provides great leverage for speeding up computation, with few thousands of degree of freedom for a standard column, taking a few seconds of computation at most.

With this limitation in mind, it is possible to use the proposed numerical method to vary virtually process parameters. The next section explores this possibility by varying process parameters to control the layer thickness, with the ultimate goal of minimising the mass of 3D-printed samples.

## 4. Mass and print time optimisation routine

### 4.1. Objective function

One optimisation algorithm application is to improve the 3DCP process by reducing material consumption and printing time. Both objectives are interesting from an economic and environmental point of view. For a given geometry, material consumption can be reduced by varying the thickness B of the layers. As a result, areas where there is a concentration of stress could benefit from larger layers, whereas low-load areas may employ thinner layers. This respects the paradigm of additive manufacturing, for which the material is applied only where it is needed. On the other hand, reducing layer thickness also reduces printing time since the total volume V is reduced, for a fixed flow rate value Q. Hence, another way to reduce printing time is to increase Q for a fixed volume V. However, the greatest mass reduction that can be achieved without failure will not occur at the same time as the greatest printing time reduction. This is so that finer structures may be printed because delayed printing increases the material's resistance to failure as a result of its structuration in time. Therefore, rather of referring to a single optimal solution, we should consider a set of them, identified in a form of a Pareto front.

The volume optimisation problem is explicated in Eq. (17). This is a constrained optimisation problem: at any time step t, the structure has

to satisfy the ULS requirements for stress and linear buckling.

minimise $V(Q, B)$

subject to

$$ 
\begin{array}{l}
\forall t, \lambda(t) > \lambda_{cr} \\
\forall t, \sigma_{cr}(t) < \sigma_{d}(t) \\
B_{min} < B < B_{max}
\end{array} 
\qquad (17) 
$$

In Eq. (17), $\sigma_d$ is the design maximal allowable Von Mises stress, defined as the product between the measured yield strength and a material safety factor taken arbitrarily equal to $\gamma = 1.3$.

$$ \sigma_d = \frac{\tau_c}{\gamma} \qquad (18) $$

and $\sigma_{cr}(t)$ is the critical stress, taken in the following as the 99th percentile for Von Mises stress, in order to avoid instabilities due to local stress concentrations. A more refined approach involving limit analysis of the shell could be employed as an alternative to an elastic-plastic modelling [32,33]. However, it is assumed here that elastic modelling provides the best compromise between accuracy and speed.

The optimisation problem from Eq. (17) is in fact a robust optimisation problem:

minimise $V(Q, B)$

subject to

$$ 
\begin{array}{l}
\min_t \lambda(t) > \lambda_{cr} \\
\max_t \sigma_{cr}(t) < \sigma_{d}(t) \\
B_{min} < B < B_{max}
\end{array} 
\qquad (19) 
$$

The optimisation problem can be relaxed thanks to a penalty method. The penalty function $g(x) = \max(0, x)^2$ is used, and the problem can be re-written as follows:

$$ 
\begin{aligned}
\min_{Q,B} V(\mathbf{B}) + \mu \cdot & \left( g \left( \lambda_{cr} - \min_t \lambda \right) + g \left( \max_t \sigma_{cr} - \sigma_d \right) \right. \\
& \left. + g \left( \max_t B - B_{max} \right) + g \left( B_{min} - \min_t B \right) \right) 
\end{aligned} 
\qquad (20) 
$$

It will be shown in the next section that the constraint $B_{min} < B < B_{max}$ can be enforced by linear box constraints on the parameters describing B. While not strictly equivalent to problem (19), Eq. (21) can be used as a relaxed optimisation problem.

$$ \min_{Q,B} V(\mathbf{B}) + \mu \cdot \left( g \left( \lambda_{cr} - \min_t \lambda \right) + g \left( \max_t \sigma_{cr} - \sigma_d \right) \right) \qquad (21) $$

The printing time optimisation problem can be formulated in similar fashion by noticing that the printing time T and total volume are related by $T = V/Q$. The rest of this paper only treats the penalised problem from Eq. (21).

### 4.2. Definition of parameters

Defining the optimisation variables is the first step into implementing the mass and print time optimisation routine. Layer thickness B and flow rate Q are two parameters of the 3DCP process that have an effect on these objectives, as previously shown. The flow rate depends on the pumping system, here a range of 0.6–1.5 L/min is chosen. For the sake of simplicity, it is assumed that the flow rate does not vary throughout the process, therefore Q is a scalar value. On contrary, the layer thickness can be varied by changing the robot speed and it assumed that B is a function. With the same objective of simplification, it is assumed that B is constant along any given lace, but can vary from one lace to the other. If the shape is sliced according to altitude, the consequence is that B is only a function of the altitude z.

In the followings, the variations of thickness are here parametrised by a Bézier spline, as shown in Fig. 8. In the image, the curve is a function from the z domain to the B-domain. The user can define a number of parameters (5 on the Figure) corresponding to the B value
<figure><img src="image_9.png" /><figcaption>Fig. 8. Parametrisation of thickness with control points.</figcaption></figure>

of the control points, whose z-coordinates are evenly spread between the extremal values $B_{min}$ and $B_{max}$. Due to the strong convex hull property of NURBS, this means that the constraint $B_{min} < B < B_{max}$ is automatically satisfied. Therefore, the optimisation problem defined in Eq. (21) is solved. The vector **B** contains the **B** values of the control points and therefore unequivocally encodes the thickness distribution. This CAD-based approach is well-established in the literature on shell optimisation [34].

The parameter space is therefore formally defined by the tuple (Q, B), but the parameter Q does not need to vary for the volume optimisation problem as we shall see in Appendix. Indeed, the optimal value for volume is always reached for the minimal value of flow rate $Q_{min}$.

Volume and printing time can be evaluated from the parameters (Q, B). The remaining task is thus to evaluate the ULS requirements. As highlighted in (19), this involves the solution of a secondary optimisation problem. The problem being a discrete problem (there is a finite number of layers), it is possible to find the worst printing stage by brute force. However, this would require hundreds of evaluations and would be too costly in practice. Alternatively, one can leverage the continuous nature of the layer thickness B, and opt to only analyse predefined steps. Additionally, the failure indicators $\lambda(z)$ and $\tau_c(z)$ also tends to exhibit continuity throughout the calculations. Hence, by choosing only a few steps to be evaluated by optimisation iteration, the computation cost can be significantly reduced, while the continuity of the problem tends to assure the remaining steps will also respect the given constraints. In the end of the optimisation phase, the user can then verify all steps for the optimal solutions selected, to assure their validity.

### 4.3. Choice of optimisation algorithm

Several optimisation procedures can be used to solve the minimisation problem. The reader should be aware that many local optima exist, and that gradient descent algorithms are not guaranteed to converge towards the global optimum. One possible alternative is to perform successive refinements of the NURBS parametrised by the **B** vector. Another alternative, chosen here, was motivated by the low number of variables: first a global optimisation was performed during 15 min using the DIViding RECTangles algorithm [35]. Then, a local optimisation was performed starting from the new configuration. We chose the COBYLA algorithm, a derivative-free method for the local optimisation step [36]. The Grasshopper plug-in GOAT, which is based on a C# wrapper of the library NLOpt was used [37]. The implementation could be vastly improved through parallel computing of sequencing steps. Indeed, the linearity assumption makes it possible to compute each step independently. The computation time for a single step takes 250 ms on a laptop with Intel(R) Core(TM) i7-8650 CPU processor with 32 GB of RAM.

The whole optimisation procedure is fully automated and does not require interaction with the designer. It is thought of as a routine that could be run before printing a batch of pieces.

<table><caption>**Table 2**
Fixed parameters during optimisation.</caption><thead><tr><th>Type</th><th>Data</th><th>Value</th></tr></thead><tbody><tr><td>Geometry</td><td>Surface</td><td>Hourglass and Barrel</td></tr><tr><td rowspan="6">Material</td><td>ρ</td><td>23 kN/m<sup>3</sup></td></tr><tr><td>ν</td><td>0.3</td></tr><tr><td>τ<sub>0</sub></td><td>1 kPa</td></tr><tr><td>A<sub>thix</sub></td><td>70 Pa/s</td></tr><tr><td>G<sub>0</sub></td><td>190 kPa</td></tr><tr><td>G<sub>thix</sub></td><td>4.5 Pa/s<sup>2</sup></td></tr><tr><td rowspan="4">Process</td><td>H</td><td>6 mm</td></tr><tr><td>β</td><td>0.8</td></tr><tr><td>D</td><td>15</td></tr><tr><td>Q<sup>*</sup></td><td>0.6-1.5 L/min</td></tr></tbody></table>

### 4.4. Design workflow

The described optimisation approach opens the door to envisioning a comprehensive design workflow for printable elements (see Fig. 9). The process starts with the choice of geometry, material and process parameters. A first optional step is to perform a sequencing of the whole printing to identify potential weak points in the structure. The result of this operation is the selection of few steps where buckling and plastic failure are evaluated. As discussed, this avoids the solution of an auxiliary optimisation problem for each value of (B, Q). Another alternative is to select by default a certain number of evaluation steps. In this paper, 5 steps were chosen. Then, the constrained optimisation routine is performed. The result satisfies the buckling and plastic failure criteria for the predefined steps, but not necessarily for all steps. Therefore, a second run of “Step by step” evaluation is performed on the optimal solution. If buckling or plastic failure is observed, there are two alternatives:

* Perform a local optimisation by running all the steps. This approach is valid only if the solution is close to be feasible;

* Add evaluation steps to the optimisation routine and re-run the procedure.

If no buckling or plastic failure is observed, then the solution is printable, and it is possible to export the toolpath. The workflow was implemented as a Grasshopper plug-in, which allows interoperability with CAD environment (Mesh or NURBS modelling), but also Finite Element Analysis or CAM software [23]. Interoperability between software is of paramount importance in free-form projects: this allows for smooth iterations between designer and manufacturers [38,39].

### 4.5. Geometries and model parameters

Two distinct shapes were selected to assess the value of the proposed modelling and optimisation - namely, the “hourglass” and “barrel” shapes (as depicted in Fig. 10). These particular geometries were chosen for their distinctive attributes, combining straightforward geometric definitions with intricate printing complexities, making them ideal subjects for the study.

The hourglass and barrel geometries are surfaces of revolution, with a vertical axis of revolution. They share common dimensions, with a height of 1.0 meter and a base diameter of 40 centimeters. The equation for the meridian $\rho(z)$ follows:

$$ \begin{cases} \rho_{\text{Barrel}}(z) = 0.4 + 0.25 \sin^2 \pi z \\ \rho_{\text{Hourglass}}(z) = 0.4 - 0.25 \sin^2 \pi z \end{cases} \qquad (22) $$

The reference case for the model parameters is shown in Table 2. It was defined on the basis of the material and setup available for this study. It should be noted that **B** and **Q** are the variables of the optimisation process yet to be performed.

The modelling methodology having been introduced, the last section discusses results of the optimisation.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Proposed design flow for process optimisation. Optional steps are dashed.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Hourglass (left) and barrel (right) geometries.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Optimal mass as a function of printing time for the barrel, with constant and variable thickness. All solutions for the constant thickness are dominated by at least one solution of the variable thickness.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Optimal mass as a function of printing time for the barrel, with constant and variable thickness. All solutions for the constant thickness are dominated by at least one solution of the variable thickness.</figcaption></figure>

# 5. Results and discussion

## 5.1. Results for the barrel geometry

Fig. 11 shows the compromise curve between printing time and mass for the barrel shown in Fig. 10. Each point plotted in the Figure is the result of a mass optimisation at prescribed flow rate $Q$, with $Q$ between 0.6L/min and 1.5L/min, and other parameters recalled in Table 2. The results are plotted in the (printing time, mass) space: the lower the point, the lighter the structure, and the more to the left, the faster the printing.

The constant thickness strategy corresponds to the case where the shell thickness is constant throughout the printing. Formally, this corresponds to the case $B = (B, B)$ with $B_{min} < B < B_{max}$. The variable thickness strategy is described by a vector $B \in \mathbb{R}^7$.

The constant thickness strategy shows an interesting trend: the optimal weight varies greatly depending on the flow rate, and thus the printing speed. The minimal mass is 199 kg, while the minimal mass for $Q = 1.5$ L/min is 296 kg. The manufacturer can make a compromise between material efficiency and printing time, and thus cost.

It can be seen that both curves create a Pareto front, where no solution dominates strictly the other. However, the Pareto front for the constant thickness is entirely dominated by the Pareto front for the variable thickness. Consider for example the grey rectangle drawn on the image, which contains three blue points. In other terms, three points are both lighter and take less time to print.

The Pareto front for the variable thickness features an abrupt change of slope at (88 min, 224 kg), corresponding to a flow rate of 1.1 L/min. This corresponds to the point where one of the variables $B$ reaches the maximal value $B_{max}$. For lower flow rates, the variable thickness uses 22% less material for the same printing time, and takes 28% less time for the same mass. For this particular example, there seems to be a strong advantage in varying shell thickness throughout printing.

## 5.2. Results for the hourglass geometry

Fig. 12 shows the Pareto front for the hourglass shown in Fig. 10. Variable thickness seems to be even more beneficial in the case of

hourglass. Some points of the “constant thickness” Pareto front are dominated by all the points from the “variable thickness” Pareto Front. The constant thickness structure uses 49% more material when printed in the same timespan.

The main difference between the two shapes is that the narrowing of the hourglass at mid-height implies a much shorter inter-layer time. In this region, the buckling load factor for shells of constant thickness is much lower than any other printing stage, as shown in Fig. 13. This kind of behaviour may be expected for shells with special features.
<figure><img src="image_14.png" /><figcaption>Fig. 13. Evolution of buckling load factor with printed height for the hourglass.</figcaption></figure>

### 5.3. Post-processing of results: compromise between productivity and environmental footprint

As seen previously, the variable thickness methodology allows significant savings in printing time and material usage. However, it is in general not possible to optimise for both objectives at the same time, since reducing printing speed results in higher structuration and improved stability of the component. This trade-off can be viewed as a compromise between economical and environmental performance. As a simplified approach, we estimate the economical cost $c_{eco}$ of a printed piece by a linear function:

$$c_{eco} = c_{cell} \cdot t + m \cdot c_{concrete} \quad (23)$$

The notation $c_{cell}$ corresponds to the operating cost of the cell in €/min and $c_{concrete}$ is the cost of concrete in €/kg.

Although Roux et al. have shown that cement alone is not sufficient to explain the environmental impact of 3D-printed concrete [18,40], it remains the main contributor to climate change of the process. The environmental cost $c_{environment}$ can thus be approximated in a crude way by comparing the binder content, which is found as the product between mass and binder intensity $b_l$ of the mix in kg of cement per kg of concrete.

$$c_{environment} \sim m \cdot b_l \quad (24)$$

Fig. 14 shows the visual representation of compromise between environmental footprint and cost minimisation. Iso-values for total binder content are horizontal lines, while iso-values for cost are inclined lines, with a slope of $-\frac{c_{cell}}{c_{concrete}}$. For the sake of representation, following values were chosen:

* The cost of operation of the cell was chosen to be equal to the average cost of labor for two people in France, assuming that two people are needed to operate a cell. $c_{cell} = 1.28€/min$

* The cost of concrete was estimated to $c_{concrete} = 0.45€/kg$.

* The cement content was estimated to 0.160 kg cement/kg concrete based on the mix B3 from [41].

The cost ratio can be modified by the user depending on the specificity of the printing system, material and labour cost. In any case, the iso-values are straight lines and it is relatively simple to identify Pareto front.

With the proposed assumptions, the gain in cost and cement content for the barrel and the hourglass are summed up in Table 3. For the hourglass, the Pareto front of the constant thickness is entirely dominated by the whole Pareto front of the variable thickness. This means that any solution resulting from variable thickness optimisation is both

<figure><img src="image_15.png" /><figcaption>Fig. 14. Compromise between cost and cement usage for the barrel. Thin grey lines correspond to iso-cost of manufacturing, dotted blue curves correspond to iso-values for cement content.</figcaption></figure>

cheaper and lighter. A subset of the optimal variable thickness subset dominates entirely the constant thickness optimal results. For example, the lightest barrel with variable thickness is both less expensive than the cheapest barrel with constant thickness and lighter than the lightest barrel with constant thickness.

It can be inferred that thickness optimisation should be investigated for highly irregular shapes, which are the main domain of application of 3D printing.

### 5.4. Role of thickness variation

Thickness variation seems to provide a significant advantage for the reduction of material. The numerical experiments highlight another phenomenon, already discussed in [4]: the change of interlayer-time through thickness variation. Indeed thickness variation is the result of a change of robot velocity as recalled in the second section of this article: for a given flow rate, increasing shell thickness also decreases the loading rate of the structure, and increases the inter-layer time.

Consider the optimisation of a dome, with base diameter of 2 m, a height of 1 m, truncated at 0.95 m to prevent collapse of the last layers. This geometry is a scale-up version of the one studied in [4]. For this example, we chose $\lambda_d = 2$. The thickness distribution resulting from the optimisation is shown in Fig. 15. It appears that the optimal distribution tends to increase the thickness at the top of the dome. This may seem counter-intuitive at first glance, since stresses are larger at the bottom of the dome. To better understand this, it is necessary to look at the time evolution of the printed material.

Fig. 16 shows the inter-layer time, assuming a layer height $h = 10$ mm. The y-axis of the graph corresponds to the altitude, and the x-axis of the graph to inter-layer time. For the constant thickness solution, the interlayer time is much larger for low altitude. However, the optimal result with variable thickness has a more even distribution of inter-layer time. For both cases, the inter-layer time converges to 36 s at the crown. The maximal inter-layer time is due to the lower length of laces, as already observed in [4].

Therefore, the role of thickness variation can be understood in two manners: as a way to increase stiffness locally, and as a way to control inter-layer time.

### 5.5. Performance of optimisation

As explained in Section 4.2, the chosen optimisation scheme starts by a global optimisation for 15 min, and refines it with a local optimisation usually carried in less than 3 min. This was done to maximise the
<table><caption>**Table 3**
Comparison of best designs for the constant and variable thickness strategy.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Cheapest solution</th><th colspan="2">Lightest solution</th></tr><tr><th>Cost</th><th>Cement</th><th>Cost</th><th>Cement</th></tr></thead><tbody><tr><td rowspan="2">Hourglass</td><td>Constant thickness</td><td>157 €</td><td>25.1 kg</td><td>170 €</td><td>19.9 kg</td></tr><tr><td>Variable thickness</td><td>129€(-18%)</td><td>18.8 kg (-25%)</td><td>137€(-20%)</td><td>15.9 kg (-20%)</td></tr><tr><td rowspan="2">Barrel</td><td>Constant thickness</td><td>243€</td><td>44.4 kg</td><td>274€</td><td>32 kg</td></tr><tr><td>Variable thickness</td><td>214 €(-12%)</td><td>36 kg (-19%)</td><td>238€(-13%)</td><td>27.8 kg (-13%)</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Visualisation of thickness gradient resulting from the proposed thickness optimisation for a dome of diameter 2 m and height 1 m: side view (top) and top view (bottom).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 17.** Mass resulting from the optimisation with COBYLA for the hourglass with $Q = 1.0L/min$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Inter-layer time for the constant and variable thickness.</figcaption></figure>

chance to converge to a global optimum of the objective function. Likewise, the number of 7 parameters for the thickness vector **B** was chosen as a reference, since further refinement does not provide significant improvement of the objective function for the treated examples. Fig. 17 displays the objective function of the hourglass for several dimensions of **B** and computation time for an optimisation with the COBYLA algorithm. The optimal value found in the previous section is represented with a thick horizontal line. Notice that the computation with 5 dimensions converges to this optimal value. With sufficient computational time, all methods outperform the single parameter optimisation, which corresponds to the constant thickness case.

This short study shows that the performance of the optimisation procedure depends on the choice of the dimensions of the vector **B**. Satisfactory results, within 3% of the best found optimum could be achieved within 5 min. Unfortunately, there is no general rule to select the appropriate parameter space, which depends on the complexity of

the input geometry. The use of iterative refinement of the **B** vector, in the spirit of [42], could result in much faster and robust convergence.

Finally, all the examples were treated with 5 printing steps, which were evaluated sequentially in each evaluation step. The last step accounts for 32% of the computational time for the linear buckling analysis. We may thus expect a factor 3 speed-up if parallel computation is enabled, although this would require to leave the Grasshopper environment.

# 6. Experimental validation

## 6.1. *Experimental protocole*

The proposed modelling methodology has been applied to the modelling of the 3D printing of components with complex geometries. The experimental equipment consists of a 6-axis robotic arm equipped with a 2K printhead developed by XTreeE. The nozzle velocity, mortar flow rate and accelerator flow rate are prescribed through RobotStudio, The printing material is similar to the B1 mix developed by de Bono et al. [41].

## 6.2. *Material parameters evaluation*

Most of the material parameters can be evaluated with tests that are now widely adopted by the 3DCP community. The initial yield stress $\tau_0$ was measured with the slug test, developed by [43]. The initial yield stress is not highly sensitive to the flow rate of accelerator, and is in the range of $\tau_0$ = 1 kPa. The slug test is also used to calibrate the flow rat, as the density of the mortar is known, with $\rho$ = 2100 kg/m$^3$. The structuration rate $A_{thix}$ was measured with the pocket shear vane introduced by Demont et al. [21]. The value for $A_{thix}$ varies between 40 Pa/s and 100 Pa/s depending on the rate of accelerator. The structuration rate is extremely high compared to the initial yield stress: at 40 Pa/s, the contribution of structuration becomes dominant after 25 s, which is close to the typical interlayer time of the printed object. Therefore, the mechanical behaviour is much more dependent to the structuration rate than to the initial yield stress.
<figure><img src="image_19.png" /><figcaption>Fig. 18. Post-buckling failure for a regular corolla and a stiffened one. Notice the orthoradial crack at mid-height for the stiffened corolla.</figcaption></figure>

The evaluation of the Young's modulus of printable materials with a single test is still an open research question. Indeed, the drastic evolution of material properties of 2K materials through time makes it difficult to apply conventional testing, such as uniaxial compression tests. The storage modulus has been evaluated with oscillatory rheometry, but it is necessary to validate this with experiments with the printed material. In the following section, an experimental technique is proposed to calibrate the $G_{thix}$. Indeed, the initial Young modulus of printable material seems to be within the range of 100 kPa, as tested by [15].

### 6.3. Full-scale buckling testing

A buckling experiment has been conducted on printed structures to have an indirect validation for the choice of Young's modulus of the printed material. It is indeed well-known that elastic buckling depends on the bending stiffness of the structure, and is independent of plastic yield stress. It is also well-known that elastic buckling is favoured by high slenderness ratio [9]. Therefore, it is possible to evaluate the bending stiffness, and thus the Young's modulus of a structure through the monitoring of the buckling load factor.

The chosen geometry is an upside down corolla with the intent to favour elastic buckling rather than plastic failure. The structure is in compression only during the printing. The outer diameter is of 1625 mm, while the inner diameter is of 525 mm. The wall thickness is constant, equal to 25 mm, while the layer height $H$ is set to 8 mm. Stiffeners can be incorporated in order to increase the buckling capacity of the shell. The cohesion between the printed layer and the printing bed was ensured by using sandpaper. The printing led to buckling failure at the 8th layer for the regular corolla (printing time of 600 s) and at the 15th layer for the reinforced corolla (printing time of 1800 s). Interestingly, for the stiffened shell, material failure occurs at mid-height. This material failure is not the main cause for failure, but rather a result of excessive deformation after initiation of buckling (see Fig. 18).

The buckling load is evaluated with the proposed Finite Element Modelling scheme. The pressing force is evaluated with Eq. (5), and is equal to 0.5N. The properties obtained with oscillatory rheometry of $G_0 = 190$ kPa and $G_{thix} = 4.45$Pa.s$^{-2}$ were used as material parameters for the Finite Element Model. The results of finite element simulations for various layers are shown in Fig. 19. The load buckling factor at failure is of 3.0 and 2.7 respectively for the regular and stiffened shells. This is in line with other experimental results for buckling of thin shell structures [44]. Notice also that buckling load factors for previous layers are higher.

The numerical model is thus in qualitative and quantitative agreement with the experiments. The testing procedure seems thus to provide satisfactory results and validates the main modelling assumptions.

<figure><img src="image_20.png" /><figcaption>Fig. 19. Linear load buckling factor $\lambda_{cr}$ with respect to layer index resulting from the proposed FEA simulation. The vertical coloured lines correspond to the measured data. Both result in similar buckling load factor.</figcaption></figure>

However, it should be reminded that buckling is heavily influenced by boundary conditions and imperfections. At the hardened state, imperfections are taken into account through the adoption of specific safety factors, and the designer can work within a well-established framework. However, choosing the appropriate safety coefficient for limiting buckling during printing is a topic that requires further exploration.

### 6.4. Thickness optimisation of a shell structure

Finally, an optimisation of shell thickness was performed to ensure printability of a complex object, shown in Fig. 20. The critical load buckling factor was chosen as $\lambda_{cr} = 5.0$. The optimisation results in a maximal layer thickness of 35 mm, while the maximal layer thickness printable with the system is of 45 mm. The total material savings between the safest printing option of constant thickness of 45 mm, which would have been adopted in an industrial context, and the optimised result exceeds 30%, illustrating potential savings offered by the customisation and optimisation of thickness in 3D printed structures. This experimentation shows that process simulation and thickness optimisation can be used in practice. However, the true limits of economy of materials with process optimisation are yet to be found. Indeed, the choice of the critical buckling factor was made with perhaps excessive precaution, and further material savings could be made with a highly accurate prediction of buckling loads including loads and material imperfections.
<figure><img src="image_21.png" /><figcaption>Fig. 20. Successful printing of a complex object with customised thickness.</figcaption></figure>

# 7. Conclusions and future research

## 7.1. Summary of intellectual contributions

This research was motivated by three questions, introduced in the end of the first Section.

**Is it possible to link process parameters to fast finite element simulation to allow for optimisation of process parameters?** Section 2 has shown that the published literature allows to draw conclusions on practical bounds for layer pressing, and that the pressing force can be evaluated with a simple analytical formula. A simplified modelling and thickness optimisation framework for stability evaluation in 3D printed structures, based on plate formulation has been proposed. The use of parametric software, such as Grasshopper, allows to automatically compute material age for any given shape and to evaluate the buckling or plastic failure at any given stage.

**Is there potential to vary thickness during printing to achieve material savings without loss of performance?** This article has discussed the potential of thickness optimisation in 3D concrete printing. Numerical case studies, with reasonable time evolution law for the printed material show potential for material savings from 13% to 25% compared to shells with constant thickness without compromising mechanical performance during printing.

**Which properties should be optimised?** Finally, we showed that, under the assumption of a constant flow rate, the optimisation for printing time and mass are equivalent. As a result, for a given flow rate, a single value optimisation provides the best compromise between economical and environmental performance. Thus, the Pareto front can be plotted by performing independent single-valued optimisation runs. An interpretation of the (time, mass) in terms of (cost, global warming potential) was proposed and resulted in a simplified way to visualise compromise between economical and environmental performance of 3D printed concrete structures.

## 7.2. Limitations and future work

One of the main unknowns is the confidence in the buckling analysis and the choice of a reasonable safety factor for linear buckling analysis $\lambda_d$. The first challenge is the evaluation of storage modulus at fresh state. There is no universally adopted methodology to evaluate it at early stage of printing. Some authors use uniaxial tests, while other evaluate a dynamic modulus from ultrasounds [29]. Oscillatory rheometry was used and we proposed arguments to explain why a quadratic evolution of storage modulus was observed in the case of ettringitic reaction. In 2K systems, it is expected that the adjuvantation rate will influence the evolution of storage modulus through time: this variable is left aside from the present study and could be integrated once more reliable measurement of $G_{thix}$ are made. The other unknown is the choice of the safety factor for linear buckling analysis $\lambda_d$, which is still an active topic of research. The calibration with respect to large-scale buckling measurements is a possible way, but at the cost of high material consumption. The use of probabilistic modelling to take into account for material and geometrical imperfections seems promising, although the heavy computational time makes this approach difficult to combine with optimisation [45]. Perhaps the most reasonable solution is to perform a sensitivity analysis to imperfection before the optimisation procedure.

With a potential of reduction up to 25%, the thickness optimisation seems to be a promising strategy that should be further explored in future research. Indeed, thickness variation does not influence the geometry of the printed component and can be leveraged by manufacturer with minimal impact on the project. Thickness optimisation can also be combined with shape optimisation as shown in [34]: the potential for mass reduction can thus also be explored by architects and structural engineers. Finally, the optimisation principles can also be extended to an optimisation considering both printability and load bearing capacity of the hardened component.

### CRediT authorship contribution statement

**Romain Mesnil:** Writing – review & editing, Writing – original draft, Validation, Supervision, Methodology, Investigation, Funding acquisition, Conceptualization. **Pedro Sarkis Rosa:** Writing – review & editing, Writing – original draft, Validation, Software, Investigation. **Léo Demont:** Writing – review & editing, Validation, Investigation.

### 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.

### Acknowledgement

This work was carried with partial funding from the French State represented Caisse des Dépôts et Consignations in the context of the Beam-O research project, itself part of Construire au Futur, Habiter le Futur project lead by Région Île-de-France. The authors thank Nicolas Ducoulombier for providing data about storage modulus at early stage. We also thank our colleague Jean-François Caron as well as the reviewers for valuable discussions.

## Appendix A. Optimisation results

Tables A.7, A.6, A.4 and A.5 contains the parameters of the optimal configurations found by optimisation. The values of linear buckling load factor $\lambda_i$ are also displayed. Here, $\lambda_i$ is evaluated at a height of 0.2(1 + i) meters. The costs for material and fabrication are calculated with the assumptions proposed in the Section 6.

It can be noticed that all $\lambda_i$ are close to 5 for the variable thickness optimisation, whereas only $\lambda_2$ is equal to 5 for constant thickness optimisation.
<table><caption>**Table A.4**

Results for the variable thickness optimisation of the hourglass.</caption><tr><td>Q (L/min)</td><td>0.6</td><td>0.65</td><td>0.7</td><td>0.75</td><td>0.8</td><td>0.85</td><td>0.9</td><td>0.95</td><td>1</td></tr><tr><td>B<sub>0</sub> (mm)</td><td>15</td><td>15</td><td>15</td><td>15</td><td>15</td><td>15</td><td>15</td><td>15</td><td>15</td></tr><tr><td>B<sub>1</sub> (mm)</td><td>22.2</td><td>23</td><td>24</td><td>25.2</td><td>26.4</td><td>27.6</td><td>29.8</td><td>32.6</td><td>35</td></tr><tr><td>B<sub>2</sub> (mm)</td><td>31.4</td><td>32.8</td><td>34</td><td>34.2</td><td>34.2</td><td>34.6</td><td>35</td><td>35</td><td>35</td></tr><tr><td>B<sub>3</sub> (mm)</td><td>21.6</td><td>23</td><td>25</td><td>26.4</td><td>27.8</td><td>28.8</td><td>30.6</td><td>32.8</td><td>35</td></tr><tr><td>B<sub>4</sub> (mm)</td><td>21</td><td>21.8</td><td>22.8</td><td>23.6</td><td>24.4</td><td>25</td><td>26</td><td>27.4</td><td>28.6</td></tr><tr><td>B<sub>5</sub> (mm)</td><td>19.6</td><td>20</td><td>20.4</td><td>21.4</td><td>22.4</td><td>22.8</td><td>24.2</td><td>25</td><td>26.4</td></tr><tr><td>B<sub>6</sub> (mm)</td><td>15</td><td>15.8</td><td>16.4</td><td>16.8</td><td>17.2</td><td>17.8</td><td>18</td><td>18.4</td><td>18.6</td></tr><tr><td>mass (kg)</td><td>99.1</td><td>103</td><td>107.2</td><td>110.6</td><td>113.9</td><td>116.7</td><td>121.7</td><td>126.9</td><td>132.1</td></tr><tr><td>Printing time (min)</td><td>71.8</td><td>68.9</td><td>66.6</td><td>64.1</td><td>61.9</td><td>59.7</td><td>58.8</td><td>58.1</td><td>57.5</td></tr><tr><td>λ<sub>1</sub> (-)</td><td>5.1</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5.1</td><td>5.6</td><td>6.3</td><td>6.8</td></tr><tr><td>λ<sub>2</sub> (-)</td><td>5.5</td><td>5.6</td><td>5.8</td><td>5.5</td><td>5.2</td><td>5</td><td>5</td><td>5</td><td>5</td></tr><tr><td>λ<sub>3</sub> (-)</td><td>5</td><td>5.3</td><td>5.8</td><td>6</td><td>6.3</td><td>6.3</td><td>6.8</td><td>7.6</td><td>8.4</td></tr><tr><td>λ<sub>4</sub> (-)</td><td>5.1</td><td>5</td><td>5</td><td>5.1</td><td>5.2</td><td>5.1</td><td>5.4</td><td>5.8</td><td>6.2</td></tr><tr><td>λ<sub>5</sub> (-)</td><td>5</td><td>5.1</td><td>5.1</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td></tr><tr><td>Operational cost (€)</td><td>92</td><td>88</td><td>85</td><td>82</td><td>79</td><td>76</td><td>75</td><td>75</td><td>74</td></tr><tr><td>Material cost (€)</td><td>45</td><td>46</td><td>48</td><td>50</td><td>51</td><td>53</td><td>55</td><td>57</td><td>59</td></tr><tr><td>Total cost (€)</td><td>137</td><td>134</td><td>133</td><td>132</td><td>131</td><td>129</td><td>130</td><td>132</td><td>133</td></tr></table>

<table><caption>**Table A.5**

Results for the constant thickness optimisation of the hourglass.</caption><tr><td>Q (L/min)</td><td>0.6</td><td>0.65</td><td>0.7</td><td>0.75</td><td>0.8</td><td>0.85</td><td>0.9</td><td>0.95</td><td>1</td></tr><tr><td>B (mm)</td><td>27.4</td><td>28.4</td><td>29.2</td><td>30.2</td><td>31.2</td><td>32</td><td>32.8</td><td>33.8</td><td>34.6</td></tr><tr><td>mass (kg)</td><td>123.8</td><td>128.3</td><td>132.3</td><td>136.4</td><td>140.9</td><td>144.6</td><td>148.7</td><td>152.7</td><td>156.3</td></tr><tr><td>Printing time (min)</td><td>89.7</td><td>85.8</td><td>82.2</td><td>79.1</td><td>76.6</td><td>74</td><td>71.8</td><td>69.9</td><td>68</td></tr><tr><td>λ<sub>1</sub> (-)</td><td>11.8</td><td>11.9</td><td>11.7</td><td>12</td><td>12.3</td><td>12.3</td><td>12.3</td><td>12.8</td><td>12.9</td></tr><tr><td>λ<sub>2</sub> (-)</td><td>5</td><td>5</td><td>4.9</td><td>5</td><td>5</td><td>5</td><td>4.9</td><td>5</td><td>5</td></tr><tr><td>λ<sub>3</sub> (-)</td><td>10.8</td><td>10.7</td><td>10.3</td><td>10.3</td><td>10.3</td><td>10.1</td><td>10</td><td>10</td><td>9.9</td></tr><tr><td>λ<sub>4</sub> (-)</td><td>14.1</td><td>13.8</td><td>13.3</td><td>13.1</td><td>13.1</td><td>12.8</td><td>12.5</td><td>12.6</td><td>12.4</td></tr><tr><td>λ<sub>5</sub> (-)</td><td>31.7</td><td>31</td><td>29.7</td><td>29.4</td><td>29.1</td><td>28.4</td><td>27.8</td><td>27.8</td><td>27.3</td></tr><tr><td>Operational cost (€)</td><td>115</td><td>110</td><td>105</td><td>101</td><td>98</td><td>95</td><td>92</td><td>89</td><td>87</td></tr><tr><td>Material cost (€)</td><td>56</td><td>58</td><td>60</td><td>61</td><td>63</td><td>65</td><td>67</td><td>69</td><td>70</td></tr><tr><td>Total cost (€)</td><td>170</td><td>168</td><td>165</td><td>163</td><td>161</td><td>160</td><td>159</td><td>158</td><td>157</td></tr></table>

<table><caption>**Table A.6**

Results for the variable thickness optimisation of the barrel.</caption><tr><td>Q (L/min)</td><td>0.6</td><td>0.7</td><td>0.8</td><td>0.9</td><td>1</td><td>1.1</td><td>1.2</td><td>1.3</td><td>1.4</td><td>1.5</td></tr><tr><td>B<sub>0</sub> (mm)</td><td>15</td><td>16</td><td>17.8</td><td>18</td><td>18.2</td><td>18.4</td><td>21</td><td>21</td><td>21.8</td><td>23.4</td></tr><tr><td>B<sub>1</sub> (mm)</td><td>15</td><td>17</td><td>18.8</td><td>20.8</td><td>22.6</td><td>24.4</td><td>24</td><td>23.6</td><td>23.4</td><td>23</td></tr><tr><td>B<sub>2</sub> (mm)</td><td>18.8</td><td>18.8</td><td>18.6</td><td>18.4</td><td>18.2</td><td>18</td><td>21.2</td><td>25.6</td><td>29.2</td><td>33</td></tr><tr><td>B<sub>3</sub> (mm)</td><td>20.6</td><td>22</td><td>23.4</td><td>24.8</td><td>26.2</td><td>27.6</td><td>28.2</td><td>29.2</td><td>31.2</td><td>33</td></tr><tr><td>B<sub>4</sub> (mm)</td><td>24.8</td><td>26.6</td><td>28.2</td><td>30</td><td>31.8</td><td>33.6</td><td>33.8</td><td>34.2</td><td>34.4</td><td>34.6</td></tr><tr><td>B<sub>5</sub> (mm)</td><td>20.4</td><td>22</td><td>23.4</td><td>24.6</td><td>25.8</td><td>26.6</td><td>28.8</td><td>31</td><td>33</td><td>35</td></tr><tr><td>B<sub>6</sub> (mm)</td><td>17.6</td><td>20.4</td><td>23.2</td><td>25.2</td><td>27</td><td>28.8</td><td>31.6</td><td>33</td><td>34.2</td><td>35</td></tr><tr><td>mass (kg)</td><td>172.9</td><td>184.3</td><td>195</td><td>204.9</td><td>214.6</td><td>223.9</td><td>235.6</td><td>249.3</td><td>263.1</td><td>276.8</td></tr><tr><td>Printing time (min)</td><td>125</td><td>114</td><td>106</td><td>99</td><td>93</td><td>88</td><td>85</td><td>83</td><td>82</td><td>80</td></tr><tr><td>λ<sub>1</sub> (-)</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td></tr><tr><td>λ<sub>2</sub> (-)</td><td>12.3</td><td>10.1</td><td>8.3</td><td>7.1</td><td>6.2</td><td>5.5</td><td>6.8</td><td>9</td><td>10.7</td><td>12.7</td></tr><tr><td>λ<sub>3</sub> (-)</td><td>23.4</td><td>22.1</td><td>21.1</td><td>20.3</td><td>19.5</td><td>19</td><td>18.5</td><td>19.3</td><td>21.2</td><td>22.9</td></tr><tr><td>λ<sub>4</sub> (-)</td><td>5</td><td>5</td><td>5.1</td><td>5.2</td><td>5.3</td><td>5.4</td><td>5</td><td>5</td><td>5.1</td><td>5.1</td></tr><tr><td>λ<sub>5</sub> (-)</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td></tr><tr><td>Operational cost (€)</td><td>160</td><td>147</td><td>136</td><td>127</td><td>119</td><td>113</td><td>109</td><td>107</td><td>105</td><td>103</td></tr><tr><td>Material cost (€)</td><td>78</td><td>83</td><td>88</td><td>92</td><td>97</td><td>101</td><td>106</td><td>112</td><td>118</td><td>125</td></tr><tr><td>Total cost (€)</td><td>238</td><td>229</td><td>223</td><td>219</td><td>216</td><td>214</td><td>215</td><td>219</td><td>223</td><td>227</td></tr></table>

<table><caption>**Table A.7**

Results for the constant thickness optimisation of the barrel.</caption><tr><td>Q (L/min)</td><td>0.6</td><td>0.7</td><td>0.8</td><td>0.9</td><td>1</td><td>1.1</td><td>1.2</td><td>1.3</td><td>1.4</td><td>1.5</td></tr><tr><td>B</td><td>23</td><td>24.6</td><td>26</td><td>27.4</td><td>28.6</td><td>29.8</td><td>31</td><td>32</td><td>33.2</td><td>34.4</td></tr><tr><td>mass (kg)</td><td>199.2</td><td>212.6</td><td>225.3</td><td>236.7</td><td>247.7</td><td>257.8</td><td>267.6</td><td>276.6</td><td>286.6</td><td>296.9</td></tr><tr><td>Printing time (min)</td><td>144.4</td><td>132.1</td><td>122.5</td><td>114.3</td><td>107.7</td><td>101.9</td><td>96.9</td><td>92.5</td><td>89</td><td>86.1</td></tr><tr><td>λ<sub>1</sub> (-)</td><td>16.2</td><td>15.6</td><td>14.9</td><td>14.6</td><td>14</td><td>13.7</td><td>13.4</td><td>13</td><td>13</td><td>13</td></tr><tr><td>λ<sub>2</sub> (-)</td><td>24</td><td>22.3</td><td>21</td><td>20.3</td><td>19.3</td><td>18.6</td><td>18.1</td><td>17.4</td><td>17.2</td><td>17.1</td></tr><tr><td>λ<sub>3</sub> (-)</td><td>29.3</td><td>28.7</td><td>27.8</td><td>27.2</td><td>26.2</td><td>25.5</td><td>25.1</td><td>24.4</td><td>24.3</td><td>24.4</td></tr><tr><td>λ<sub>4</sub> (-)</td><td>4.9</td><td>4.9</td><td>4.9</td><td>5</td><td>4.9</td><td>4.9</td><td>5</td><td>4.9</td><td>5.1</td><td>5.2</td></tr><tr><td>λ<sub>5</sub> (-)</td><td>6.4</td><td>6.2</td><td>5.8</td><td>5.7</td><td>5.4</td><td>5.3</td><td>5.2</td><td>5</td><td>5</td><td>5</td></tr><tr><td>Operational cost (€)</td><td>185</td><td>169</td><td>157</td><td>146</td><td>138</td><td>130</td><td>124</td><td>118</td><td>114</td><td>110</td></tr><tr><td>Material cost (€)</td><td>90</td><td>96</td><td>101</td><td>107</td><td>111</td><td>116</td><td>120</td><td>124</td><td>129</td><td>134</td></tr><tr><td>Total cost (€)</td><td>274</td><td>265</td><td>258</td><td>253</td><td>249</td><td>246</td><td>244</td><td>243</td><td>243</td><td>244</td></tr></table>
## Appendix B. On the dependency of objective function to flow rate

Consider two configurations ($Q_1$, $B_1$) and ($Q_2$, $B_1$), with $Q_1 < Q_2$. When any given altitude $z$ is reached, the printing time $t_1(z)$ for the flow rate $Q_1$ will always be higher than the printing time $t_2(z)$ corresponding to $Q_2$. This is due to the fact that at any altitude, the volume of the structures $V(z)$ is identical, and because $V = Q \cdot t$.

Since the mass does not vary, and since the pressing force only depends on $B$, it can be concluded that the structures are subject to the same load. Additionally, it is clear that the following relation holds:

$$ \begin{cases} \lambda_1(t_1) > \lambda_2(t_2) \\ \sigma_d(t_1) > \sigma_d(t_2) \end{cases} \qquad (B.1) $$

Therefore, any solution that satisfies the constraints for the flow rate $Q_2$ also satisfies them for the flow rate $Q_1$. This shows that the set of feasible solutions $\mathcal{F}(Q_2)$ is included in $\mathcal{F}(Q_1)$.

$$ (Q_1 < Q_2) \Rightarrow \mathcal{F}(Q_2) \in \mathcal{F}(Q_1) \qquad (B.2) $$

The order relation between the optima is a consequence of this inclusion, since the volume only depends on $B$.

$$ (Q_1 < Q_2) \Rightarrow \min_{B \in \mathcal{F}(Q_1)} V(B) \le \min_{B \in \mathcal{F}(Q_2)} V(B) \qquad (B.3) $$

## Appendix C. On the Pareto front

The Pareto front is a concept from multi-objective optimisation: it is the set of all the Pareto optimal points. Consider two performance metrics $f_1(Q, B)$ and $f_2(Q, B)$, then $(Q^*, B^*)$ is a Pareto optimum if it is not strictly dominated by another point.

We consider the case where $f_1$ is the mass and $f_2$ is the printing time, both to minimise. As already noticed, the two values are linked by the following relation:

$$ f_1(Q, B) = Q \cdot f_2(Q, B) \qquad (C.1) $$

At any given $Q = Q_0$, the minimiser of $f_1$ is also the minimiser of $f_2$. This minimiser is therefore the only candidate in the hyper-plane $(Q_0, B)$ to be a Pareto optimum. Therefore, the Pareto front can be plotted by solving a series of single-valued optimisation problem. These optimisation routines can be launched in parallel, with little cost overhead.

## Data availability

Data will be made available on request.

## References

[1] F. Craveiro, J.P. Duarte, H. Bartolo, P.J. Bartolo, Additive manufacturing as an enabling technology for digital construction: A perspective on construction 4.0, Autom. Constr. (2019) http://dx.doi.org/10.1016/j.autcon.2019.03.011.

[2] M. Batikha, R. Jotangia, M.Y. Baaj, I. Mousleh, 3D concrete printing for sustainable and economical construction: A comparative study, Autom. Constr. 134 (2022) 104087, http://dx.doi.org/10.1016/j.autcon.2021.104087.

[3] H. Van Damme, Concrete material science: Past, present, and future innovations, Cem. Concr. Res. 112 (2018) 5–24, http://dx.doi.org/10.1016/j.cemconres.2018.05.002, SI: Digital concrete 2018.

[4] 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.

[5] V. Nguyen-Van, S. Li, J. Liu, K. Nguyen, P. Tran, Modelling of 3D concrete printing process: A perspective on material and structural simulations, Additive Manuf. 61 (2023) 103333, http://dx.doi.org/10.1016/j.addma.2022.103333.

[6] S.A. Khan, M. Koç, Numerical modelling and simulation for extrusion-based 3D concrete printing: The underlying physics, potential, and challenges, Results Mater. 16 (2022) 100337, http://dx.doi.org/10.1016/j.rinma.2022.100337.

[7] A. Perrot, A. Pierre, V. Nerella, R. Wolfs, E. Keita, S. Nair, N. Neithalath, N. Roussel, V. Mechtcherine, From analytical methods to numerical simulations: A process engineering toolbox for 3D concrete printing, Cem. Concr. Compos. 122 (2021) 104164, http://dx.doi.org/10.1016/j.cemconcomp.2021.104164.

[8] 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, SI: Digital concrete 2018.

[9] A. Suiker, Mechanical performance of wall structures in 3D printing processes: Theory, design tools and experiments, Int. J. Mech. Sci. 137 (2018) 145–170, http://dx.doi.org/10.1016/j.ijmecsci.2018.01.010.

[10] J. Kruger, S. Zeranka, G. van Zijl, 3D concrete printing: A lower bound analytical model for buildability performance quantification, Autom. Constr. 106 (2019) 102904, http://dx.doi.org/10.1016/j.autcon.2019.102904.

[11] V. Nguyen-Van, H. Nguyen-Xuan, B. Panda, P. Tran, 3D concrete printing modelling of thin-walled structures, Structures 39 (2022) 496–511, http://dx.doi.org/10.1016/j.istruc.2022.03.049.

[12] G. Vantyghem, T. Ooms, W.D. Corte, FEM modelling techniques for simulation of 3D concrete printing, 2020, http://dx.doi.org/10.48550/arxiv.2009.06907, arXiv:2009.06907.

[13] 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.

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

[15] 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.

[16] R. Mesnil, V. Poussard, K. Sab, J.-F. Caron, On the geometrical origin of the anisotropy in extrusion-based 3d printed structures, Eng. Struct. 275 (2023) 115082, http://dx.doi.org/10.1016/j.engstruct.2022.115082.

[17] V. Nguyen-Van, B. Panda, G. Zhang, H. Nguyen-Xuan, P. Tran, Digital design computing and modelling for 3-D concrete printing, Autom. Constr. 123 (2021) 103529, http://dx.doi.org/10.1016/j.autcon.2020.103529.

[18] C. Roux, K. Kuzmenko, N. Roussel, R. Mesnil, A. Feraille, Life cycle assessment of a concrete 3D printing process, Int. J. Life Cycle Assess. 28 (1) (2023) 1–15, http://dx.doi.org/10.1007/s11367-022-02111-3.

[19] 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.

[20] I. Dressler, N. Freund, D. Lowke, The effect of accelerator dosage on fresh concrete properties and on interlayer strength in shotcrete 3D printing, Materials 13 (2) (2020) http://dx.doi.org/10.3390/ma13020374.

[21] L. Demont, R. Mesnil, N. Ducoulombier, J.-F. Caron, Affordable online structure 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.

[22] R. Comminal, M.P. Serdeczny, D.B. Pedersen, J. Spangenberg, Numerical modeling of the strand deposition flow in extrusion-based additive manufacturing, Addit. Manuf. 20 (2018) 68–76, http://dx.doi.org/10.1016/j.addma.2017.12.013.

[23] T. Gobin, S. Andraos, T. Schwartz, R. Vriet, A framework for modelling, simulating and parametrically programming heterogeneous industrial machines, Constr. Robot. 6 (2) (2022) 101–108, http://dx.doi.org/10.1007/s41693-022-00078-z.

[24] A. Suiker, R. Wolfs, S. Lucas, T. Salet, Elastic buckling and plastic collapse during 3D concrete printing, Cem. Concr. Res. 135 (2020) 106016, http://dx.doi.org/10.1016/j.cemconres.2020.106016.

[25] N. Roussel, A thixotropy model for fresh fluid concretes: Theory, validation and applications, Cem. Concr. Res. 36 (10) (2006) 1797–1806, http://dx.doi.org/10.1016/j.cemconres.2006.05.025.

[26] X. Chateau, G. Ovarlez, K.L. Trung, Homogenization approach to the behavior of suspensions of noncolloidal particles in yield stress fluids, J. Rheol. 52 (2) (2008) 489–506, http://dx.doi.org/10.1122/1.2838254.

[27] A. Cody, H. Lee, R. Cody, P. Spry, The effects of chemical environment on the nucleation, growth, and stability of ettringite [Ca3Al(OH)6]2(SO4)3 ·26H2O, Cem. Concr. Res. 34 (5) (2004) 869–881, http://dx.doi.org/10.1016/j.cemconres.2003.10.023.

[28] M. Motamedi, R. Mesnil, A.-M. Tang, J.-M. Pereira, O. Baverel, Structural build-up of 3D printed earth by drying, Addit. Manuf. 95 (2024) 104492, http://dx.doi.org/10.1016/j.addma.2024.104492.

[29] R. Nicolas, B. Richard, D. Nicolas, I. Irina, K.J. Temitope, L. Dirk, M. Viktor, M. Romain, P. Arnaud, P. Ursula, R. Lex, S. Dietmar, W. Timothy, W. Rob, Z. Wenqiang, Assessing the fresh properties of printable cement-based materials: High potential tests for quality control, Cem. Concr. Res. 158 (2022) 106836, http://dx.doi.org/10.1016/j.cemconres.2022.106836.

[30] C. Preisinger, M. Heimrath, Karamba—a toolkit for parametric structural design, Struct. Eng. Int. 24 (2) (2014) 217–221, http://dx.doi.org/10.2749/101686614X13830790993483.

[31] J. Argyris, M. Papadrakakis, C. Apostolopoulou, S. Koutsourelakis, The TRIC shell element: theoretical and numerical investigation, Comput. Methods Appl. Mech. Engrg. 182 (1) (2000) 217–245, http://dx.doi.org/10.1016/S0045-7825(99)00094-8.
[32] L. Mourad, J. Bleyer, R. Mesnil, J. Nseir, K. Sab, W. Raphael, Topology optimization of load-bearing capacity, Struct. Multidiscip. Optim. 64 (3) (2021) 1367-1383, http://dx.doi.org/10.1007/s00158-021-02923-1.

[33] J. Bleyer, P. De Buhan, Lower bound static approach for the yield design of thick plates, Internat. J. Numer. Methods Engrg. 100 (11) (2014) 814-833, http://dx.doi.org/10.1002/nme.4776.

[34] M. Firl, K.-U. Bletzinger, Shape optimization of thin walled structures governed by geometrically nonlinear mechanics, Comput. Methods Appl. Mech. Engrg. 237 (2012) 107-117, http://dx.doi.org/10.1016/j.cma.2012.05.016.

[35] D.R. Jones, C.D. Perttunen, B.E. Stuckman, Lipschitzian optimization without the Lipschitz constant, J. Optim. Theory Appl. 79 (1993) 157-181, http://dx.doi.org/10.1007/BF00941892.

[36] M.J.D. Powell, A direct search optimization method that models the objective and constraint functions by linear interpolation, in: Advances in Optimization and Numerical Analysis, Springer Netherlands, Dordrecht, 1994, pp. 51-67, http://dx.doi.org/10.1007/978-94-015-8330-5_4.

[37] S.G. Johnson, The NLopt nonlinear-optimization package, 2007, https://github.com/stevengj/nlopt.

[38] H. Stehling, F. Scheurer, J. Roulier, F. Gramazio, M. Kohler, S. Langenberg, Bridging the gap from CAD to CAM: concepts, caveats and a new grasshopper plug-in, Fabricate: Negot. Des. Mak. (2014) 52-59, http://dx.doi.org/10.2307/j.ctt1tp3c5w10.

[39] R. Mesnil, T. Gobin, L. Demont, P. Margerit, N. Ducoulombier, C. Douthe, J.-F. Caron, Flexible digital manufacturing of timber construction: the design and fabrication of a free-form nexorade, Constr. Robot. 7 (2) (2023) 193-212, http://dx.doi.org/10.1007/s41693-023-00105-7.

[40] K. Kuzmenko, N. Ducoulombier, A. de Feraille, N. Roussel, Environmental impact of extrusion-based additive manufacturing: generic model, power measurements and influence of printing resolution, Cem. Concr. Res. 157 (2022) 106807, http://dx.doi.org/10.1016/j.cemconres.2022.106807.

[41] V. De Bono, N. Ducoulombier, R. Mesnil, J.F. Caron, Methodology for formulating low-carbon printable mortar through particles packing optimization, Cem. Concr. Res. 176 (2024) 107403, http://dx.doi.org/10.1016/j.cemconres.2023.107403.

[42] J.-A. Désidéri, B. Abou El Majd, A. Janka, Nested and self-adaptive bézier parameterizations for shape optimization, J. Comput. Phys. 224 (1) (2007) 117-131, http://dx.doi.org/10.1016/j.jcp.2006.12.016, Special Issue Dedicated to Professor Piet Wesseling on the occasion of his retirement from Delft University of Technology.

[43] 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.

[44] C. Calladine, A shell-buckling paradox resolved, in: Advances in the Mechanics of Plates and Shells: The Avinoam Libai Anniversary Volume, Springer Netherlands, Dordrecht, 2002, pp. 119-134, http://dx.doi.org/10.1007/0-306-46954-5_8.

[45] Z. Diab, D.-P. Do, S. Rémond, D. Hoxha, Probabilistic prediction of structural failure during 3D concrete printing processes, Mater. Struct. 56 (4) (2023) 73, http://dx.doi.org/10.1617/s11527-023-02167-z.