Numerical simulation of three dimensional concrete printing based on a unified fluid and solid mechanics formulation

Abstract

Deformation control constitutes one of the main technological challenges in three dimensional (3D) concrete printing, and it presents a challenge that must be addressed to achieve a precise and reliable construction process. Model-based information of the expected deformations and stresses is required to optimize the construction process in association with the specific properties of the concrete mix. In this work, a novel thermodynamically consistent finite strain constitutive model for fresh and early-age 3D-printable concrete is proposed. The model is then used to simulate the 3D concrete printing process to assess layer shapes, deformations, forces acting on substrate layers and prognoses of possible structural collapse during the layer-by-layer buildup. The constitutive formulation is based on a multiplicative split of the deformation gradient into elastic, aging and viscoplastic parts, in combination with a hyperelastic potential and considering evolving material properties to account for structural buildup or aging. One advantage of this model is the stress-update-scheme, which is similar to that of small strain plasticity and therefore enables an efficient integration with existing material routines. The constitutive model uses the particle finite element method, which serves as the simulation framework, allowing for modeling of the evolving free surfaces during the extrusion process. Computational analyses of three printed layers are used to create deformation plots, which can then be used to control the deformations during 3D concrete printing. This study offers further investigations, on the structural level, focusing on the potential structural collapse of a 3D printed concrete wall. The capability of the proposed model to simulate 3D concrete printing processes across the scales—from a few printed layers to the scale of the whole printed structure—in a unified fashion with one constitutive formulation, is demonstrated.

Graphical abstract

Keywords

particle finite element method / 3D concrete printing / multiplicative split / additive manufacturing / elasto-viscoplasticity

1 Introduction

Additive manufacturing techniques using fresh concrete materials, such as extrusion-based three dimensional (3D) concrete printing, have emerged in the construction sector in recent decades [1–4]. These techniques have the potential to increase the productivity of construction processes through automation and to improve the environmental impact through material-efficient, topologically optimized designs [5,6]. On the other hand, many technological challenges, such as material characterization, interlayer bond, integration of reinforcement or deformation control, need to be addressed to obtain a reliable and efficient production technique [7]. Deformation control constitutes the main challenge in finding the optimal process parameters for constructing a structure without collapse or severe deformations. In this context and in order to develop suitable models for 3D concrete printing, the 3D concrete printing process can be characterized using four production steps: pumping, extrusion, deposition and layer-by-layer buildup [8]. The analysis of the extrusion and deposition steps at the scale of a few layers can be modeled using computational approaches, which capture the complex flow behavior of the material and the large deformations involved [9]. Such approaches allow for an in-depth analysis of the deformation behavior of the layers within a time frame of seconds when the material is deposited. Understanding and taking full control of the deformations and of the printed layer shape is of high significance for enabling a more precise control of the printing process. In contrast, layer-by-layer buildup is defined as the deposition process at the structural level within a time frame of minutes to hours, which can be modeled using quasi-static analysis to assess structural collapse during printing [10].

In previous research, the numerical simulation of extrusion and deposition processes has predominantly been carried out using computational fluid dynamics (CFD) models adopting viscoplastic or elasto-viscoplastic constitutive formulations [9]. The shapes of one printed layer were analyzed based either on CFD approaches in Refs. [11–13] or based on the particle finite element method (PFEM) in Refs. [14–16]. Numerical simulations to analyze the material distribution at corners are shown in Ref. [17]. Those investigations revealed that the pressures generated under the extrusion nozzle can be a few times larger than the self-weight of one layer and can lead to deformations in substrate layers or influence the interlayer bond. When printing multiple layers, the shapes, the deformations, and the contact forces acting on the substrate, were previously analyzed in Refs. [15,18–22]. According to the computational and experimental study in Ref. [23], the pressures acting on the substrate affect the final interlayer bond strength of the structure. As outlined above, great effort has been made in previous studies to analyze and model the extrusion and deposition processes during 3D concrete printing based on fluid dynamics approaches.

For simulations of the layer-by-layer buildup of 3D concrete printing, a quasi-static solid mechanics finite element model with an element activation technique to mimic the printing process has typically been adopted [9]. For example, structural collapse of a 3D-printed cylinder was investigated in Ref. [24] using a solid mechanics finite element model and the same model was compared in Ref. [10] with an analytical solution for structural collapse of a 3D-printed wall element. A similar element activation technique using a voxel-based approach for the layer-by-layer buildup of 3D concrete printing is presented in Ref. [25]. More recently, a lattice approach for modeling 3D concrete printing, including damage has been proposed in Refs. [26,27]. The computational approaches for the layer-by-layer buildup discussed above assumed an idealized rectangular cross section and were based on a solid mechanics constitutive formulation of the fresh concrete. Influences from the extrusion and deposition stage, such as the real layer shapes, deformations, or dynamic forces, were neglected in these models. A more elegant modeling approach would aim for a unified fluid and solid mechanics constitutive formulation for the freshly printed concrete material, allowing for seamless simulations across the different processing steps: extrusion, deposition and layer-by-layer buildup.

This work presents a first step towards a unified fluid and solid mechanics constitutive formulation, which can be used for any concrete age and any process step during 3D concrete printing. Section 2 provides a detailed overview of the proposed constitutive model, based on an elasto-viscoplastic description with time evolving material properties. The constitutive model relies on a finite strain formulation, adopting a multiplicative split of the deformation gradient into elastic, aging and viscoplastic parts, and constitutes a combination of the elasto-viscoplastic and aging models as developed in Refs. [16,28]. A particular feature of this model is the stress-update-algorithm, which is similar to that in small strain plasticity schemes. The proposed constitutive model is implemented in a PFEM model to solve large deformation problems. In Section 3, a brief overview of the adopted face-based smoothed PFEM model based on Ref. [16] is given. To model the contact condition between the printed layers more accurately, a novel PFEM implementation of a tying algorithm using a Dirichlet–Neumann coupling and a transmission matrix [29,30] is shown.

The computational investigations in this work predominantly focus on the extrusion and deposition process by analyzing three layers of printed concrete, as reported in Section 4. The deformations, layer shapes and stresses, for several printing scenarios with varying process parameters, such as printing speed, standoff distance or nozzle geometry, are analyzed in detail. Deformation plots are proposed as a design tool to find optimal process parameters. Section 5 demonstrates the model’s capability of solving a structural collapse problem during 3D concrete printing based on a printing process of a rectangular hollow box. In this context, a simplified modeling strategy to account for realistic layer shapes obtained during the extrusion and deposition process is proposed.

2 Constitutive model

The analysis in the main parts of this work is restricted to the deposition of a few layers within a time frame of a few minutes at maximum. Within this time frame, irreversible creep strains and shrinkage strains are negligible [31,32] and the deformation processes are dominated by elastic and viscoplastic strains and by evolving material properties due to structural build-up or aging. This section provides a framework for a hyperelastic-viscoplastic constitutive model, incorporating evolving material properties for 3D-printed fresh concrete at finite strains.

2.1 Kinematics

To account for the large deformations of the printed layer during extrusion and deposition, a thermodynamically consistent finite strain constitutive model is required. Instead of using a hypoelastic model, as typically used for study of 3D-printed concrete [33,34], a hyperelastic formulation based on a multiplicative split of the deformation gradient is adopted [35,36]. This model can easily be shown to account for thermodynamic consistency and results in a convenient stress-update-algorithm similar to that for small strain plasticity. The deformation gradient is multiplicatively split into elastic, aging, and viscoplastic parts as described by:

(1)\mathit{F}={\mathit{F}}^{\mathrm{e}}{\mathit{F}}^{\mathrm{a}}{\mathit{F}}^{\mathrm{p}},

with F_{ij}=\mathrm{\partial }{x}_i/\mathrm{\partial }{X}_j denoting the deformation gradient of the motion at X, {\mathit{F}}^{\mathrm{e}} denoting the elastic deformation gradient, {\mathit{F}}^{\mathrm{a}} denoting the aging deformation gradient and {\mathit{F}}^{\mathrm{p}} denoting the viscoplastic deformation gradient. The superscripts ‘e’, ‘a’, and ‘p’ will be used for elastic, aging, and plastic quantities, respectively. This multiplicative split implies the assumption of the existence of a local unstressed intermediate configuration denoted, as in Fig.1, as the aging intermediate configuration. Fig.1 illustrates the local multiplicative decomposition of the deformation gradient including the coordinates X and domains \mathrm{\Omega } of the different configurations. In the following and in Fig.1, quantities of the reference configuration are labeled with capital letters, quantities of the plastic intermediate configuration are labeled with a line over a capital letter, quantities of the aging intermediate configuration are labeled with a tilde over a capital letter and quantities of the current configuration are labeled with lowercase letters. Based on the definition of the multiplicative split in Eq. (1), the velocity gradient tensor l can be defined using the definitions of various elastic, aging and viscoplastic velocity gradients with respect to the different configurations (Fig.1) as:

(2)\mathit{l}=\dot{\mathit{F}}{\mathit{F}}^{-1}=\underset{={\mathit{l}}^{\mathrm{e}}}{\underbrace{{\dot{\mathit{F}}}^{\mathrm{e}}{\mathit{F}}^{{\mathrm{e}}^{-1}}}}+\underset{={\mathit{l}}^{\mathrm{a}}}{\underbrace{{\mathit{F}}^{\mathrm{e}}\underset{={\tilde{\mathit{L}}}^{\mathrm{a}}}{\underbrace{{\dot{\mathit{F}}}^{\mathrm{a}}{\mathit{F}}^{{\mathrm{a}}^{-1}}}}{\mathit{F}}^{{\mathrm{e}}^{-1}}}}+\underset{={\mathit{l}}^{\mathrm{p}}}{\underbrace{{\mathit{F}}^{\mathrm{e}}\underset{={\tilde{\mathit{L}}}^{\mathrm{p}}}{\underbrace{{\mathit{F}}^{\mathrm{a}}\underset{={\overline{\mathit{L}}}^{\mathrm{p}}}{\underbrace{{\dot{\mathit{F}}}^{\mathrm{p}}{\mathit{F}}^{{\mathrm{p}}^{-1}}}}{\mathit{F}}^{{\mathrm{a}}^{-1}}}}{\mathit{F}}^{{\mathrm{e}}^{-1}}}}.

Each velocity gradient can be decomposed into its symmetric rate of deformation tensor d and its skew-symmetric spin tensor w with respect to the different configurations defined as:

(3)\begin{array}{lll}\mathit{l}=\mathit{d}+\mathit{w},& \mathit{d}={\displaystyle \frac{1}{2}}\left(\mathit{l}+{\mathit{l}}^{\mathrm{T}}\right),& \mathit{w}={\displaystyle \frac{1}{2}}\left(\mathit{l}-{\mathit{l}}^{\mathrm{T}}\right),\\ {\mathit{l}}^{\mathrm{a}}={\mathit{d}}^{\mathrm{a}}+{\mathit{w}}^{\mathrm{a}},& {\mathit{d}}^{\mathrm{a}}={\displaystyle \frac{1}{2}}\left({\mathit{l}}^{\mathrm{a}}+{\mathit{l}}^{{\mathrm{a}}^{\mathrm{T}}}\right),& {\mathit{w}}^a={\displaystyle \frac{1}{2}}\left({\mathit{l}}^{\mathrm{a}}-{\mathit{l}}^{{\mathrm{a}}^{\mathrm{T}}}\right),\\ {\tilde{\mathit{L}}}^{\mathrm{a}}={\tilde{\mathit{D}}}^{\mathrm{a}}+{\tilde{\mathit{W}}}^{\mathrm{a}},& {\tilde{\mathit{D}}}^{\mathrm{a}}={\displaystyle \frac{1}{2}}\left({\tilde{\mathit{L}}}^{\mathrm{a}}+{\tilde{\mathit{L}}}^{{\mathrm{a}}^{\mathrm{T}}}\right),& {\tilde{\mathit{W}}}^a={\displaystyle \frac{1}{2}}\left({\tilde{\mathit{L}}}^{\mathrm{a}}-{\tilde{\mathit{L}}}^{{\mathrm{a}}^{\mathrm{T}}}\right),\\ {\mathit{l}}^{\mathrm{p}}={\mathit{d}}^{\mathrm{p}}+{\mathit{w}}^{\mathrm{p}},& {\mathit{d}}^{\mathrm{p}}={\displaystyle \frac{1}{2}}\left({\mathit{l}}^{\mathrm{p}}+{\mathit{l}}^{{\mathrm{p}}^{\mathrm{T}}}\right),& {\mathit{w}}^p={\displaystyle \frac{1}{2}}\left({\mathit{l}}^{\mathrm{p}}-{\mathit{l}}^{{\mathrm{p}}^{\mathrm{T}}}\right),\\ {\tilde{\mathit{L}}}^{\mathrm{p}}={\tilde{\mathit{D}}}^{\mathrm{p}}+{\tilde{\mathit{W}}}^{\mathrm{p}},& {\tilde{\mathit{D}}}^{\mathrm{p}}={\displaystyle \frac{1}{2}}\left({\tilde{\mathit{L}}}^{\mathrm{p}}+{\tilde{\mathit{L}}}^{{\mathrm{p}}^{\mathrm{T}}}\right),& {\tilde{\mathit{W}}}^p={\displaystyle \frac{1}{2}}\left({\tilde{\mathit{L}}}^{\mathrm{p}}-{\tilde{\mathit{L}}}^{{\mathrm{p}}^{\mathrm{T}}}\right),\\ {\overline{\mathit{L}}}^{\mathrm{p}}={\overline{\mathit{D}}}^{\mathrm{p}}+{\overline{\mathit{W}}}^{\mathrm{p}},& {\overline{\mathit{D}}}^{\mathrm{p}}={\displaystyle \frac{1}{2}}\left({\overline{\mathit{L}}}^{\mathrm{p}}+{\overline{\mathit{L}}}^{{\mathrm{p}}^{\mathrm{T}}}\right),& {\overline{\mathit{W}}}^p={\displaystyle \frac{1}{2}}\left({\overline{\mathit{L}}}^{\mathrm{p}}-{\overline{\mathit{L}}}^{{\mathrm{p}}^{\mathrm{T}}}\right).\end{array}

