1 Introduction

In the past decades, additive manufacturing technologies such as 3D printing of concrete have become increasingly important for the industrial precast construction of supporting structures such as steel-reinforced beams or columns, stairs or benches but also in-situ construction of complete buildings or cyclist/pedestrian bridges [1,2,3]. In contrast to conventionally cast concrete, 3D printing of concrete has the potential of a fast and economical construction method for building structures with complex geometries without the use of complicated and cost-intensive formwork [4].

For guaranteeing a sufficient quality of the hardened 3D printed structure, 3D printed concrete (3DPC) needs to conform special requirements in both, the fresh and the hardened state [5, 6]. In the fresh state, rheological properties as well as a controllable hydration time [7,8,9] play an important role, e.g., a rapid hydration time prevents undesired material failure or buckling during the printing process [10,11,12,13]. In the hardened state, a sufficient material stiffness and material strength need to be ensured, which are significantly influenced by, e.g., the printing conditions, the width of the printed layers, the interfacial bond properties between the hardened filaments and the porosity of the hardened 3DPC.

Typical numerical modeling approaches focus on either the printing process or the hardened state of 3DPC.

The first category is concerned with obtaining a deeper understanding of the flow properties, the layer deformations and the occurring stresses during the printing process. To this end, numerical analyses of the printing process based on the Lagrangian-based particle finite element method [14] or Eulerian-based finite volume method combined with volume-of-fluid interface tracking [15] were performed. Since the focus of this paper is on the second category, i.e., the representation of the mechanical behavior of hardened 3D printed concrete, an overview on existing modeling approaches for concrete and 3DPC is provided in the following.

The mechanical behavior of hardened concrete—either conventionally cast or considered within the filaments of 3DPC—exhibits the following characteristics: (i) irreversible deformations associated with strain hardening in the pre-peak regime of the stress–strain curve and (ii) strain softening in the post-peak regime accompanied by the degradation of the material stiffness due to damaging processes. For modeling the highly nonlinear, isotropic mechanical behavior of concrete, several isotropic constitutive models were proposed in the literature. Among standard continuum-based single-phase models, plasticity models [16, 17] are suitable for describing irreversible strains in the pre-peak regime whereas damage models [18] are capable of capturing the degradation of material stiffness in the post-peak regime. For considering both, the combination of plasticity and damage models, the damage-plastic model for concrete proposed by Grassl and Jirásek [19] is a popular representative and provides the fundamental basis for the orthotropic material model for 3DPC to be presented in this publication. It is formulated by a stress-based plasticity formulation in the effective stress space, considering strain hardening by an evolving yield function in the pre-peak regime as well as strain softening accompanied by the degradation of material stiffness in the post-peak regime by a scalar damage variable. For ensuring a robust application of the model in finite element simulations with mesh-insensitive results in the softening regime, enriched formulations of the model have been proposed, such as an over-nonlocal implicit gradient-enhanced formulation in [20] or a formulation within the gradient-enhanced micropolar continuum in [21]. In addition, the model was extended in [22] for capturing the time-dependent behavior of shotcrete and even further extended for representing nonlinear creep as well as failure due to creep in [23]. Another category of models, particularly popular for considering the damage-induced anisotropic material behavior, are micro-plane models [24,25,26].

Lastly, a potentially attractive modeling category for describing the time-dependent material behavior of the hardening 3D printed concrete throughout the complete construction process could be offered by multi-phase models, e.g., by extending respective models for conventionally cast concrete [27,28,29,30,31]. They enable the time-dependent modeling of the hydration process incorporating the evolution of material strength and stiffness, temperature or creep and shrinkage by considering physical processes on the microstructure.

So far, the literature review highlighted advanced models for conventionally cast concrete. In the following, essential modeling components for 3DPC are reviewed: Due to the printing process of 3DPC, such as the extrusion-based 3D printing by applying fresh concrete layer-upon-layer [32] or 3D printing by binder jetting of cementitious powders [33], interfacial bond properties between the layered filaments significantly influence the mechanical behavior of 3DPC [34,35,36]. This results in an inherent orthotropic mechanical behavior as reported by many experimental investigations [10, 32, 33, 37, 38].

For considering the influence of interfacial bonds on the mechanical behavior of 3DPC in finite element simulations, typically two different categories of modeling approaches can be distinguished: In the first category, an isotropic constitutive model for concrete is employed for the filaments. The orthotropic behavior is considered as a structural effect by modeling the interfaces between the individual filaments by incorporating planes of weaknesses [35] or by using contact mechanics [8, 39]. They allow for detailed analyses of the interface regions at the expense of requiring a high computational cost, making them unsuitable for large-scale finite element simulations of 3DPC structures. In contrast, in the second category, the interfaces are smeared such that the layered material is substituted by a homogeneous anisotropic continuum, which is also the scope of the current publication. The latter models offer an approximate but computationally more efficient approach for modeling 3DPC in particular for larger 3DPC structures with multiple layers, as inherent interfaces and concrete are modeled in a smeared manner as a quasi-continuum. So far, only a few research groups made first steps to model the orthotropic, nonlinear material behavior of 3DPC via such models, e.g., proposing an uniaxial material law in total form [33] for limiting the principal stresses, neglecting irreversible deformations as well as damaging processes.

However, in the past, several approaches for continuum-based modeling of other inherently layered materials, such as timber or rocks, have been proposed and may be considered for modeling 3DPC. For example, in [40] a hybrid multi-phase-field model was proposed capable for capturing the orthotropic failure mechanisms of wood. For modeling the direction-dependent behavior of rock, typical approaches comprise the ubiquitous joint concept [41], the introduction of scalar anisotropy parameters [42, 43] or the projection into a fictitious isotropic configuration [44,45,46,47,48,49]. For a detailed comparison of these mentioned modeling approaches for rock, interested readers are referred to [50].

To the best knowledge of the authors, an advanced continuum-based model of the orthotropic nonlinear mechanical behavior of 3DPC is not available so far. Hence, an orthotropic damage plasticity model for 3DPC will be presented, denoted as 3DPCDP model. Thereby, the focus lies on the combined consideration of inherent orthotropic mechanical behavior due to the printing process and the highly nonlinear quasi-brittle material behavior. For this purpose, the isotropic concrete damage plasticity (CDP) model by [19] will be extended to orthotropic material behavior by mapping of the effective stress tensor into a fictitious isotropic configuration following the approach of [44, 45]. The latter was used for the extension of a rock damage plasticity model [51] to transversely isotropic material behavior in [50, 52] and applied to 2D finite element simulations of tunnel excavation in layered rock mass [53]. For obtaining mesh-insensitive results in the softening regime, the formulation of the CDP model is formulated within the framework of a gradient-enhanced continuum proposed in [20]. The capabilities of the model will be assessed by comparing the results of finite element simulations of uniaxial compression tests and three-point bending tests on extrusion-based 3DPC, performed for different loading directions with respect to the principal directions of the material, with the experimental results documented in [32]. In addition, the applicability of the model to replicate the mechanical behavior of 3DPC manufactured by an alternative 3D printing process by binder jetting of a cementitious powder [33] will be demonstrated by 3D finite element simulations of an arch structure with varying orientations of the loading direction relative to the layering.

The remainder of the paper is structured as follows: In Sect. 2, the constitutive equations of the 3DPCDP model for modeling the orthotropic nonlinear material behavior of 3DPC will be presented. In Sect. 3, the influence of newly introduced parameters for considering orthotropic plastic material behavior on the yield surface and the evolution of plastic strain will be explained. Section 4 comprises a numerical study on extrusion-based 3DPC considering the results of the experimental campaign of [32] with varying orientations of the loading direction relative to the layering. The material parameters are calibrated by considering data from uniaxial compression tests. Subsequently, the proposed model is validated by finite element simulations of three-point bending tests. In Sect. 5, based on a proof-of-concept study, the application of the 3DPCDP model to replicate the behavior of 3D printed concrete structures built by the alternative printing technique binder jetting of cementitious powder is demonstrated by 3D finite element simulations of an arch structure. Finally, in Sect. 6 the paper is closed with a summary and an outlook for recommended future research activities.

2 A damage plasticity model for 3D printed concrete

In this section, a new model for describing the orthotropic and highly nonlinear material behavior of 3DPC is proposed, denoted as 3D printed conrete damage plasticity (3DPCDP) model. It is formulated by the extension of the well-established, isotropic damage-plastic model for concrete (CDP) by Grassl and Jirásek [19] to orthotropic material behavior. In addition, it is enriched by a gradient-enhanced continuum formulation by Poh and Swaddiwudhipong [20] for obtaining mesh-insensitive results in the softening regime. Its mathematical framework is very similar to phase-field models [54] and compared to other existing regularization techniques [55] and fracture models [56, 57], it provides a framework that (i) is applicable for representing various types of failure modes, (ii) can be calibrated from lab tests and (iii) is computationally efficient to realize large-scale 3D structural simulations. To this end, the original transversely isotropic mapping tensor in [44, 45] is used as a basis for formulating the orthotropic mapping tensor. The isotropic CDP model is characterized by coupling the theories of hardening plasticity and damage mechanics. Thus, it is perfectly suited to represent the highly nonlinear material behavior of 3DPC.

2.1 Stress–strain relation

The stress–strain relation of the orthotropic 3DPCDP model, expressed in the global coordinate system, reads as

$$\begin{aligned} \varvec{\sigma }= \left( 1-\omega \right) \cdot \mathbb {C}^{(\text {o})}:\underbrace{\left( \varvec{\varepsilon }- \varvec{\varepsilon }^{(p)} \right) }_{\varvec{\varepsilon }^{(e)}}. \end{aligned}$$
(1)

Here, \(\varvec{\sigma }\) represents the nominal stress tensor and \(\varvec{\varepsilon }\) the total strain tensor, which is decomposed into its plastic part \(\varvec{\varepsilon }^{(p)}\) and elastic part \(\varvec{\varepsilon }^{(e)}\) by an additive split. \(\omega \) denotes a scalar damage variable ranging between 0 (undamaged material) and 1 (completely damaged material). According to continuum damage mechanics, the nominal stress tensor \(\varvec{\sigma }\) (force per total area) is linked to the effective stress tensor \(\bar{\varvec{\sigma }}\) (force per undamaged area) by \(\varvec{\sigma }= (1 - \omega )\,\bar{\varvec{\sigma }}\).

