A curing model for the numerical simulation within additive manufacturing of soft polymers using peridynamics

Within this paper, the modelling and simulation of extrusion-based Additive Manufacturing (AM) processes of curing polymers is presented. The challenge of the AM is the adjustment of processing parameters. This includes the application of laser radiation to locally accelerate the curing in order to control the final geometry of the implant. Since complex multi-physical coupling effects are hardly predictable by operator experience, numerical simulations are beneficial. When the underlying physical effects of the AM processes are captured correctly within the simulations, a realistic representation of the process is possible. To model the material behaviour during the process, a process-dependent large strain curing model is formulated, considering the stress free curing behaviour of the material. State-of-the-art models are not able to model the fluid-like behaviour of low cured polymers. This needs a formulation that takes into account finite deformations. Hence, the current model is extended to finite plasticity using a process-dependent yield function. This allows the modelling of material spreading in the fluid-like state by simultaneously reducing the accumulation of elastic stored energy, which would lead to an unintentional and non-physical bounce-off behaviour otherwise. For the numerical simulations, an enhanced version of the peridynamic correspondence formulation using fractional subfamilies with associated volume weighting factors is introduced and implemented. Besides the specific laser modelling as a volumetric heat source, a local–non-local coupling of the arising thermo-chemo-mechanical coupled equations is introduced within the peridynamic framework. Within the simulations, the applicability of the plasticity-based approach to model material spreading in the fluid-like state is presented. Finally, the software for extrusion-based printing processes is developed and the complete thermo-chemo-mechanical coupled AM process is simulated. It is shown that higher geometrical precision is obtainable in terms of a reduced material spreading by the application of a laser radiation. The model constitutes the first step of the virtual implant development regarding the optimisation possibilities during the AM process.


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  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 highspeed 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 devel-opment, 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 stressfree 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, threedimensional, weakly compressible, large strain viscoelasticplastic 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.

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 Here, ρ 0 denotes the volume in the reference configuration, u the displacements, P the first Piola-Kirchhoff stresses, b the body force, D int the dissipation, 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 is introduced. It defines the fraction of accumulated released heat at time t, H (t), divided by the total accumulated released heat, H (∞), 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.

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 is adopted from [19]. For the chemical as well as for the thermal deformation gradient, isotropic volumetric approaches are introduced where the volumetric change is governed by Within the functions J C and J Θ Θ 0 is the reference temperature, β c the shrinkage and β Θ 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. 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 fluidlike behaviour at the beginning of the process, the yield stress σ y is formulated with respect to the degree of cure. Up to the gelation point α < α 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. σ y (α ≥ α gel ) = ∞. 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

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] with 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 F p M . Therefore, an associated plasticity model using the von Mises yield criterion Therein τ M denotes the mechanical Kirchhoff stresses and τ M its deviatoric part. In the current configuration, it yields for the associated flow rule where d p M is the mechanical plastic rate of the deformation tensor. This evolution equation has to be solved ensuring the classical Karush-Kuhn-Tucker conditions (λ ≥ 0, Φ ≥ 0, λΦ = 0). In contrast to classical formulations, the yield function depends on the degree of cure and constitutes a new modelling approach of curing polymers. The last evolution equation is related to the elasticinelastic split of the viscoelastic deformation gradient within 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 The temperature-dependent part is modelled by where p μ Θ is the exponent of the temperature dependency. Furthermore, the shear modulus is modelled as a function of cure via to exhibit linear growth, cf. [8]. The parameters p μ α 1 and p μ α 2 represent the temperature-independent shear moduli for the uncured and fully cured material.

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 Herein, referring to [15], the specific thermo-chemical free energy function follows the experiment-based ansatz 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 To model weakly compressible material behaviour, the equilibrium part of the mechanical free energy is additively decomposed by 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) 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] w(J ve M ) = 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 thoughC ve M =C 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

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 bỹ The stressS M is defined by the differentiation of the free energy function with respect to the viscoelastic mechanical right Cauchy-Green tensor 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. For the free energy at hand, the mechanical second Piola-Kirchhoff stresses in the intermediate configuration are then where is a fourth-order projection tensor. The first Piola-Kirchhoff stresses required for the simulation are then obtained by Besides the stress definition, the evaluation of the Clausius-Duhem inequality leads to the entropy definition (25) and to the inequality Utilising the first law of thermodynamics the energy equation is derived where the mechanical Mandel stress is defined by M = C ve MS M and is the Laplace operator resulting by the utelisation of Fourier's law of heat conduction with a constant heat conductivity κ. Moreover, a process-dependent specific heat capacity 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 ΘC .

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 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 δ. Referring to [20], a horizon size of 3 x , where x 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 ξ = X − X are considered, which are defined as the vectors between an initial point X I and all points X within its family H X . The resulting equation of motion for a state-based 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 T[X, t], that has to be determined by a constitutive model.

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 for the force vector state. The correspondence formulation includes the scalar state w ξ , which represents an weighting influence function. By definition it has to be zero outside the family, and it depends only on the magnitude of ξ . Furthermore, it requires the inverse of the so-called shape tensor K, containing information about the particle distribution in the initial configuration. Moreover, the non-local deformation gradient F is defined based on deformed bonds ζ by 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 H = n s a=1 H a (33) and to compute the related non-local deformation gradient F a and shape tensor 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 of the related correspondence formulation. Here, the influence function w a is defined as where Within the formulation, the partial free energy of a subdomain is assumed to be directly related to its partial volume where n b is the number of particles in the respective subdomain. Consequently, the volume of the entire family is defined by V fam = n s a=1 V a sup 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.