The constitutive model, given in the following section, requires the definition of the elastic left Cauchy–Green deformation tensor as:

(4){\mathit{b}}^{\mathrm{e}}={\mathit{F}}^{\mathrm{e}}{\mathit{F}}^{{\mathrm{e}}^{\mathrm{T}}}.

The time derivative of b^e is further obtained, using Eqs. (2)–(4), which yields:

(5){\dot{\mathit{b}}}^{\mathrm{e}}=\mathit{l}{\mathit{b}}^{\mathrm{e}}+{\mathit{b}}^{\mathrm{e}}{\mathit{l}}^{\mathrm{T}}-2{\mathit{F}}^{\mathrm{e}}{\tilde{\mathit{D}}}^{\mathrm{a}}{\mathit{F}}^{{\mathrm{e}}^{\mathrm{T}}}-2{\mathit{F}}^{\mathrm{e}}{\tilde{\mathit{D}}}^{\mathrm{p}}{\mathit{F}}^{{\mathrm{e}}^{\mathrm{T}}}.

Equation (5) is later used to evaluate the dissipation inequality and to show thermodynamic consistency of the constitutive model. The following sections focus on the constitutive assumptions required to describe the kinematics as provided in this section (above). These assumptions include the definition of a hyperelastic potential and the specification of the evolution equations for the aging and plastic deformation rate tensors {\tilde{\mathit{D}}}^{\mathrm{a}} and {\tilde{\mathit{D}}}^{\mathrm{p}} of the aging intermediate configuration in Eq. (5). The evolution equations presented in the following sections underly the assumptions that the aging spin tensor and the plastic spin tensor are zero ({\mathit{w}}^{\mathrm{a}}={\mathit{w}}^{\mathrm{p}}=0). As a consequence of these assumptions, it follows that {\tilde{\mathit{W}}}^{\mathrm{a}}={\tilde{\mathit{W}}}^{\mathrm{p}}=0, but {\overline{\mathit{W}}}^{\mathrm{p}}\ne 0. In fact, the order of F^{\mathrm{a}} and F^{\mathrm{p}} is arbitrary in Eq. (1) and can be changed. In this case, the obtained model becomes equivalent to the present one as long as {\mathit{w}}^{\mathrm{a}}={\mathit{w}}^{\mathrm{p}}=0 and the deformation rate tensors {\mathit{d}}^{\mathrm{a}} and {\mathit{d}}^{\mathrm{p}} of the current configuration are the same as those presented in the following sections.

2.2 Hyperelastic potential

The elastic material response of freshly printed concrete is assumed to be isotropic with time evolving material properties due to the thixotropic nature of the material [37]. The hyperelastic potential is given as a function of the left Cauchy–Green deformation tensor b^e and an internal parameter \alpha, and it describes the current state of aging, accounting for time evolving elastic constants. The internal parameter \alpha can be interpreted as the degree of flocculation or hydration of the aging fresh concrete. A suitable physically or experimentally motivated evolution law for \alpha must be carefully selected to ensure thermodynamic consistency of the model. Further details on \alpha and a specific evolution law are discussed in the following sections. According to Ref. [28], a hyperelastic potential based on logarithmic strains results in a convenient stress-update-scheme. The potential is formulated in terms of principal logarithmic strains as:

(6)\begin{array}{rl}{\rho }_0\Psi =& \alpha G_0\left({\left({\epsilon }_1^{\mathrm{e}}\right)}^2+{\left({\epsilon }_2^{\mathrm{e}}\right)}^2+{\left({\epsilon }_3^{\mathrm{e}}\right)}^2\right)\\ & +\left(\frac{\alpha K_0}{2}-\frac{\alpha G_0}{3}\right){\left({\epsilon }_1^{\mathrm{e}}+{\epsilon }_2^{\mathrm{e}}+{\epsilon }_3^{\mathrm{e}}\right)}^2,\end{array}

with the reference bulk modulus K_0, the reference shear modulus G_0 and the Ath elastic logarithmic principal eigenstrain {\epsilon }_{\mathrm{A}}^{\mathrm{e}}=\mathrm{l}\mathrm{n}{\lambda }_{\mathrm{A}}^{\mathrm{e}} (A=1,2,\mathrm{a}\mathrm{n}\mathrm{d}3), where {\lambda }_{\mathrm{A}}^{\mathrm{e}} denotes the Ath elastic principal eigenstretch. {\lambda }_{\mathrm{A}}^{\mathrm{e}} is obtained from the spectral decomposition of the elastic left Cauchy–Green deformation tensor {\mathit{b}}^{\mathrm{e}} as:

(7){\mathit{b}}^{\mathrm{e}}={\sum }_{A=1}^3{\left({\lambda }_A^{\mathrm{e}}\right)}^2{\mathit{n}}_A\otimes {\mathit{n}}_A,

where {\mathit{n}}_A denotes the Ath eigenvector. Equation (6) is based on the assumption that the shear and bulk moduli evolve equivalently with respect to \alpha, which is reasonable for aging fresh concrete [9,24]. A more general formulation for the evolution of the shear and bulk moduli can be found in Ref. [28]. The Ath principal Kirchhoff stress {\tau }_A can be obtained from Eq. (6) via the relation {\tau }_A={\rho }_0{\displaystyle \frac{\mathrm{\partial }\Psi }{\mathrm{\partial }{\epsilon }_A^{\mathrm{e}}}}, yielding:

(8){\tau }_A=2\alpha G_0{\epsilon }_A^{\mathrm{e}}+\left(\alpha K_0-\frac{2}{3}\alpha G_0\right)\left({\epsilon }_1^{\mathrm{e}}+{\epsilon }_2^{\mathrm{e}}+{\epsilon }_3^{\mathrm{e}}\right).

The Kirchhoff stresses are related to the Cauchy stresses via \mathit{\tau }=\mathit{\sigma }J, where J denotes the determinant of the deformation gradient.

2.3 Aging model

The evolution of elastic constants due to aging can be interpreted macroscopically through the creation of new unloaded springs in the current configuration. For example, the increase of the material stiffness over time can be discretely modeled by inserting an unloaded elastic spring at the beginning of every time step during the analysis. As the creation of new springs in each time step is unfeasible in computational models due to high memory requirement, a single spring element, with evolving material properties as introduced in Subsection 2.2, is used in this work. The aging deformation gradient must be further introduced to mimic the behavior of irreversible age-related strains, as Eq. (6) alone would result in physically impossible strain jumps.

An evolution equation for the age-related strains can be derived from the condition that the stresses must remain constant while the material properties evolve; the deformation gradient is also held constant. This can be formulated as in Ref. [28]:

(9)\frac{\mathrm{\partial }\mathit{\tau }}{\mathrm{\partial }\alpha }\stackrel{!}{=}0,\text{while}\frac{\mathrm{\partial }\mathit{F}}{\mathrm{\partial }\alpha }\stackrel{!}{=}0,\frac{\mathrm{\partial }{\mathit{F}}^{\mathrm{e}}}{\mathrm{\partial }\alpha }\ne 0,\frac{\mathrm{\partial }{\mathit{F}}^{\mathrm{a}}}{\mathrm{\partial }\alpha }\ne 0.

The evolution equation for the age-related strains was derived in Ref. [28], and relates to the hyperelastic model introduced in Subsection 2.2. The aging rate of deformation tensor is written under the assumption that {\mathit{w}}^{\mathrm{a}}={\tilde{\mathit{W}}}^{\mathrm{a}}=0, so that:

(10){\mathit{d}}^{\mathrm{a}}={\mathit{F}}^{\mathrm{e}}{\tilde{\mathit{D}}}^{\mathrm{a}}{\mathit{F}}^{{\mathrm{e}}^{-1}}=\frac{1}{2}\frac{\dot{\alpha }}{\alpha }\mathrm{l}\mathrm{n}\left[{\mathit{b}}^{\mathrm{e}}\right],

where \mathrm{l}\mathrm{n}\left[{\mathit{b}}^{\mathrm{e}}\right] denotes the tensor logarithm of {\mathit{b}}^{\mathrm{e}}. Note that Eq. (10) is only valid for an increasing internal parameter, \alpha \ne 0, and for the hyperelastic potential introduced in Subsection 2.2. In this work, \alpha is described as a function of time, but generally any thermodynamically consistent evolution law can be postulated for \alpha. When another hyperelastic model is to be used, a new evolution equation for the age-related strains must be derived based on the aging condition in Eq. (9).

2.4 Viscoplastic model

The flow of fresh concrete is characterized by a viscoplastic response, which can be described by Bingham or Herschel–Bulkley models [38]. The stress state in a viscoplastic material must exceed the yield stress for the material to undergo viscous deformation. Implementations of viscoplastic fluids in computational models are often based on regularized approaches, which result in a convenient algorithmic treatment of the material response for stress states below the yield stress based on a pseudo-rigid behavior [39]. To accurately track the deformations in substrate layers during 3D-printing of fresh concrete, the material response for stress states below the yield stress cannot be approximated using such pseudo-rigid regularized approaches. Therefore, in the proposed model, the hyperelastic model, as presented in Subsection 2.2, is employed. While the yield criterion in regularized viscoplastic models is smooth, the adoption of a hyperelastic-viscoplastic model requires a viscoplastic sub-model with a yield criterion as is summarized below.

A Bingham model equivalent to the one presented in Ref. [16] is used for describing the viscoplastic material behavior with a von Mises yield function given as:

(11)\varphi \left(\mathit{\tau }\right)={\tau }_{\mathrm{e}\mathrm{q}}-{\tau }_0,

where {\tau }_0 is the yield stress and {\tau }_{\mathrm{e}\mathrm{q}}=\sqrt{{\displaystyle \frac{1}{2}}\mathrm{d}\mathrm{e}\mathrm{v}\left[\mathit{\tau }\right]:\mathrm{d}\mathrm{e}\mathrm{v}\left[\mathit{\tau }\right]} is the equivalent Kirchhoff stress, with the deviatoric Kirchhoff stresses \mathrm{d}\mathrm{e}\mathrm{v}\left[\mathit{\tau }\right]. Note, that the yield stress {\tau }_0 is evolving over time due to aging. The flow rule for the viscoplastic part is assumed to be associative and is formulated under the assumption that {\mathit{w}}^{\mathrm{p}}={\tilde{\mathit{W}}}^{\mathrm{p}}=\mathbf{0} as:

(12){\mathit{d}}^{\mathrm{p}}={\mathit{F}}^{\mathrm{e}}{\tilde{\mathit{D}}}^{\mathrm{p}}{\mathit{F}}^{{\mathrm{e}}^{-1}}=\dot{\gamma }\frac{\mathrm{\partial }\text{ϕ}\left(\mathit{\tau }\right)}{\mathrm{\partial }\mathit{\tau }}=\dot{\gamma }\frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[\mathit{\tau }\right]}{2{\tau }_{\mathrm{e}\mathrm{q}}},

with \dot{\gamma } denoting the plastic multiplier and {\displaystyle \frac{\delta \mathrm{\Phi }\left(\mathit{\tau }\right)}{\delta \mathit{\tau }}} denoting the direction of plastic flow. To obtain a Bingham viscoplastic behavior, the plastic multiplier can be described based on an overstress function [40] as:

(13)

where \mu denotes the plastic viscosity and , which ensures that viscoplastic flow only occurs when the stress state exceeds the yield stress. The plasticity model employs a standard finite strain plasticity formulation [41,42], which can be easily converted into a Drucker–Prager model to account for the influence of pressure on the yield criterion. This becomes important as the fresh concrete ages [24].

2.5 Thermodynamic consistency

In this section, thermodynamic consistency of the constitutive model is shown. Assuming an isothermal process, the Clausius–Duhem inequality is written as:

(14)\mathit{\tau }:\mathit{d}-{\rho }_0\dot{\mathrm{\Psi }}\ge 0,

with the definition of the Kirchhoff stresses for isotropic constitutive models as \mathit{\tau }=2{\rho }_0{\displaystyle \frac{\mathrm{\partial }\mathrm{\Psi }}{\mathrm{\partial }{\mathit{b}}_{\mathrm{e}}}}{\mathit{b}}_{\mathrm{e}}, the time evolution of the hyperelastic potential yields:

(15){\rho }_0\dot{\mathrm{\Psi }}=\mathit{\tau }:\mathit{d}-\mathit{\tau }:{\mathit{d}}^{\mathrm{a}}-\mathit{\tau }:{\mathit{d}}^{\mathrm{p}}+{\rho }_0\frac{\mathrm{\partial }\dot{\mathrm{\Psi }}}{\mathrm{\partial }\alpha }\dot{\alpha }.

After a lengthy derivation, which is given in Ref. [28], it holds that:

(16){\rho }_0\frac{\mathrm{\partial }\dot{\mathrm{\Psi }}}{\mathrm{\partial }\alpha }\dot{\alpha }=-\frac{1}{2}\mathit{\tau }:{\mathit{d}}^{\mathrm{a}}.

Substitution of Eqs. (15) and (16) into the Clausius–Duhem inequality, Eq. (14), yields:

(17)\frac{3}{2}\mathit{\tau }:{\mathit{d}}^{\mathrm{a}}+\mathit{\tau }:{\mathit{d}}^{\mathrm{p}}\ge 0.