The orthogonal transformation matrix \(a_{ij} = \cos \angle \left( {\varvec{e}}^{\prime \,(i)},{\varvec{e}}^{(j)}\right) \) consists of the direction cosines between the base vectors of the local and the global coordinate system \({\varvec{e}}^{\prime \,(i)}\) and \({\varvec{e}}^{(i)},i\in \{1,2,3\}\), respectively. The notation of a superscript prime \(\left( \bullet \right) ^\prime \) is used to refer to variables formulated in the local coordinate system aligned with the principal material directions (cf. Fig. 1).

Fig. 1
figure 1

Sketch of a 3D printed orthotropic material with \({\varvec{e}}^{\prime (i)}, i \in \{1,2,3\}\), denoting the orthogonal base vectors of the local coordinate system aligned with the principal material directions, i.e., the three axis of material symmetry, and \({\varvec{e}}^{(i)}\) the orthogonal base vectors of the global coordinate system

Full size image

By means of the transformation matrix, the fourth-order orthotropic elastic stiffness tensor, formulated with respect to the local coordinate system, reads in Voigt notation as

$$\begin{aligned} \mathbb {C}^{\prime (\text {o})}= {\begin{bmatrix} \frac{1}{E_1} &{} -\frac{\nu _{12}}{E_2} &{} -\frac{\nu _{13}}{E_3} &{} 0 &{} 0 &{} 0\\ -\frac{\nu _{12}}{E_2} &{} \frac{1}{E_2} &{} -\frac{\nu _{23}}{E_3} &{} 0 &{} 0 &{} 0\\ -\frac{\nu _{13}}{E_3} &{} -\frac{\nu _{23}}{E_3} &{} \frac{1}{E_3} &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} \frac{1}{G_{12}} &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} \frac{1}{G_{13}} &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \frac{1}{G_{23}}\end{bmatrix}}^{-1}. \end{aligned}$$
(2)

It is transformed into the global coordinate system by

$$\begin{aligned} \mathbb {C}^{(\text {o})}_{ijkl}=a_{mi}\, a_{nj}\, a_{pk}\, a_{ql}\, \mathbb {C}^{\prime (\text {o})}_{mnpq}. \end{aligned}$$
(3)

The orthotropic elastic stiffness tensor is characterized by nine independent elastic material parameters, i.e., the three Young’s moduli \(E_i\in \{E_1, E_2, E_3\}\) measured along the corresponding principal material directions \({\varvec{e}}^{\prime \,(i)}\), the three Poisson’s ratios \(\nu _{ij}\in \{\nu _{12}, \nu _{13}, \nu _{23}\}\), corresponding to a contraction in direction \({\varvec{e}}^{\prime \,(j)}\) due to an extension in direction \({\varvec{e}}^{\prime \,(i)}\), and the three shear moduli \(G_{ij}\in \{G_{12}, G_{13}, G_{23}\}\) in the planes containing the principal directions \({\varvec{e}}^{\prime \,(i)}\) and \({\varvec{e}}^{\prime \,(j)}, i\ne j\). Due to the symmetry of \(\mathbb {C}^{\prime (\text {o})}\), the Poisson’s ratios \(\nu _{21}, \nu _{31}\) and \(\nu _{32}\) follow from the relations

$$\begin{aligned} \frac{\nu _{21}}{E_1} = \frac{\nu _{12}}{E_2}, \hspace{.5cm} \frac{\nu _{31}}{E_1} = \frac{\nu _{13}}{E_3}, \hspace{.5cm} \frac{\nu _{23}}{E_2} = \frac{\nu _{32}}{E_3}. \end{aligned}$$
(4)

2.2 Orthotropic hardening plasticity

Orthotropic plastic material behavior is incorporated based on the notion of a fictitious isotropic configuration, denoted by the superscript \((\bullet )^*\). The effective stress tensor \(\bar{\varvec{\sigma }}\) is projected to a mapped effective stress tensor

$$\begin{aligned} \bar{\varvec{\sigma }}^*= \mathbb {O}:\bar{\varvec{\sigma }}\end{aligned}$$
(5)

in the fictitious isotropic configuration based on a fourth-order mapping tensor \(\mathbb {O}\). Inspired by the mapping tensor for transversely isotropic behavior presented in [44, 45, 50] and applied in [53] for modeling the transversely isotropic material behavior of layered rock, a novel mapping tensor \(\mathbb {O}\) for rendering orthotropic material behavior is proposed:

$$\begin{aligned} \mathbb {O} = \alpha \,\mathbb {E}^{(11)} + \beta \, \mathbb {E}^{(22)} + \gamma \, \mathbb {E}^{(33)} + \zeta \, \mathbb {S}^{(12)} + \xi \, \mathbb {S}^{(13)} + \eta \,\mathbb {S}^{(23)}. \end{aligned}$$
(6)

The orthotropic mapping tensor \(\mathbb {O}\) represents a fourth-order linear tensor with major and minor symmetry and is constituted from the six sub-tensors:

$$\begin{aligned} \begin{aligned} \mathbb {E}^{(11)}_{ijkl}&= \varphi _{ij}\,\varphi _{kl},\hspace{.5cm} \mathbb {E}^{(22)}_{ijkl} = \chi _{ij}\,\chi _{kl},\hspace{.5cm} \mathbb {E}^{(33)}_{ijkl} = \psi _{ij}\,\psi _{kl}, \\ \mathbb {S}^{(12)}_{ijkl}&= \frac{1}{2}\,\left( \varphi _{ik}\,\chi _{jl} + \varphi _{il}\,\chi _{jk} + \varphi _{jl}\,\chi _{ik} + \varphi _{jk}\,\chi _{il}\right) , \\ \mathbb {S}^{(13)}_{ijkl}&= \frac{1}{2}\,\left( \varphi _{ik}\,\psi _{jl} + \varphi _{il}\,\psi _{jk} + \varphi _{jl}\,\psi _{ik} + \varphi _{jk}\,\psi _{il}\right) , \\ \mathbb {S}^{(23)}_{ijkl}&= \frac{1}{2}\,\left( \chi _{ik}\,\psi _{jl} + \chi _{il}\,\psi _{jk} + \chi _{jl}\,\psi _{ik} + \chi _{jk}\,\psi _{il}\right) . \end{aligned} \end{aligned}$$
(7)

The latter are formulated in terms of the three second-order microstructure tensors

$$\begin{aligned} \begin{aligned} \varphi _{ij} = e^{\prime (1)}_i\,e^{\prime (1)}_j,\hspace{.5cm} \chi _{ij} = e^{\prime (2)}_i\,e^{\prime (2)}_j,\hspace{.5cm} \psi _{ij} = e^{\prime (3)}_i\,e^{\prime (3)}_j, \end{aligned} \end{aligned}$$
(8)

defined by the three base vectors of the principal directions of the material \({\varvec{e}}^{\prime \,(i)}, i \in \{1, 2, 3\}\) (cf. Fig. 1). If an effective stress tensor is transformed into the local coordinate system, the superscripts \((\bullet )^{(ij)}\) of the tensors \(\mathbb {E}\) and \(\mathbb {S}\) indicate that they only influence the mapping of the ijth (and correspondingly the jith) component of the effective stress tensor into the isotropic configuration. This can be shown by writing the mapping of an effective stress tensor aligned with the principal material directions into its fictitious isotropic configuration in Voigt notation

$$\begin{aligned} \begin{bmatrix} {\bar{\sigma }}^{\prime \,*}_{11}\\ {\bar{\sigma }}^{\prime \,*}_{22}\\ {\bar{\sigma }}^{\prime \,*}_{33}\\ {\bar{\sigma }}^{\prime \,*}_{12}\\ {\bar{\sigma }}^{\prime \,*}_{13}\\ {\bar{\sigma }}^{\prime \,*}_{23} \end{bmatrix} = \begin{bmatrix} \alpha &{}0&{}0&{}0&{}0&{}0\\ 0&{}\beta &{}0&{}0&{}0&{}0\\ 0&{}0&{}\gamma &{}0&{}0&{}0\\ 0&{}0&{}0&{}\zeta &{}0&{}0\\ 0&{}0&{}0&{}0&{}\xi &{}0\\ 0&{}0&{}0&{}0&{}0&{}\eta \end{bmatrix}\,\begin{bmatrix} {\bar{\sigma }}^{\prime }_{11}\\ {\bar{\sigma }}^{\prime }_{22}\\ {\bar{\sigma }}^{\prime }_{33}\\ {\bar{\sigma }}^{\prime }_{12}\\ {\bar{\sigma }}^{\prime }_{13}\\ {\bar{\sigma }}^{\prime }_{23} \end{bmatrix}, \end{aligned}$$
(9)

for which the mapping tensor \(\mathbb {O}\) boils down to a diagonal matrix with the six orthotropy parameters \(\alpha , \beta , \gamma , \zeta , \xi , \eta \) representing the diagonal entries. Thus, the latter represent proportionality factors between the respective effective stress component and the mapped effective stress component. Following a similar calibration approach as proposed in [50] by considering limit states of uniaxial compression/simple shearing for different loading directions, the calibration of the orthotropy parameters is proposed as:

$$\begin{aligned} \alpha = \frac{f_\text {cu}^{*}}{f_\text {cu}^{(1)}}, \hspace{.5cm} \beta = \frac{f_\text {cu}^{*}}{f_\text {cu}^{(2)}}, \hspace{.5cm} \gamma = \frac{f_\text {cu}^{*}}{f_\text {cu}^{(3)}} \end{aligned}$$
(10)

and

$$\begin{aligned} \zeta = \frac{f_\text {sh}^{*}}{f_\text {sh}^{(12)}}, \hspace{.5cm} \xi = \frac{f_\text {sh}^{*}}{f_\text {sh}^{(13)}}, \hspace{.5cm} \eta = \frac{f_\text {sh}^{*}}{f_\text {sh}^{(23)}}. \end{aligned}$$
(11)