Numerical treatment
Similar to [24,26], the meshfree discretised equation of motion is obtained for the formulation at hand by 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 δ = m x 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 for positions and velocities.

Non-local energy equation
In the local energy Eq. (27), the partial derivative of the diffusion term κ Θ 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.

Thermal diffusion
In the derived energy Eq. (27), thermal diffusion is modelled in terms of the Laplacian by c pΘ = 1 ρ 0 κ Θ, 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 h(X, t) X − X , the local thermal conduction part is replaced by the non-local thermal diffusion term For an applied bond-based approach, i.e. the heat flow scalar state within a bond ξ only depends on the associated temperatures of X and X; the following mathematical definition holds. Consequently, the non-local diffusion is formulated with respect to the pairwise heat flow density function For its definition, the temperature scalar state is introduced as Postulating the peridynamic micro-potential where κ m is the so-called thermal micro-conductivity, leads then to the definition for the pairwise heat flow density. Taking the peridynamic thermal potential within a complete family into account and comparing the result with its local counterpart leads to the relationship 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 H X I . The correction factors g I are then defined by g I = Z ∞ Z I , where Z ∞ 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 familyV = H X dV X and defining the correction factors as g I =V max 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.

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 laser . Taking into account a radial intensity function I rad (r ), the radial dependent power per area is defined, whereby the intensity function I rad (r ) has to satisfy the normalisation constraint 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 μ = r i + r a −r i 2 . The standard derivation is now defined with respect to the ±3σ -interval which is supposed to lie within the inner and outer radii. Thus, it yields 3σ = r a −r i 2 , respectively σ = r a −r i 6 . Consequently, the Gaussian radial intensity distribution is given by where the peak intensity I 0 is obtained, using the normalisation constraint (48), by The modelled normalised intensity distribution is illustrated for the parameters r i = 1 mm and r a = 3 mm in Fig. 4. On the left, the intensity is plotted in a Cartesian frame and on the right over the radius. 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 x . This leads to the definition of the radial power per volume P e laser (r ) dv = P e laser I rad (r )I pen , where the normalised intensity in the irradiation direction is defined by I pen = 1 d 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 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 x having a higher z-position than the considered particle.

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-chemomechanical 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 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 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 gradienṫ and the exponential map integrator is used for implicit time integration. In this framework, the plastic evolution equation is solved in terms of the mechanical inverse plastic deformation gradient. Furthermore, the viscous evolution equations are approximated by an Euler backwards scheme C ve,n+1 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 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 sup stresses are defined on nodal level. This is why averaged mechanical quantities, denoted by a hat, computed with the volume weighting factors

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 ρ 0 = 1100 kg m 3 is used for the simulations.

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 Pa and a Poisson's ratio of ν = 0.499. To neglect the viscous influence within the viscoelastic-plastic approach, only a single Maxwellelement is taken into account and the related shear modulus μ 1 = 0.05 Pa and viscosity η 1 = 0.035 Pa s are chosen. Additionally, a constant yield stress of σ y = 0.25 Pa is selected. Further, due to the small domain and the contact constraint a time step of t = 10 −7 s is used.
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 s almost to the initial cylinder (a), while the viscoelastic- 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 associ-ated 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 π 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 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 viscoelasticplastic approach together with the introduction of plastic deformations allows the simulation of spreading of soft materials.

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 t = 10 −7 s. The applied temperatures are Θ = 20, 50 and 100 • C. For a better comparison, the shrinkage and thermal expansion parameters β C and β Θ are set to zero. Furthermore, a degree of cure and temperature-dependent shear modulus is used and the gelation point is set to α gel = 0.1. The other applied material parameters are related to the process-dependent material functions and are depicted in Table 1.
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 ν(α = 0) = 0.4999 and to a maximal Poisson's ratio of ν(α = 1) = 0.4851. Furthermore, the parameters of evolution Eq. (6) related to curing are shown in Table 2.
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  (16), the parameters of Table 3 are applied.
In the subfigures of Fig. 7, the deformed cylinders in the final time step at t = 0.02 s are depicted for the three mentioned isothermal temperatures. With an increasing temperature from 20 (a) to 50 (b) to 100 • 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.
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 • C, where 76.93 % of the potential energy remains after a time of t = 0.006 s. For a temperature of 50 • C the minimum potential energy is obtained in a reduced time of t = 0.0056 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 • C for which a minimum potential energy of 86.41 % in t = 0.0044 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.
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 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 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.

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δ. 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 Fig. 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 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.

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.
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 v m and the extrusion velocity v 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 m s in the x-direction is applied and the extrusion velocity is varied from v e = − 0.0025 m s to v e = − 0.005 m s and to v e = − 0.01 m s . In all cases, a total time of t = 0.5 s is simulated with a time step of t = 0.5 × 10 −7 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 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 pn = 0.3375 mm, whereby the total number of initial particles is 13,280. From a peridynamic point of view, a horizon of δ = 3 x is used for all simulations and a constant influence function of w ξ = 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, β C = − 0.03, and thermal expansion, β Θ = 10 −5 1 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 κ = 0.27 W 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.
For the intermediate extrusion speed of v e = − 0.005 m 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 e = − 0.01 m s the ratio is 1.810 at the minimum x-position and takes a maximum value of 1.824. Up to x = 0.2773 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 cir- At the minimum xposition, 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.

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 laser = 1.5 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 mm and r a = 0.35 mm are chosen, such that the initial particles are not irradiated. Moreover, the fast extrusion velocity of v e = − 0.01 m 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 s and at the final time at t = 0.05 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 s, cf. Fig. 11a, a portion of the extruded material has already been irradiated and the highest temperature of Θ = 200.46 • 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 s is depicted. Again the highest temperature coincides with the current position of the intensity peak of the laser and takes a value of Θ = 208.32 • C. Altogether, a moving temperature profile is obtained during the AM process due to the moving laser. Consequently, the thermo-chemomechanical coupled behaviour is continuously induced for newly extruded material. Hence, the temperature is increased Fig. 13 Comparison of ratios of final to initial diameters with and without laser radiation using the fast extrusion velocity 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 xposition 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 %.
In conclusion, the postulated effect of a reduced material spreading during the AM process as a consequence of a highspeed 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.

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.
The mechanical second Piola-Kirchoff stress tensor is defined by a pull-back of the Cauchy stresses with the mechanical deformation gradient and lives in the thermochemical intermediate configuration. With this the second Piola-Kirchhoff stresses are given by From the split of the deformation gradient (3) follows the Green-Lagrange strain tensor with the definition of the mechanical right Cauchy-Green tensor and the mechanical Green-Lagrange tensor The rate of the Green-Lagrange strain tensor is defined bẏ Further, the chemical and thermal velocity gradients are defined by ∂ g ∂αα 1. (64) Taking the definition of the stress power P int = S :Ė and inserting (60)  and thermal stress power