Based on the definition for {\mathit{d}}^{\mathrm{a}} given in Subsection 2.3, the first term of Eq. (17) is always larger than zero when \alpha is an increasing function in time, which is assumed a priori in this work. Note that when destructuration processes in fresh concrete are to be considered, a new evolution law for {\mathit{d}}^{\mathrm{a}} must be derived. The second term in Eq. (17) constitutes the contribution from the plasticity model, which has been shown to satisfy thermodynamic consistency in previous works [41,42].

2.6 Time marching scheme and stress update algorithm

The numerical simulation of large deformation problems requires an adequate time marching scheme with an incremental solution procedure to solve the governing equations. In finite strain algorithms [41], the deformation gradient is multiplicatively updated from time step t_n to t_{n+1} via:

(18)\begin{array}{c}{\mathit{F}}_{n+1}={\mathit{f}}_{n+1}{\mathit{F}}_n,\end{array}

with the incremental deformation gradient {\mathit{f}}_{n+1}=\mathit{I}+ {\mathrm{\nabla }}_n\left(\mathrm{\Delta }{\mathit{u}}_{n+1}\right), the displacement increment in one time step \mathrm{\Delta }{\mathit{u}}_{n+1}, I_{ij}={\delta }_{ij}, the Kronecker delta {\delta }_{ij} and the Nabla operator {\mathrm{\nabla }}_n=\mathrm{\partial }◻/\mathrm{\partial }{\mathit{x}}_n evaluated at the current configuration at time step n. Recall that the solution of a first order linear problem \dot{\mathit{A}}\left(t\right)=\mathit{B}\left(t\right)\mathit{A}\left(t\right) in [t_n,t_{n+1}] is the exponential integration rule, as {\mathit{A}}_{n+1}=\mathrm{e}\mathrm{x}\mathrm{p}\left[{\int }_{t_{n+1}}^{t_n}\mathit{B}\left(t\right)\mathrm{d}t\right] {\mathit{A}}_n [41]. Using the exponential integration rule, the evolution equations of the aging and viscoplastic parts, Eqs. (10) and (12), and the definition of the velocity gradients {\tilde{\mathit{L}}}^{\mathrm{a}}={\dot{\mathit{F}}}^{\mathrm{a}}{{\mathit{F}}^{\mathrm{a}}}^{-1} and {\overline{\mathit{L}}}^{\mathrm{p}}={\dot{\mathit{F}}}^{\mathrm{p}}{{\mathit{F}}^{\mathrm{p}}}^{-1} (Eq. (2)) yield, after some reformulations:

(19){\mathit{F}}_{n+1}^{\mathrm{a}}=({\mathit{F}}_{n+1}^{\mathrm{e}}{)}^{-1}\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{{\alpha }_{n+1}-{\alpha }_n}{{\alpha }_n}\frac{1}{2}\mathrm{l}\mathrm{n}\left[{\mathit{b}}_{n+1}^{\mathrm{e}}\right]\right]{\mathit{F}}_{n+1}^{\mathrm{e}}{\mathit{F}}_n^{\mathrm{a}},

(20){\mathit{F}}_{n+1}^{\mathrm{p}}={\left({\mathit{F}}_{n+1}^{\mathrm{a}}\right)}^{-1}{\left({\mathit{F}}_{n+1}^{\mathrm{e}}\right)}^{-1}\mathrm{e}\mathrm{x}\mathrm{p}\left[\mathrm{\Delta }{\gamma }_{n+1}\frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\mathit{\tau }}_{n+1}\right]}{2{\tau }_{\mathrm{e}\mathrm{q},n+1}}\right]{\mathit{F}}_{n+1}^{\mathrm{e}}{\mathit{F}}_{n+1}^{\mathrm{a}}{\mathit{F}}_n^{\mathrm{p}},

where \mathrm{e}\mathrm{x}\mathrm{p}\left[◻\right] is the tensor exponential and \mathrm{\Delta }{\gamma }_{n+1} is the viscoplastic increment within one time step. Note that in Eqs. (19) and (20) we made use of the expression \mathrm{e}\mathrm{x}\mathrm{p}\left[{\mathit{A}}^{-1}\mathit{B}\mathit{A}\right]={\mathit{A}}^{-1}\mathrm{e}\mathrm{x}\mathrm{p}\left[\mathit{B}\right]\mathit{A}, which holds for any invertible tensor A. The updated elastic deformation gradient and the updated elastic left Cauchy–Green deformation tensor are found from Eqs. (1), (4), (18)–(20), leading to:

(21)\begin{array}{rl}{\mathit{F}}_{n+1}^{\mathrm{e}}=& \mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\alpha }_{n+1}-{\alpha }_n}{{\alpha }_n}\frac{1}{2}\mathrm{l}\mathrm{n}\left[{\mathit{b}}_{n+1}^{\mathrm{e}}\right]\right]\\ & \mathrm{e}\mathrm{x}\mathrm{p}\left[-\mathrm{\Delta }{\gamma }_{n+1}\frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\mathit{\tau }}_{n+1}\right]}{2{\tau }_{\mathrm{e}\mathrm{q},n+1}}\right]{\mathit{F}}_{n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}},\end{array}

(22)\begin{array}{rl}{\mathit{b}}_{n+1}^{\mathrm{e}}=& \mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\alpha }_{n+1}-{\alpha }_n}{{\alpha }_n}\mathrm{l}\mathrm{n}\left[{\mathit{b}}_{n+1}^{\mathrm{e}}\right]\right]\\ & \mathrm{e}\mathrm{x}\mathrm{p}\left[-2\mathrm{\Delta }{\gamma }_{n+1}\frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\mathit{\tau }}_{n+1}\right]}{2{\tau }_{\mathrm{e}\mathrm{q},n+1}}\right]{\mathit{b}}_{n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}},\end{array}

with the trial elastic deformation gradient {\mathit{F}}_{n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}={\mathit{f}}_{n+1}{\mathit{F}}_n^{\mathrm{e}} and the trial left Cauchy–Green deformation tensor {\mathit{b}}_{n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}={\mathit{F}}_{n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}{{\mathit{F}}_{n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}}^{\mathrm{T}}. Note that these reformulations are only possible because \mathrm{d}\mathrm{e}\mathrm{v}\left[{\mathit{\tau }}_{n+1}\right], {\mathit{b}}_{n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}, {\mathit{b}}_{n+1}^{\mathrm{e}}, and their tensor exponentials and logarithms commute. The proposed time integration scheme is a standard procedure in multiplicative plasticity [41] and further details on the derivation of the time integration of the aging part can be found in Ref. [28]. The novelty of the proposed scheme lies in the combination of the aging and viscoplastic models, which is shown below to result in a convenient formulation in the principal direction. When applying the tensor logarithm to Eq. (22) and making use of the fact that all tensors in Eq. (22) commute, an expression for the Ath updated elastic logarithmic eigenstretch is given as:

(23){\epsilon }_{A,n+1}^{\mathrm{e}}=\underset{={\hat{\epsilon }}_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}}{\underbrace{\frac{{\alpha }_n}{{\alpha }_{n+1}}\mathrm{l}\mathrm{n}\left[{\lambda }_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}\right]}}-\frac{{\alpha }_n}{{\alpha }_{n+1}}\mathrm{\Delta }{\gamma }_{n+1}\frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\tau }_{A,n+1}^{\mathrm{t}\mathrm{r}}\right]}{2{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}},

with the trial Kirchhoff stresses {\mathit{\tau }}_{n+1}^{\mathrm{t}\mathrm{r}} (see Box 1), which share the same eigenvectors with {\mathit{\tau }}_{n+1}. Based on Eq. (23) a standard return mapping algorithm can be formulated to yield a convenient stress update procedure in principal direction similar to that obtained from small strain plasticity models [42,43]. The contribution to account for the time dependent elastic constants constitutes only small modifications of a standard plasticity algorithm.

Box 1: Stress updating algorithm

Given {\mathit{b}}_n^{\mathrm{e}},{\alpha }_n,{\alpha }_{n+1},\mathrm{\Delta }{\mathit{u}}_{n+1}

1. Calculation of the incremental deformation gradient:

{\mathit{f}}_{n+1}=\mathit{I}+{\displaystyle \frac{\mathrm{\partial }\left(\mathrm{\Delta }{\mathit{u}}_{n+1}\right)}{\mathrm{\partial }{\mathit{x}}_n}}.

2. Calculation of the trial elastic left Cauchy–Green 　　 deformation tensor and the spectral decomposition:

{\mathit{b}}_{n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}={\mathit{f}}_{n+1}{\mathit{b}}_n^{\mathrm{e}}({\mathit{f}}_{n+1}{)}^{\mathrm{T}}=\sum _{A=1}^3 ({\lambda }_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}{)}^2{\mathit{n}}_A\otimes {\mathit{n}}_A.

3. Calculation of the modified trial elastic principal　　 logarithmic eigenstretches:

{\hat{\epsilon }}_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}={\displaystyle \frac{{\alpha }_n}{{\alpha }_{n+1}}}\mathrm{l}\mathrm{n}\left[{\lambda }_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}\right].

4. Calculation of the trial principal Kirchhoff stresses:

{\tau }_{A,n+1}^{\mathrm{t}\mathrm{r}}={\alpha }_{n+1}\left[2G_0{\hat{\epsilon }}_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}+\left(K_0-{\displaystyle \frac{2}{3}}{G}_0\right)\left({\hat{\epsilon }}_{1,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}+{\hat{\epsilon }}_{2,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}+{\hat{\epsilon }}_{3,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}\right)\right].

5. Check of the yield criterion and calculation of the 　　 viscoplastic increment, if required:

\begin{array}{c}\text{i}\text{f}{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}-{\tau }_{0,n+1}>0\to \mathrm{\Delta }{\gamma }_{n+1}={\displaystyle \frac{{\alpha }_{n+1}}{{\alpha }_n}}\left({\displaystyle \frac{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}-{\tau }_{0,n+1}}{{\alpha }_{n+1}{G}_0+\mu /\mathrm{\Delta }t}}\right)\\ \text{e}\text{l}\text{s}\text{e}{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}-{\tau }_{0,n+1}\le 0\to \mathrm{\Delta }{\gamma }_{n+1}=0.\end{array},

6. Calculation of the updated elastic principal 　　 logarithmic eigenstretches and principal Kirchhoff 　　 stresses:

\begin{array}{l}{\epsilon }_{A,n+1}^{\mathrm{e}}={\hat{\epsilon }}_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}-{\displaystyle \frac{{\alpha }_n}{{\alpha }_{n+1}}}\mathrm{\Delta }{\gamma }_{n+1}{\displaystyle \frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\tau }_{A,n+1}^{\mathrm{t}\mathrm{r}}\right]}{2{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}}},\\ {\tau }_{A,n+1}={\tau }_{A,n+1}^{\mathrm{t}\mathrm{r}}-{\displaystyle \frac{{\alpha }_n}{{\alpha }_{n+1}}}\mathrm{\Delta }{\gamma }_{n+1}{\displaystyle \frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\tau }_{A,n+1}^{\mathrm{t}\mathrm{r}}\right]}{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}}}{\alpha }_{n+1}{G}_0.\end{array}

7. Calculation of the updated elastic left Cauchy–　　 Green deformation tensor and the updated Kirchhoff 　　 stresses:

\begin{array}{rl}{\mathit{b}}_{n+1}^{\mathrm{e}}=& {\sum }_{A=1}^3{\left(\mathrm{e}\mathrm{x}\mathrm{p}\left[{\epsilon }_{A,n+1}^{\mathrm{e}}\right]\right)}^2{\mathit{n}}_A\otimes {\mathit{n}}_A\mathrm{a}\mathrm{n}\mathrm{d}{\mathit{\tau }}_{n+1}\\ =& {\sum }_{A=1}^3 {\tau }_{A,n+1}{\mathit{n}}_A\otimes {\mathit{n}}_A.\end{array}

The stress update algorithm of the proposed constitutive model for printed fresh concrete is summarized in Box 1, which is prepared so that it can be readily implemented in standard finite element codes for finite strain problems. Note that the stress update algorithm must be repeated within each global iteration i of the finite element solution procedure in each time step t_n to t_{n+1}. The incremental displacements \mathrm{\Delta }{\mathit{u}}_{n+1}={\mathit{u}}_{n+1}-{\mathit{u}}_n within one time step are obtained from the global solution procedure of the finite element model.

2.7 Spatial algorithmic tangent modulus

The algorithmic tangent modulus is required during the global solution procedure in implicit computational finite strain models. To derive the algorithmic tangent modulus used in a Newton–Raphson algorithm, the stresses must be linearized in direction of the trial strains [42], and this process yields the proposed constitutive model in principal direction the spatial algorithmic tangent moduli as [44–46]:

(24)\begin{array}{rl}J\mathit{c}=& {\sum }_{A,B=1}^3\frac{\mathrm{\partial }{\tau }_{A,n+1}}{\mathrm{\partial }{\epsilon }_{B,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}}{\mathit{n}}_{AABB}-{\sum }_{A=1}^{3}2{\tau }_{A,n+1}{\mathit{n}}_{AAAA}\\ & +{\sum }_{\begin{array}{c}A,B=1,\\ A\ne B\end{array}}^3\frac{{\tau }_{A,n+1}{\left({\lambda }_{B,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}\right)}^2-{\tau }_{B,n+1}{\left({\lambda }_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}\right)}^2}{{\left({\lambda }_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}\right)}^2-{\left({\lambda }_{B,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}\right)}^2}\left({\mathit{n}}_{ABAB}+{\mathit{n}}_{ABBA}\right),\end{array}