Therein, the parameters \(\alpha , \beta \) and \(\gamma \) represent the ratios between the equivalent isotropic uniaxial compressive strength \(f_\text {cu}^*\), i.e., in the case of 3DPC corresponding to the uniaxial compressive strength of the cast concrete used for 3D printing, and the values for the uniaxial compressive strength \(f_\text {cu}^{(i)}\in \{f_\text {cu}^{(1)},f_\text {cu}^{(2)},f_\text {cu}^{(3)}\}\) of the hardened 3D printed concrete, loaded in the three principal directions of the material. Likewise, the parameters \(\zeta , \xi \) and \(\eta \) describe the ratios between the equivalent isotropic shear strength \(f_\text {sh}^*\) and the shear strength \(f_\text {sh}^{(ij)}\in \{f_\text {sh}^{(12)},f_\text {sh}^{(13)},f_\text {sh}^{(23)}\}\) of the 3D printed concrete loaded in the \({\varvec{e}}^{\prime \,(i)} {\varvec{e}}^{\prime \,(j)}\)-planes.

For delimiting the elastic domain by an orthotropic yield function, the Haigh-Westergaard coordinates of the mapped effective stress tensor \(\bar{\varvec{\sigma }}^*\), i.e., the mean stress \({\bar{\sigma }}_\text {m}^*\), the deviatoric radius \({\bar{\rho }}^*\) and the lode angle \({\bar{\theta }}^*\)

$$\begin{aligned} {\bar{\sigma }}_\text {m}^*= {\bar{I}}^{*}_1/3,\hspace{.5cm} {\bar{\rho }}^*= \sqrt{2\,{\bar{J}}^{*}_2} \hspace{.5cm} \text {and} \hspace{.5cm} {\bar{\theta }}^*= \frac{1}{3}\arccos \left( \frac{3\sqrt{3}}{2}\frac{{\bar{J}}^{*}_3}{\sqrt{\left( {\bar{J}}^{*}_2\right) ^3}}\right) , \end{aligned}$$
(12)

are incorporated into the original isotropic yield function of the CDP model [19] by replacing the respective quantities. \({\bar{I}}^{*}_1\) denotes the first invariant of \(\bar{\varvec{\sigma }}^*, and {\bar{J}}^{*}_2\) and \({\bar{J}}^{*}_3\) the second and third invariants of the deviatoric part of \(\bar{\varvec{\sigma }}^*\). Thus, the orthotropic yield function in the fictitious isotropic configuration reads as

$$\begin{aligned} \begin{aligned} f\left( {\bar{\sigma }}_\text {m}^*,{\bar{\rho }}^*,{\bar{\theta }}^*,q_\text {h}\left( \alpha _\text {p}\right) \right)&= \left( \left( 1-q_\text {h}\left( \alpha _\text {p}\right) \right) \left( \frac{{\bar{\rho }}^*}{\sqrt{6}\,f_\text {cu}^{*}} + \frac{{\bar{\sigma }}_\text {m}^*}{f_\text {cu}^{*}}\right) ^2 + \sqrt{\frac{3}{2}}\frac{{\bar{\rho }}^*}{f_\text {cu}^{*}} \right) ^2 \\&\quad + m_0^*\,q_\text {h}^2\left( \alpha _\text {p}\right) \left( \frac{{\bar{\rho }}^*}{\sqrt{6}\,f_\text {cu}^{*}}\,r\left( {\bar{\theta }}^*, e\right) + \frac{{\bar{\sigma }}_\text {m}^*}{f_\text {cu}^{*}}\right) - q_\text {h}^2\left( \alpha _\text {p}\right) , \end{aligned} \end{aligned}$$
(13)

with \(f_\text {cu}^{*}\) denoting the equivalent isotropic uniaxial compressive strength and \(m_0^*\) the equivalent isotropic friction parameter, to be measured on cast samples of the concrete used for 3D printing. The yield function evolves depending on the stress-like hardening variable \(\dfrac{f_\text {cy}^{*}}{f_\text {cu}^{*}}\le q_\text {h}\le 1\). In Appendix A, the orthotropic yield function in the fictitious isotropic configuration is specified for simple loading scenarios, i.e., uniaxial compression with loading in the different principal directions of the material and simple shearing in planes perpendicular to the planes of material symmetry. The shape of the yield function in deviatoric planes is controlled by the Willam–Warnke function [16]

$$\begin{aligned} r\left( {\bar{\theta }}^*, e\right) = \frac{4\,\left( 1 - e^2\right) \,\cos ^2\left( {\bar{\theta }}^*\right) + \left( 2\,e - 1\right) ^2}{2\,\left( 1-e^2\right) \,\cos \left( {\bar{\theta }}^*\right) + \left( 2\,e - 1\right) \,\sqrt{4\,\left( 1 - e^2\right) \,\cos ^2\left( {\bar{\theta }}^*\right) + 5e^2-4e}}, \end{aligned}$$
(14)

with the eccentricity parameter

$$\begin{aligned} e = \frac{1 + \epsilon }{2 - \epsilon }, \end{aligned}$$
(15)

restricted to values between 0.5 (triangular shape) and 1 (circular shape) for ensuring a convex shape of the yield function. The auxiliary parameter

$$\begin{aligned} \epsilon = \frac{f_\text {tu}^*\,\left( f_\text {cb}^{*\,2} - f_\text {cu}^{*\,2}\right) }{f_\text {cb}^*\,\left( f_\text {cu}^{*\,2} - f_\text {tu}^{*\,2}\right) } \end{aligned}$$
(16)

is calculated by the equivalent isotropic uniaxial tensile strength \(f_\text {tu}^*\), the equivalent isotropic uniaxial compressive strength \(f_\text {cu}^*\) and the equivalent isotropic equi-biaxial compressive strength \(f_\text {cb}^*\), all referring to parameters determined from lab tests on cast concrete specimens used for 3D printing. Based on the eccentricity parameter (15), the equivalent isotropic friction parameter \(m_0^*\) is determined as

$$\begin{aligned} m_0^* = 3\,\frac{f_\text {cb}^{*\,2} - f_\text {tu}^{*\,2}}{f_\text {cu}^*\,f_\text {tu}^*}\frac{e}{e+1}. \end{aligned}$$
(17)

The evolution of the plastic strain \(\varvec{\varepsilon }^{(p)}\) is described by the non-associated flow rule

$$\begin{aligned} \dot{\varvec{\varepsilon }}^{(p)} = {\dot{\lambda }} \frac{\partial \,g\left( {\bar{\sigma }}_\text {m}^*,{\bar{\rho }}^*,{\bar{\theta }}^*,q_\text {h}\left( \alpha _\text {p}\right) \right) }{\partial \,\bar{\varvec{\sigma }}^*}:\frac{\partial \,\bar{\varvec{\sigma }}^*}{\partial \,\bar{\varvec{\sigma }}}= {\dot{\lambda }}\,\mathbb {O}:\frac{\partial \,g\left( {\bar{\sigma }}_\text {m}^*,{\bar{\rho }}^*,{\bar{\theta }}^*,q_\text {h}\left( \alpha _\text {p}\right) \right) }{\partial \,\bar{\varvec{\sigma }}^*}, \end{aligned}$$
(18)

with the plastic potential function

$$\begin{aligned} \begin{aligned} g\left( {\bar{\sigma }}_\text {m}^*,{\bar{\rho }}^*,{\bar{\theta }}^*,q_\text {h}\left( \alpha _\text {p}\right) \right)&= \left( \left( 1-q_\text {h}\left( \alpha _\text {p}\right) \right) \left( \frac{{\bar{\rho }}^*}{\sqrt{6}\,f_\text {cu}^{*}} + \frac{{\bar{\sigma }}_\text {m}^*}{f_\text {cu}^{*}}\right) ^2 + \sqrt{\frac{3}{2}}\frac{{\bar{\rho }}^*}{f_\text {cu}^{*}} \right) ^2 \\&\quad + q_\text {h}^2\left( \alpha _\text {p}\right) \left( \frac{m_0^*\,{\bar{\rho }}^*}{\sqrt{6}\,f_\text {cu}^{*}} + \frac{m_\text {g}^*({\bar{\sigma }}_\text {m}^*)}{f_\text {cu}^{*}}\right) . \end{aligned} \end{aligned}$$
(19)

By analogy to the yield function, the latter is formulated in terms of the mapped effective stress tensor \(\bar{\varvec{\sigma }}^*\) for considering the influence of orthotropic material behavior on \(\varvec{\varepsilon }^{(p)}\). The evolution of the volumetric plastic strain is controlled by the pressure-dependent dilatancy function

$$\begin{aligned} m_\text {g}^*({\bar{\sigma }}_\text {m}^*) = A_\text {g}\,B_\text {g}\,f_\text {cu}^{*}\,\exp \left( \frac{{\bar{\sigma }}_\text {m}^*- f_\text {tu}^{*}/3}{B_\text {g}\,f_\text {cu}^{*}}\right) , \end{aligned}$$
(20)

with

$$\begin{aligned} A_\text {g}= 3\,\frac{f_\text {tu}^{*}}{f_\text {cu}^{*}} + \frac{m_0^*}{2} \end{aligned}$$
(21)

and

$$\begin{aligned} B_\text {g}= \frac{1}{3}\,\frac{1 + f_\text {tu}^*/f_\text {cu}^*}{\ln \left( A_\text {g}\right) - \ln \left( 2D_\text {f}- 1\right) - \ln \left( 3 + m_0^*/2\right) + \ln \left( D_\text {f}+ 1\right) }, \end{aligned}$$
(22)

depending on the mapped effective mean stress \({\bar{\sigma }}_\text {m}^*\). \(D_\text {f}\) denotes an independent model parameter with a suggested default value of 0.85 according to [19].

Hardening plastic material behavior is controlled by the stress-like hardening variable

