1 Introduction

AM is a rising technology, increasingly used in industrial applications. Particularly, selective laser beam melting and filament extrusion are widely applied for the production of prototypes, small batches or otherwise not producible pieces. However, AM could also be applied for soft polymers to produce patient-specific implants, as shown in [29], where it is the goal to improve the medical functionality of neuroimplants by patient specific manufacturing.

This requires the processing of room temperature vulcanisation medical grade silicone. In the beginning, the material can be characterised as a viscous fluid, which transforms to a viscoelastic solid due to the building of crosslinks during curing. The challenge during the printing process is to predict the material behaviour correctly that depend on the processing parameters with emphasis on the material spreading during the process. Important processing parameters are the nozzle diameter, the extrusion velocity and the translational velocity of the extruder. To control the material spreading, a high-speed curing is further induced by the application of high local temperatures using an infrared laser.

Since the thermo-chemo-mechanical coupled behaviour during the printing process is not fully understood, simulations should be performed to support the implant development, respectively, the optimisation of their production by AM. In this paper, a large strain material model applicable for material spreading is developed. The aim is to reproduce the phenomenological observation of a reduced material spreading by the application of laser radiation.

To describe the material behaviour of the medical grade silicone, a suitable material model has to be applied. Small strain curing models considering a process-dependent viscoelastic approach are found among others in [1, 8, 17, 33]. The key point of the formulations is to ensure the stress-free curing behaviour, i.e. new crosslinks are formed in an undeformed configuration and thus do not lead to additional stresses in the predeformed configuration. Hereby the material parameters within the formulations depend on an internal chemical variable. An extension to small strain viscoplasticity is exemplarily shown by [18], considering additionally plastic deformations. In contrast to the small strain regiment, large strain curing models are still rare. In [9], a hypoelastic approach and in [10] its extension to viscoelasticity is used, where a superimposed shrinkage function is applied for the right Cauchy–Green tensor. In the formulation of [21], the logarithmic Henky strain is additively decomposed into a mechanical, thermal and chemical part and an exemplary FE-simulation of a deep drawing process is presented. In this paper, the main idea of the multiplicative decomposition of the deformation gradient into its mechanical, thermal and chemical part, as firstly introduced by [19], is used. A feasibility of this approach is shown by the associated implementation of [23]. A similar approach with the extension to weakly compressible material can be found in [16]. In the mentioned contributions, highly dynamic processes during curing, like the underlying extrusion and subsequent material spreading during the AM process, are not considered and firstly investigated in this paper. Furthermore, the curing models have only been applied in classical FEM implementations, which are not suitable for the simulation of extrusion processes. This is why a suitable meshfree solution scheme has to be applied.

Since multiple problems are tackled in this paper, it is divided into five parts. Firstly, the multi-physical coupling and the related formulation of a process-dependent, three-dimensional, weakly compressible, large strain viscoelastic–plastic curing model is presented. Afterwards, the numerical stable enhanced peridynamic correspondence solution scheme is introduced. In the third part, the AM-specific laser is modelled as a volumetric heat source with the penetration depth of the particle spacing and the classical peridynamic bond-based approach for thermal diffusion is presented. In the next step, the developed material model is applied within the meshfree framework using a local–non-local coupling and the applicability of the plasticity-based approach to model material spreading is shown. At the end, the results of first AM simulations are presented, whereby the complex thermo-chemo-mechanical coupled behaviour is correctly captured and a reduced material spreading during the AM process is obtained by the application of laser radiation.

2 Multi-physical coupling

The simulation of the AM process requires the solution of two fundamental continuum mechanical equations. On the one hand, the balance of momentum has to be solved to describe the deformation. On the other hand, the first law of thermodynamics has to be formulated to simulate the temperature evolution. With respect to the initial configuration, the equations are given by

$$\begin{aligned} \begin{aligned} \rho _0 \ddot{\mathbf {u}}&= \mathsf {Div} \mathbf {P} + \rho _0 \mathbf { b}\\ \rho _0\varTheta \dot{s}&= D_{\mathrm{int}} - \mathsf {Div} \mathbf {Q} + \rho _0 r. \end{aligned} \end{aligned}$$
(1)

Here, \(\rho _0\) denotes the volume in the reference configuration, \(\mathbf {u}\) the displacements, \(\mathbf {P}\) the first Piola–Kirchhoff stresses, \(\mathbf {b}\) the body force, \(D_{\mathrm{int}}\) the dissipation, \(\mathbf {Q}\) the Cauchy heat flux and s the entropy. In general, the stress and entropy are material dependent and have to be determined using a suitable model. In this case, the modelling of curing polymers requires the consideration of an underlying exothermic chemical reaction leading to the formation of crosslinks which solidifies the material. Therefore, the chemical variable

$$\begin{aligned} \alpha (t) = \frac{H(t)}{H(\infty )} \in [0,1] \end{aligned}$$
(2)

is introduced. It defines the fraction of accumulated released heat at time t, H(t), divided by the total accumulated released heat, \(H(\infty )\), of the fully cured material. For the considered stoichiometric curing reaction, the fraction is equivalent to the mass fraction of already cured material, and this is why the chemical variable is denoted as the degree of cure. Due to the modelling of thermo-chemo-mechanical coupled behaviour  the stresses as well as the entropy depend on the deformation, the temperature and the degree of cure.

2.1 Kinematics

To formulate a suitable large strain curing model for the AM process, the idea of a multiplicative split of the deformation gradient into a mechanical, a chemical and thermal part

$$\begin{aligned} \mathbf {F} = \mathbf {F}_{\mathrm{M}} \mathbf {F}_{C} \mathbf {F}_{\varTheta } \end{aligned}$$
(3)

is adopted from [19]. For the chemical as well as for the thermal deformation gradient, isotropic volumetric approaches are introduced

$$\begin{aligned} \mathbf {F}_{C} = J_C^{1/3} \mathbf {1}, \quad \mathbf {F}_{\varTheta } = J_{\varTheta }^{1/3} \mathbf {1}, \end{aligned}$$
(4)

where the volumetric change is governed by

$$\begin{aligned} \begin{aligned} J_C&= 1+ \beta _{c} \alpha \quad \text {with }\quad \beta _{c} \le 0\\ J_{\varTheta }&= 1+ \beta _{\varTheta }(\varTheta - \varTheta _0) \quad \text {with }\quad \beta _{\varTheta }\ge 0. \end{aligned} \end{aligned}$$
(5)

Within the functions \(J_C\) and \( J_{\varTheta }\) \(\varTheta _0\) is the reference temperature, \(\beta _{c}\) the shrinkage and \(\beta _{\varTheta }\) the thermal expansion parameter. Thus, the effects of chemical shrinkage and thermal expansion are induced by the related deformation gradients. For the sake of simplicity, it is assumed that the chemical shrinkage and the thermal expansion are not directly coupled.

Fig. 1
figure 1

Kinematics of large strain curing model

Full size image

The additional kinematics of the underlying mechanics are explained in Fig. 1, where the consecutive split of the deformation gradient (a) and a related rheological model (b) are depicted. The rheological model consists of a process-dependent Maxwell model, serially connected with a process-dependent frictional element. The shear modulus of the equilibrium spring depends on the degree of cure and temperature to model the solidification and temperature-dependent stiffness of the material. To model process-dependent viscous behaviour, the relaxation times of the individual Maxwell-branches are assumed to depend on the degree of cure and on the temperature. To allow fluid-like behaviour at the beginning of the process, the yield stress \(\sigma _y\) is formulated with respect to the degree of cure. Up to the gelation point \(\alpha < \alpha _{\mathrm{gel}}\), the material is supposed to undergo mainly plastic deformations so that fluid-like behaviour is induced. From a micro-mechanical point of view, cf. [6], these plastic deformations can be motivated by a sliding of uncured polymer chains in the polymer melt. Thus, the macro-mechanical deformation shows a plastic behaviour until the gelation point is reached and the polymer chains cannot move freely anymore. After the gelation point is reached, the behaviour is solely characterised as viscoelastic, i.e. \(\sigma _y(\alpha \ge \alpha _{\mathrm{gel}}) = \infty \). Consequently, the transition from fluid to solid-like behaviour is achieved by an increasing yield stress which reduces the rate of plastic deformations to zero until the gelation point is reached as well as by an increasing shear modulus.

In the consecutive split of the deformation gradient, the resulting kinematics for the large strain model are depicted in Fig. 1b and in Fig. 1a an associated rheological model is depicted. As described in Eq. (3), the deformation gradient is firstly decomposed into a mechanical, chemical and thermal part. Due to the serial connection of the generalised Maxwell model with the frictional element, \(\mathbf {F}_{\mathrm{M}}\) is decomposed into a viscoelastic part \(\mathbf {F}_{\mathrm{M}}^{\mathrm{ve}}\) and a dissipative plastic part \(\mathbf {F}_{\mathrm{M}}^{p}\).

The remaining viscoelastic deformation is now split in two different ways related to the viscoelastic approach at large strains and the consideration of weakly compressible material. The first split is related to a volumetric isochoric split \(\mathbf {F}_{\mathrm{M}}^{\mathrm{ve}} = \mathbf {F}_{\mathrm{M}}^{\mathrm{ve, \,vol}} {\bar{\mathbf {F}}}_{\mathrm{M}}^{\mathrm{ve}}\) of the deformation within the equilibrium spring to model weakly compressible material. Thus, both deformations share the elastic stored energy. The second split belongs to the n Maxwell-branches. In each branch, the viscoelastic deformation is split into an elastic \(\mathbf {F}_{M_{el_i}}^{\mathrm{ve}}\) and an inelastic part \(\mathbf {F}_{M_{in_i}}^{\mathrm{ve}}\), where the inelastic part is related to the viscous deformation in the dashpot.