where {\mathit{n}}_{ABCD}={\mathit{n}}_A\otimes {\mathit{n}}_B\otimes {\mathit{n}}_C\otimes {\mathit{n}}_D. A detailed derivation of the algorithmic tangent modulus can be found in Refs. [44–46]. In case of equal eigenvalues, L’Hôpital rule is used with respect to the third part on the right-hand side of Eq. (24) [45]. The first part on the right-hand side can be derived using the relations of the updated quantities in Box 1 leading to similar expressions as in analysis of small strain plasticity [42]:

(25)\begin{array}{rl}\frac{\mathrm{\partial }{\tau }_{A,n+1}}{\mathrm{\partial }{\epsilon }_{A,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}}=& {\alpha }_n\left({\displaystyle \frac{4}{3}}{G}_0+K_0\right)-{\displaystyle \frac{{\displaystyle \frac{4}{3}}{\alpha }_n^2{G}_0^2\mathrm{\Delta }{\gamma }_{n+1}}{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}}}\\ & +\left[\frac{{\alpha }_n^2{G}_0^2\mathrm{\Delta }{\gamma }_{n+1}}{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}}-\frac{{\alpha }_n{\alpha }_{n+1}{G}_0^2}{{\alpha }_{n+1}{G}_0+\mu /\mathrm{\Delta }t}\right]\\ & \frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\tau }_{A,n+1}^{\mathrm{t}\mathrm{r}}\right]}{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}}\frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\tau }_{A,n+1}^{\mathrm{t}\mathrm{r}}\right]}{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}},\end{array}

(26)\begin{array}{rl}\frac{\mathrm{\partial }{\tau }_{A,n+1}}{\mathrm{\partial }{\epsilon }_{B,n+1}^{\mathrm{e},\mathrm{t}\mathrm{r}}}=& {\alpha }_n\left(K_0-\frac{2}{3}{G}_0\right)+\frac{{\displaystyle \frac{2}{3}}{\alpha }_n^2{G}_0^2\mathrm{\Delta }{\gamma }_{n+1}}{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}}\\ & +\left[\frac{{\alpha }_n^2{G}_0^2\mathrm{\Delta }{\gamma }_{n+1}}{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}}-\frac{{\alpha }_n{\alpha }_{n+1}{G}_0^2}{{\alpha }_{n+1}{G}_0+\mu /\mathrm{\Delta }t}\right]\\ & \frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\tau }_{A,n+1}^{\mathrm{t}\mathrm{r}}\right]}{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}}\frac{\mathrm{d}\mathrm{e}\mathrm{v}\left[{\tau }_{B,n+1}^{\mathrm{t}\mathrm{r}}\right]}{{\tau }_{\mathrm{e}\mathrm{q},n+1}^{\mathrm{t}\mathrm{r}}}.\end{array}

2.8 Evolution of the material properties of fresh concrete

This section describes the time evolution of the material properties of printed fresh concrete, also referred to as thixotropy or structural build-up [49]. Strongly sheared and fresh concrete mixtures exhibit the so-called dynamic yield stress {\tau }_{\mathrm{d}}, which characterizes the stress level at which viscous material flow stops. As soon as the material is at rest, flocculation based on colloidal interactions is the main phenomenon and drives the time evolution of the yield stress. After a few minutes, the so-called static yield stress {\tau }_{\mathrm{s}} is reached, which characterizes the stress level that must be overcome to (re)initiate viscous material flow. Over longer resting times of the material, the time evolution of the yield stress is dominated by early hydration products. While flocculation is considered reversible, early hydration is considered irreversible when the material is sheared [49]. In 3D printing of concrete, we assume that the material is sheared heavily after extrusion and remains at rest after deposition. Destructuration processes are therefore not considered in the proposed computational model and the yield stress is formulated as a function of time after deposition. To the best of the authors’ knowledge only the bilinear empirical model developed in Refs. [47,50] takes into account the evolution from a dynamic to a static yield stress for 3D-printed concrete. The computational studies presented in this work are based on the experimental values of the printable high-performance concrete, as summarized in Ref. [47]. The time evolution of the yield stress is illustrated in Fig.2(a). The elastic constants of fresh concrete were not measured in Ref. [47]; they are taken from the experimental results in Ref. [48], as depicted in Fig.2(b) and Fig.2(c). By choosing the reference shear modulus and bulk modulus as G_0=15.18\mathrm{k}\mathrm{P}\mathrm{a} and K_0=32.90\mathrm{k}\mathrm{P}\mathrm{a}, the evolution of \alpha can be defined, with t in seconds, as:

(27)\alpha \left(t\right)=1+7.09\cdot {10}^{-4}\left[\frac{1}{s}\right]\cdot t.

The plastic viscosity is chosen as \mu =5\mathrm{P}\mathrm{a}\cdot \mathrm{s} based on typical values found in Refs. [12,15,51]. The effect of aging on plastic viscosity is neglected in this work, as it is assumed to have little effect once the material is deposited and no longer flows.

3 Simulation framework

The proposed constitutive model from Section 2 is implemented in a numerical model based on the PFEM to deal with large deformation problems during 3D concrete printing. This section summarizes the simulation framework, which is based on the PFEM implementation proposed in Ref. [16]. First, the balance equations of momentum and mass are introduced, which is followed by a summary of the PFEM model. Finally, a novel PFEM implementation of a domain composition algorithm is proposed to tie non-matching meshes together. The domain composition algorithm is required to accurately model the contact between two layers during 3D concrete printing.

3.1 Balance of momentum and mass

The strong form of the balance of momentum and mass equations solved in the PFEM model are summarized as follows. The balance of momentum of a domain Ω during a time interval [0,T] is given as:

(28)\rho \mathit{a}-\mathrm{\nabla }\cdot \mathit{\sigma }-\rho \mathit{b}=0,\mathrm{i}\mathrm{n}\mathrm{\Omega }\times (0,T),

where the density is \rho, the external body force is \rho \mathit{b} and the acceleration is \mathit{a}. Equation (28) requires Dirichlet and Neumann boundary conditions defined on the Dirichlet boundary {\mathrm{\Gamma }}^{\mathrm{u}} and on the Neumann boundary {\mathrm{\Gamma }}^{\mathrm{t}} as:

(29)\mathit{u}={\mathit{u}}^{\ast },\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}^{\mathrm{u}}\times (0,T),

(30)\mathit{\sigma }\cdot \mathit{n}={\mathit{t}}^{\ast },\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}^{\mathrm{t}}\times (0,T),

with the prescribed displacements {\mathit{u}}^{\ast }, the tractions {\mathit{t}}^{\ast } and the normal vector n. To deal with plastic incompressibility during viscoplastic flow of the material, a mixed displacement-pressure formulation is adopted. The Cauchy stresses are decomposed as \mathit{\sigma }=\mathrm{d}\mathrm{e}\mathrm{v}\left[\overline{\mathit{\sigma }}\right]+p\mathit{I}, where \overline{\mathit{\sigma }} are the Cauchy stresses of an integration point obtained from the elastic left Cauchy–Green deformation tensor, I_{ij}={\delta }_{ij}, and p is the pressure at an integration point obtained after discretization from the nodal values. Based on the constitutive model, the pressure p is related to the determinant of the deformation gradient J via the balance of mass as:

(31)\frac{\mathrm{l}\mathrm{n}\left(J\right)}{J}-\frac{p}{\alpha K_0}=0,\text{i}\text{n}\mathrm{\Omega }\times (0,T).

Equations (28) and (31) are the balance equations, which are discretized and solved using finite elements. Details of the derivation of the weak form of Eqs. (28) and (31) can be found in Ref. [16].

3.2 Particle finite element method

Numerical simulations of 3D concrete printing extrusion and deposition processes require a numerical method capable of dealing with large deformations and evolving free boundaries [9]. In previous studies, the PFEM was found to be an excellent method for dealing with the large deformations by using an updated Lagrangian description of motion and a finite element discretization with triangular and tetrahedral elements [15,16]. To avoid severe mesh distortions, PFEM uses re-meshing techniques to guarantee a good mesh quality. The Lagrangian description of motion also leads to an easier integration of constitutive models compared to Eulerian-based methods. This allows the adoption of, for example, standard return mapping algorithms in plasticity and the adoption of the stress update scheme of the proposed constitutive model for 3D-printed fresh concrete summarized in Box 1.

The adopted PFEM model is based on the implementation proposed in Ref. [16], which differs from the original PFEM implementation [52,53] in several aspects. In Ref. [16], a re-meshing procedure is based on a constrained Delaunay triangulation in conjunction with re-triangulation and refinement of the free surface discretization when the mesh becomes too distorted; see Fig.3 and Ref. [54]. This procedure is advantageous over the original alpha shape approach with respect to the conservation of mass and the definition of the free surface. However, the free surface must remain continuous without self-contact or excessive motions, which is the case when simulating 3D concrete printing extrusion and deposition processes. Note, that the initial free surface triangulation is still found from an alpha shape re-meshing approach (Fig.3). After re-meshing, the internal variables must be projected from the old to the new discretization. As proposed in Ref. [16] the isochoric elastic left Cauchy–Green tensor is projected based on a nearest neighbor interpolation and the volumetric part of the elastic left Cauchy–Green tensor is recovered from an incremental superconvergent patch recovery projection to obtain a more robust algorithm.

The balance of momentum and mass equations introduced in Subsection 3.1 are solved implicitly in time using a face-based smoothed finite element discretization of the weak form. The solution procedure of the PFEM algorithm is summarized in Box 2, for one iteration i, with the matrices and vectors summarized in Box 3, all of which are defined in the updated configuration. Details about the derivation of the individual matrices and properties of the smoothed PFEM implementation can be found in Ref. [16]. Note that the face-based smoothed finite element formulation is the 3D counterpart of the 2D edge-based smoothed finite element formulation with stress points located on the faces.

Box 2: Iterative solution scheme of the PFEM algorithm [16]

Initialize the solution variables: {\overline{\mathit{x}}}_{n+1}^i={\overline{\mathit{x}}}_n; \mathrm{\Delta }{\overline{\mathit{u}}}_{n+1}^i=0; {\overline{\mathit{p}}}_{n+1}^i={\overline{\mathit{p}}}_n; {\mathit{b}}_{n+1}^{\mathrm{e},i}={\mathit{b}}_n^{\mathrm{e}}; {\overline{\mathit{v}}}_{n+1}^i={\overline{\mathit{v}}}_n; {\overline{\mathit{a}}}_{n+1}^i={\overline{\mathit{a}}}_n

For each nonlinear iteration i:

1. Performance of the stress update algorithm given in 　　Box 1.

2. Calculation of the residuals:

\begin{array}{l}{\mathit{R}}_{\mathrm{u},n+1}^i={\mathit{F}}_{\mathrm{u},\mathrm{i}\mathrm{n}\mathrm{t},n+1}^i-{\mathit{F}}_{\mathrm{u},\mathrm{e}\mathrm{x}\mathrm{t},n+1}^i+{\mathit{M}}_{\mathrm{u}\mathrm{u},n+1}^i{\overline{\mathit{a}}}_{n+1}^i,\\ {\mathit{R}}_{\mathrm{p},n+1}^i={\mathit{F}}_{\mathrm{p},n+1}^i-{\mathit{K}}_{\mathrm{p}\mathrm{p},\mathrm{K},n+1}^i{\overline{\mathit{p}}}_{n+1}^i-{\mathit{K}}_{\mathrm{p}\mathrm{p},\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b},n+1}^i{\overline{\mathit{p}}}_{n+1}^i.\end{array}

3. Calculation of the tangent matrices:

\mathit{K}=\left[\begin{array}{cc}{\mathit{K}}_{\mathrm{u}\mathrm{u}}& {\mathit{K}}_{\mathrm{u}\mathrm{p}}\\ {\mathit{K}}_{\mathrm{p}\mathrm{u}}& {\mathit{K}}_{\mathrm{p}\mathrm{p}}\end{array}\right]=\left[\begin{array}{cc}{\displaystyle \frac{1-{\alpha }_{\mathrm{\Delta }t}}{\mathrm{\Delta }{t}^2{\beta }_{\mathrm{\Delta }t}}}{\mathit{M}}_{\mathrm{u}\mathrm{u}}+{\mathit{K}}_{\mathrm{u}\mathrm{u},\mathrm{g}\mathrm{e}\mathrm{o}}+{\mathit{K}}_{\mathrm{u}\mathrm{u},\mathrm{m}\mathrm{a}\mathrm{t}}& {\mathit{K}}_{\mathrm{u}\mathrm{p}}\\ {\mathit{K}}_{\mathrm{p}\mathrm{u}}& -{\mathit{K}}_{\mathrm{p}\mathrm{p},\mathrm{K}}-{\mathit{K}}_{\mathrm{p}\mathrm{p},\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}\end{array}\right].

4. Calculation of the displacement and pressure 　　increments:

\left[\begin{array}{c}\mathrm{\Delta }\hat{\mathit{u}}\\ \mathrm{\Delta }\overline{\mathit{p}}\end{array}\right]=-{\mathit{K}}^{-1}\left[\begin{array}{c}{\mathit{R}}_{\mathrm{u},n+1}^i\\ {\mathit{R}}_{\mathrm{p},n+1}^i\end{array}\right].

5. Update the solution variables:

\begin{array}{l}\mathrm{\Delta }{\overline{\mathit{u}}}_{n+1}^{i+1}=\mathrm{\Delta }{\overline{\mathit{u}}}_{n+1}^i+\mathrm{\Delta }\hat{\mathit{u}},\\ {\overline{\mathit{p}}}_{n+1}^{i+1}={\overline{\mathit{p}}}_{n+1}^i+\mathrm{\Delta }\overline{\mathit{p}},\\ {\overline{\mathit{x}}}_{n+1}^{i+1}={\overline{\mathit{x}}}_n+\mathrm{\Delta }{\overline{\mathit{u}}}_{n+1}^{i+1},\\ {\overline{\mathit{a}}}_{n+1}^{i+1}={\displaystyle \frac{1}{\mathrm{\Delta }{t}^2{\beta }_{\mathrm{\Delta }t}}}\mathrm{\Delta }{\overline{\mathit{u}}}_{n+1}^{i+1}-{\displaystyle \frac{1}{\mathrm{\Delta }t{\beta }_{\mathrm{\Delta }t}}}{\overline{\mathit{v}}}_n+\left(1-{\displaystyle \frac{1}{2{\beta }_{\mathrm{\Delta }t}}}\right){\overline{\mathit{a}}}_n,\\ {\overline{\mathit{v}}}_{n+1}^{i+1}={\overline{\mathit{v}}}_n+\mathrm{\Delta }t(1-{\gamma }_{\mathrm{\Delta }t}){\overline{\mathit{a}}}_n+\mathrm{\Delta }t{\gamma }_{\mathrm{\Delta }t}{\overline{\mathit{a}}}_{n+1}^{i+1}.\end{array}

6. Checking for convergence:

.

When not fulfilled, action returns to step 1, and i\leftarrow i+1.

The individual matrices and vectors are summarized in Box 3. Note that {\alpha }_{\mathrm{\Delta }t},{\beta }_{\mathrm{\Delta }t},\mathrm{a}\mathrm{n}\mathrm{d}{\gamma }_{\mathrm{\Delta }t} are time integration constants from the α-Bossak time integration scheme [55]. The convergence tolerances e_{\mathrm{u}}\mathrm{a}\mathrm{n}\mathrm{d}{e}_{\mathrm{p}} are taken as 10⁻⁸.

Box 3: Matrices and vectors used in the iterative solution scheme of the PFEM algorithm [16]

The matrices and vectors used in Box 2 are summarized in 3D as:

\begin{array}{rl}& {\mathit{F}}_{\mathrm{u},\mathrm{i}\mathrm{n}\mathrm{t}}=\stackrel{NF}{\underset{k=1}{A}}{\tilde{\mathit{B}}}^{k^{\mathrm{T}}}[{\mathit{\sigma }}^k{]}^{\mathrm{V}}{V}^k,\\ & {\mathit{F}}_{\mathrm{u},\mathrm{e}\mathrm{x}\mathrm{t}}=\stackrel{NE}{\underset{e=1}{A}}{\int }_{{\mathrm{\Gamma }}^{e,t}}{\mathit{N}}_u^{\mathrm{T}}{\mathit{t}}^{\ast }\mathrm{d}{\mathrm{\Gamma }}^{\mathrm{t}}+\stackrel{NE}{\underset{e=1}{A}}{\int }_{{\mathrm{\Omega }}^e}{\mathit{N}}_u^{\mathrm{T}}\rho \mathit{b}\mathrm{d}\mathrm{\Omega },\\ & {\mathit{M}}_{\mathrm{u}\mathrm{u}}=\stackrel{NE}{\underset{e=1}{A}}{\int }_{\mathrm{\Omega }e}{\mathit{N}}_{\mathrm{u}}^{\mathrm{T}}\rho {\mathit{N}}_{\mathrm{u}}\mathrm{d}\mathrm{\Omega },\\ & {\mathit{F}}_{\mathrm{p}}=\stackrel{NF}{\underset{k=1}{A}}{\mathit{N}}_{\mathrm{p}}^{k^{\mathrm{T}}}{\displaystyle \frac{\mathrm{l}\mathrm{n}\left(J\right)}{J^2}}{V}^k,\\ & {\mathit{K}}_{\mathrm{u}\mathrm{u},\mathrm{g}\mathrm{e}\mathrm{o}}=\stackrel{NF}{\underset{k=1}{A}}{\tilde{\mathit{G}}}^{k^{\mathrm{T}}}\left[{\mathit{\sigma }}^k\right]{\tilde{\mathit{G}}}^k{V}^k,\end{array}

\begin{array}{rl}& {\mathit{K}}_{\mathrm{u}\mathrm{u},\mathrm{m}\mathrm{a}\mathrm{t}}=\stackrel{NF}{\underset{k=1}{A}}{\tilde{\mathit{B}}}^{k^{\mathrm{T}}}{\mathit{c}}^k{\tilde{\mathit{B}}}^k{V}^k,\\ & {\mathit{K}}_{\mathrm{u}\mathrm{p}}=\stackrel{NF}{\underset{k=1}{A}}{\tilde{\mathit{B}}}^{k^{\mathrm{T}}}\mathit{m}{\mathit{N}}_{\mathrm{p}}^k{V}^k,\\ & {\mathit{K}}_{\mathrm{p}\mathrm{u}}=\stackrel{NF}{\underset{k=1}{A}}{\mathit{N}}_{\mathrm{p}}^{k^{\mathrm{T}}}{\displaystyle \frac{1-\mathrm{l}\mathrm{n}\left(J\right)}{J^2}}{\mathit{m}}^{\mathrm{T}}{\tilde{\mathit{B}}}^k{V}^k,\\ & {\mathit{K}}_{\mathrm{p}\mathrm{p},\mathrm{K}}=\stackrel{NE}{\underset{e=1}{A}}{\int }_{{\mathrm{\Omega }}^e}{\mathit{N}}_{\mathrm{p}}^T{\displaystyle \frac{1}{J}}{\displaystyle \frac{1}{\alpha K_0}}{\mathit{N}}_{\mathrm{p}}\mathrm{d}\mathrm{\Omega },\\ & {\mathit{K}}_{\mathrm{p}\mathrm{p},\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}=\stackrel{NE}{\underset{e=1}{A}}{\int }_{{\mathrm{\Omega }}^e}{\displaystyle \frac{{\alpha }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}}}{JG}}\left({\mathit{N}}_{\mathrm{p}}^{\mathrm{T}}{\mathit{N}}_{\mathrm{p}}-{\stackrel{⌣}{\mathit{N}}}_{\mathrm{p}}^{\mathrm{T}}{\stackrel{⌣}{\mathit{N}}}_{\mathrm{p}}\right)\mathrm{d}\mathrm{\Omega },\end{array}

where {\mathit{N}}_{\mathbf{u}} and {\mathit{N}}_{\mathrm{p}} are the linear shape functions of the displacement and pressure fields, A is the finite element assembly operator, NE and NF are the number of elements and faces, {\left[{\mathit{\sigma }}^k\right]}^{\mathrm{V}} are the Cauchy stresses in Voigt notation, {\alpha }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}} is a stabilization parameter [56], {\stackrel{⌣}{\mathit{N}}}_{\mathrm{p}}={\left[1/4,1/4,1/4,1/4\right]}^{\mathrm{T}} are the discontinuous polynomials and \mathit{m}={\left[1,1,1,0,0,0\right]}^{\mathrm{T}}. Note that the lumped mass matrix {\mathit{M}}_{\mathrm{u}\mathrm{u}} is used. The matrices {\tilde{\mathit{B}}}^k,{\tilde{\mathit{G}}}_I^k, and \left[{\mathit{\sigma }}^k\right] are defined as:

\begin{array}{c}{\tilde{\mathit{B}}}^k=\left[{\tilde{\mathit{B}}}_1^k,\cdots ,{\tilde{\mathit{B}}}_I^k,\cdots ,{\tilde{\mathit{B}}}_{N^k}^k\right],{\tilde{\mathit{G}}}^k=\left[{\tilde{\mathit{G}}}_1^k,\cdots ,{\tilde{\mathit{G}}}_I^k,\cdots ,{\tilde{\mathit{G}}}_{N^k}^k\right],\\ {\tilde{\mathit{B}}}_I^k=\left[\begin{array}{ccc}{\tilde{b}}_{I1}^k& 0& 0\\ 0& {\tilde{b}}_{I2}^k& 0\\ 0& 0& {\tilde{b}}_{I3}^k\\ {\tilde{b}}_{I2}^k& {\tilde{b}}_{I1}^k& 0\\ 0& {\tilde{b}}_{I3}^k& {\tilde{b}}_{I2}^k\\ {\tilde{b}}_{I3}^k& 0& {\tilde{b}}_{I1}^k\end{array}\right],{\tilde{\mathit{G}}}_I^k=\left[\begin{array}{ccc}{\tilde{b}}_{I1}^k& 0& 0\\ {\tilde{b}}_{I2}^k& 0& 0\\ {\tilde{b}}_{I3}^k& 0& 0\\ 0& {\tilde{b}}_{I1}^k& 0\\ 0& {\tilde{b}}_{I2}^k& 0\\ 0& {\tilde{b}}_{I3}^k& 0\\ 0& 0& {\tilde{b}}_{I1}^k\\ 0& 0& {\tilde{b}}_{I2}^k\\ 0& 0& {\tilde{b}}_{I3}^k\end{array}\right],\\ \left[{\mathit{\sigma }}^k\right]=\left[\begin{array}{ccc}{\mathit{\sigma }}^k& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathit{\sigma }}^k& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathit{\sigma }}^k\end{array}\right],\end{array}

with the number of nodes N^k connected with a face k (N^k=5 for an internal face and N^k=4 for a free surface face) and with the face-based smoothed derivatives of the shape functions {\tilde{b}}_{I1}^k [16], which is the volume average of the shape function derivatives of the elements sharing the same face when using linear shape functions [57].

3.3 Domain composition algorithm for contact surfaces

In PFEM, the contact between two domains is naturally accounted for by newly formed elements between the domains, which are obtained by the alpha-shape-based re-meshing algorithm when the domains are close enough [52]. A drawback is that these elements create a non-physical gap, between the domains, of the size of an element. In contrast, when using a constrained Delaunay triangulation algorithm, the domains cannot even make contact and a contact algorithm is required. Oliver et al. developed the contact domain method [58,59] to simulate contact problems when using a PFEM algorithm based on a constrained Delaunay triangulation algorithm. However, due to the fresh nature of the printed material in 3D concrete printing, it is reasonable to assume an ideal bond between the layers during the deposition process, and that no complex contact condition between the layers is necessary.

Instead of using the contact domain method, an algebraic-based tying algorithm is used to bond the layers together, which also has the advantage of condensing some degrees of freedom during the solution procedure. The adopted algorithm is based on a domain composition method with a Dirichlet–Neumann coupling of the surfaces using a transmission matrix obtained from a projection scheme, as described in Refs. [29,30]. By arbitrarily defining one target surface and one source surface of the two surfaces to be linked, the continuity condition is weakly defined over the target surface as:

(32){\int }_{\mathrm{\Gamma }} \delta {\mathit{u}}_{\mathrm{t}}{\mathit{u}}_{\mathrm{t}}\mathrm{d}\mathrm{\Gamma }={\int }_{\mathrm{\Gamma }} \delta {\mathit{u}}_{\mathrm{t}}{\mathit{u}}_{\mathrm{s}}\mathrm{d}\mathrm{\Gamma },

with the target surface \overline{\mathrm{\Gamma }}, the target test function \delta {\mathit{u}}_{\mathrm{t}}, the displacement on the target surface {\mathit{u}}_{\mathrm{t}} and the displacement on the source surface {\mathit{u}}_{\mathrm{s}}. In the spirit of PFEM, a finite element discretization of the surfaces with linear shape functions is used for the surface displacements, and these are defined for an element e, as:

(33){\mathit{u}}_{\mathrm{t}}={\mathit{N}}_{\mathrm{t}}{\overline{\mathit{u}}}_{\mathrm{t}}^{\mathrm{e}}\mathrm{a}\mathrm{n}\mathrm{d}{\mathit{u}}_{\mathrm{s}}={\mathit{N}}_{\mathrm{s}}{\overline{\mathit{u}}}_{\mathrm{s}}^{\mathrm{e}},

where {\mathit{N}}_{\mathrm{t}} denotes the shape functions of the target surface element, {\overline{\mathit{u}}}_{\mathrm{t}}^{\mathrm{e}} denotes the nodal displacements of the target surface element, {\mathit{N}}_{\mathrm{s}} denotes the shape functions of the source surface element and {\overline{\mathit{u}}}_{\mathrm{s}}^{\mathrm{e}} denotes the nodal displacements of the source surface element. Using the finite element discretization in Eq. (33), Eq. (32) can be recast as:

(34)\underset{{\mathit{A}}_{\mathrm{t}\mathrm{t}}}{\underbrace{\stackrel{NT}{\underset{e=1}{A}}\left({\int }_{{\overline{\mathrm{\Gamma }}}^e}{\mathit{N}}_{\mathrm{t}}^{\mathrm{T}}{\mathit{N}}_{\mathrm{t}}\mathrm{d}\mathrm{\Gamma }\right)}}{\overline{\mathit{u}}}_{\mathrm{t}}=\underset{{\mathit{A}}_{\mathrm{t}\mathrm{s}}}{\underbrace{\stackrel{NT}{\underset{e=1}{A}}\left({\int }_{{\overline{\mathrm{\Gamma }}}^e}{\mathit{N}}_{\mathrm{t}}^{\mathrm{T}}{\mathit{N}}_{\mathrm{s}}\mathrm{d}\mathrm{\Gamma }\right)}}{\overline{\mathit{u}}}_{\mathrm{s}},

with the finite element assembly operator A and the number of target surface elements NT. Equation (34) can be re-written as:

(35){\overline{\mathit{u}}}_{\mathrm{t}}=({\mathit{A}}_{\mathrm{t}\mathrm{t}}{)}^{-1}{\mathit{A}}_{\mathrm{t}\mathrm{s}}{\overline{\mathit{u}}}_{\mathrm{s}}=\mathit{T}{\overline{\mathit{u}}}_{\mathrm{s}},