$$\begin{aligned} q_\text {h}(\alpha _\text {p}) = {\left\{ \begin{array}{ll} \dfrac{f_\text {cy}^{*}}{f_\text {cu}^{*}} + \left( 1-\dfrac{f_\text {cy}^{*}}{f_\text {cu}^{*}}\right) \,\alpha _p\,\left( \alpha ^2_p - 3\,\alpha _p + 3 \right) &{} \alpha _p < 1\\ 1 &{} \alpha _p \ge 1 \end{array}\right. }, \end{aligned}$$
(23)

ranging between \(\dfrac{f_\text {cy}^{*}}{f_\text {cu}^{*}}\) (initial state) and 1 (fully hardened state), with \(f_\text {cy}^{*}\) denoting the equivalent isotropic uniaxial compressive yield stress determined from the cast concrete used for 3D printing. The stress-like hardening variable \(q_\text {h}\) is driven by the strain-like internal variable \(\alpha _\text {p}\), according to the smooth hardening function

$$\begin{aligned} {\dot{\alpha }}_\text {p}= \frac{\left| \left| \dot{\varvec{\varepsilon }}^{(p)}\right| \right| }{x_\text {h}\left( {\bar{\sigma }}_\text {m}^*\right) }\,\left( 1 + 3\,\frac{{\bar{\rho }}^{*2}}{{\bar{\rho }}^{*2}+ f_\text {cu}^{*\,2} \cdot 10^{-8}}\,\cos ^2\left( \frac{3\,{\bar{\theta }}^*}{2}\right) \right) , \end{aligned}$$
(24)

proposed in [58]. The hardening ductility measure

$$\begin{aligned} x_\text {h}\left( {\bar{\sigma }}_\text {m}^*\right) = {\left\{ \begin{array}{ll} \left( B_\text {h}- D_\text {h}\right) \,\exp \left( \dfrac{R_\text {h}\left( {\bar{\sigma }}_\text {m}^*\right) \left( A_\text {h}- B_\text {h}\right) }{\left( B_\text {h}- D_\text {h}\right) \,C_\text {h}}\right) &{} R_\text {h}\left( {\bar{\sigma }}_\text {m}^*\right) <0\\ A_\text {h}- \left( A_\text {h}- B_\text {h}\right) \,\exp \left( \dfrac{-R_\text {h}\left( {\bar{\sigma }}_\text {m}^*\right) }{C_\text {h}}\right) &{} R_\text {h}\left( {\bar{\sigma }}_\text {m}^*\right) \ge 0 \end{array}\right. }, \end{aligned}$$
(25)

controls the ductility behavior in the hardening regime, characterized by four independent material parameters \(A_\text {h}, B_\text {h}\), \(C_\text {h}\) and \(D_\text {h}\) and the orthotropic transition point \(R_\text {h}\left( {\bar{\sigma }}_\text {m}^*\right) = -{\bar{\sigma }}_\text {m}^*/f_\text {cu}^{*} - 1/3\), the latter ensuring a more ductile pre-peak response in compression than in tension. According to [19], \(A_\text {h}= 0.08, B_\text {h}= 0.003, C_\text {h}= 2\) and \(D_\text {h}= 10^{-6}\) are suggested as default values for the hardening parameters.

Fig. 2
figure 2

Evolving yield surface of the 3DPCDP model in the mapped principal effective stress space (A) and the principal effective stress space (B) for \(q_\text {h}= 0\) (orange), \(q_\text {h}= 1/2\) (gray) and \(q_\text {h}= 1\) (blue) for principal stress directions aligned with the principal directions of the material; the orthotropy parameters are set to \(\alpha = 1.0, \beta = 1.5, \gamma = 0.5\), \(\zeta = \xi = \eta = 1.0\) (colour figure online)

Full size image

Figure 2 shows the evolution of the yield surface (13) of the 3DPCDP model during strain hardening, illustrated in (A) in terms of the mapped principal effective stress space \({\bar{\sigma }}^*_\text {I}\) and (B) the principal effective stress space \({\bar{\sigma }}_\text {I}\) with \(\text {I} \in \{1,2,3\}\), evaluated for principal stress directions aligned with the principal directions of the material (cf. Fig. 1). From Fig. 2A follows, that the elastic domain within the mapped principal effective stress space, i.e., corresponding to the one of the isotropic reference configuration, is bounded by the isotropic yield surface. The latter is characterized by equal values for the uniaxial compressive/tensile strength along the \({\bar{\sigma }}^*_1\)-, \({\bar{\sigma }}^*_2\)- and \({\bar{\sigma }}^*_3\)-axes, independent from the loading direction. In contrast, Fig. 2 B shows that the yield surface of the 3DPCDP model plotted with respect to the principal effective stress space is clearly orthotropic, visible by the dependence of, e.g., the values for the uniaxial compressive/tensile strength along the \({\bar{\sigma }}_1\)-, \({\bar{\sigma }}_2\)- and \({\bar{\sigma }}_3\)-axes on the loading direction. Further, it can be seen that due to the orthotropic material behavior, the apex of the orthotropic yield function is not located on the hydrostatic axis anymore.

2.3 Damage formulation for inherent orthotropic material behavior

Once the yield function (13) attains its fully hardened state, i.e., at \(q_\text {h} = \alpha _\text {p}= 1\), damage is initiated. Physically, this transition is manifested by the coalescence of microscopic cracks to a marcoscopic crack, resulting in a degradation of the material stiffness and strength. The softening behavior is described by a scalar damage variable

$$\begin{aligned} \begin{aligned} 0\ \le \ \omega ({\hat{\alpha }}_\text {d}) = 1 - \exp \left( \frac{-{\hat{\alpha }}_\text {d}\left( {\bar{\alpha }}_d,\alpha _\text {d}\right) }{\varepsilon ^{*}_{\text {f}}}\right) \ \le \ 1, \end{aligned} \end{aligned}$$
(26)

which is controlled by the equivalent isotropic softening modulus \(\varepsilon ^{*}_{\text {f}}\), referring to the value of the cast concrete used for 3D printing and calibrated from the specific mode I fracture energy \(G^{*}_\text {fI}\) of the latter. In order to remedy the pathological mesh-sensitivity of finite element simulations in the softening regime, an over-nonlocal implicit gradient-enhanced formulation [20] is used. To this end, the damage law (26) is formulated in terms of the over-nonlocal softening variable

$$\begin{aligned} {\hat{\alpha }}_\text {d} = m\,{\bar{\alpha }}_\text {d} + (1 - m)\,\alpha _\text {d}, \end{aligned}$$
(27)

with \({\bar{\alpha }}_d\) denoting the nonlocal softening variable, \(\alpha _\text {d}\) the local softening variable and m a weighting parameter for controlling the degree of nonlocality typically chosen to be slightly larger than one, e.g., 1.05. The nonlocal softening variable \({\bar{\alpha }}_\text {d}\) is implicitly determined by the screened Poisson equation

$$\begin{aligned} \alpha _\text {d} = {\bar{\alpha }}_\text {d} - l^2\,\nabla ^2\,{\bar{\alpha }}_\text {d} \end{aligned}$$
(28)

from the local softening variable \(\alpha _\text {d}\) by smearing the latter over a finite but small domain determined by the internal length parameter l. The weak form of (28) combined with the one of the equilibrium equations (cf. Appendix B) represent, after discretization in space and time, a coupled system of algebraic equations to be solved for the displacement field and the nonlocal softening variable. The evolution of the local softening variable \(\alpha _\text {d}\) is determined by the rate equation

$$\begin{aligned} \begin{aligned} {\dot{\alpha }}_{\text {d}} = {\left\{ \begin{array}{ll} 0 &{} \alpha _{\text {p}} < 1\\ \dfrac{{\dot{\varepsilon }}^{(p)}_{\text {vol}}}{x_\text {s}\left( \dot{\varvec{\varepsilon }}^{(p)}\right) } &{} \alpha _{\text {p}} \ge 1 \end{array}\right. }. \end{aligned} \end{aligned}$$
(29)

It depends on the volumetric plastic strain rate \({\dot{\varepsilon }}^{(p)}_{\text {vol}} = \text {tr}\left( \dot{\varvec{\varepsilon }}^{(p)}\right) \) and the softening ductility measure

$$\begin{aligned} x_\text {s}\left( \dot{\varvec{\varepsilon }}^{(p)}\right) = {\left\{ \begin{array}{ll} 1 + A_\text {s}\,R_\text {s}\left( \dot{\varvec{\varepsilon }}^{(p)}\right) ^2 &{} R_\text {s}\left( \dot{\varvec{\varepsilon }}^{(p)}\right) < 1\\ 1 - 3A_\text {s}+ 4A_\text {s}\sqrt{R_\text {s}\left( \dot{\varvec{\varepsilon }}^{(p)}\right) } &{} R_\text {s}\left( \dot{\varvec{\varepsilon }}^{(p)}\right) \ge 1 \end{array}\right. }. \end{aligned}$$
(30)

Therein, \(A_\text {s}\) denotes a softening ductility parameter in the post-peak regime to be calibrated from uniaxial compression tests on the isotropic cast concrete used for 3D printing or alternatively assumed as the recommended default value of \(A_\text {s}= 15\) according to [19]. Based on the normalized compressive part of the volumetric plastic strain rate, represented by

$$\begin{aligned} R_\text {s}\left( \dot{\varvec{\varepsilon }}^{(p)}\right) = \frac{1}{{\dot{\varepsilon }}^{(p)}_{\text {vol}}}\,\sum ^3_{\text {I} = 1}\left\langle - {\dot{\varepsilon }}^{(p)}_\text {I}\right\rangle \end{aligned}$$
(31)

with \({\dot{\varepsilon }}^{(p)}_\text {I}\) denoting the principal plastic strain rates and \(\left\langle \bullet \right\rangle \) the Macaulay brackets, the influence of multiaxial stress states on the evolution of the softening variable is taken into account. The influence of orthotropic material behavior on the damage behavior is inherently considered by the dependence of \(\alpha _\text {d}\) (29) on the rate of the plastic strain tensor (18).

The presented constitutive equations of the 3DPCDP model are integrated numerically by employing an implicit (backward) Euler algorithm and were implemented into the framework of the Marmot material modeling toolbox library [59], providing a C++ library within the in-house finite element codes EdelweissFE [60] and mpFEM [61] but also into commercial finite element programs such as Abaqus [62], Plaxis [63] and OpenSees [64].

3 Influence of the orthotropy parameters on the yield function and the plastic potential function

This section provides an illustrative explanation of the orthotropy parameters, introduced in (10) and (11), by demonstrating their influence on the yield function and the plastic potential function.

For evaluating the influence of the orthotropy parameters \(\alpha \), \(\beta \) and \(\gamma \), deviatoric sections through the ultimate strength surface at \(q_\text {h}= 1\) are illustrated in the principal effective stress space in Fig. 3. The principal stress directions are aligned with the principal directions of the material. The orthotropy parameters are set to representative values of \(\alpha = 1.0 , \beta = 1.5, \gamma = 0.5\). Since the orthotropy parameters \(\zeta , \xi \) and \(\eta \) do not influence the deviatoric sections of the yield surface for principal stress directions aligned with the principal material directions, they are set to \(\zeta = \xi = \eta = 1.0\). Figure 3A shows the section of the yield surface in the mapped principal effective stress space, i.e., in the isotropic reference configuration for \({\bar{\sigma }}_\text {m}^*= f_\text {cu}^*/3\). The isotropic yield surface is characterized by mirror-symmetry with respect to the principal stress axes. In Fig. 3B, the section of the orthotropic yield surface in the principal effective stress space is shown for \({\bar{\sigma }}_\text {m}= f_\text {cu}^{(1)}/3\) (gray), \({\bar{\sigma }}_\text {m}= f_\text {cu}^{(2)}/3\) (orange) and \({\bar{\sigma }}_\text {m}= f_\text {cu}^{(3)}/3\) (blue), representing the three uniaxial limit stress states \({\bar{\sigma }}_1=f_\text {cu}^{(1)}, {\bar{\sigma }}_2=f_\text {cu}^{(2)}\) and \({\bar{\sigma }}_3=f_\text {cu}^{(3)}\). For uniaxial compression in \({\varvec{e}}^{\prime \,(1)}\)-direction, as indicated by the gray arrow stress path, the mapped principal effective stress component and the corresponding principal effective stress component are proportional. Hence, \({\bar{\sigma }}^*_1 = \alpha \,{\bar{\sigma }}_1\). If the ultimate stress \(f_\text {cu}^{(1)}\) is attained in the principal effective stress space, its relation to the equivalent isotropic uniaxial compressive strength is obtained from \(f_\text {cu}^{(1)} = f_\text {cu}^*/\alpha \), as introduced by (10)\(_1\). Similarly, it can be seen that for uniaxial loading in \({\varvec{e}}^{\prime \,(2)}\)-direction (indicated by the orange arrow) and in \({\varvec{e}}^{\prime \,(3)}\)-direction (indicated by the blue arrow), the orthotropy parameters \(\beta \) and \(\gamma \) represent the proportionality coefficients between the mapped effective stress component and the effective stress component, i.e., \({\bar{\sigma }}^*_2 = \beta \,{\bar{\sigma }}_2\) and \({\bar{\sigma }}^*_3 = \gamma \,{\bar{\sigma }}_3\).

Fig. 3
figure 3

Deviatoric sections through A the equivalent isotropic yield surface and B the orthotropic yield surface for principal stress directions aligned with the principal directions of the material (orthotropy parameters chosen as \(\alpha = 1.0, \beta = 1.5, \gamma = 0.5, \zeta = \xi = \eta = 1.0\))

Full size image

In Fig. 4, the influence of \(\alpha \), \(\beta \) and \(\gamma \) on the yield surface is shown for plane stress conditions again for stress paths aligned with the principal directions of the material for the special case \(\alpha = 1.0, \beta = 1.5\) and \(\gamma = 0.5\). For this purpose, sections through the orthotropic (orange curve) and the equivalent isotropic yield surface (blue curve) are illustrated in the \({\bar{\sigma }}_1{\bar{\sigma }}_2\)-plane at \({\bar{\sigma }}_3 = 0\) (Fig. 4 (left)), in the \({\bar{\sigma }}_1{\bar{\sigma }}_3\)-plane at \({\bar{\sigma }}_2 = 0\) (Fig. 4 (middle)) and in the \({\bar{\sigma }}_2{\bar{\sigma }}_3\)-plane at \({\bar{\sigma }}_1 = 0\) (Fig. 4 (right)).

Fig. 4
figure 4

Plane stress sections through the orthotropic yield surface (orange line) and the equivalent isotropic yield surface (blue line) for principal stress directions aligned with the principal directions of the material: \({\bar{\sigma }}_1{\bar{\sigma }}_2\)-plane (left), \({\bar{\sigma }}_1{\bar{\sigma }}_3\)-plane (middle) and \({\bar{\sigma }}_2{\bar{\sigma }}_3\)-plane (right); the orthotropy parameters are set to \(\alpha = 1.0, \beta = 1.5, \gamma = 0.5\), \(\zeta = \xi = \eta = 1.0\) (colour figure online)

Full size image

Again, for uniaxial loading along the principal stress axes, the proportionality between the mapped principal effective stress and the corresponding principal effective stress component through the orthotropy parameters \(\alpha , \beta \) and \(\gamma \) becomes apparent. Furthermore, considering the different planes of plane stress conditions, it is seen that \(\alpha = f_\text {cu}^{*}/f_\text {cu}^{(1)}\) according to (10)\(_1\) only affects the section of the yield surface in the \({\bar{\sigma }}_1{\bar{\sigma }}_2\)-plane (cf. Fig. 4 (left)) and \({\bar{\sigma }}_1{\bar{\sigma }}_3\)-plane (cf. Fig. 4 (middle)). By analogy, \(\beta \) and \(\gamma \) have only an influence on the yield surface in the \({\bar{\sigma }}_1{\bar{\sigma }}_2\)-plane and \({\bar{\sigma }}_2{\bar{\sigma }}_3\)-plane and \({\bar{\sigma }}_1{\bar{\sigma }}_3\)-plane and \({\bar{\sigma }}_2{\bar{\sigma }}_3\)-plane, respectively.

For evaluating the influence of the orthotropy parameters \(\zeta \), \(\xi \) and \(\eta \) (11) on the yield surface (13), sections of plane stress conditions, similar to the ones of Fig. 4 are considered. However, for illustrating simple shear stress states in the \({\varvec{e}}^{\prime (i)}{\varvec{e}}^{\prime (j)}\)-planes, \(i\ne j\), of the principal directions of the material, the latter are rotated by \(45^\circ \) around the normal vector of the corresponding plane, i.e., the \({\varvec{e}}^{\prime (1)}{\varvec{e}}^{\prime (2)}\)-plane is rotated around the \({\varvec{e}}^{\prime (3)}\)-axis (Fig. 5 (left)), the \({\varvec{e}}^{\prime (1)}{\varvec{e}}^{\prime (3)}\)-plane around the \({\varvec{e}}^{\prime (2)}\)-axis (Fig. 5 (middle)) and the \({\varvec{e}}^{\prime (2)}{\varvec{e}}^{\prime (3)}\)-plane around the \({\varvec{e}}^{\prime (1)}\)-axis (Fig. 5 (right)). For the sake of illustration, the orthotropy parameters influencing the shear behavior are set to \(\zeta = 1.0, \xi = 1.2\) and \(\eta = 0.8\) and the parameters \(\alpha , \beta \) and \(\gamma \) set to 1.0.

Fig. 5
figure 5

Plane stress sections through the orthotropic yield surface (orange line) and the equivalent isotropic yield surface (blue line) in the \({\bar{\sigma }}_1{\bar{\sigma }}_2\)-plane (left), \({\bar{\sigma }}_1{\bar{\sigma }}_3\)-plane (middle) and \({\bar{\sigma }}_2{\bar{\sigma }}_3\)-plane (right) of the global coordinate system for principal material directions \(45^\circ \) rotated around the respective normal vector of the corresponding plane; the orthotropy parameters are set to \(\zeta = 1.0, \xi = 1.2, \eta = 0.8, \alpha = \beta = \gamma = 1.0\) (colour figure online)

Full size image

From Fig. 5 (left) follows that the yield function in the \({\bar{\sigma }}_1{\bar{\sigma }}_2\)-plane is only affected by the parameter \(\zeta \) with a resulting shear strength \(f_\text {sh}^{(12)} = f_\text {sh}^*\) because of \(\zeta = 1.0\). Similarly, Fig. 5 (middle) shows that the parameter \(\xi \) only influences the yield function in the \({\bar{\sigma }}_1{\bar{\sigma }}_3\)-plane by the shear strength \(f_\text {sh}^{(13)} = f_\text {sh}^*/\xi \) for simple shearing in the \({\bar{\sigma }}_1{\bar{\sigma }}_3\)-plane. For the sake of completeness, the influence of the parameter \(\eta \) on the yield function is illustrated in Fig. 5 (right), showing that the shear strength \(f_\text {sh}^{(23)} = f_\text {sh}^*/\eta \) for simple shearing in the \({\bar{\sigma }}_2{\bar{\sigma }}_3\)-plane is controlled by the orthotropy parameter \(\eta \).

In Fig. 6, deviatoric sections through the fully hardened plastic potential function (19) at \(g=0\) are shown in the mapped principal effective stress space (cf. Fig. 6A) and the principal effective stress space (cf. Fig. 6B, C) for principal stress directions aligned with the principal material directions. Figure 6B shows the special case of transversely isotropic behavior with the axis of material symmetry in \({\varvec{e}}^{\prime \,(3)}\)-direction by setting the orthotropy parameters to \(\alpha = \beta = 1.5, \gamma = 1.0\), \(\zeta = \xi = \eta = 1.0\). In Fig. 6C, an orthotropic parameter configuration, i.e., \(\alpha = 1.0\), \(\beta = 0.5, \gamma = 1.5, \zeta = \xi = \eta = 1.0\) is shown.

The directions of the plastic strain rate in the isotropic reference configuration \(\partial \,g/\partial \,\bar{\varvec{\sigma }}^*\) (cf. Fig. 6A) and in the real configuration \(\partial \,g/\partial \,\bar{\varvec{\sigma }}= \mathbb {O}:\partial \,g/\partial \,\bar{\varvec{\sigma }}^*\) (cf. Fig. 6B, C) are illustrated for uniaxial compressive stress states in the \({\varvec{e}}^{\prime \,(1)}\)-direction (indicated by the orange arrow), \({\varvec{e}}^{\prime \,(2)}\)-direction (indicated by the gray arrow) and \({\varvec{e}}^{\prime \,(3)}\)-direction (indicated by the blue arrow).

Fig. 6
figure 6

Deviatoric sections of the fully hardened plastic potential function \(g\left( {\bar{\sigma }}_\text {m}^*,{\bar{\rho }}^*,{\bar{\theta }}^*,q_\text {h}= 1\right) =0\) in the mapped principal effective stress space (A) and the principal effective stress space for the special case of transversely isotropic behavior (B) and for orthotropic behavior (C). Principal stress directions are aligned with the principal directions of the material

Full size image

Figure 6A shows that in the isotropic reference configuration the directions of plastic strain rates for the analyzed stress states are aligned with the projected axes of the mapped principal effective stress space. This is due to the fact that for isotropic materials the plastic strain components in directions perpendicular to the loading direction are equal, independent from the direction of axial loading. Considering, e.g., uniaxial compression in \({\varvec{e}}^{\prime \,(3)}\)-direction, negative plastic strains are induced in \({\varvec{e}}^{\prime \,(3)}\)-direction, i.e., \({\dot{\varepsilon }}^{(p)\,*}_3 = {\dot{\lambda }}\,\partial \,g/\partial \,{\bar{\sigma }}^*_3 < 0\) accompanied by equal positive values of the plastic strains in \({\varvec{e}}^{\prime \,(1)}\)- and \({\varvec{e}}^{\prime \,(2)}\)-direction, i.e., \({\dot{\varepsilon }}^{(p)\,*}_1 = {\dot{\varepsilon }}^{(p)\,*}_2 = {\dot{\lambda }}\,\partial \,g/\partial \,{\bar{\sigma }}^*_1 = {\dot{\lambda }}\,\partial \,g/\partial \,{\bar{\sigma }}^*_2 > 0\).

For transversely isotropic behavior with the axis of material symmetry in \({\varvec{e}}^{\prime \,(3)}\)-direction (cf. Fig. 6B), for uniaxial compression in \({\varvec{e}}^{\prime \,(3)}\)-direction (indicated by the blue arrow) the same behavior as for the isotropic reference material is obtained, i.e., the direction of the plastic strain rate is still aligned with the projected axis of the principal stress space. This follows from the fact that the plastic strain components in transverse directions, i.e., in \({\varvec{e}}^{\prime \,(1)}\)- and \({\varvec{e}}^{\prime \,(2)}\)-direction, are in the isotropic plane, and thus are equal, i.e., \({\dot{\varepsilon }}^{(p)}_1 = {\dot{\varepsilon }}^{(p)}_2 = {\dot{\lambda }}\,\alpha \cdot \partial \,g/\partial \,{\bar{\sigma }}^*_1 = {\dot{\lambda }}\,\beta \cdot \partial \,g/\partial \,{\bar{\sigma }}^*_2 > 0\). However, for uniaxial compression in \({\varvec{e}}^{\prime \,(1)}\)- (indicated by the orange arrow) or \({\varvec{e}}^{\prime \,(2)}\)-direction (indicated by the gray arrow), the directions of plastic strain rates change with respect to the corresponding stress axis due to the different material behavior in \({\varvec{e}}^{\prime \,(3)}\)-direction. Thus, since the orthotropy parameters in the isotropic plane are chosen as \(\alpha = \beta = 1.5\) and \(\gamma = 1.0\), the plastic strain component in \({\varvec{e}}^{\prime \,(3)}\)-direction is smaller than the respective transverse plastic strain components, i.e., for uniaxial compression in \({\varvec{e}}^{\prime \,(1)}\)-direction, \({\dot{\varepsilon }}^{(p)}_2 = {\dot{\lambda }}\,\beta \cdot \partial \,g/\partial \,{\bar{\sigma }}^*_2> {\dot{\varepsilon }}^{(p)}_3 = {\dot{\lambda }}\,\gamma \cdot \partial \,g/\partial \,{\bar{\sigma }}^*_3 > 0\) and for uniaxial compression in \({\varvec{e}}^{\prime \,(2)}\)-direction, \({\dot{\varepsilon }}^{(p)}_1 = {\dot{\lambda }}\,\alpha \cdot \partial \,g/\partial \,{\bar{\sigma }}^*_1> {\dot{\varepsilon }}^{(p)}_3 = {\dot{\lambda }}\,\gamma \cdot \partial \,g/\partial \,{\bar{\sigma }}^*_3 > 0\).

For orthotropic plastic behavior (cf. Fig. 6C), it is seen that the directions of plastic strain rates for uniaxial loading in all three considered directions are inclined with respect to the corresponding projected axis of the principal stress space due to the different material behavior in each principal direction of the material.

Summarizing, the introduction of the orthotropy parameters to the mapping tensor (5) allows to control the orthotropic material behavior with respect to the uniaxial compressive strengths as well as the shear strengths and thus to limit the elastic domain by an orthotropic yield surface (cf. Figs. 3, 4, 5). Furthermore, it has been shown that the orthotropy parameters significantly influence the plastic potential function (cf. Fig. 6) and consequently lead to a direction-dependent evolution of the plastic strain.

4 Calibration and validation of the 3DPCDP model for extrusion-based 3D printed concrete

In this section, the capability of the 3DPCDP model to reproduce the orthotropic, nonlinear mechanical behavior of 3DPC is evaluated in a numerical study based on the experimental campaign documented in [32]. Therein, 3DPC specimens were fabricated by extrusion-based 3D printing, applying the fresh concrete layer-by-layer through a printing nozzle. First, the parameters of the 3DPCDP model are determined from the available experimental data on cast and 3D printed concrete cube specimens. Second, for validation of the 3DPCDP model, 2D finite element simulations of three-point bending tests with different orientations of the principal directions of the material with respect to the loading direction are performed. The numerical simulations presented within this section and the subsequent section were performed using the open-source finite element code EdelweissFE [60].

4.1 Parameter identification

The material parameters of the 3DPCDP model are identified by data from uniaxial compression tests performed on cube specimens with a side length of 40 mm as documented in [32] and illustrated in Fig. 7. The concrete cubes, fabricated from a standard 3D printing mortar C35/45—CEM I with a maximum aggregate size of 1 mm, were produced by cutting from layered extrusion-based 3D printed specimens and by conventional casting.

Fig. 7
figure 7

Illustration of the uniaxial compression tests, documented in [32], performed for different orientations \(\alpha _\text {lay}\) of the printing layers with respect to the direction of axial loading

Full size image

For the layered 3DPC, the tests were performed for uniaxial compression in \({\varvec{e}}^{\prime (1)}\)-, \({\varvec{e}}^{\prime (2)}\)- and \({\varvec{e}}^{\prime (3)}\)-direction as well as inclined by \(45^\circ \) with respect to the \({\varvec{e}}^{\prime (3)}\)-direction, i.e., for the orientation of the printing layers \(\alpha _\text {lay} = 45^\circ \) as illustrated in Fig. 7 and documented in [32].

The average value of the measured uniaxial compressive strength of the cast concrete is \(f_\text {cu}^* = 51.03\,\)MPa. Combined with the averaged values \(f_\text {cu}^{(1)} = 47.26\,\)MPa, \(f_\text {cu}^{(2)} = 42.52\,\)MPa and \(f_\text {cu}^{(3)} = 49.16\,\)MPa of the 3D printed concrete, the orthotropy parameters \(\alpha , \beta \) and \(\gamma \) are calculated following (10) as \(\alpha = 1.08, \beta = 1.2\) and \(\gamma = 1.04\). According to [65], the equivalent isotropic uniaxial tensile strength \(f_\text {tu}^*\) is assumed as \(f_\text {cu}^*/10\), the equivalent isotropic uniaxial compressive yield stress \(f_\text {cy}^*\) as \(f_\text {cu}^*/3\) and the equivalent isotropic equi-biaxial compressive strength \(f_\text {cb}^*\) is set to \(1.16\,f_\text {cu}^*\). For the hardening parameters \(A_\text {h}, B_\text {h}, C_\text {h}\) and \(D_\text {h}\) in (25), the damage parameter \(A_\text {s}\) in (30), and the model parameter \(D_\text {f}\) in (22), the default values according to [19] are employed.

The orthotropy parameter \(\eta \), controlling the shear behavior in the \({\varvec{e}}^{\prime \,(2)}{\varvec{e}}^{\prime \,(3)}\)-plane, is determined from the uniaxial compression tests for the loading direction \(45^{\circ }\) inclined to the \({\varvec{e}}^{\prime \,(3)}\)-direction, i.e., \(\alpha _\text {lay} = 45^\circ \) (cf. Fig. 7). As documented in [32], the corresponding uniaxial compressive strength was determined to 44.48 MPa. Thus, following the alternative calibration procedure for the orthotropy parameters \(\zeta , \xi \) and \(\eta \) in Appendix C, the orthotropy parameter \(\eta \) is calibrated to 1.13. Combined with the equivalent isotropic shear strength \(f_\text {sh}^* = 5.06\,\)MPa according to (A4), this results in the shear strength in the \({\varvec{e}}^{\prime \,(2)}{\varvec{e}}^{\prime \,(3)}\)-plane of \(f_\text {sh}^{(23)} = 4.48\,\)MPa. Due to the limited amount of test results and the resulting lack of information on the shear strength in the \({\varvec{e}}^{\prime \,(1)}{\varvec{e}}^{\prime \,(3)}\)-plane, the orthotropy parameter \(\xi \) is assumed to be equal to \(\eta \), i.e., \(\xi = 1.13\) and the shear strength in the \({\varvec{e}}^{\prime \,(1)}{\varvec{e}}^{\prime \,(2)}\)-plane is assumed to be identical to the equivalent isotropic shear strength, resulting in \(\zeta = 1.0\).

4.2 2D finite element simulations of three-point bending tests on 3DPC specimens

The 3DPCDP model is validated by 2D finite element simulations of three-point bending tests (cf. Fig. 8 (left)), documented in [32], on both, cast concrete specimens and extrusion-based 3DPC specimens for loading directions parallel (\(\alpha _\text {lay} = 90^\circ \)), perpendicular (\(\alpha _\text {lay} = 0^\circ \)) and inclined by \(\alpha _\text {lay} = 45^\circ \) to the printing layers.

Fig. 8
figure 8

Test setup of the three-point bending tests [32] (left) and the associated 2D finite element model (right)

Full size image

The dimensions of the specimens were 40 \(\times \) 40 \(\times \) 160 mm and the same concrete mixture as for the uniaxial compressive tests, described in Sect. 4.1, was used. The distance between the supports was \(110\,\)mm. The tests were performed displacement-controlled by applying the displacement \(u_3\) at the middle of the top surface of the specimens with a speed of 0.25 mm/min (cf. Fig. 8).

The 2D finite element model (Fig. 8 (right)) is characterized by a structured mesh consisting of reduced integrated 8-node quadrilateral plane stress finite elements (CPS8R). The maximum element size is 8 mm, and for capturing the steep gradients of the nonlocal softening variable \({\bar{\alpha }}_\text {d}\) in (28), the mesh is refined to a minimum element size of 0.625 mm in the center region of the specimens, corresponding to the zone where the crack is expected to propagate.

The additional parameters for the over-nonlocal implicit gradient-enhanced regularization of the softening behavior in (26), (27) and (28) are chosen as \(\varepsilon _\text {f}^* = 3\cdot 10^{-2}, m = 1.05\) and \(l = 1.25\,\)mm, preserving the specific mode I fracture energy \(G_\text {fI}^* \approx 0.11\,\)N/mm. The latter is estimated for the concrete used for 3D printing of the strength class C 35/45 according to the relation of the CEB-FIB model code [66] by

$$\begin{aligned} G_\text {fI}^* = G_\text {fI,0}^*\left( \frac{f^*_\text {cm}}{f^*_\text {cm0}}\right) ^{0.7}, \end{aligned}$$
(32)

with the base value of the specific mode I fracture energy \(G_\text {fI,0}^* = 0.04\) N/mm, the mean value of the uniaxial compressive strength \(f^*_\text {cm} = 43\) N/mm\(^2\) and \(f^*_\text {cm0} = 10\) N/mm\(^2\).

The experimentally measured load–displacement curves, documented in [32], exhibit a very ductile behavior accompanied by an increasing stiffness at the beginning of loading. This may be caused by (i) influences of the compliance of the frame structure of the testing machine on the measured displacements, (ii) settlements in the porous structure of the concrete or (iii) a non-planar contact area between the loading stamp and the specimens. Thus, in order to obtain comparable results in the numerical simulations, the Young’s moduli of the 3DPCDP model are estimated from the experimental load–displacement curves to \(E_3 = 1800\,\)MPa, \(E_1 = E_2 = 2400\,\)MPa. According to [67] the Poisson’s ratios are assumed as \(\nu _{12} = 0.21, \nu _{13} = \nu _{23} = 0.24\) and the shear moduli \(G_{ij} \in \{G_{12}, G_{13}, G_{23}\}\) are approximated by extending Saint Venant’s formula [45, 68], originally formulated for transversely isotropic materials, to orthotropic material behavior:

$$\begin{aligned} \frac{1}{G_{ij}} \approx \frac{1}{E_i} + \frac{1}{E_j} + 2\,\frac{\nu _{ij}}{E_j} \hspace{1cm}\text {for} \ i \ne j. \end{aligned}$$
(33)

For modeling the isotropic material behavior of the specimens made of cast concrete, the elastic properties and the orthotropy parameters are set to \(E_1 = E_2 = E_3 = 1800\,\)MPa, \(\nu _{12} = \nu _{13}=\nu _{23} = 0.25\) and \(\alpha = \beta = \gamma = \zeta = \xi = \eta = 1.0\), respectively.

In Fig. 9, the computed load–displacement curves (orange curves) are compared to the experimental data (discrete points) [32] for the specimens made of cast concrete used for 3D printing (A) and the specimens with layered 3DPC for loading directions parallel (B), perpendicular (C) and \(45^\circ \) inclined (D) to the printing layers. Due to the aforementioned unrealistic result of a gradual stiffness increase of the experimentally determined curves at the beginning of loading, the latter are considered only after this initial loading phase.

Fig. 9
figure 9

Computed load–displacement curves (orange curves) of the three-point bending tests compared to the experimental data (discrete points) [32] for the specimens made of cast concrete used for 3D printing (A) and the layered 3DPC specimens for loading directions parallel (B), perpendicular (C) and \(45^\circ \) inclined (D) to the printing layers (colour figure online)

Full size image

It can be seen that in the pre-peak regime the numerically determined load–displacement curves are in good agreement with the experimental data for the cast concrete specimens and the layered 3DPC specimens for all investigated loading directions. In addition, the load-bearing capacities \(P_{\max }\) as well as the corresponding displacement \(u_\text {p}\) of the 3DPC specimens including the direction-dependent variation are predicted satisfactorily with a maximum deviation of \(\approx \,8\%\). Unfortunately, the experimental data do not provide detailed information on the softening response, indicated by the brittle decrease of the load after the load-bearing capacity was reached.

In Fig. 10 (right), the evolution of the stress component \(\sigma _{22}\) and the damage variable \(\omega \) are shown for the symmetry plane with increasing loading. Correspondingly, in Fig. 10 (left) fully damaged zones of one half of the specimen are illustrated. Since the structural behavior of the cast concrete specimens and the 3DPC specimens is identical, only the respective distributions for the former are shown.

Fig. 10
figure 10

Fully damaged zones in the right half of the specimen (left) and evolution of the stress component \(\sigma _{22}\) and the damage variable \(\omega \) in the symmetry plane (right) illustrated from crack initiation (\(u_3 = u_\text {p} = 0.35\,\)mm, blue curves) up to failure at \(u_3 = 0.85\,\)mm (orange curves) (colour figure online)

Full size image

The blue curves represent the stress and damage distributions when the load-bearing capacity is attained, i.e., at the displacement \(u_3 = u_\text {p} = 0.35\,\)mm. Until this stage of loading, damage has not been initiated as the uniaxial tensile strength of \(f_\text {cu}^*/10\) is only reached at the bottom of the cross section in the plane of symmetry. Upon further increase in the top displacement \(u_3\), damage is initiated at the bottom face and the crack propagates toward the upper face. Concomitantly, this results in a progressively reduced load-bearing area and, consequently, a decreasing load (cf. Fig. 9A). Structural failure was indicated in the simulation at a computed displacement of \(u_3 = 0.85\,\)mm (orange curves) with a remaining intact cross sectional area of \(\approx 10 \times 40 = 400\,\) mm\(^2\), by divergence of the Newton solution scheme.

Fig. 11
figure 11

Damage distribution of the deformed 3DPC specimens in the final stage (left) and deformed 3DPC specimens after the three-point bending tests taken from [32] (right)

Full size image

Figure 11 shows the deformed 3DPC specimens in the final stage of the simulation. In the experiments, formation of a macro-crack in the central domain was found to be the dominant failure mode, as indicated by the corresponding photos [32], which is also reflected by the plotted distribution of the damage variable \(\omega \). For loading perpendicular (\(\alpha _\text {lay} = 0^\circ \)) and parallel (\(\alpha _\text {lay} = 90^\circ \)) to the printing layers a symmetric damage distribution is obtained, whereas for loading inclined to the printing layers (\(\alpha _\text {lay} = 45^\circ \)) a slightly asymmetric damage distribution is computed. Hence, the predicted location and orientation of the macro-crack of the three specimens is in good agreement with the experimental results.

For verifying the mesh insensitivity of the results obtained by the 3DPCDP model, a mesh-refinement study is performed considering three different element sizes, i.e., a minimum edge length of 2.5 mm (coarse mesh), of 1.25 mm (medium mesh), and of 0.625 mm (fine mesh). In Fig. 12, the computed damage distribution (left) and the corresponding load–displacement curves (right) are shown for the three investigated element sizes for the loading direction \(45^\circ \) inclined to the printing layers. From Fig. 12 (left) follows that the damaged zone is equal for the medium and the fine mesh. This becomes also evident from the corresponding load–displacement curves in Fig. 12 (right), which demonstrates converging structural load-bearing behavior upon mesh refinement.

Fig. 12
figure 12

Damage distribution (left) and load–displacement curves (right) in the final stage for three different element sizes: coarse (2.5 mm), medium (1.25 mm) and fine (0.625 mm)

Full size image

5 Application of the 3DPCDP model to the structural analysis of an arch printed by binder jetting of cementitious powders

The 3DPCDP model, described in Sect. 2, is applied to 3D finite element simulations of an arch printed by binder jetting of cementitious powders as documented in [33]. In contrast to the extrusion-based printing method, described in Sect. 4, the printing layers are formed by spreading cementitious powder. In the next step, the printing nozzle applies the binder to the powder for constructing the geometry of the layer. This procedure is repeated until the complete structure is printed. Finally, the excess powder is removed. For a detailed explanation of the printing process with cementitious powders, interested readers are referred to [33].

5.1 Calibration of the model parameters

The Young’s moduli \(E_1, E_2, E_3\), the Poisson’s ratios \(\nu _{12}, \nu _{13}, \nu _{23}\) as well as the uniaxial compressive strengths \(f_\text {cu}^{(1)}, f_\text {cu}^{(2)}, f_\text {cu}^{(3)}\) of the orthotropic material are calibrated by uniaxial compressive tests documented in [33]. Due to limitations of the printing capacity of the 3D printer for printing cylindrical specimens [33], the tests were performed on cube specimens with a side length of 50 mm. The specimens were loaded in the principal directions of the material, i.e., the \({\varvec{e}}^{\prime (1)}\)-, \({\varvec{e}}^{\prime (2)}\)- and \({\varvec{e}}^{\prime (3)}\)-direction. For measuring the axial and transverse strains, two pairs of strain gauges were applied at the center of the two opposite faces of each specimen, as described in detail in [33]. The measured material properties are listed in Table 1. Due to lack of information, the shear moduli \(G_{ij} \in \{G_{12}, G_{13}\), \(G_{23}\}\), listed in Table 1, are approximated again by Saint Venant’s formula (33).

Table 1 Measured and identified material parameters for concrete printed with cementitious powders according to the experimental data documented in [33]
Full size table

The uniaxial compressive strengths were obtained by dividing the maximum load by the face area of the cubes of \(2500\,\text {mm}^2\), resulting in \(f_\text {cu}^{(1)} = 16.8\) MPa, \(f_\text {cu}^{(2)} = 11.6\) MPa and \(f_\text {cu}^{(3)} = 13.2\) MPa. In contrast to the extrusion-based 3D printing method presented in the previous Sect. 4, no experimental data of the isotropic reference material are provided in [33]. Thus, for calibrating the orthotropy parameters \(\alpha , \beta \) and \(\gamma \), the \({\varvec{e}}^{\prime \,(1)}\)-direction is defined as reference direction for uniaxial compression, i.e., \(f_\text {cu}^{(1)} = f_\text {cu}^* \Rightarrow \alpha = 1.0\). The resulting orthotropy parameters \(\beta \) and \(\gamma \) are determined according to (10)\(_{2}\) and (10)\(_{3}\) as listed in Table 1. According to [65], the uniaxial compressive yield stress is set to \(f_\text {cy}^* = f_\text {cu}^*/3\), the uniaxial tensile strength to \(f_\text {tu}^* = f_\text {cu}^*/10\) and the equi-biaxial compressive strength to \(f_\text {cb}^*\) = \(1.16\,f_\text {cu}^*\).

The shear strengths \(f_\text {sh}^{(12)}, f_\text {sh}^{(13)}\) and \(f_\text {sh}^{(23)}\) were measured in direct double shear tests as \(f_\text {sh}^{(12)} = 0.63\) MPa, \(f_\text {sh}^{(13)} = 1.26\) MPa and \(f_\text {sh}^{(23)} = 4.63\) MPa. For calibrating the orthotropy parameters \(\zeta , \xi \) and \(\eta \) based on the latter, the fictitious isotropic reference shear strength \(f_\text {sh}^*\) is determined according to (A4) to \(f_\text {sh}^* = 1.67\) MPa. Thus, \(\zeta , \xi \) and \(\eta \) can be calculated according to (11).

As for the extrusion-based 3D printing method, presented in Sect. 4, default values according to [19] are employed for the remaining parameters, i.e., the hardening parameters \(A_\text {h}, B_\text {h}, C_\text {h}\) and \(D_\text {h}\), the damage parameter \(A_\text {s}\) and the model parameter \(D_\text {f}\).

Figure 13 shows a comparison of the experimentally determined stress–strain curves of the uniaxial compression tests, loaded in the principal directions of the material, with the computed stress–strain curves using the calibrated 3DPCDP model. The latter are obtained by numerical simulations at a single integration point.

Fig. 13
figure 13

Comparison of experimental results by Feng et.al. [33] (discrete data points) with numerical results (continuous lines) of stress–strain curves of uniaxial compression tests on specimens printed with cementitious powders for loading directions aligned with the principal axes of the material: \({\varvec{e}}^{\prime (1)}\)- (left), \({\varvec{e}}^{\prime (2)}\)- (middle) and \({\varvec{e}}^{\prime (3)}\)-direction (right)

Full size image

It can be seen that the computed stress–strain curves are in good agreement with the experimental data for each direction of axial loading. Since the formation of a shear band resulting in failure represents a structural phenomenon, the stress–strain curves are only shown for the pre-peak regime.

5.2 3D finite element simulations of an arch structure

For demonstrating the capabilities of the 3DPCDP model, calibrated in Sect. 5.1, the latter is applied to 3D finite element simulations of a fictitious example of a 3D arch structure, printed by binder jetting of cementitious powders. This example is inspired from the numerical study by Feng et al. [33].

The geometry of the finite element model together with the prescribed boundary conditions is illustrated in Fig. 14.

Fig. 14
figure 14

Sketch of the geometry and boundary conditions of the finite element model with printing layers perpendicular to the base vector \({\varvec{e}}^{(1)}\) (A), \({\varvec{e}}^{(3)}\) (B) and \({\varvec{e}}^{(2)}\) (C) of the global coordinate system

Full size image

The arch is simply supported at the lower end and loaded at the top of the arch by prescribing the vertical displacement \(u_3\). For reducing the number of degrees of freedom, symmetry with respect to the global \({\varvec{e}}^{(1)}{\varvec{e}}^{(3)}-\)plane is exploited by applying the respective boundary conditions at the top of the arch. The geometry is discretized by a structured mesh, consisting of reduced integrated 20-node 3D quadrilateral finite elements with a maximum element size of \(\approx 40\) mm. For capturing the steep gradients of the nonlocal softening variable, the mesh is refined at the potential cracking locations, i.e., in a small region at the top and at an azimuth angle of \(\phi \approx 45^\circ \) of the arch, to a minimum element size of \(\approx 6\) mm. The parameters for regularizing the softening behavior are chosen as \(l = 6\,\)mm, \(m = 1.05\) and \(\varepsilon _\text {f}^* = 9\cdot 10^{-2}\), preserving the assumed mode I fracture energy \(G_\text {fI}^* \approx 130\,\)N/m.

For investigating the influence of the orthotropic material behavior on the load-bearing capacity of the structure, three different orientations of the printed layers are considered according to Fig. 14.

In Fig. 15, the computed displacement vectors at the front face and the distribution of the damage variable \(\omega \) are shown for the directions (A) and (B) of the printed layers in the final stage of the simulation. Since for case (C) a similar failure mode is obtained as for case (B), only the latter is illustrated.

Fig. 15
figure 15

Predicted displacement vectors of the front face and distribution of the damage variable in the final stage of the simulation of the 3DPC arch (deformation scale factor 5) for the printing layers perpendicular to the base vectors \({\varvec{e}}^{\prime \,(1)}\) (A) and \({\varvec{e}}^{\prime \,(3)}\) (B)

Full size image

It can be seen that for the two loading directions (A) and (B) failure is manifested by the accumulation of damage in two fully cracked regions. The location of the second crack, i.e., the one along the sidewalls, depends on the orientation of the printed layers.

This behavior is also reflected in the corresponding load–displacement curves and contour plots of the damage variable \(\omega \) in slices through the 3D FE-model, illustrated in Fig. 16 for the three investigated cases (A), (B) and (C).

Fig. 16
figure 16

Load–displacement curves of the arch for different loading directions with respect to the principal material directions (left) and contour plots of the damage variable \(\omega \) in the final stage (right)

Full size image

From Fig. 16 follows that for the three cases the first crack opens at the top of the arch marked as \(\text {CA}_1, \text {CB}_1\) and \(\text {CC}_1\). However, the location of the second crack strongly depends on the orientation of the printed layers. For case (A), the second crack \(\text {CA}_2\) is predicted at an azimuth angle of \(\phi \approx 60^\circ \) with a crack orientation slightly inclined to the global \({\varvec{e}}^{(2)}\)-direction. In contrast, the locations of the second cracks \(\text {CB}_2\) and \(\text {CC}_2\) for cases (B) and (C) are predicted at azimuth angles of \(\phi \approx 30^\circ \) and \(\phi \approx 40^\circ \), respectively.

The load–displacement curves in Fig. 16 (left) show that for case (A), the load-bearing capacity is attained at \(P_3 = 63.4\,\)kN, corresponding to the propagation of the first crack \(\text {CA}_1\) at the top of the arch. Further increase of the displacement \(u_3\) yields a slightly decreasing load up to a displacement of \(\approx 13\,\)mm, until the second crack \(\text {CA}_2\) has propagated through the entire cross section, initiating a very steep softening behavior until structural failure. The load–displacement curves of cases (B) and (C) are characterized by a similar behavior. In contrast to case (A), the load still increases after the first crack (\(\text {CB}_1\) and \(\text {CC}_1\), respectively) has formed at the top of the arch at \(u_3 \approx 5\,\)mm. The load-bearing capacities of \(63.7\,\)kN for case (B) and \(66.4\,\)kN for case (C) are attained at \(u_3 \approx 18\,\)mm, when the second crack (\(\text {CB}_2\) and \(\text {CC}_2\), respectively) has propagated through the cross section. Further increase in \(u_3\) again yields a steep decrease in the load up to structural failure.

6 Conclusions and outlook

A new 3D printed concrete damage plasticity (3DPCDP) model was presented for modeling the inherently orthotropic and highly nonlinear constitutive behavior of hardened 3D printed concrete. To this end, the well-established isotropic damage-plastic model for conventionally cast concrete (CDP model) by Grassl and Jirásek [19], formulated within the framework of the gradient-enhanced continuum according to Poh and Swaddiwudhipong [20] for obtaining mesh-independent results, is used as basis. The latter was extended to orthotropic damage-plastic material behavior by mapping of the effective stress tensor into a fictitious isotropic configuration following the approach by Boehler [44] based on the introduction of a novel orthotropic linear mapping tensor.

A straightforward calibration procedure for identifying the material parameters governing the orthotropic damage-plastic material behavior of the 3DPCDP model was proposed: The approach for identifying the parameters for controlling the damage-plastic behavior in the fictitious isotropic configuration can be adopted from [19], considering results from standard laboratory tests on cast concrete used for 3D printing. The additionally introduced six orthotropy parameters can be determined from (i) uniaxial compression tests for loading directions in the principal directions of the material and (ii) simple shear tests in the principal planes of the material. Alternatively, in absence of data from (ii), also uniaxial compression tests for loading directions inclined by \(45^\circ \) to the principal material directions can be considered.

The 3DPCDP model was validated by 2D finite element simulations of three-point bending tests considering the experimental campaign documented in [32] for both, cast concrete specimens and 3DPC specimens. Good agreement of the computed load–displacement curves and the predicted failure modes with the experimental results, reported in [32], was obtained. For 3DPC this holds for the three investigated orientations of the printed layers parallel, perpendicular and inclined by \(45^\circ \) to the printing layers.

For demonstrating the versatile applicability of the 3DPCDP model, numerical simulations of 3DPC specimens printed by binder jetting of cementitious powders were performed. The parameters were identified considering the experimental results reported in [33] for uniaxial compression tests and direct double shear tests performed for different loading directions. The capabilities of the 3DPCDP model at the structural level were demonstrated by 3D finite element simulations of an arch structure for different orientations of the printed layers. It was shown that the 3DPCDP model is capable of predicting the orthotropic nonlinear mechanical behavior up to failure in a realistic manner.

Summarizing, the 3D printed concrete damage plasticity model is well suited for predicting the inherent orthotropic nonlinear mechanical behavior of 3DPC structures made by extrusion-based 3D printing or alternative techniques such as binder jetting of cementitious powders. Due to the formulation within the framework of a gradient-enhanced continuum, the model is suitable for large-scale 3D finite element simulations of 3DPC structures by ensuring mesh-independent results in the softening regime. This opens the door to further model extensions such as considering different orthotropic behavior of 3DPC in compression and tension or crack-induced anisotropic behavior.