2.2 Evolution equations and process dependencies

For the underlying model, three different kinds of evolution equations have to be solved. The first one is related to the degree of cure, cf. Eq. (2), for which different possibilities exist to formulate associated evolution equations, cf. [12]. Within this work, the dual Arrhenius approach of [11], respectively, [28]

$$\begin{aligned} {\dot{\alpha }}(\varTheta ) = [A_1(\varTheta ) + A_2(\varTheta ) \alpha ^m](1-\alpha )^n \end{aligned}$$
(6)

with

$$\begin{aligned} A_1(\varTheta )= A_{c_1} e^{-\frac{B_1}{\varTheta }} \quad \text {and} \quad A_2(\varTheta )= A_{c_2} e^{-\frac{B_2}{\varTheta }} \end{aligned}$$
(7)

is applied. Herein \(A_{c_1}, B_1, A_{c_1}\) and \(B_2\) are material parameters influencing the curing kinematics.

The second evolution equation considers the rate of the mechanical plastic deformation gradient \(\mathbf {F}_{\mathrm{M}}^p\). Therefore, an associated plasticity model using the von Mises yield criterion

$$\begin{aligned} \varPhi = \sqrt{\frac{3}{2}} \Vert \varvec{\tau }_{\mathrm{M}}'\Vert - \varvec{\sigma }_y(\alpha ) \quad \text {with} \quad \varvec{\tau }_{\mathrm{M}}'= \varvec{\tau }_{\mathrm{M}} - \frac{1}{3} \mathrm {tr}(\varvec{\tau }_{\mathrm{M}}) \mathbf { 1}. \end{aligned}$$
(8)