where \mathit{T}={\left({\mathit{A}}_{\mathrm{t}\mathrm{t}}\right)}^{-1}{\mathit{A}}_{\mathrm{t}\mathrm{s}} and denotes the transmission matrix, which discretely links the target surface displacements with the source surface displacements. To simplify the inversion of {\mathit{A}}_{\mathrm{t}\mathrm{t}}, a lumped version of this matrix is used. The calculation of {\mathit{A}}_{\mathrm{t}\mathrm{s}} is more challenging as it requires an approximation of the source shape functions on the target surface. Various projection or interpolation schemes, summarized in Ref. [30], can be used to calculate {\mathit{A}}_{\mathrm{t}\mathrm{s}}. In the context of the PFEM, it was found that Gauss integration in conjunction with an orthogonal projection of the Gauss points onto the source surface was the most robust scheme [29]. The matrix {\mathit{A}}_{\mathrm{t}\mathrm{s}} can then be calculated as:

(36){\mathit{A}}_{\mathrm{t}\mathrm{s}}=\stackrel{NT}{\underset{e=1}{A}}\left({\sum }_{i=1}^{n_{\mathrm{g}\mathrm{p}}} w^i{\mathit{N}}_{\mathrm{t}}^{\mathrm{T}}\left({\mathit{\xi }}^i\right){\mathit{N}}_{\mathrm{s}}\left({\mathrm{\Pi }}_{\mathrm{s}}\right({\mathit{\xi }}^i\left)\right)\left|{\mathit{J}}^{\mathrm{e}}\right({\mathit{\xi }}^i\left)\right|\right),

with the number n_{\mathrm{g}\mathrm{p}} of Gauss points, the natural Gauss point coordinates \mathit{\xi }, the Gauss point weight w, the Jacobian {\mathit{J}}^e and the natural coordinates of the Gauss point projection from the target onto the source surface {\mathrm{\Pi }}_{\mathrm{s}}\left({\mathit{\xi }}^i\right); see Fig.4. To calculate Eq. (36), the active target surface elements must be identified. A target surface element is considered active, when an adjacent node has slightly penetrated into the source domain, which is checked in each time step. For all active target surface elements, the Gauss points are projected onto the source surface to find the source surface shape function values for each Gauss point. The number of source surface elements tied with a target surface is generally unknown and depends on the discretization (nodal spacing, element shapes, etc.) and the number of Gauss points, which should be large enough to allow an accurate integration of Eq. (36). In this work, three Gauss points for 3D triangular surface elements were found to yield sufficiently good results. The nodal spacing of the target surface discretization should also be approximately the same as or finer than that of the source surface discretization to yield accurate results.

The target surface degrees of freedom can be condensed from the solution system using the transmission matrix \mathit{T}, which can be calculated using Eqs. (34)–(36). In the PFEM, a modified transmission matrix \tilde{\mathit{T}} is calculated, which links the condensed solution system with the full solution system. The modified transmission matrix \tilde{\mathit{T}} is of size \left(n_{\mathrm{d}\mathrm{o}\mathrm{f}}-n_{\mathrm{d}\mathrm{o}\mathrm{f},\mathrm{t}}\right) times n_{\mathrm{d}\mathrm{o}\mathrm{f}}, where n_{\mathrm{d}\mathrm{o}\mathrm{f}} is the number of all degrees freedom and n_{\mathrm{d}\mathrm{o}\mathrm{f},\mathrm{t}} is the number of target surface degrees of freedom. An identity relation in \tilde{\mathit{T}} is used to guarantee that the degrees of freedom not in contact remain unchanged. The solution algorithm summarized in Box 2 is modified by using the condensed stiffness matrices {\tilde{\mathit{K}}}_{\mathrm{u}\mathrm{u}}={\tilde{\mathit{T}}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{u}\mathrm{u}}\tilde{\mathit{T}}, {\tilde{\mathit{K}}}_{\mathrm{u}\mathrm{p}}={\tilde{\mathit{T}}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{u}\mathrm{p}} and the condensed displacement residual {\tilde{\mathit{R}}}_{\mathrm{u}}={\tilde{\mathit{T}}}^{\mathrm{T}}{\mathit{R}}_{\mathrm{u}}. After solving the condensed solution system, the displacement increments \mathrm{\Delta }\hat{u} are obtained from the condensed displacement increments \mathrm{\Delta }\tilde{u} via \mathrm{\Delta }\hat{u}=\tilde{\mathit{T}}\mathrm{\Delta }\tilde{u}.

The PFEM solution procedure of one time step can be summarized as.

1) Free surface refinement and constrained Delaunay triangulation.

2) Calculation of the modified transmission matrix \tilde{\mathit{T}}.

3) Calculation of {\alpha }_{n+1}.

4) Iterative solution of the balance of momentum and mass (Box 3).

To obtain a more robust and efficient algorithm, the modified transmission matrix is only calculated once during each time step and remains constant during the iterative solution scheme. Due to small time steps used in PFEM algorithms, an accurate representation of the tied surfaces is still guaranteed.

4 Computational investigation of 3D concrete printing on the scale of three printed layers

This section describes the computational investigation of 3D concrete printing extrusion and deposition processes of multiple printed layers. The presented investigation includes the analysis of the layer shapes of three printed layers, the forces acting on substrate layers and the resulting deformations in substrate layers for different process parameters, such as printing speed, time gap between the layers, standoff distance, flow rate and nozzle diameter. As the correlations between process parameters and deformations or forces are highly complex and nonlinear, design tools are required to control the deformations during 3D concrete printing. Therefore, one aim is to generate deformation plots, which can be used to find optimal process parameters. In addition, the simulations demonstrate the effectiveness and applicability of the proposed numerical model for 3D concrete printing processes on the scale of a few layers. In Section 5, the applicability of the proposed model is also demonstrated to simulate 3D concrete printing processes on the structural scale based on the analysis of a 3D-printed box.

4.1 Simulation outline

The computational investigation is based on the simulation of three printed concrete layers with varying process parameters. Fig.5 shows the geometry of a layer, of length l, during the printing process with the definitions of the process parameters, including the diameter of the round extrusion nozzle d, the standoff distance between the substrate and the nozzle outlet h, the printing speed v_{\mathrm{p}}, and the extrusion rate Q. The numerical simulations are conducted under the following assumptions.

1) A fresh on fresh contact between the layers is assumed, which is modeled with ideal no-slip tying conditions between the layers and the substrate.

2) The printed material is assumed as homogeneously sheared after extrusion and the dynamic yield stress is taken as the initial yield stress value after extrusion.

3) Destructuration processes are assumed negligible during deposition processes.

4) Plastic shrinkage, cracking, creep and desiccation are assumed irrelevant during the studied time period.

5) Flow processes within the extruder are neglected and a constant inlet velocity is assumed at the nozzle outlet.

In previous studies [15,16], the constant inlet velocity assumption was found to be suitable for 3D concrete printing, which resulted in accurate predictions of the layer shape. As in those previous studies, the inlet flow is realized using frequently new created elements, which are pushed into the domain. Accordingly, the nodes at the nozzle outlet are fixed by the inlet velocity plus the printing speed, and a row of new elements is created every couple of time steps when the old row of elements has been completely pushed into the domain. This procedure leads to an irregular periodical behavior at the nozzle outlet, which, however, has no influence on the printed cross section or the steady-state forces obtained from the averaged values over multiple time steps; this point was analyzed in preliminary studies. As depicted in Fig.5, the flow rate Q is scaled with the printing speed, the standoff distance and the nozzle diameter via Q=qhr{v}_{\mathrm{p}}, with the flow rate factor q. When q>1 it is guaranteed that the available space under the extrusion nozzle is always completely filled. In most of the conducted simulations q is taken as 1.4, which theoretically results in a layer width of 1.4·d, where d is the nozzle diameter. According to the definition in Ref. [60], the printing process lies in the controlled pressing regime. An overview of the parametric studies is provided in Tab.1. The process parameters under investigation are the nozzle diameter d, the standoff distance h, the flow rate factor q, the printing speed v_{\mathrm{p}} and the time gap t_{\mathrm{g}}, which is the time gap between printing of each layer. These process parameters are systematically analyzed in the next sections. The simulations are carried out with the material properties, as provided in Subsection 2.8, using a density of {\rho }_0=2078\mathrm{k}\mathrm{g}/{\mathrm{m}}^3. In all simulations, symmetry across the xy-plane in the layer center was taken into account and a time step size of \mathrm{\Delta }t=0.0002·5/v_{\mathrm{p}} (with v_{\mathrm{p}} in cm/s) was used. A nodal spacing of around d/13 resulted in approximately {10}^5 elements when all layers were printed. The typical simulation time of one scenario was around 3 d on a standard computer.

4.2 Post-processing strategy

In this section, the post-processing strategies regarding how the results are processed, presented and analyzed are outlined. The analysis of the results includes the investigation of the printed cross sectional layer shapes at the time instances when one two and three layers were printed. This is required due to deformations in lower layers as a result of the successively increasing loads of subsequently printed layers. The deformations are evaluated in the lowest layer as the relative horizontal deformations via:

(37){\mathit{\text{ϵ}}}_{\mathrm{w}}=\frac{w_j-w_1}{w_1},

where w_1 is the width of the lowest layer at the time instances when one layer was printed and w_j is the width of the lowest layer at the time instances when two (j=2) or three (j=3) layers were printed. The deformations in substrate layers are driven by the forces of subsequently printed layers in additive manufacturing deposition processes. In this work, the impact of these forces was measured using the vertical contact tractions (force per area) between two layers or between the layer and the substrate. Note that the vertical contact tractions can only be interpreted as contact pressure when they act normal to the contact surface, which is not necessarily guaranteed. It was shown in Ref. [15] that these tractions can be decomposed into self weight and the tractions acting on the substrate layer generated by the material extrusion under the extrusion nozzle. In Ref. [15], these tractions were mainly analyzed in 2D for one layer printing tests with a purely viscoplastic model. In the following, computational analyses of 3D multi-layer printing are discussed based on the proposed elasto-viscoplastic material formulation including evolving material properties, which allows a more realistic assessment of the deformations. To simplify the interpretation of the results from the 3D analysis and to draw meaningful conclusions, a new post-processing strategy is required. Fig.6 illustrates the extraction of the vertical contact tractions between the layers in a three-layer printing process. These tractions are further averaged over the layer width and over multiple time steps and plotted over the distance to the extrusion nozzle in Fig.6 right, to allow a more convenient interpretation of the results. Under the assumption that vertical tractions dominate the stress state in substrate layers, viscoplastic deformations are expected when the yield criterion is satisfied:

(38)\mathrm{m}\mathrm{a}\mathrm{x}\left(\text{Tractions}\right)>\sqrt{3}{\tau }_{0,\mathrm{s}\mathrm{u}\mathrm{b}},

where {\tau }_{0,\mathrm{s}\mathrm{u}\mathrm{b}} is the yield stress in the substrate layer. Based on the simplified yield criterion in Eq. (38) the deformations in substrate layers can be characterized as purely elastic or elasto-viscoplastic directly below the extrusion nozzle. As depicted in Fig.6 right, the self-weight contribution to the tractions is constant along the layer and given as \rho gh. The contribution of the extrusion tractions can be obtained by subtracting the self weight contribution from the total tractions. The contribution of the extrusion tractions decreases (and ultimately becomes zero) with increasing distance from the extrusion nozzle and has its peak (maximum negative value) right under the nozzle, as shown in Fig.6, right. The negative traction values imply that they act as a compressive traction on the lower layers.

4.3 Analysis

The analyses of the parametric studies are organized in separate subsections focusing on the individual process parameters. The simulations associated with the analysis of a certain process parameter are grouped as indicated in Tab.1. First, simulation 1 is analyzed, which serves as a reference simulation in this computational study. Fig.7(a) shows the side view on the printing process during the time when the third layer is printed. Based on the equivalent Kirchhoff stresses on the symmetry plane depicted in Fig.7(b), a stress concentration under the extrusion nozzle in substrate layers can be observed due to the tractions from self-weight and extrusion. These tractions are large enough for the substrate layer to yield locally under the extrusion nozzle, leading to viscoplastic deformations as indicated by the visualization of the yield criterion in Fig.7(c) (yielded region in red color). The observation that viscoplastic deformations only occur in the currently printed layer and the substrate layer is supported by analyzing the printed cross section after printing one, two and three layers, as shown in Fig.8(a). Larger deformations can be observed in the first layer when the second layer is printed, indicating viscoplastic deformations. Only slight deformations are found in the first layer when the third layer is printed indicating that only elastic deformations occur in the first layer.

The observations above are consistent with the simplified yield criterion in Eq. (38), using the maximum value of the averaged contact tractions taken from Fig.8(b). Fig.8(b) shows the averaged tractions relative to the distance to the extrusion nozzle, which exceed the yield point (\sqrt{3}{\tau }_{0,\mathrm{d}}) right under the extrusion nozzle. Note that the influence of aging on the yield stress is negligible when the time gap t_{\mathrm{g}} between the layers is zero (simulations 1–19) and a substrate yield stress of {\tau }_{\mathrm{s}\mathrm{u}\mathrm{b}}\approx {\tau }_{\mathrm{d}}\approx 0.56\mathrm{k}\mathrm{P}\mathrm{a} can be assumed. Fig.8(b) also depicts the averaged tractions when printing the first, second and third layer. The tractions decrease with the number of layers, which is due to the stiffer properties of the substrate when fewer substrate layers exist. The analysis of the reference simulation 1 is concluded by the fact that viscoplastic deformations can occur locally under the extrusion nozzle in the substrate due to tractions from self-weight and extrusion. The simplified yield criterion in Eq. (38) is found to be useful in assessing whether these viscoplastic deformations occur or not.

4.3.1 Influence of the nozzle diameter d

The printed cross sections using different nozzle diameters d at the time instances when one (blue), two (red) and three layers (green) were printed are provided in Fig.8(a) and in Fig.9(a)–Fig.9(d) obtained from simulations 1–5. As can be seen, layers with a larger width are printed for larger extrusion nozzle diameters due to more available space under the nozzle that must be filled. In addition, larger deformations occurred in the substrate layers when using larger extrusion nozzles, which is confirmed by evaluation of the relative deformations in width direction based on Eq. (37), shown in Fig.10(a). Comparably small differences are observed between the deformations after printing two and three layers, which indicates that only elastic deformations are generated in the first layers when the third layers were printed. The increase in tractions for larger extrusion nozzles is illustrated in Fig.10(b), which is attributed to the more constrained situation under the extrusion nozzle when using larger nozzle diameters. Note that the time interval of a material point in the substrate layer being exposed to the extrusion tractions is also larger for larger nozzle diameters, which supports the build-up of viscoplastic deformations. The simplified yield criterion as illustrated in Fig.10(b) is exceeded for all nozzle diameters under the extrusion nozzle indicating that viscoplastic deformations occur in substrate layers. However, the simplified yield criterion is only slightly exceeded for the two lowest nozzle diameters, resulting in only small deformations; see Fig.10(a). Note that the averaged tractions over the width are shown and that the tractions in the center of the layers can also be larger. The kinky shape of the blue line (d = 1.5 cm) in Fig.10(b) is attributed to the mesh discretization and the interpolation of the tractions during post-processing. Possibly, the use of adaptive refinement schemes could yield smoother curves. Such schemes, however, are not yet implemented in the present PFEM version. Based on Fig.10(a), the expected deformations in the printed layer can be assessed when choosing a suitable nozzle diameter for a printing process. Such a plot can be used as a design tool for 3D concrete printing.

4.3.2 Influence of the standoff distance h

The results of the printed cross sections with varying standoff distances h after one, two and three layers were printed are shown in Fig.11(a)–Fig.11(c) obtained from simulations 6–8. The deformations show no clear correlation with the standoff distance, indicating that the smallest deformations occur for the medium standoff distance h = 2 cm. The investigation of the relative deformations in Fig.12(a) reveals a decrease of the deformations with increasing standoff distance after two layers were printed, but a sudden increase of the deformations for the largest standoff distance after three layers were printed. The apparent paradoxical situation can be explained based on the average tractions. While the contribution to the tractions from self weight decreases with decreasing standoff distances, the contribution from extrusion increases and dominates the overall tractions under the extrusion nozzle (Fig.12(b)). The increase in tractions under the extrusion nozzle explains the increase in deformations with decreasing standoff distances when only two layers were printed. The sudden increase in deformation for the largest standoff distance when three layers were printed is attributed to a global plastic collapse of the lowest layer due to the large dead weight of the upper layers when using a larger standoff distance. Global plastic collapse cannot be identified from the averaged contact tractions and must be checked based on the total dead weight of the upper layers [10].

4.3.3 Influence of the flow rate factor q and the printing speed v_{\mathrm{p}}

The relative horizontal deformations when using different flow rate factors q are depicted in Fig.13(a) obtained from simulations 1, 9–11 (Tab.1). As could be expected, the deformations increase with increasing flow rate, and this is attributed to the increasing contact tractions acting on the substrate layer under the extrusion nozzle (Fig.13(b)). The increase in flow rate leads to larger extrusion velocities, which manifests itself in larger tractions acting on the substrate layer. Similarly, larger printing speeds also result in larger extrusion velocities, which lead to larger tractions acting on the substrate layer and to larger deformations (Fig.13(c) and Fig.13(d)). The results for different printing speeds v_{\mathrm{p}} were obtained from simulations 1, 12, 13, 14, and 15 according to Tab.1. While the increase in the peak values of the average tractions is in the same range for the increasing flow rate and the increasing printing velocity, as shown in Fig.13(b) and Fig.13(d), the resulting deformations are smaller for increasing printing velocities (Fig.13(a) and Fig.13(c)). The reason for the smaller displacements is the smaller time interval during which a material point in the substrate layer is exposed to the extrusion tractions, for larger printing speeds. The observed viscoplastic deformations result from opposing influences of the maximum contact tractions and the effective time interval during which a material point is exposed to the tractions. In the investigation above, the maximum contact tractions were the dominating effect, which is not necessarily always the case. As only few parametric studies have been carried out here, so that an ad hoc statement on the dominating effect is not possible.

4.3.4 Influence of the nozzle diameter d and standoff distance h with h/d = constant

The influence of the nozzle diameter d and standoff distance h are analyzed for the case when the ratio h/d = 0.5 and is constant based on simulations 1, 16–19. Fig.14(a) shows the relative deformations, which increase with an increasing nozzle diameter and standoff distance. The average tractions in Fig.14(b) reveal an increasing trend with an increasing nozzle diameter and standoff distance resulting in the larger deformations. As is the case for the study on different standoff distances based on Fig.12(a) and Fig.12(b), the contributions of the tractions from self-weight and extrusion vary inversely with respect to the nozzle diameter and standoff distance for h/d constant at 0.5. While the self weight increases with an increasing nozzle diameter and standoff distance, the contribution from extrusion decreases. Fig.14(c) depicts the extrusion tractions without the contribution from self weight over the scaled distance x/d to the extrusion nozzle. The extrusion tractions increase with decreasing nozzle diameter and standoff distance due to larger stresses in smaller extrusion nozzles [15]. As is observed in Fig.12(b), the influence of self weight dominates the surface tractions under the extrusion nozzle. This is in contrast to the study on the influence of standoff distances, where the contribution from extrusion dominated; see Fig.12(b). The above investigations show that the resulting deformations and tractions in 3D concrete printing are the result of various partially opposing effects. Which of these effects dominates the overall response is a challenging question that in our opinion can only be answered based on numerical simulations.

4.3.5 Influence of the time gap t_{\mathrm{g}}

To control the deformations in 3D concrete printing, substrate layers must be strong enough to carry the loads from extrusion and self-weight. Time gaps between the layers are required so that the printed concrete can age and get stronger. Fig.15(a) shows the relative deformations over different time gaps based on simulations 1, 20–24. As can be observed, the deformations decrease with an increasing time gap. At a time gap of around 3 min, the deformations stay almost constant leading to the conclusion that only elastic deformations occur in substrate layers. The simplified yield criterion in Eq. (38) can also be used to check whether only elastic deformations occur. Using a time gap of 3 min, the substrate yield stress is {\tau }_{\mathrm{s}\mathrm{u}\mathrm{b}}\approx 0.93\mathrm{k}\mathrm{P}\mathrm{a} resulting in \sqrt{3}{\tau }_{\mathrm{s}\mathrm{u}\mathrm{b}}\approx 1.61\mathrm{k}\mathrm{P}\mathrm{a}, which easily exceeds the maximum average tractions under the extrusion nozzle (Fig.15(b)). Based on the simplified yield criterion, no viscoplastic deformations occur in substrate layers. Note that only the average tractions over the width are shown and that the actual tractions are slightly larger. As can be seen in Fig.15(b), the tractions slightly increase with an increasing time gap, which is attributed to the stiffer substrate material properties.

4.3.6 Discussion

The previous analyses reveal that the computational model can give adequate results for individual printing scenarios for three layers. In particular, viscoplastic and elastic deformations can be analyzed and the effect of aging on the structural response is accounted for. The analysis of the results shows that the interactions between the process parameters and the resulting deformations and layer shapes are highly complex and nonlinear, so that ad hoc design rules cannot be specified. Numerical analyses can only provide answers for individual printing scenarios. The deformation plots based on the various process parameters as given in Fig.10(a), Fig.12(a), Fig.13(a), Fig.13(c), Fig.14(a), and Fig.15(a) have the potential of becoming suitable design tools to find optimal process parameters. However, a larger database for several combinations of printing cases and parameters would be helpful to fully cover the space of relevant parameters in 3D printing processes.

5 Computational investigation of 3D concrete printing on the structural scale

Structural collapse during printing is a major issue in 3D concrete printing due to typically very thin structures in combination with the soft nature of the freshly printed material [10,37]. Computational simulations can be used to study the mechanisms of structural collapse during 3D concrete printing. In this section, computational simulations are presented to demonstrate the capability of the proposed model to simulate structural collapse during 3D concrete printing processes on the structural scale. These studies also reveal that using a realistic (in this case simulated) instead of an idealized rectangular shape of the printed layer can influence the stability of the structure during printing. Such influences can only be captured by knowing the actual shape of the layer, and such influences are directly assessed in the simulations discussed in Section 4.

A set of simulations on the structural scale was performed with the aim of numerically determining the effect of the printing speed on the maximum number of layers that can be printed before failure occurs. A simple hollow structure with a rectangular cross section with rounded corners (shown in Fig.16(a)) is considered. A 3D rendering of the final form is shown in Fig.16(b). The same material and process parameters as in the reference model in Section 4 are adopted; however, several printing speeds are examined (as shown in Tab.2).

Using the same solution framework presented in Section 4 to run this set of simulations is unfeasible due to the immense computational time required. Therefore, it is common in practice to use the element-activation strategy to run such simulations on this scale, where the finite element mesh for the entire model is a priori generated. At the beginning of the simulation all elements are deactivated, then at each new time-step a new set of elements, which represents a part that is being printed in this time-step, is activated [24,25]. In this strategy an ideal and uniform shape of the printed layers is usually assumed to generate the mesh of the model. Evidently, the effect of the deposition process currently is not taken into account in this modeling strategy. However, as was shown in Section 4, the printing of new layers significantly affects the shape of the substrate layers.

Hence, we propose the following simulation scheme: We divide the model into different components related to different printing paths. In this case, the structure is characterized by a straight and a curved path (as shown in red and green in Fig.16(c), respectively). Then we run a multiple-layers deposition simulation of a model that has two unique paths (Fig.16(d)) to capture the deformations in the substrate layers. For this deposition simulation, we use the same computational framework and parameters as described in Section 4. We observe, in this case, that the shape of the most substrate layer converges, i.e., stops changing due to the effects of extrusion and deposition, after printing of 4 layers on top of it. Therefore, we run this simulation for 5 layers first. Then, the shapes of the layers obtained from this simulation are used as the input geometry to generate the mesh for the full-scale model which subsequently will be simulated using the element-activation strategy. By adopting this approach, the impact of elastic deformations on the printed layer shape is neglected when more than 5 layers are activated. Nonetheless, this influence is very small and can be considered negligible.

To demonstrate the effect of the actual shape of the layers on the global behavior of the structure during printing, the same element-activation simulations were performed with an idealized uniform shape of the layers assuming that they have a rectangular cross section. The results of our simulations are also compared with the analytical solution of the elastic buckling case of a straight 3D-printed wall element proposed in Ref. [61]. Tab.2 summarizes the obtained results.

Fig.17 compares the maximum number of printed layers before collapse for the cases of numerical models with layers of idealized shape and simulated shape and for the analytical solution. An overestimation in the buildability is observed for the case of idealized-shaped layers and the analytical solution in comparison to the results from the computational simulations based on realistic shapes of the layers. The analytical solution assumes a rectangular cross section of the layers and is limited to failure only due to linear buckling. The numerical model of layers assuming an ideal shape also assumes a rectangular cross section; however, it does include an elastoplastic material model that accounts for the aging effect. Therefore, good agreement with the analytical solution is observed except for the case of a relatively high printing speed (0.1 cm/s), where the structure fails due to plastic strains, which is not accounted for in the analytical solution.

The numerical models with the simulated realistic cross section, although having the same material model as the models based on the rectangular cross section, show less capacity in terms of number of printable layers. This is a result of the geometry, as the quasi-oval shape of the layers results in less contact area between the layers in comparison to rectangular-shaped layers, and therefore, earlier buckling is to be expected. Fig.18 shows the buckling of model 1 (printing speed 0.005 cm/s) upon printing the 15th layer.

The results discussed indicate that the proposed model can adequately model structural failure of 3D concrete printing on the structural scale. The presented approach shows a first attempt on how simulations on the layer and structural scale can be combined by using the proposed unified fluid and solid mechanics material formulation. The relevance of such an approach is supported by the results, revealing an influence of the actual or realistic layer shape on the structural stability during 3D concrete printing.

6 Conclusions

In this work, a novel elasto-viscoplastic constitutive model for fresh concrete, during aging, was proposed to simulate 3D concrete printing processes using the PFEM. Computational investigations reveal the capability of the model in simulating extrusion and deposition processes on the scale of a layer, and also structural collapse during printing on the scale of the structure. The proposed model represents the first step towards a unified fluid and solid mechanics formulation that can be used for any age of the printed concrete and for every production step in 3D concrete printing.

The material formulation is based on a multiplicative split of the deformation gradient into elastic, aging and viscoplastic parts. The aging deformation gradient adopted a hyperelastic potential with evolving material properties. A Bingham viscoplastic formulation was used to simulate viscoplastic flow during the material extrusion in 3D concrete printing. Due to the existence of a hyperelastic potential, thermodynamic consistency of the material model was easily accounted for. One advantage of the model is the stress update algorithm, which is similar to that in small strain plasticity due to the adoption of a hyperelastic potential with logarithmic strains. The stress update algorithm was then implemented in a face-based smoothed PFEM model to simulate 3D concrete printing processes of multiple layers. To allow a more accurate representation of the contact condition between the layers, a novel PFEM tying algorithm was introduced to tie non-matching layer discretizations during 3D concrete printing together. The tying algorithm is based on a Dirichlet–Neumann formulation using a transmission matrix to condense the tied degrees of freedom from the solution system.