Therein \(\varvec{\tau }_{\mathrm{M}}\) denotes the mechanical Kirchhoff stresses and \(\varvec{\tau }_{\mathrm{M}}'\) its deviatoric part. In the current configuration, it yields for the associated flow rule

$$\begin{aligned} \mathbf {d}_{\mathrm{M}}^{p} = {\dot{\lambda }} \mathbf {n} \quad \text {with} \quad \mathbf {n} = \frac{\partial \varPhi }{\partial \varvec{\tau }_{\mathrm{M}} }, \end{aligned}$$
(9)

where \(\mathbf {d}_{\mathrm{M}}^{p}\) is the mechanical plastic rate of the deformation tensor. This evolution equation has to be solved ensuring the classical Karush–Kuhn–Tucker conditions (\(\lambda \ge 0, \ \varPhi \ge 0, \ \lambda \varPhi = 0\)). In contrast to classical formulations, the yield function depends on the degree of cure

$$\begin{aligned} \sigma _y(\alpha ) = {\left\{ \begin{array}{ll} \sigma _0 +(\sigma _{\mathrm{gel}}- \sigma _0) \frac{\alpha }{\alpha _{\mathrm{gel}}}&{} \quad \text {if} \ \alpha \le \alpha _{\mathrm{gel}}\\ \infty &{} \quad \text {otherwise} \end{array}\right. } \end{aligned}$$
(10)

and constitutes a new modelling approach of curing polymers. The last evolution equation is related to the elastic-inelastic split of the viscoelastic deformation gradient within each Maxwell-branch i. Within this paper, the formulation of finite linear viscosity

$$\begin{aligned} \dot{\mathbf {C}}_{M_{in_i}}^{\mathrm{ve}} = \frac{1}{\tau (\alpha ,\varTheta )}\left( \mathbf {C}_{\mathrm{M}}^{\mathrm{ve}}-\mathbf {C}_{M_{in_i}}^{\mathrm{ve}}\right) \quad \text {with} \quad \tau (\alpha ,\varTheta ) = \tau = \text {const.} \end{aligned}$$
(11)

is applied. Initially, the process dependency of the viscosities is neglected, such that a constant relaxation time will be applied for the simulations.

The degree of cure and temperature-dependent behaviour of the shear modulus is expressed by two independent functions

$$\begin{aligned} \mu (t) = \mu (\varTheta (t), \alpha ) = f^\mu _{\varTheta } \mu _\alpha . \end{aligned}$$
(12)

The temperature-dependent part is modelled by

$$\begin{aligned} f^\mu _{\varTheta } = \left( \frac{\varTheta _0}{\varTheta (t)}\right) ^{p^\mu _\varTheta }, \end{aligned}$$
(13)

where \(p^\mu _\varTheta \) is the exponent of the temperature dependency. Furthermore, the shear modulus is modelled as a function of cure via

$$\begin{aligned} \mu _\alpha (\alpha ) = p^\mu _{\alpha _1}(1-\alpha ) + p^\mu _{\alpha _2} \alpha \end{aligned}$$
(14)

to exhibit linear growth, cf. [8]. The parameters \(p^\mu _{\alpha _1}\) and \(p^\mu _{\alpha _2}\) represent the temperature-independent shear moduli for the uncured and fully cured material.

2.3 Free energy function

Within the thermo-chemo-mechanical coupled framework, the specific Helmholtz free energy is additively decomposed into a mechanical and a thermo-chemical part

$$\begin{aligned} \varPsi :=\varPsi _{\mathrm{M}} + \varPsi _{\varTheta C}. \end{aligned}$$
(15)

Herein, referring to [15], the specific thermo-chemical free energy function follows the experiment-based ansatz

$$\begin{aligned} \varPsi _{\varTheta C}(\alpha , \varTheta )= & {} h_{F0} + \Delta h_{FS} \alpha - \left( \frac{1}{2} a_F \varTheta ^2+\frac{1}{6} b_F \varTheta ^3\right) (1-\alpha )\nonumber \\&- \left( \frac{1}{2} a_S \varTheta ^2+\frac{1}{6} b_S \varTheta ^3\right) \alpha . \end{aligned}$$
(16)

The parameters can be fitted by the means of DSC experiments. The specific mechanical free energy is now deduced from the one-dimensional rheological model. It is computed by the equivalent specific free energy of the equilibrium spring and the sum of specific free energies from the Maxwell-branches

$$\begin{aligned} \varPsi _{\mathrm{M}} = \varPsi _{M_\infty } + \sum _{i=1}^{n} \varPsi _{M_i}. \end{aligned}$$
(17)

To model weakly compressible material behaviour, the equilibrium part of the mechanical free energy is additively decomposed by

$$\begin{aligned} \rho _0 \varPsi _{M_\infty }= & {} \rho _0 \varPsi _{M_\infty ^{\mathrm{iso}}}({\bar{\mathbf {C}}}_{\mathrm{M}}^{\mathrm{ve}}, \varTheta )+\rho _0 \varPsi _{M_\infty ^{\mathrm{vol}}}(J_{\mathrm{M}}^{\mathrm{ve}}) \nonumber \\= & {} \frac{1}{2} \left[ - f_{\varTheta } \int _{-\infty }^{t} \mu _\alpha (s)\left( \frac{\mathrm {d}}{\mathrm {d}s} \bar{\mathbf {C}}_{\mathrm{M}}^{{ve}^-1}(s) \right) \mathrm {d}s \right] :\bar{\mathbf {C}}_{\mathrm{M}}^{\mathrm{ve}} + K w(J_{\mathrm{M}}^{\mathrm{ve}})\nonumber \\ \end{aligned}$$
(18)

into its isochoric and volumetric part. The integral formulation of the isochoric part is ascribed to a upper-convected Maxwell model considering a time-dependent shear modulus by simultaneously neglecting the classical time-dependent viscosity part. By this formulation, the stress-free curing behaviour is ensured, i.e. no stresses are induced by an increasing stiffness when the deformation is kept constant. For a detailed description of the derivation, the reader is referred to [15]. Within Eq. (18) \(f_{\varTheta } = \left( \frac{\varTheta _0}{\varTheta }\right) ^{p_\varTheta }\) denotes a linear temperature function describing the temperature influence on the shear modulus with respect to a reference temperature, K is the bulk modulus and w is the compressible extension term, for which the ansatz of [5]

$$\begin{aligned} w(J_{\mathrm{M}}^{\mathrm{ve}})= \frac{1}{4}\left( {J_{\mathrm{M}}^{\mathrm{ve}}}^2-1\right) -\frac{1}{2} \ln {J_{\mathrm{M}}^{\mathrm{ve}}} \end{aligned}$$
(19)

is used. Due to the definition as mechanical stored free energy, they depend solely on mechanical kinematic quantities or more precisely on kinematic quantities defined in the thermo-chemical plastic intermediate configuration. Even though \({\bar{\mathbf {C}}}_{\mathrm{M}}^{\mathrm{ve}} = {\bar{\mathbf {C}}}^{\mathrm{ve}}\), the isochoric part is defined with respect to the mechanical intermediate configuration to allow a derivation of stresses with respect to the mechanical configuration later on. In contrast, the viscoelastic subscript of the mechanical Jacobean is further neglected due to the incompressibility of the plastic deformation gradient. The still missing free energies of the Maxwell-branches are modelled with the elastic part of the Neo-Hookean free energy by

$$\begin{aligned} \rho _0 \varPsi _{M_i}(\mathbf {C}_{M_{el_i}}^{\mathrm{ve}}) =\frac{\mu _i}{2}\left( I_{\mathbf { C}_{M_{el_i}}^{\mathrm{ve}}} -3\right) . \end{aligned}$$
(20)

2.4 Constitutive equations

Similar to [19], the Clausius–Duhem inequality is formulated with respect to the thermo-chemo-mechanical split (3), whereby more detailed information is provided in the “Appendix”. The mechanical second Piola–Kirchhoff stresses in the plastic intermediate configuration are introduced as pull back of the Cauchy stresses in terms of the viscoelastic mechanical deformation gradient by

$$\begin{aligned} {\tilde{\mathbf {S}}}_{\mathrm{M}} = J_{\mathrm{M}}\mathbf {F}_{\mathrm{M}}^{ve^{-1}}\varvec{\sigma } \mathbf {F}_{\mathrm{M}}^{ve^{-T}}. \end{aligned}$$
(21)

The stress \({\tilde{\mathbf {S}}}_{\mathrm{M}}\) is defined by the differentiation of the free energy function with respect to the viscoelastic mechanical right Cauchy–Green tensor

$$\begin{aligned} {\tilde{\mathbf {S}}}_{\mathrm{M}} = 2 \rho _0 \frac{\partial \varPsi }{\partial \mathbf {C}_{\mathrm{M}}^{\mathrm{ve}}}. \end{aligned}$$
(22)

An overview about the underlying configurations, the stresses in the different configurations and the required deformation gradients for pull-back and push-forward operations are depicted in Fig. 2.

Fig. 2
figure 2

Push-forward and pull-back operators for the introduced intermediate configurations

Full size image

For the free energy at hand, the mechanical second Piola–Kirchhoff stresses in the intermediate configuration are then obtained with respect to the viscoelastic kinematic quantities. From the Clausius–Duhem inequality results, the relationship

$$\begin{aligned} {\tilde{\mathbf {S}}}_{\mathrm{M}}= & {} \frac{2}{J_C J_{\varTheta }} \left[ \left( \frac{\partial \rho _0 \varPsi _{M_\infty }^{\mathrm{iso}}}{\partial {\bar{\mathbf {C}}}_{\mathrm{M}}^{\mathrm{ve}}} {\bar{\mathbf {C}}}^{\mathrm{ve}}_{\mathrm{M}}\right) _{\mathrm{dev}} \mathbf {C}_{\mathrm{M}}^{ve^{-1}}\right. \nonumber \\&\left. +\frac{ \partial \rho _0 \varPsi _{M_\infty }^{\mathrm{vol}}}{\partial J_{\mathrm{M}}} \frac{J_{\mathrm{M}}}{2} \mathbf {C}_{\mathrm{M}}^{ve^{-1}} \right. \nonumber \\&\left. + \sum _{i=1}^{n} \mathbf {F}_{M_{in_i}}^{ve^{-1}} \frac{\partial \rho _0 \varPsi _{M_i}}{\partial \mathbf {C}_{M_{el_i}}^{\mathrm{ve}}} :\mathbb {P}^{\mathrm{T}} \mathbf {F}_{M_{in_i}}^{ve^{-T}} \right] , \end{aligned}$$
(23)

where \(\mathbb {P}^{\mathrm{T}}= J_{\mathrm{M}}^{-\frac{2}{3}}(\mathbb {I}-\frac{1}{3} \mathbf {C}_{\mathrm{M}}^{\mathrm{ve}} \otimes \mathbf {C}_{\mathrm{M}}^{{ve}^{-1}})\) is a fourth-order projection tensor. The first Piola–Kirchhoff stresses required for the simulation are then obtained by

$$\begin{aligned} \mathbf {P} = \mathbf {F}_{\mathrm{M}}^{\mathrm{ve}} {\tilde{\mathbf {S}}}_{\mathrm{M}} \mathbf {F}_{\mathrm{M}}^{\mathrm{ve}} \mathbf {F}^{-T}. \end{aligned}$$
(24)

Besides the stress definition, the evaluation of the Clausius–Duhem inequality leads to the entropy definition

$$\begin{aligned} s = \frac{J_C J_{\mathrm{M}} \mathrm {tr}(\mathbf {S})}{3 \rho _0} \frac{\partial J_{\varTheta }}{\partial \varTheta } -\frac{\partial \varPsi _{M_\infty }^{\mathrm{iso}}}{\partial \varTheta } - \frac{\partial \varPsi _{\varTheta C}}{\partial \varTheta } \end{aligned}$$
(25)

and to the inequality

$$\begin{aligned} \left( \frac{J_{\varTheta } J_{\mathrm{M}} \mathrm {tr}(\mathbf {S})}{3} \frac{\partial J_C}{\partial \alpha } - \frac{\partial \rho _0 \varPsi _{\varTheta C}}{\partial \alpha }\right) {\dot{\alpha }} \ge 0. \end{aligned}$$
(26)

This inequality cannot be fulfilled a priori, but similar to [15] it can be shown that the thermodynamic consistency is ensured within the range of application.

Utilising the first law of thermodynamics the energy equation is derived

$$\begin{aligned} \begin{aligned} c(\alpha ,\varTheta ) {\dot{\varTheta }}&= \left( \frac{J_{\varTheta } J_{\mathrm{M}} \mathrm {tr}(\mathbf {S})}{3 \rho _0} \frac{\partial J_C}{\partial \alpha } - \frac{\partial \varPsi _{\varTheta C}}{\partial \alpha }\right. \\&\quad \left. - \varTheta \frac{J_{\mathrm{M}} \mathrm {tr}(\mathbf {S})}{3 \rho _{0}} \frac{\partial J_{\varTheta }}{\partial \varTheta }\frac{\partial J_C}{\partial \alpha } + \frac{\partial ^2 \varPsi _{\varTheta C}}{\partial \varTheta \partial \alpha } \right) {\dot{\alpha }} \\&\quad + \varTheta \left( \underbrace{\frac{\partial ^2 \varPsi ^{\mathrm{iso}}_{M_\infty }}{\partial \varTheta \partial \mathbf {C}_{\mathrm{M}}^{\mathrm{ve}}} :\dot{\mathbf {C}}_{\mathrm{M}}^{\mathrm{ve}}}_{\text {Gough Joule}} -\frac{J_C J_{\mathrm{M}} \mathrm {tr}(\dot{\mathbf {S}})}{3 \rho _{0}} \frac{\partial J_{\varTheta }}{\partial \varTheta }\right) \\&\quad + \varvec{\Sigma }_{\mathrm{M}} :\mathbf {L}_{\mathrm{M}}^p + \frac{1}{\rho _0} \kappa \Delta \varTheta + r, \end{aligned} \end{aligned}$$
(27)

where the mechanical Mandel stress is defined by \(\varvec{\Sigma }_{\mathrm{M}} = \mathbf {C}_{\mathrm{M}}^{\mathrm{ve}} {\tilde{\mathbf {S}}}_{\mathrm{M}}\) and \(\Delta \) is the Laplace operator resulting by the utelisation of Fourier’s law of heat conduction with a constant heat conductivity \(\kappa \). Moreover, a process-dependent specific heat capacity

$$\begin{aligned} c(\alpha ,\varTheta )= \varTheta \frac{ \partial s}{\partial \varTheta } = \underbrace{-\varTheta \frac{\partial ^2 \varPsi _{\varTheta C }}{\partial \varTheta ^2}}_{c_{\varTheta C}} + \frac{J_C J_{\mathrm{M}} \mathrm {tr}(\mathbf {S})}{3 \rho _{0}} \underbrace{\frac{\partial ^2 J_{\varTheta }}{\partial \varTheta ^2}}_{=0} \underbrace{-\frac{\partial ^2 \varPsi ^{\mathrm{iso}}_{M_\infty }}{\partial \varTheta ^2}}_{c_{\mathrm{M}}} \end{aligned}$$
(28)

is taken into account. In the following, the mechanical part of the specific heat capacity will be neglected. Consequently, the specific heat capacity is equal to the thermo-chemical part \(c_{\varTheta C}\).

3 Peridynamic framework

In the next step, a suitable numerical scheme has to be chosen. Due to large deformations during the extrusion and the subsequent material spreading, mesh-based solution schemes do not appear to be optimal for the simulation of the AM process. Hence, a meshfree solution scheme is advantageous, and peridynamics, firstly introduced by [25], is selected. It is based on non-local particle interactions over a specific radius and is written in terms of an integro-differential equations over a family \(\mathcal {H}\) without spatial derivatives. As depicted in Fig. 3, for a master particle I this family represents the domain of influence and contains all neighbouring particles J which distances are less equal than the horizon size \(\delta \). Referring to [20], a horizon size of \(3\Vert {\Delta x}\Vert \), where \(\Vert {\Delta x}\Vert \) is the nodal spacing, leads to accurate results. Since the framework of peridynamics is of the Lagrangian type, particle interactions do not change during the computation, but the individual families are deforming. In general, deformations of bonds \(\varvec{\xi } = \mathbf {X}'- \mathbf {X}\) are considered, which are defined as the vectors between an initial point \(\mathbf {X}_I\) and all points \(\mathbf {X}\) within its family \(\mathcal {H}_{\mathbf {X}}\). The resulting equation of motion for a state-based model

$$\begin{aligned} \rho _0(\mathbf {X}) \ddot{\mathbf {u}}(\mathbf {X},t)= & {} \int _{\mathcal {H}_{\mathbf {X}}} {\underline{\mathbf {T}}}[\mathbf {X},t] \langle \mathbf {X'}- \mathbf {X} \rangle \nonumber \\&- \mathbf {\underline{T}}[\mathbf {X'},t] \langle \mathbf {X}- \mathbf {X'} \rangle \mathrm {d} V_{\mathbf {X'}} + \rho _0(\mathbf {X}) \mathbf {b}(\mathbf {X},t)\nonumber \\ \end{aligned}$$
(29)

is an integro-differential equation, in which the integral term replaces the divergence of stresses of classical continuum mechanics. It contains the so-called force vector state \({\underline{\mathbf {T}}}[\mathbf {X},t]\), that has to be determined by a constitutive model.

Fig. 3
figure 3

Family \(\mathcal {H}\) of particle I with interacting neighbouring particles J

Full size image

3.1 Correspondence formulation

In this work, an enhanced version of the correspondence formulation is applied for the computation of force states. The original correspondence formulation, firstly introduced in [27], is based on the specific free energy function. Postulating the correspondence of virtual work of a peridynamic and a related continuum mechanical model leads to the definition of

$$\begin{aligned} {\underline{\mathbf {T}}} = {\underline{w}} \langle \varvec{\xi }\rangle \mathbf {P}(\mathcal {F}) \mathbf {K}^{-1} \varvec{\xi } \end{aligned}$$
(30)

for the force vector state. The correspondence formulation includes the scalar state \({\underline{w}} \langle \varvec{\xi }\rangle \), which represents an weighting influence function. By definition it has to be zero outside the family, and it depends only on the magnitude of \(\varvec{\xi }\). Furthermore, it requires the inverse of the so-called shape tensor \(\mathbf {K}\),

$$\begin{aligned} \mathbf {K} = \int _{\mathcal {H}_{\mathbf {X}}} {\underline{w}} \langle \varvec{\xi }\rangle \varvec{\xi } \otimes \varvec{\xi } \mathrm {d} V_{\mathbf {X'}}, \end{aligned}$$
(31)

containing information about the particle distribution in the initial configuration. Moreover, the non-local deformation gradient \(\mathcal {F}\) is defined based on deformed bonds \(\varvec{\zeta }\) by

$$\begin{aligned} \mathcal {F} = \int _{\mathcal {H}_{\mathbf {X}}} {\underline{w}} \langle \varvec{\xi }\rangle \varvec{\zeta } \otimes \varvec{\xi } \mathrm {d} V_{\mathbf {X'}} \mathbf {K}^{-1}. \end{aligned}$$
(32)

Due to the formulation with respect to the non-local deformation gradient, an arbitrary material model from continuum mechanics can be applied within the correspondence formulation.

Referring to [3, 31], among others, the derived meshfree correspondence framework is susceptible to instabilities similar to the appearance of zero-energy-modes, especially in regions of high gradients. The reason behind the instabilities is the definition of the non-local deformation gradient. An overview about different approaches to prevent the instabilities is presented in [7]. A specific approach was introduced in [4] and is further denoted as the enhanced correspondence formulation. The basic idea is to divide the horizon into \(n_s\) subhorizons

$$\begin{aligned} \mathcal {H} = \bigcup \limits _{a=1}^{n_s} \mathcal {H}_a \end{aligned}$$
(33)

and to compute the related non-local deformation gradient \(\mathcal {F}_a\) and shape tensor \(\mathbf {K}_a\) for each subhorizon. As proposed in [7], non-uniformly distributed free energies among the subdivisions are taken into account, leading under preserving linear and angular momentum to the force state

$$\begin{aligned} {\underline{\mathbf {T}}} = \sum _{a=1}^{n_s} \frac{V_{\mathrm{sup}}^a}{V_{\mathrm{fam}}} {\underline{w}}_a \langle \varvec{\xi }\rangle \mathbf {P}_a(\mathcal {F}_a) \mathbf {K}_a^{-1} \varvec{\xi } \end{aligned}$$
(34)

of the related correspondence formulation. Here, the influence function \({\underline{w}}_a\) is defined as

$$\begin{aligned} {\underline{w}}_a\langle \varvec{\xi }\rangle = I_a\langle \varvec{\xi } \rangle {\underline{w}}\langle \varvec{\xi }\rangle , \end{aligned}$$
(35)

where

$$\begin{aligned} I_a\langle \varvec{\xi } \rangle = {\left\{ \begin{array}{ll} 1 &{} \quad \text {if} \ \mathbf {X}' \in \mathcal {H}_{X_a}\\ 0 &{} \quad \text {otherwise}. \end{array}\right. } \end{aligned}$$
(36)

Within the formulation, the partial free energy of a subdomain is assumed to be directly related to its partial volume \(V_{\mathrm{sup}} = \sum _{j=1}^{n_b} V_J\), where \(n_b\) is the number of particles in the respective subdomain. Consequently, the volume of the entire family is defined by \(V_{\mathrm{fam}} = \sum _{a=1}^{n_s} V_{\mathrm{sup}}^a \) and a volume weighting is performed for the computation of the force vector state.

As shown in [7], the application of the enhanced correspondence formulation prevents the appearances of instabilities in the deformation field but is not completely free of non-physical behaviour.

3.1.1 Numerical treatment

Similar to [24, 26], the meshfree discretised equation of motion is obtained for the formulation at hand by

$$\begin{aligned} \begin{aligned} \rho _{0_I} \ddot{\mathbf {u}}_I&= \sum _{J \in \mathcal {H}_{\mathbf {X}_I}} ({\underline{\mathbf {T}}}_I \langle \varvec{\xi }_{IJ}\rangle - {\underline{\mathbf {T}}}_J \langle \varvec{\xi }_{JI}\rangle ) V_J + \rho _{0_I} \mathbf {b}_I\\ =&\sum _{J \in \mathcal {H}_{\mathbf {X}_I}} {\underline{w}}_{IJ} \left( \sum _{a=1}^{n_{s_I}} \frac{V_{\mathrm{sup}}^a}{V_{\mathrm{fam}}} I_{a} \mathbf {P}_{a} \mathbf {K}_{a}^{-1} \right. \\&\quad \left. + \sum _{b=1}^{n_{s_J}} \frac{V_{\mathrm{sup}}^b}{V_{\mathrm{fam}}} I_{b} \mathbf {P}_{b} \mathbf {K}_{b}^{-1}\right) \varvec{\xi }_{IJ} V_J + \rho _{0_I} \mathbf {b}_I. \end{aligned} \end{aligned}$$
(37)

Due to the meshfree particle discretisation, kinematic boundaries cannot be directly imposed since the initial particle positions are defined as centroids of volume elements. Therefore, ghost layers of particles are used. For kinematic boundaries, the number of layers is chosen to be proportional to the horizon size and for a size of \(\delta =m\Vert {\Delta x}\Vert \ \) m layers are applied.

For the time integration of the equation of motion, the Velocity–Stoermer–Verlet scheme is further used. This results into the updates

$$\begin{aligned} \begin{aligned}&\mathbf {x}(t+\Delta t)= \mathbf {x}(t)+ \mathbf {v}(t) \Delta t + \frac{1}{2} \mathbf {a}(t) \Delta t^2\\&\mathbf {v}(t+\Delta t)= \mathbf {v}(t) + \frac{\mathbf {a}(t) + \mathbf {a}(t+ \Delta t)}{2} \Delta t \end{aligned} \end{aligned}$$
(38)

for positions and velocities.

3.2 Non-local energy equation

In the local energy Eq. (27), the partial derivative of the diffusion term \( \kappa \Delta \varTheta \) is now replaced by a non-local counterpart. Additionally, a process-specific laser radiation is modelled as volumetric heat source term r. In the following, the respective parts of the thermal effects within the energy equation are derived separately. In the subsequent Section, the full energy equation is then defined and discretised.

3.2.1 Thermal diffusion

In the derived energy Eq. (27), thermal diffusion is modelled in terms of the Laplacian by \(c_p {\dot{\varTheta }} = \frac{1}{\rho _0} \kappa \Delta \varTheta \), representing a partial differential equation. Motivated by a non-local peridynamic approach, the Laplacian is replaced by an integral equation as described by [2]. In the following, the formulation introduced by [22] is applied. Based on the heat flow scalar state \({\underline{h}}(\mathbf {X},t) \langle \mathbf {X'}- \mathbf {X} \rangle \), the local thermal conduction part is replaced by the non-local thermal diffusion term

$$\begin{aligned} \kappa \Delta \varTheta = \int _{\mathcal {H}_{\mathbf {X}}} \underline{h}[\mathbf {X},t] \langle \mathbf {X'}- \mathbf {X} \rangle - \underline{h}[\mathbf {X'},t] \langle \mathbf {X}- \mathbf {X'} \rangle \mathrm {d} V_{\mathbf {X'}}.\nonumber \\ \end{aligned}$$
(39)

For an applied bond-based approach, i.e. the heat flow scalar state within a bond \(\varvec{\xi }\) only depends on the associated temperatures of \(\mathbf {X'}\) and \(\mathbf {X}\); the following mathematical definition

$$\begin{aligned} {\underline{h}}(\mathbf {X},t) \langle \mathbf {X'}- \mathbf {X} \rangle = -\,{\underline{h}}(\mathbf {X'},t) \langle \mathbf {X}- \mathbf {X'} \rangle \end{aligned}$$
(40)

holds. Consequently, the non-local diffusion is formulated with respect to the pairwise heat flow density function

$$\begin{aligned} f^h(\mathbf {X},\mathbf {X'},t)= & {} {\underline{h}}(\mathbf {X},t) \langle \mathbf {X'}- \mathbf {X} \rangle -{\underline{h}}(\mathbf {X'},t) \langle \mathbf {X}- \mathbf {X'} \rangle \nonumber \\= & {} 2 \underline{h}(\mathbf {X},t) \langle \mathbf {X'}- \mathbf {X} \rangle . \end{aligned}$$
(41)

For its definition, the temperature scalar state is introduced as

$$\begin{aligned} {\underline{\tau }}(\mathbf {X},t)\langle \mathbf {X'}- \mathbf {X} \rangle = \varTheta (\mathbf {X'},t)-\varTheta (\mathbf {X},t). \end{aligned}$$
(42)

Postulating the peridynamic micro-potential

$$\begin{aligned} z = \kappa _m \frac{{\underline{\tau }}^2}{2||\varvec{\xi }||}, \end{aligned}$$
(43)

where \(\kappa _m\) is the so-called thermal micro-conductivity, leads then to the definition

$$\begin{aligned} f^h(\mathbf {X},\mathbf {X'},t) = \kappa _m \frac{\underline{\tau }(\mathbf {X},\mathbf {X'},t)}{||\varvec{\xi }||} \end{aligned}$$
(44)

for the pairwise heat flow density. Taking the peridynamic thermal potential

$$\begin{aligned} Z(\mathbf {X},t) = \frac{1}{2} \int _{\mathcal {H}_{\mathbf {X}}} z(\mathbf {X'},\mathbf {X},t) \mathrm {d} V_{\mathbf {X'}} \end{aligned}$$
(45)

within a complete family into account and comparing the result with its local counterpart leads to the relationship

$$\begin{aligned} \kappa _m = \frac{6 \kappa }{\delta ^4 \pi } \end{aligned}$$
(46)

between the thermal micro-conductivity and the thermal conductivity k. It is obvious that for an incomplete support, a correction is required since the potential is not complete within the family. A common approach is to apply a linear temperature field and to evaluate the peridynamic potential (45) within each family \(\mathcal {H}_{\mathbf {X}_I}\). The correction factors \(g_I\) are then defined by \(g_I = \frac{Z_\infty }{Z_I}\), where \(Z_\infty \) is the potential of a complete family. Within this work, a simpler yet equally effective approach is formulated, considering only the integrated volume over the family \({\bar{V}} = \int _{\mathcal {H}_{\mathbf {X}}} \mathrm {d} V_{\mathbf {X'}}\) and defining the correction factors as \(g_I = \frac{{\bar{V}}_{\mathrm{max}}}{\bar{V}_I}\). This approach is reasonable and more intuitive since the computed correction factors are directly approximated by volume fractions, which in turn are directly related to the fractional thermal potential.

3.2.2 Laser radiation

A model for the laser radiation has to be introduced. Here, a general laser modelling is used, that is not restricted to the application within a peridynamic framework and could be similarly applied within a different meshfree method.

In the present case, an intensity distribution similar to [32] is used for the laser modelling as a volumetric heat source r. The starting point of the formulation is the rate of absorbed heat being equal to the effective, i.e. absorbed laser power \(P^e_{\mathrm{laser}}\). Taking into account a radial intensity function \(I_{\mathrm{rad}}(r)\), the radial dependent power per area

$$\begin{aligned} \frac{P_{\mathrm{laser}}^e(r)}{\mathrm {d}a} = P_{\mathrm{laser}}^e I_{\mathrm{rad}}(r) \end{aligned}$$
(47)

is defined, whereby the intensity function \(I_{\mathrm{rad}}(r)\) has to satisfy the normalisation constraint

$$\begin{aligned} \int _{0}^{\infty } I_{\mathrm{rad}}(r) \ \mathrm {d}a =\int _{r_i}^{r_a} \int _{0}^{2 \pi r } I_{\mathrm{rad}}(r) \ \mathrm {d}\phi \ \mathrm {d}r \overset{!}{=} 1. \end{aligned}$$
(48)

With respect to the specific AM process a radial Gaussian distribution within an annulus with inner radius \(r_i\) and outer radius \(r_a\) is modelled. Thus, the expected value of the associated normal distribution is \(\mu = r_i + \frac{r_a-r_i}{2}\). The standard derivation is now defined with respect to the \(\pm 3 \sigma \)-interval which is supposed to lie within the inner and outer radii. Thus, it yields \(3 \sigma = \frac{r_a-r_i}{2}\), respectively \(\sigma = \frac{r_a-r_i}{6}\). Consequently, the Gaussian radial intensity distribution is given by

$$\begin{aligned} I_{\mathrm{rad}}(r) = I_0 e^{-\frac{(r-\mu )^2}{2 \sigma ^2}}= I_0 e^{-\frac{18 \left( r-\frac{r_i}{2}- \frac{r_a}{2}\right) ^2}{(\frac{r_a}{2}- \frac{r_i}{2})^2}}, \end{aligned}$$
(49)

where the peak intensity \(I_0\) is obtained, using the normalisation constraint (48), by

$$\begin{aligned} I_0 = \frac{3 \sqrt{2}}{\pi ^{3/2}(r_a-r_i)(r_a+r_i) {{\,\mathrm{erf}\,}}\left( {\frac{3}{\sqrt{2}}}\right) }. \end{aligned}$$
(50)

The modelled normalised intensity distribution is illustrated for the parameters \(r_i=1\,{\mathrm{mm}}\) and \(r_a=3\,{\mathrm{mm}}\) in Fig. 4. On the left, the intensity is plotted in a Cartesian frame and on the right over the radius.

Fig. 4
figure 4

Modelled radial intensity distribution for inner and outer radii of \(r_i=1\,{\mathrm{mm}}\) and \(r_a=3\,{\mathrm{mm}}\)

Full size image

Due to the fast decay of the intensity over the thickness, cf. [30], the laser irradiation is modelled as a surface effect, i.e. only particles on the surface are affected. Therefore, it is assumed that a constant intensity is observed within a penetration depth being equal to the nodal spacing \(d_{\mathrm{pen}} = \Vert {\Delta x}\Vert \). This leads to the definition of the radial power per volume

$$\begin{aligned} \frac{P_{\mathrm{laser}}^e(r)}{\mathrm {d}v} = P_{\mathrm{laser}}^e I_{\mathrm{rad}}(r) I_{\mathrm{pen}}, \end{aligned}$$
(51)

where the normalised intensity in the irradiation direction is defined by \(I_{\mathrm{pen}}= \frac{1}{d_{\mathrm{pen}}}\). Consequently, it is possible to model the laser impact as a volumetric heat source term. In the energy Eq. (27), the volumetric heat source term resulting from laser irradiation is defined by

$$\begin{aligned} r = \frac{1}{\rho _0}\frac{P^e_{\mathrm{laser}} I_{\mathrm{rad}}(r)}{J_I d_{\mathrm{pen}}}. \end{aligned}$$
(52)

Since only the upper layer is supposed to be irradiated, a detection of surface particles is required. Therefore, a geometrical algorithm is developed incorporating all particles intersecting the beam area in the x–y plane. A particle is defined as a surface particle if there is no other particle in the projected x–y plane within the distance \(0.5 \Vert {\Delta x}\Vert \) having a higher z-position than the considered particle.

4 Application of large strain curing framework within Peridynamics

4.1 Local–non-local coupling

In the next step, the developed multi-physical large strain curing model is applied within the enhanced correspondence framework. Therefore, based on the local thermo-chemo-mechanical split of the deformation gradient (3) and the non-local nature of deformation gradients within Peridynamics it has to be distinguished between local and non-local tensors. This includes the definition of all multi-physical quantities in the present approach. From the chemical point of view, all related quantities are computed and stored at particles characterising a local level. The same holds for all thermal quantities. Thus, from the split of the non-local deformation gradient

$$\begin{aligned} {\mathcal {F}} = \mathbf {F}_{\mathrm{M}} {\mathbf {F}_{C} \mathbf {F}_{\varTheta }} \end{aligned}$$
(53)

follows a local–non-local coupled computation of mechanical deformation gradients. Additionally, due to the split of families non-local sub-deformation gradients have to be considered, where the chemical and thermal deformation gradients are constant within a family. Consequently, all mechanical related evolution equations have to be solved \(n_{\mathrm{sup}}\) times per particle, while the evolution equation for the degree of cure is only solved once and used for all subfamilies.

Within a time step, the degree of cure is updated first using an Euler forward scheme. The viscoelastic–plastic formulation requires a coupled solution for plastic as well as viscous internal variables Therefore, referring to [13], the flow rule is rewritten in terms of the mechanical plastic deformation gradient

$$\begin{aligned} \dot{\mathbf {F}}_{\mathrm{M}}^{p} = {\dot{\lambda }} \mathbf {F}_{\mathrm{M}}^{ve^{-1}} \mathbf {n} \mathbf {F}_{\mathrm{M}}^{\mathrm{ve}} \mathbf {F}_{\mathrm{M}}^{p} \end{aligned}$$
(54)

and the exponential map integrator is used for implicit time integration. In this framework, the plastic evolution equation

$$\begin{aligned} \underbrace{\mathbf {F}_{\mathrm{M}}^{n+1} \mathbf {F}_{\mathrm{M}}^{p^{-1},n+1}}_{\mathbf {F}_{\mathrm{M}}^{ve,n+1}} =\exp {\Delta \lambda \mathbf {n}} \mathbf {F}_{\mathrm{M}}^{ve,n+1} \end{aligned}$$
(55)

is solved in terms of the mechanical inverse plastic deformation gradient. Furthermore, the viscous evolution equations are approximated by an Euler backwards scheme

$$\begin{aligned} \mathbf {C}_{M_{in_i}}^{ve,n+1} = \frac{\mathbf {C}_{M_{in_i}}^{ve,n}+ \frac{\Delta t}{\tau ^{n+1}} \mathbf {C}_{\mathrm{M}}^{ve,n+1} }{1 + \frac{\Delta t}{\tau ^{n+1}}}. \end{aligned}$$
(56)

As a result, the nine components of the mechanical inverse plastic deformation gradient, the plastic increment and six components of the mechanical viscoelastic-inelastic right Cauchy Green tensor for each Maxwell branch are taken simultaneously into account for the implicit solution. Therefore, the residual is formulated with respect to the implicit formulated Eqs. (55), (56) and the yield function. Due to the complexity of the associated model, the automatic code generator AceGen, cf. [14], is applied to compute the derivatives of the residual with respect to all unknown parameters within a Newton–Raphson algorithm.

Additionally, the energy equation has to be solved within the developed framework. Therefore, the locally derived Eq. (27) is reformulated under the negligence of the Gough Joule effect, replacing the local diffusion part by its local counterpart and adding the volumetric heat source term for the laser radiation into

$$\begin{aligned} \begin{aligned} c_{\varTheta C} {\dot{\varTheta }}_I&= \underbrace{\left( \frac{\varphi {\hat{J}}_{\mathrm{M}} \mathrm {tr}\hat{(\varvec{\sigma })}}{3 \rho _0} \frac{\partial g}{\partial \alpha } - \frac{\partial \varPsi _{\varTheta C}}{\partial \alpha } - \varTheta _I \frac{J_{\mathrm{M}} \mathrm {tr}\hat{(\varvec{\sigma })}}{3 \rho _{0}} \frac{\partial \varphi }{\partial \varTheta }\frac{\partial g}{\partial \alpha } + \frac{\partial ^2 \varPsi _{\varTheta C}}{\partial \varTheta \partial \alpha } \right) {\dot{\alpha }}}_{\text {curing}} \\&\quad \underbrace{-\frac{g \hat{J}_{\mathrm{M}} \mathrm {tr}\hat{\dot{(\varvec{\sigma })}}}{3 \rho _{0}} \frac{\partial \varphi }{\partial \varTheta }}_{\text {multi-physical split}} + \underbrace{(\hat{\varvec{\Sigma }_{\mathrm{M}} :\mathbf {L}_{\mathrm{M}}^p})}_{\text {plasticity}}\\&\quad + \underbrace{\frac{1}{\rho _0} \sum _{J \in \mathcal {H}_{\mathbf {x}_I}} g_{IJ} \kappa _p \frac{\varTheta _J - \varTheta _I}{||\varvec{\xi }_{IJ}||} V_J}_{\text {diffusion}}\\&\quad +\underbrace{\frac{1}{{\hat{J}} \rho _0}\frac{P^e_{\mathrm{laser}} I_{\mathrm{rad}}(r)}{d_{\mathrm{pen}}}}_{\text {laser}}. \end{aligned} \end{aligned}$$
(57)

With the exception of temperature, the subscript I has been omitted for the individual terms in Eq. (57) above. It is clear that all quantities have to be evaluated with respect to particle I on the basis of the meshfree discretisation. Moreover, the distinct origins of the resulting energy equation are indicated with an under-brace. Due to nodal integration and the subdomain decomposition, \(n_{\mathrm{sup}}\) stresses are defined on nodal level. This is why averaged mechanical quantities, denoted by a hat, computed with the volume weighting factors \(\frac{V_{\mathrm{sup}}}{V_{\mathrm{fam}}}\) are used. The averaged mechanical quantities are defined by

$$\begin{aligned} {\hat{\mathbf {A}}} = \sum _{a=1}^{n_s} \frac{V_{\mathrm{sup}}^a}{V_{\mathrm{fam}}} \mathbf {A}^a, \quad \text {respectively} \quad {\hat{A}} = \sum _{a=1}^{n_s} \frac{V_{\mathrm{sup}}^a}{V_{\mathrm{fam}}} A^a. \end{aligned}$$
(58)

4.2 Simulation of material spreading

This section deals with the investigation of the material spreading behaviour of a horizontally orientated cylinder with diameter and length 0.25 mm subjected to gravitational force. As depicted in Fig. 5, the cylinder is discretised with 6640 particles, whereby the cylinder is subjected to a stick contact constraint at \(z=0\). Consequently, particles which z-position comes below \(z=0\) are perpendicular projected to the printing plate at \(z=0\) and are treated as fixed boundary particles. An initial density of \(\rho _0 = 1100\) \(\frac{{\text {kg}}}{{\text {m}}^3}\)is used for the simulations.

Fig. 5
figure 5

Discretised cylinder for spreading simulations

Full size image

4.2.1 Evolution of energies

At first, a classical Neo-Hookean model is compared with the developed large strain curing model without process dependencies and not considering the stress free curing behaviour. Thus, the model is reduced to a viscoelastic–plastic material model with a Neo-Hookean equilibrium free energy within the Maxwell model. The simulations are performed with a Young’s modulus of \(E = 15\,\hbox {Pa}\) and a Poisson’s ratio of \(\nu = 0.499\). To neglect the viscous influence within the viscoelastic–plastic approach, only a single Maxwell-element is taken into account and the related shear modulus \(\mu _1 = 0.05\,\hbox {Pa}\) and viscosity \(\eta _1= 0.035\) Pa s are chosen. Additionally, a constant yield stress of \(\sigma _y = 0.25\,\hbox {Pa}\) is selected. Further, due to the small domain and the contact constraint a time step of \(\Delta t = 10^{-7}\) s is used.

Fig. 6
figure 6

Normalised energy and shape comparison between a Neo-Hookean and a viscoelastic–plastic approach for the simulation of a cylinder under gravity loading

Full size image

In the subfigures of Fig. 6, the resulting shapes (a, b) and normalised energies (c, d) are depicted for the Neo-Hookean and the viscoelastic–plastic approaches. Utilising the Neo-Hookean approach leads after the simulated time of \(0.02\,\hbox {s}\) almost to the initial cylinder (a), while the viscoelastic–plastic approach leads to a spreading of the cylinder with a maximum height of \(0.074\,{\mathrm{mm}}\) (b).

For the further evaluation, the evolution of the total, the potential, the elastic stored and the kinematic energy, is considered. In case of the Neo-Hookean approach, the total energy of the system remains constant overall, whereby minor non-physical oscillations are observable. It is assumed that these are neither related to the material model nor to the enhanced peridynamic correspondence formulation but to the application of the stick contact constraint. Since a simple penalty approach is utilised and the equation of motion is solved in an explicit manner, such deviations in total energy are to be expected. Since the deviations only have a minor influence on the global solution and do not lead to a steady increase or decrease in energy, they are acceptable.

Furthermore, the individual energies of the system are closely connected. Based on the gravitational force and the application of a nearly incompressible material model, the cylinder is compressed in the z-direction and associated deformations in the x-direction and the y-direction are induced. Thus, the potential energy of the system transforms into kinematic and elastic stored energy until the cylinder is not compressed any more. At this point, the kinematic energy is zero, the elastic stored energy is at its maximum, and the potential energy is at its minimum. Thus, the potential energy is purely transformed into elastic stored energy. Afterwards, the behaviour reverses and the elastic stored energy transforms back into potential energy. This typical dynamic behaviour of elastic materials is known as bounce-off and repeats such that the oscillation of the potential energy has a phase-shift of \(\pi \) compared to the oscillation of the elastic stored energy. As a consequence, it is not possible to induce a physically meaningful dynamic spreading behaviour using the Neo-Hookean approach.

In contrast, the desired dissipative spreading behaviour is induced by the utilisation of the viscoelastic–plastic approach. As illustrated in Fig. 6d, the total energy of the system is reduced until an almost steady state is reached at \(t=0.013\,\hbox {s}\). Equivalent to the Neo-Hookean approach, the potential energy of the system is reduced by the gravitational load, whereby the induced deformation is mainly plastic and hence has a dissipative character. Consequently, only a small fraction is transformed into elastic stored energy which is smaller by a magnitude than the potential and kinematic energy and is thus not visible in the plot. Furthermore, a monotonic decrease of potential energy is at hand before an almost unnoticeable oscillation of potential, elastic stored and kinematic energy takes place, being related to the bounce-off behaviour of the remaining elastic stored energy. In summary, the application of the viscoelastic–plastic approach together with the introduction of plastic deformations allows the simulation of spreading of soft materials.

4.2.2 Isothermal temperatures

Within the next simulations, the process-dependent material parameters are applied for the thermo-chemo-mechanical coupled viscoelastic–plastic large strain curing model. The simulations are oriented on a printing approach and consider the spreading kinematics during curing which depends on temperature. Therefore, the curing of the horizontally orientated cylinder under gravity loading is simulated for different isothermal temperatures and a time step size of \(\Delta t = 10^{-7}\) s. The applied temperatures are \(\varTheta = 20, 50\) and \(100\,^{\circ }\hbox {C}\). For a better comparison, the shrinkage and thermal expansion parameters \(\beta _C\) and \(\beta _{\varTheta }\) are set to zero. Furthermore, a degree of cure and temperature-dependent shear modulus is used and the gelation point is set to \(\alpha _{\mathrm{gel}}=0.1\). The other applied material parameters are related to the process-dependent material functions and are depicted in Table 1.

Table 1 Material paremeters for isothermal material spreading simulations
Full size table

Note that the utilisation of a constant compression modulus by simultaneously applying an increasing shear modulus automatically leads to a decrease in Poisson’s ratio. This is why a high compression modulus is selected for the simulations. At the reference temperature, the applied parameters lead to a starting Poisson’s ratio of \(\nu (\alpha = 0) = 0.4999\) and to a maximal Poisson’s ratio of \(\nu (\alpha = 1) = 0.4851\). Furthermore, the parameters of evolution Eq. (6) related to curing are shown in Table 2.

Table 2 Material parameters for the evolution equation of the degree of cure
Full size table

They are selected such that fast curing takes place and the model characteristics are verifiable without acknowledging the exact material behaviour. The target is only to show a temperature-dependent spreading behaviour without capturing extensive material spreading. For the thermo-chemical free energy function (16), the parameters of Table 3 are applied.

Table 3 Fitted material parameters for the definition of thermo-chemical quantities based on DSC experiment for acrylic bone cement, cf. [15], whereby \(b_F\) and \(b_S\) have been adjusted
Full size table

In the subfigures of Fig. 7, the deformed cylinders in the final time step at \(t = 0.02\,\hbox {s}\) are depicted for the three mentioned isothermal temperatures. With an increasing temperature from 20 (a) to 50 (b) to \(100\,^{\circ }\hbox {C}\) (c) the spreading of the initial cylinder decreases, i.e. smaller deformations in all directions are observable. Since a higher temperature is associated with a faster curing, the material stiffens and transforms faster to a solid.

Fig. 7
figure 7

Cylinder spreading for three different isothermal temperatures in the final time step at \(t = 0.02\,\hbox {s}\)

Full size image

To quantify the spreading effect, the potential energy of the cylinder and its evolution over time is taken as the initial measurement. This measurement is directly coupled with the final height of the cylinder. The results for the three different isothermal temperatures are depicted in Fig. 8. Here, the potential energy is normalised to allow a better comparison between the different isothermal temperatures.

In alignment with the results of Fig. 7, the highest decrease in potential energy is observed for the lowest temperature of \(20\,^{\circ }\hbox {C}\), where 76.93 % of the potential energy remains after a time of \(t = 0.006\,\hbox {s}\). For a temperature of \(50\,^{\circ }\hbox {C}\) the minimum potential energy is obtained in a reduced time of \(t = 0.0056\,\hbox {s}\) and accounts for 79.38 % of the initial value. An even shorter time and increase in the remaining potential energy is obtained for an isothermal temperature of \(100\,^{\circ }\hbox {C}\) for which a minimum potential energy of 86.41 % in \(t = 0.0044\,\hbox {s}\) is reached. Furthermore, after reaching their minima the potential energies oscillate over time, whereby the period length decreases from low to high temperature and the magnitudes seem to be slightly damped.

Fig. 8
figure 8

Evolution of normalised potential energies over time for different isothermal temperatures

Full size image

The observed oscillations are related to the bounce-off effect of the material. Even though most of the energy is dissipated by plastic deformations during the spreading, elastic energy is still stored to a certain extent. The higher the temperature, the shorter is the period length of the oscillations because of a smaller spreading velocity. Moreover, the damping in the oscillations is associated with the increasing stiffness and by that with increasing stresses during the deformations within the oscillations. An explanation for the slight damping effect are then additional plastic deformations due to a higher von Mises stress leading to additional energy dissipation. This only takes place if the increase in von Mises stress is higher than the increase of the process-dependent yield stress during the oscillations. Thus, it depends on the modelling functions for the shear modulus as well as for the yield stress.

Altogether, the developed material model and its application within the enhanced peridynamic correspondence formulation is suitable to exhibit the required spreading dynamics of uncured material as well as reduced spreading at higher temperatures. Consequently, the developed material model is applicable for the desired simulation of the AM process.

5 Simulation of AM processes

5.1 Software development

In comparison with classical computations in solid mechanics, the simulation of the underlying AM process is related to an extrusion process. Thus, the initial configuration is defined by the material inside the extruder. In order to accelerate the computation, an automatic material addition inside the extruder is implemented instead of modelling the whole dispenser. It is based on the discretisation algorithm developed for the initial material at the beginning of the computation. Due to the circular shape of the nozzle, a circle with its diameter is discretised and then extruded in a third dimension to generate the initial three-dimensional discretisation, further denoted as initial cylinder.

This cylinder is then extruded and as soon as half of it passes the nozzle the initial cylinder is duplicated and put on top of the initial one with a consecutive particle numeration. The initial particle positions of the added cylinder are not defined by their actual position while adding but under the appearance that they have already been initialised in the beginning of the computation. Due to the non-locality of the model, this would usually require an additional family search and the computation of related subshape tensors in the transition zone of \(2 \delta \). The need for this is prevented by the special discretisation, respectively, by the duplication of the two-dimensional discretisation into the third dimension. The required information of family memberships, associated bond vectors and related subshape tensors is obtained for each particle in the transition zone by a perpendicular projection of the individual particles in the extrusion direction and using the information of the first particle which has a compact support in the extrusion direction. Since the discretisation is layer-based, it is sufficient to store the particle numbers of two specific layers from which the required information can be extracted.

5.2 Isothermal processes

This section deals with the application of the developed peridynamic framework for the simulation of the AM process. Isothermal extrusion processes are considered at first applying the developed large strain curing model. The simulations are based on the same initial body, which is depicted in Fig. 9.

Fig. 9
figure 9

Illustration of the initial body for the simulation of AM processes, whereby particles inside the nozzle are depicted in blue, the printing plate in black, the nozzle in light red and the remaining particles are colour coded with respect to their normalised distance to the z-position of the nozzle

Full size image

Since the z-position of the nozzle does not vary, the normalised colour scale remains constant during the simulations, and therefore, the scale in the following figures is omitted. When referring to the scale the phrase ‘normalised distance to the printing plate’ is used to indicate the z-position of the extruded material.

The initial body is generated by a duplication of an initial cylinder and an applied bending, such that one end of the initial body is parallel and the other end perpendicular to the printing plate. Moreover, a stick contact boundary condition is applied as soon as a particle reaches the printing plate. Additionally, the displacement of the most left layer of particles is constrained in the x-direction and a kinematic extruder boundary condition is applied, i.e. the motion of particles inside the nozzle is determined by the horizontal extruder movement speed \(\mathbf {v}_m\) and the extrusion velocity \(v_{\mathrm{e}}\). As soon as a particle leaves the extruder, its boundary condition is removed and its deformation is determined by the equation of motion.

In the following, a fixed translational velocity of \(v_m=0.01\,\frac{{\mathrm{m}}}{\hbox {s}}\) in the x-direction is applied and the extrusion velocity is varied from \(v_{\mathrm{e}}=-\,0.0025\,\frac{{\mathrm{m}}}{\hbox {s}}\) to \(v_{\mathrm{e}}=-\,0.005\,\frac{{\mathrm{m}}}{\hbox {s}}\) and to \(v_{\mathrm{e}}=-\,0.01\,\frac{{\mathrm{m}}}{\hbox {s}}\). In all cases, a total time of \(t = 0.5\,\hbox {s}\) is simulated with a time step of \(\Delta t = 0.5 \times 10^{-7}\,\hbox {s}\). Furthermore, an automatic material addition inside the extruder is formulated. As soon as the amount of material inside the extruder is below a predefined threshold, new material is automatically added.

The diameter of the nozzle and thus the diameter of the initial cylinder equals the cylinder diameter of the previous simulations and amounts to \(0.25\,{\mathrm{mm}}\). Furthermore, the distance between the nozzle and the extrusion plate is chosen to be 1.5 times the discretised cylinder diameter plus the particle spacing. For the discretisation with 332 particles per layer, the distance is \(d_{\mathrm{pn}} = 0.3375\,{\mathrm{mm}}\), whereby the total number of initial particles is 13,280. From a peridynamic point of view, a horizon of \(\delta = 3 \Vert {\Delta x}\Vert \) is used for all simulations and a constant influence function of \({\underline{w}} \langle \varvec{\xi }\rangle = 1\) is applied.

For the large strain curing model, the material parameters of the previous section are applied, cf. Sect. 4.2. An exception is the application of nonzero chemical shrinkage, \(\beta _C = -\,0.03\), and thermal expansion, \(\beta _{\varTheta } = 10^{-5}\) \(\frac{{\text {1}}}{{\text {K}}}\), parameters to include the full thermo-chemo-mechanical coupling of the material model. As a model parameter for the thermal diffusion, a constant thermal conductivity of \(\kappa = 0.27\) \(\frac{{\text {W}}}{{\text {m K}}}\)is applied.

In Figure 10a–c, the final shapes using the distinct extrusion velocities are depicted. To evaluate the results, a quantitative comparison of the final shapes between the distinct extrusion speeds is performed and illustrated in Fig. 10d. Therefore, the deformed diameters of the initially circular particle layers are considered, whereby the two particles with the largest distance in the y-direction are taken as a reference in each layer. As a measurement, the ratio of their distance in the final time step to their initial distance, which is equal to the discretised cylinder diameter, is used. Consequently, values higher than 1 are related to an increased diameter and thus to a certain extent, to material spreading in the y-direction. On the other hand, values smaller than 1 are related to a decreased diameter and thus to a certain extent, to material necking.

Fig. 10
figure 10

Final shapes of 3D-printed material for different extrusion velocities and quantitative comparison of initial to final diameters across the extruded material

Full size image

For the intermediate extrusion speed of \(v_{\mathrm{e}} = -\,0.005\,\frac{{\mathrm{m}}}{\hbox {s}}\), the ratio assumes a value of 1.096 at the minimum x-position and remains then almost constant before approaching the nozzle, whereby a maximum value of 1.108 is reached. Consequently, a slight spreading of the extruded material occurs.

For the fast extrusion speed of \(v_{\mathrm{e}} = -\,0.01\,\frac{{\mathrm{m}}}{\hbox {s}}\) the ratio is 1.810 at the minimum x-position and takes a maximum value of 1.824. Up to \(x= 0.2773\,{\mathrm{mm}}\) the ratio remains above 1.7 before dropping to 1 near the nozzle. This result confirms that doubling the extrusion velocity has a high influence of the deposited material. In the present case, an extensive material spreading is perceived since the diameter of the initially circular body in the y-direction increases by a factor of up to 1.824.

The application of the slow extrusion speed of \(v_{\mathrm{e}} = -\,0.0025\,\frac{{\mathrm{m}}}{\hbox {s}}\) leads to the inverted result. At the minimum x-position, the ratio has a minimum value of 0.7362 before the ratio monotonically increases up to the nozzle. Thus, the previously described necking is quantitatively captured. However, contrary to the intermediate and the fast extrusion speed no plateau of approximately constant ratios is observable. It is assumed that the already-extruded material has not reached its equilibrium state and that a plateau would be reached for the ratios after a longer simulation time. The reason for this assumption is that up to the considered simulation time only a small fraction of the extruded material has reached the printing plate. Thus, the material is still under tension in the x-direction leading to an additional necking over time.

In summary, the simulations of the AM process reveal a high dependency of the manufactured shapes and of the dynamics during the AM process from the extrusion velocity. This supports the initial statement that the results of the considered AM process are highly dependent on operator experience in terms of choosing proper printing parameters. Within this work, only the extrusion speed is varied. However, the distance from the nozzle to the printing plate and the translational velocity are crucial parameters as well. Together with the laser radiation, their combination in close connection with the resulting impact of gravity is the crucial point when it comes to the optimisation of extrusion-based AM processes.

In general, it is possible to capture the high influence of applied printing parameters, in terms of varying only one parameter, within the designed, peridynamic-based AM framework and the application of the developed large strain curing model within.

Fig. 11
figure 11

Temperature distribution as a consequence of laser radiation and thermal diffusion during the AM process

Full size image
Fig. 12
figure 12

Comparison of final shapes of 3D-printed material using the fast extrusion velocity with and without laser radiation

Full size image

5.3 Fully multi-physical coupled extrusion process

In the last example, a fully thermo-chemo-mechanical coupled AM process including laser radiation is considered. Here the goal is to show that the material spreading is reduced in terms of an applied laser radiation. Therefore, a laser power of \(P_{\mathrm{laser}}= 1.5\,\hbox {W}\) with an absorbed fraction of 0.99 is applied. The laser speed \(v_p\) is further defined by the translational velocity of the extruder, and the laser centre is identical with the geometrical centre of the nozzle during the entire simulation. For the inner and outer radii of the laser \(r_i= 0.25\,\hbox {mm}\) and \(r_a= 0.35\,\hbox {mm}\) are chosen, such that the initial particles are not irradiated. Moreover, the fast extrusion velocity of \(v_{\mathrm{e}}= -\,0.01\,\frac{{\mathrm{m}}}{\hbox {s}}\) is applied.

Before the impact of the laser radiation on the dynamics is investigated, a short explanation about its influence on the temperature distribution during the AM process is given. For this purpose, the temperature distributions after half of the simulation time at \(t = 0.025\,\hbox {s}\) and at the final time at \(t = 0.05\,\hbox {s}\) are plotted in the subfigures of Fig. 11. By a progressing nozzle movement, the mechanically linked laser moves with the same velocity, such that the extruded material is continuously irradiated by the trailing edge of the laser. Thus, thermal energy is transferred into the body and the temperature of the irradiated surface particles increases.

At \(t = 0.025\,\hbox {s}\), cf. Fig. 11a, a portion of the extruded material has already been irradiated and the highest temperature of \(\varTheta = 200.46\,^{\circ }\hbox {C}\) is obtained at the intensity peak of the laser. Moreover, the temperature of the previously irradiated particles has decreased due to thermal diffusion. With progressing time, more material is extruded and the temperature diffuses further for the initially extruded particles. This is seen from Fig. 11b, where the temperature distribution at the final time of \(t = 0.05\,\hbox {s}\) is depicted. Again the highest temperature coincides with the current position of the intensity peak of the laser and takes a value of \(\varTheta = 208.32\,^{\circ }\hbox {C}\). Altogether, a moving temperature profile is obtained during the AM process due to the moving laser. Consequently, the thermo-chemo-mechanical coupled behaviour is continuously induced for newly extruded material. Hence, the temperature is increased due to high-speed laser radiation and the curing is accelerated. Due to that the material stiffens faster and as a result the spreading is reduced. In the next step, the final shapes with and without laser radiation are compared to investigate the impact of the applied laser radiation during the AM. The normalised distances to the printing plate of the final shape for the simulation with laser radiation is depicted in Fig. 12a and the result without laser radiation in Fig. 12b. On the basis of the shape and the colour scale, it can be clearly seen that the distances to the printing plate for the extruded material are larger when the laser treatment is applied. Thus, the deformation in the z-direction and by that the spreading is reduced.

To quantify the effect the ratio of final to initial diameter is again used as measurement. In Fig. 13, the ratios are plotted against the associated averaged x-positions of the individual particle layers for the simulation with a laser (red) and without a laser (black). By the application of laser radiation, the ratio is reduced from 1.810 to 1.242 at the minimum x-position and the maximum ratio reduces from 1.824 to 1.315. Thus, a decrease in 56.8 % in material spreading is obtained at the minimum x-position and the maximum obtained diameter is decreased by 50.9 %.

Fig. 13
figure 13

Comparison of ratios of final to initial diameters with and without laser radiation using the fast extrusion velocity

Full size image

In conclusion, the postulated effect of a reduced material spreading during the AM process as a consequence of a high-speed IR laser radiation is captured within the performed simulation. For this very reason, a successful proof of concept for the simulation of a laser-induced high speed curing with a resulting reduction of material spreading is achieved. Consequently, the framework is well-suited for the simulation of the AM process for RTV medical grade silicones. However, one has to properly fit the developed material model as soon as sophisticated experimental data is available.

6 Conclusion

In this paper, the formulation of a large strain curing model for silicones inhibiting viscoelastic–plastic behaviour is shown. Moreover, the successful coupling of the local curing model within the non-local framework of Peridynamics is developed.

It is shown that inclusion of process-dependent plasticity within the curing model allows to represent a material spreading within simulations. In addition, a reduced spreading of a horizontally orientated cylinder by the application of higher isothermal temperatures is demonstrated. Thus, the thermo-chemo-mechanical coupled material behaviour is successfully captured.

At the end, the developed framework has been applied for the simulation of isothermal AM processes and finally for the fully thermo-chemo-mechanical coupled problem including a high-speed laser curing. Here, the ultimate goal was to represent the reduced material spreading during the AM process as a consequence of laser radiation. This phenomenological observation is captured correctly in the simulation, such that a successful proof of concept has been performed.

When it comes to high-fidelity simulations, it is an absolute requirement to exploit the full power of more sophisticated AM process simulations later on. Possible extensions are the consideration of surface tension effects, modelling the dynamics inside the nozzle and a virtual material development. In [30], the impact of surface wetting on the droplet spreading of medical grade silicones is investigated. It turns out that the interfacial tension between the surface and the silicone can be increased or reduced by the application of specific substrates. Here, the contact angle of the droplet and by that the shape of the finally cured material is directly influenced. Simulations could help to optimise the wetting influence, such that derivations in the final shape of the additive manufactured body with respect to the virtually generated shape are minimised. A separate topic on its own is the simulation of dynamic extrusion effects inside the nozzle. It is desirable to capture the dynamics including friction and adhesion effects between the uncured material and the nozzle wall since the extrusion is only possible up to a minimal diameter of the nozzle. A prediction of these effects with respect to the nozzle size and material parameters will help to understand and optimise the extrusion process further. Moreover, the processabilty of the material in a pre-cured state could be virtually tested with the goal to print the material in a less fluid like state which would reduce the material spreading.