The computational model was then used to simulate 3D concrete printing extrusion and deposition processes of three layers using experimentally calibrated models from the literature for time evolving material properties. Parametric studies on the process parameters revealed the influence of various process parameters on the layer shapes, deformations, and the contact tractions between the layers. It was shown that viscoplastic deformations may occur immediately below the extrusion nozzle in substrate layers due to tractions from self weight and extrusion. The average contact tractions over the layer width under the extrusion nozzle acting on the substrate layer were then used to formulate a simplified yield criterion to characterize whether a substrate layer deforms plastically. In this context, the contact tractions were also decomposed into contributions from self-weight and extrusion, which were found to inversely correlate with the size of the layer. Smaller layers have a smaller contribution from self-weight, but a higher contribution from material extrusion due to the more constrained situation under the nozzle. The resulting deformations result from opposing tendencies of the maximum contact tractions and the time interval of a material point being exposed to the tractions. This opposition varies with the nozzle size and the printing speed. Which of the effects above dominates the overall response cannot be stated ad hoc. The deformations and contact tractions in 3D concrete printing are the result of a complex interaction between the influences of the material and process parameters, which can only be analyzed using numerical models. From several parametric studies of the process parameters, deformation plots were generated that can be used as design tools for engineers to control the deformations generated during extrusion and deposition.

In a final computational investigation, structural collapse was analyzed based on a reference case of a 3D printed rectangular hollow box. As simulations of the whole structure are computationally too expensive, an element activation-based simulation strategy was proposed, which took a realistic printed layer shape into account. The realistic printed layer shape was obtained from extrusion and deposition process simulations of a few layers. Results from analyses with a varying printing speed indicated that the structures simulated using the realistic printed layer shape had a lower capacity with respect to structural collapse than the structures simulated using an idealized rectangular shape. The lower capacity was attributed to the smaller contact surface between the layers when using the realistic printed layer shape.

The proposed numerical model using the novel constitutive formulation has been shown to be capable of adequately simulating multiple layers in 3D concrete printing extrusion and deposition processes, as well as simulations of structural collapse during 3D concrete printing. In the future, the constitutive model will be extended to account for creep and shrinkage phenomena using a physical description of the concrete aging. Since the time dependent material response of printed fresh concrete is quickly dominated by hydration after deposition, a formulation based on the degree of hydration needs to be considered

References

[1] Khoshnevis B. Automated construction by contour crafting-related robotics and information technologies. Automation in Construction, 2004, 13(1): 5–19

[2] Lim S, Buswell R A, Le T T, Austin S A, Gibb A G F, Thorpe T. Developments in construction-scale additive manufacturing processes. Automation in Construction, 2012, 21: 262–268

[3] Bos F P, Wolfs R J M, Ahmed Z Y, Salet T A M. Additive manufacturing of concrete in construction: Potentials and challenges of 3D concrete printing. Virtual and Physical Prototyping, 2016, 11(3): 209–225

[4] MechtcherineVNerellaV N. 3D-Druck mit Beton: Sachstand, entwicklungstendenzen, herausforderungen. Bautechnik, 2018, 95(4): 275−287 (in German)

[5] de Schutter G, Lesage K, Mechtcherine V, Nerella V N, Habert G, Agusti-Juan I. Vision of 3D printing with concrete-technical, economic and environmental potentials. Cement and Concrete Research, 2018, 112: 25–36

[6] Menna C, Mata-Falcon J, Bos F, Vantyghem G, Ferrara L, Asprone D, Salet T, Kaufmann W. Opportunities and challenges for structural engineering of digitally fabricated concrete. Cement and Concrete Research, 2020, 133: 106079

[7] El-Sayegh S, Romdhane L, Manjikian S. A critical review of 3D printing in construction: Benefits, challenges, and risks. Archives of Civil and Mechanical Engineering, 2020, 20(2): 34

[8] Mechtcherine V, Bos F P, Perrot A, Leal da Silva W R, Nerella V N, Fataei S, Wolfs R J, Sonebi M, Roussel M N. Extrusion-based additive manufacturing with cement-based materials—Production steps, processes, and their underlying physics: A review. Cement and Concrete Research, 2020, 132: 106037

[9] Roussel N, Spangenberg J, Wallevik J, Wolfs R J M. Numerical simulations of concrete processing: From standard formative casting to additive manufacturing. Cement and Concrete Research, 2020, 135: 106075

[10] Wolfs R J M, Suiker A S J. Structural failure during extrusion-based 3D printing processes. International Journal of Advanced Manufacturing Technology, 2019, 104(1−4): 565–584

[11] Liu Z, Li M, Weng Y, Qian Y, Wong T N, Tan M J. Modelling and parameter optimization for filament deformation in 3D cementitious material printing using support vector machine. Composites. Part B, Engineering, 2020, 193: 108018

[12] Comminal R, Leal da Silva W R, Andersen T J, Stang H, Spangenberg J. Modelling of 3D concrete printing based on computational fluid dynamics. Cement and Concrete Research, 2020, 138: 106256

[13] Wolfs R J M. Salet T A M, Roussel N. Filament geometry control in extrusion-based additive manufacturing of concrete: The good, the bad and the ugly. Cement and Concrete Research, 2021, 150: 106615

[14] Reinold J, Meschke G. Particle finite element simulation of fresh cement paste-inspired by additive manufacturing techniques. Proceedings in Applied Mathematics and Mechanics, 2019, 19(1): e201900198

[15] Reinold J, Nerella V N, Mechtcherine V, Meschke G. Extrusion process simulation and layer shape prediction during 3D-concrete-printing using the particle finite element method. Automation in Construction, 2022, 136: 104173

[16] Reinold J, Meschke G. A mixed U-P edge-based smoothed particle finite element formulation for viscous flow simulations. Computational Mechanics, 2022, 69(4): 891–910

[17] Liu Z, Li M, Tay J W D, Weng Y, Wong T N, Tan M J. Rotation nozzle and numerical simulation of mass distribution at corners in 3D cementitious material printing. Additive Manufacturing, 2020, 34: 101190

[18] He L, Chow W T, Li H. Effects of interlayer notch and shear stress on interlayer strength of 3D printed cement paste. Additive Manufacturing, 2020, 36: 101390

[19] Mollah M T, Comminal R, Serdeczny M P, Pedersen D B, Spangenberg J. Stability and deformations of deposited layers in material extrusion additive manufacturing. Additive Manufacturing, 2021, 46: 102193

[20] Spangenberg J, Leal da Silva W R, Comminal R, Mollah M T, Andersen T J, Stang H. Numerical simulation of multi-layer 3D concrete printing. RILEM Technical Letters, 2021, 6: 119–123

[21] Pan T, Teng H, Liao H, Jiang Y, Qian C, Wang Y. Effect of shaping plate apparatus on mechanical properties of 3D printed cement-based materials: Experimental and numerical studies. Cement and Concrete Research, 2022, 155: 106785

[22] Mollah M, Comminal R, Serdeczny M P, Seta B, Spangenberg J. Computational analysis of yield stress buildup and stability of deposited layers in material extrusion additive manufacturing. Additive Manufacturing, 2023, 71: 103605

[23] PanTGuoRJiangYJiX. How do the contact surface forces affect the interlayer bond strength of 3D printed mortar? Cement and Concrete Composites, 2022, 133: 104675

[24] Wolfs R J M, Bos F P, Salet T A M. Early age mechanical behaviour of 3D printed concrete: Numerical modelling and experimental testing. Cement and Concrete Research, 2018, 106: 103–116

[25] Vantyghem G, Ooms T, de Corte W. VoxelPrint: A Grasshopper plug-in for voxel-based numerical simulation of concrete printing. Automation in Construction, 2021, 122: 103469

[26] Chang Z, Zhang H, Liang M, Schlangen E, Savija B. Numerical simulation of elastic buckling in 3D concrete printing using the lattice model with geometric nonlinearity. Automation in Construction, 2022, 142: 104485

[27] Chang Z, Liang M, Xu Y, Schlangen E, Savija B. 3D concrete printing: Lattice modeling of structural failure considering damage and deformed geometry. Cement and Concrete Composites, 2022, 133: 104719

[28] Reinold J, Meschke G. Algorithm for aging materials with evolving stiffness based on a multiplicative split. Computer Methods in Applied Mechanics and Engineering, 2022, 397: 115080

[29] de Boer A, van Zuijlen A, Bijl H. Review of coupling methods for non-matching meshes. Computer Methods in Applied Mechanics and Engineering, 2007, 196(8): 1515–1525

[30] Houzeaux G, Cajas J C, Discacciati M, Eguzkitza B, Gargallo-Peiró A, Rivero M, Vázquez M. Domain decomposition methods for domain composition purpose: Chimera, overset, gluing and sliding mesh methods. Archives of Computational Methods in Engineering, 2017, 24(4): 1033–1070

[31] Esposito L, Casagrande L, Menna C, Asprone D, Auricchio F. Early-age creep behaviour of 3D printable mortars: Experimental characterisation and analytical modelling. Materials and Structures, 2021, 54(6): 207

[32] Moelich G, Kruger J, Combrinck R. Plastic shrinkage cracking in 3D printed concrete. Composites. Part B, Engineering, 2020, 200: 108313

[33] Nedjar B. On a geometrically nonlinear incremental formulation for the modeling of 3D concrete printing. Mechanics Research Communications, 2021, 116: 103748

[34] Nedjar B. Incremental viscoelasticity at finite strains for the modelling of 3D concrete printing. Computational Mechanics, 2022, 69(1): 233–243

[35] Lee E H. Elastic-plastic deformation at finite strains. Journal of Applied Mechanics, 1969, 36(1): 1–6

[36] Reese S, Govindjee S. A theory of finite viscoelasticity and numerical aspects. International Journal of Solids and Structures, 1998, 35(26−27): 3455–3482

[37] Roussel N. Rheological requirements for printable concretes. Cement and Concrete Research, 2018, 112: 76–85

[38] RousselNGramA. Simulation of Fresh Concrete Flow. Technical Report. Dodrecht: Springer Netherlands. 2014

[39] Papanastasiou T C. Flows of materials with yield. Journal of Rheology, 1987, 31(5): 385–404

[40] Perzyna P. Fundamental problems in viscoplasticity. Advances in Applied Mechanics, 1966, 10(2): 243

[41] Simo J C. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Computer Methods in Applied Mechanics and Engineering, 1992, 99(1): 61–112

[42] deSouza Neto E APerićDOwenD R J. Computational Methods for Plasticity: Theory and Applications. Hoboken: John Wiley & Sons, Ltd., 2008

[43] SimoJ CHughesT J R. Computational Inelasticity. Berlin: Springer, 1998

[44] OgdenR W. Non-linear Elastic Deformations. New York: Dover Publications, Inc., 1997

[45] BonetJWoodR D. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge: Cambridge University Press, 1997

[46] Simo J C. Numerical analysis and simulation of plasticity. Handbook of Numerical Analysis, 1998, 6: 183–499

[47] van den Heever M, Bester F, Kruger J, van Zijl G. Mechanical characterisation for numerical simulation of extrusionbased 3D concrete printing. Journal of Building Engineering, 2021, 44: 102944

[48] Wolfs R, Bos F, Salet T. Triaxial compression testing on early age concrete for numerical analysis of 3D concrete printing. Cement and Concrete Composites, 2019, 104: 103344

[49] Roussel N, Ovarlez G, Garrault S, Brumaud C. The origins of thixotropy of fresh cement pastes. Cement and Concrete Research, 2012, 42(1): 148–157

[50] Kruger J, Zeranka S, van Zijl G. 3D concrete printing: A lower bound analytical model for buildability performance quantification. Automation in Construction, 2019, 106: 102904

[51] Mohan M K, Rahul A V, van Tittelboom K, de Schutter G. Rheological and pumping behaviour of 3D printable cementitious materials with varying aggregate content. Cement and Concrete Research, 2021, 139: 106258

[52] Oñate E, Idelsohn S R, Del Pin F, Aubry R. The particle finite element method: An overview. International Journal of Computational Methods, 2004, 1(2): 267–307

[53] Idelsohn S R, Oñate E, Del Pin F. The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. International Journal for Numerical Methods in Engineering, 2004, 61(7): 964–989

[54] Carbonell J M, Rodriguez J M, Oñate E. Modelling 3D metal cutting problems with the particle finite element method. Computational Mechanics, 2020, 66(3): 603–624

[55] Wood W L, Bossak M, Zienkiewicz O C. An alpha modification of Newmark’s method. International Journal for Numerical Methods in Engineering, 1980, 15(10): 1562–1566

[56] Dohrmann C R, Bochev P. A stabilized finite element method for the stokes problem based on polynomial pressure projections. International Journal for Numerical Methods in Fluids, 2004, 46(2): 183–201

[57] Zeng W, Liu G R. Smoothed finite element methods (S-FEM): An overview and recent developments. Archives of Computational Methods in Engineering, 2018, 25(2): 397–435

[58] Oliver J, Hartmann S, Cante J, Weyler R, Hernández J. A contact domain method for large deformation frictional contact problems. Part 1: Theoretical basis. Computer Methods in Applied Mechanics and Engineering, 2009, 198(33−36): 2591–2606

[59] Hartmann S, Oliver J, Weyler R, Cante J, Hernández J. A contact domain method for large deformation frictional contact problems. Part 2: Numerical aspects. Computer Methods in Applied Mechanics and Engineering, 2009, 198(33−36): 2607–2631

[60] Carneau P, Mesnil R, Baverel O, Roussel N. Layer pressing in concrete extrusion-based 3D-printing: Experiments and analysis. Cement and Concrete Research, 2022, 155: 106741

[61] Suiker A S J. Mechanical performance of wall structures in 3D printing processes: Theory, design tools and experiments. International Journal of Mechanical Sciences, 2018, 137: 145–170

RIGHTS & PERMISSIONS

The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn