A computational phase transformation model for selective laser melting processes

Selective laser melting (SLM) has gained large interest due to advanced manufacturing possibilities. However, the growing potential also necessitates reliable predictions of structures in particular regarding their long-term behaviour. The constitutive and structural response is thereby challenging to reproduce, due to the complex material behaviour. This motivates the aims of this contribution: To establish a material model that accounts for the behaviour of the different phases occurring during SLM but that still allows the use of (basic) process simulations. In particular, the present modelling framework explicitly takes into account the mass fractions of the different phases, their mass densities, and specific inelastic strain contributions. The thermomechanically fully coupled framework is implemented into the software Abaqus. The numerical examples emphasise the capabilities of the framework to predict, e.g., the residual stresses occurring in the final part. Furthermore, a postprocessing of averaged inelastic strains is presented yielding a micromechanics-based motivation for inherent strains.


Introduction
Additive manufacturing (AM) refers to processes where a part is fabricated by sequential addition of material rather than by subtractive manufacturing or moulding. Thus, the final part has little geometrical and material restrictions. Customised parts can be produced within a reasonably low cost range and comparative short amount of time. AM is often called a near-net-shaping technique, as the part is directly built based on a computer aided design (CAD) model. The CAD model is then converted into a standard triangulation language (STL) format which is widely used for all AM machines. It represents the geometry of the CAD model by a simple mesh. This file is digitally sliced into thin crosssectional layers and can then be constructed by the chosen AM process. The layering is carried out through deposition and bonding of these layers. Thus, no tools are required to initiate the AM process (except for necessary post-processing).
B Isabelle Noll isabelle.noll@tu-dortmund.de 1 Nowadays, there are various technologies, differing in the used material (mostly powder or wire), in the heat source and in the procedure in general.
In this contribution, the focus lies on the selective laser melting (SLM) process, also referred to as laser powder-bed fusion as summarised in, e.g. [14], where a work piece is manufactured in a-in this case-metal powder bed. The process uses a laser beam as the heat source to selectively melt the metallic powder particles of the layer according to the CAD model. The material then re-solidifies after a cooling period. Afterwards, the building platform is lowered and a new layer of powder is applied. This process is repeated until the final geometry of the work piece is achieved. A detailed view of the three specific material states-powder, melting pool, and re-solidified material-is shown in Fig. 1 for the SLM process. Hence, it is possible to manufacture parts by melting metallic powder in a layer-upon-layer technique. It also improves lightweight construction due to the possibility of arbitrarily shaped parts. Especially the aerospace and automotive sector, as well as the biomedical sector investigate the SLM process.
So far, the prediction of the influence of the additive manufacturing process on the properties of the final workpiece is still challenging. The temperature gradients and cooling rates present will vary significantly depending on the posi- tion of the laser beam and the considered particle. In addition, the three different states of the material exhibit significantly different mass densities so that the phase changes evoke relevant process-induced eigenstrains. Furthermore, inhomogeneously distributed thermal strains are present. Overall, these inelastic strains cause eigenstresses which may have a major impact on the properties of the final workpiece, in particular with respect to the long-term behaviour and its lifetime. In addition, these eigenstresses of the part are directly affected by the various process parameters. So far, experiments have mostly been used to identify the optimal parameter configuration on a purely empirical level. Therefore, simulations with physically well-motivated models are necessary in order to gain a deeper understanding of the relation between process parameters and the quality of the SLM part. Over the last twenty years, starting with [44], different approaches for simulations using the Finite Element (FE) method have been available in literature regarding the modelling of the SLM process on a small scale: purely thermal models, see for example [12,15,26,29,30,38,42,45] and more sophisticated thermomechanical small strain models, i.e. [11,16,22,36,37,43,49,51]. Thermomechanical models enable the prediction of residual stresses and the geometry, as well as the influence on scanning strategies. Within these contributions, there are different approaches to model the phase changes, i.e. [22] enhances the Stefan-Neumann equation and [26] uses a mathematical phase change function calibrated by experiments. In [39], a purely thermal model is used which explicitly incorporates two state variables for both the phase and the porosity. This enables the capturing of the consolidation of the material. However, most researchers use the melting point temperature to indicate the phase changes, so that these are purely temperature driven.
The effect of volume shrinkage of the powder bed due to porosity has been investigated in [12] by using a purely thermal model. In [16,37] a thermomechanical model is used, where a pseudo-thermal expansion is implemented to take into account not only the thermal expansion of the material, but also the shrinkage of the powder bed. The advantage of a thermomechanical model and the difference in the temperature evolution is also stated in [16]. In [48], a thermomechanical finite deformation framework is developed to predict the melting and solidification process on the powder scale.
For industrial usage, the influence of process parameters, product quality and the prediction of deformation and residual stresses still remains challenging for a physically real component. Researchers are presently focusing on advanced models for the numerical investigation of rather large parts in contrary to a single-line melt tracks. On the one hand, high performing FE-algorithms are developed to speed up computation time, for example in [27,33,34]. On the other hand, the models are simplified at the macro scale: Multiscale approaches where the micro, meso, and macro scale are taken into account, as for example described in [9,24,28]. In particular, simulating each layer at the same time is used to speed up the process, see for example [17]. A similar approach is the inherent strain method which was first introduced for metal welding applications in [21,50]. This method is now applied to additive manufacturing processes, see for instance [24,31,40,41]. The challenge in using this method lies in the determination of the inherent strains, either through simulation or experiments. Altogether, there are different formulations available for the simulation of the micro and macro scale with help of the FE method. However, from our point of view, no frameworks exist which take into account the mechanical material model and phase change in a physically well-motivated approach.
Based on the framework for solid-solid phase transformations in shape memory alloys, see e.g. [3], this contribution focuses on the advancement of the model presented in [2]. The model is based on the explicit consideration of phase energy densities, where each state of the material-powder, molten and re-solidified-is a single phase (see Fig. 1). During the transformation process, the material first melts and then solidifies. These phase energy densities are combined by using a specific homogenisation assumption which yields the averaged energy density-and with this the constitutive model-for the possible phase mixture. Several inelastic strain contributions are taken into account within the material models of the single phases. More precisely speaking, viscous strains and transformation strains are considered in the molten phase as well as thermal strains, plastic strains, and transformation strains in the solid phase. The trans-formation strains are introduced as material constants and capture the significant change of the mass densities during the phase changes. The internal state variables including the mass fractions of the different phases are calculated by using thermomechanically consistent evolution equations solved at the material point level. The formulation of the material model as well as the FE model consider full thermomechanical coupling. The present investigations aim at the accurate prediction of essential process-induced quantities such as eigenstresses. Furthermore, the present modelling framework could be used to calculate effective inherent strains which could then be used in simpler yet faster process simulations.
This article is set up as follows: A brief summary of the thermomechanical field equations is given Sect. 2, where a focus is set on the thermomechanical coupling and mechanical dissipation contributions. In Sect. 3, the phase energy densities and the irreversible strains are specified to introduce the constitutive material model which is based on mass fractions. We provide deeper insight into the local algorithm and Abaqus implementation of the framework at hand in Sect. 4. Calculations at the material point level designed to yield a proof of concept as well as representative three-dimensional boundary value problems are studied in Sect. 5, before we finally conclude our findings in Sect. 6.

Thermomechanical framework
This section provides general information on the kinematics and the thermodynamically fully coupled balance equations. At this stage, the small strain theory is considered to be appropriate. Thus, the linearised strain measure with the displacement field denoted as u is applied. In general, the total strains can be decomposed additively as follows into an elastic part ε el and an irreversible part ε irr . The irreversible strains describe inelastic effects e.g. viscous, thermal or plastic strains and the combination thereof, and transformation strains, i.e. eigenstrains, capturing the phase transformation of the material, as explained in subsequent sections.
To simulate the underlying SLM-process, a fully coupled thermomechanical model is proposed. The problem at hand is based on the balance of linear momentum and on the energy equation which follows from the first law of thermodynamics, i.e.
respectively. In these equations, σ denotes the stress tensor, b describes the body forces, q represents the heat flux vector, r ext denotes the external heat supply, ρ is the mass density and e symbolises the specific internal energy. The absolute temperature is denoted by θ . In general, notation• is used to describe the time derivative of the respective quantity •. As this work proceeds, all body forces shall be neglected, i.e. b ≡ 0.
In the following, the second law of thermodynamics, respectively the entropy inequality, is used to introduce the internal dissipation. As proposed by [10], the specific Helmholtz free energy Ψ shall be used for a thermodynamic consistent derivation. In order to reformulate the problem depending on the temperature θ as state variable rather than the entropy s, the specific Helmholtz free energy Ψ is obtained via the Legendre transformation e = Ψ + s θ, In general, the Helmholtz free energy Ψ (ε, θ, V) depends on the total strains ε, the temperature θ and-at this point not further specified-on internal variables V. With this at hand, the time derivative of the Helmholtz free energy results iṅ Following the Colemann-Noll procedure, cf. [10], the entropy inequality is used and, combining (8) with the energy equation (4), the second law of thermodynamics reads This can be reformulated in terms of the dissipation symbolised with D and split into a mechanical and thermal part, i.e.
Moreover, each contribution is assumed to be non negative. This results in the Fourier inequality for the thermal contri-bution on the one hand, and in the Clausius-Planck inequality on the other, By recalling (4) and using (6), it is possible to reformulate (12) such that Inserting (7) into (13) and rearranging terms for the purpose of simplification results in where • represents a generalised scalar product. Following the lines of the Coleman-Noll procedure, the constitutive equations based on the Helmholtz free energy density can be directly extracted, i.e.
From this, the definition of entropy and stresses can be determined in a straightforward manner. As a further abbreviation, the driving forces of the internal variables are introduced. Substituting these results into (14) finally leads to the reduced form of the mechanical dissipation as With these derivations at hand, (4) can be rewritten in the more convenient form where the specific heat capacity is defined by In the following sections, it is assumed that the (unspecific) Helmholtz free energy density is defined as ψ = ρ Ψ .

Constitutive framework
The thermomechanical material model based on phase energy densities is introduced in this section. The overall material behaviour is determined according to a phase transformation model. Furthermore, the specifications of the thermal problem are introduced. Altogether, both models are thermomechanically fully coupled.

Phase energy densities
As already indicated in Sect. 1, each state of the material is represented by an energy density. Altogether, the material response of the complete model is determined with the help of a homogenisation approach, rather than by relying on solely temperature-dependent material parameters. At this point, three distinct phases representing the different states of the material, namely powder, molten, and re-solidified, are explicitly taken into account. All of these phases are modelled as a solid continuum, which may be arguable in particular for the molten phase with respect to the (viscous) solid approach and for the powder phase regarding the continuum ansatz. Motivations for these simplifications are, however, provided in the subsequent sections.
The material model is based on the averaged Helmholtz free energy density ψ, where additively decomposed energy densities ψ i each of which consists of two parts, are chosen for each phase denoted by the index i. The first part describes the mechanical energy density ψ mech i , while the second part ψ cal i is a purely caloric contribution due to the temperature dependence of the material model. This results in an energy density where the index i refers to all possible phases. The mechanical part of the energy density is defined as whereas the caloric part is a function of temperature and specific material parameters and which is to be defined later. While ε i denote the total strains of each phase, ε irr i describe the irreversible strains and E i represents an fourth-order elasticity tensor of the respective phase. The influence of the caloric part and the irreversible strains are visualised in Fig. 2. Each energy density represents a potential well which is shifted due to ψ cal i and ε irr i . In the following paragraphs, the corresponding phase energy densities ψ i will be discussed in detail. Powder particles are assumed to move freely due to the high degree of porosity between the powder grains. Therefore, the material behaves solely elastic. It has to be kept in mind that only compression strains are admissible within the powder particles. In the framework at hand, the powder is modelled as a continuum with a significantly higher compliance than the solid material. Overall, the powder shall not undergo states of high stress or strain level. In the initial state, only the powder phase is present. In addition, due to the high porosity of the powder phase, any thermal strains are neglected for this state of the material. Thus, the powder is considered as the parent phase so that no irreversible strains are incorporated for this phase, i.e. ε irr pow = 0. The energy density of the powder is approximated by where c pow = ρ pow c is the weighted heat capacity of the powder, and where L pow = ρ pow L describes the weighted latent heat of the powder at a constant reference temperature θ ref pow . The weighted values are calculated by the mass density of the powder ρ pow and the specific value of the heat capacity c and latent heat L, respectively.

Molten phase
The molten phase shall also be approximated as a solid type phase. This behaviour is assumed to be appropriate, as the molten pool is only present for a brief time span. Fluid effects like the Marangoni flow are therefore neglected. However, to include a fluid-like behaviour of the molten pool in the material model, a viscous strain contribution ε vis mel is included, which shall be used to enable full stress relaxation within the molten phase. Furthermore, the mass density of the three phases varies considerably. To take into account the shrinkage from the powder material to the molten phase, a transformation strain ε trans mel is incorporated. Both irreversible strains will be further defined in what follows. Overall, an additive split of the irreversible strains is used for the molten phase, i.e. ε irr mel = ε trans mel + ε vis mel .
As the molten phase is the high temperature phase, the influence of the latent heat is only integrated into the powder and re-solidified phase. With this at hand, the specific energy density of the molten phase is defined as with the heat capacity defined as c mel = ρ mel c and the mass density of the molten phase denoted as ρ mel .

Re-solidified phase
An additive decomposition of the irreversible strains into three parts is chosen according to for the re-solidified phase. The transformation strains ε trans sol take into account the volume changes of the material during the phase transition from the molten phase to the re-solidified phase due to the changing mass densities. In contrast to the previously introduced phases, two further inelastic strain contributions are considered. Plastic deformation can arise due to high eigenstrains after the phase transformation and due to high temperature gradients present. A plastic strain tensor ε pl sol is used to capture this behaviour whose evolution will be further defined in Sect. 3.4. Furthermore, thermal strains ε th sol are considered in the re-solidified phase to take into account inelastic strains stemming from the high temperature changes during the process. These strains will be further defined in Sect. 3.5. Finally, the energy density of the re-solidified phase is defined by analogy with the previously introduced energies as The heat capacity and the latent heat are defined in an analogous manner as c sol = ρ sol c and L sol = ρ sol L. The mass density of the solid phase is denoted by ρ sol . In addition, H sol defines the hardening modulus of the solid phase, whereas k hard sol indicates the accumulated equivalent plastic strain related to isotropic hardening.

Specification of transformation strains
The transformation strains ε trans i of the molten and resolidified phases shall be introduced next, as they are solely dependent on the transformation process and on material parameters previously introduced. These strains shall capture the volume changes of the material during phase transitions due to the different mass densities. In fact, the validity of the small strain theory needs to be carefully scrutinised in view of physically plausible transformation strains. Initially, the material solely consists of powder with the mass contribution dm 0 , thus Hypothetically assuming that this infinitesimal volume element of the parent powder phase fully transforms into one of the other phases, i.e. molten and re-solidified material, the volume of the new phase is given by Here it is assumed that the mass contribution dm 0 is constant (conservation of mass). Using (28) and (29) together with the physically sound assumption of ε trans • being purely volumetric, these transformation strains are specifically given by where I corresponds to the second order identity tensor. This results in the transformation strains for the molten phase and for the re-solidified phase.

Homogenisation via convexification
In the context of additive manufacturing and the related phase changes of the underlying material, the mass densities of the powder, molten phase and re-solidified material differ significantly. Thus, conservation of mass does not coincide with conservation of volume. The volume consequently cannot be used as a conserved quantity and the homogenisation algorithm needs to be reformulated accordingly.
Material models for solid-solid phase transformation, see, e.g. [3,4,23,25], mainly use the respective volume fractions as state variables and as conserved quantity. This is motivated by the fact that the mass densities of, e.g., austenite and martensite can be considered identical. Therefore, the algorithm has in general been introduced with respect to volume fractions where dV 0 denotes the (infinitesimal) initial volume and where dV i is the corresponding volume of the respective phase.
First of all, the averaged energy density ψ is introduced via a volume averaging as subject to 0 ≤ ξ i ≤ 1 and i ξ i = 1 as constraints for the volume fraction. The number of phases is given by n ph . The volume-also the mass for equal mass densities-is hereby used as conserved quantity. Furthermore, an additional constraint in terms of the compatibility condition has to be considered. For the application at hand, the algorithm needs to be enhanced and re-formulated with respect to the mass fractions These mass fractions can be related to the volume fractions via For the sake of consistency, the algorithm shall now be developed based on mass fractions and on the averaged mass specific energy Ψ . While working with mass fractions ζ i , some adaptions to the aforementioned approach have to be made. The overall energy Ψ is analogously calculated via a linear mixture rule of the mass specific phases Ψ , thus where the relation between the volume and mass specific energy density ψ i = ρ i Ψ i still holds. The use of mass fractions in particular affects the compatibility condition as well as as constraints for the mass fractions. One can see the similarities between (38) and (34). This is advantageous for the implementation, as will be explained in Sect. 4.2.
For the specific material model at hand, three phases are present so that n ph = 3, namely powder, molten pool and solid. The averaged or, in other words, effective energy density explicitly used in the present modelling framework is given by This averaged energy density has to be minimised subject to the constraints regarding the domain of feasible mass frac- and the domain of the admissible strain distributions ε • ∈ E with This constrained minimisation problem results in the so called convexification of Ψ , i.e.
Thus, CΨ is also known as the convex hull, which is qualitatively visualised in Fig. 3.
In contrast to standard solid-solid phase transformation, an additional constraint is necessary regarding the evolution has to be fulfilled in addition to (42). At this point, only the determination of the optimal strains shall be discussed. The evolution equations of the mass fractions will be considered in detail in Sect. 3.4. The strain states in each phase shall minimise the averaged specific energy density while considering (43). The aforementioned minimisation problem subject to the respective equality constraint can be solved with help of a Lagrangian L = Ψ + λ : r ε (46) with the Lagrange multipliers contained in λ. Accordingly, the necessary conditions for the minimisation are determined as It is possible to derive an explicit expression for the optimal strain distributions by using (47) and by solving these equations with respect to the Lagrange multipliers After taking into account (43) and (48), the following analytical results are obtained for the respective optimal strains of each phase with the abbreviation Finally, the overall stress can be determined in a compact form by using the abbreviation E = E sol : E mel : E sol . The consolidated notation for (49)- (51) and (54) can only be derived in case E • is an isotropic tensor. This assumption seems to be appropriate for the material at hand, as defined in Sect. 5.

Evolution equations
In line with, e.g., [7], it is possible to derive evolution equations for state variables V from a dissipation function C, i.e.
with the consistency parameter Γ and the generalised inequality constraint r ≤ 0. For specific cases, it may be more convenient to use the dual dissipation function C instead. This quantity depends on the driving force and can be computed by applying the Legendre-Fenchel transformation As an alternative to (55), the evolution equation is then given bẏ This concept is now applied in order to derive the evolution equations for the model-specific independent internal state variables The mass fraction of the powder is substituted by according to conservation of mass. As the volume and mass fractions can be converted into one another, compare (37), the following framework is derived based on volume fractions ξ • . This enables the enforcement of the minimisation based on the Helmholtz free energy ψ without weighting the equation itself with the respective mass density ρ • . The stresses can then be calculated in a straightforward manner via (54).

Volume fractions
In view of establishing evolution equations for the volume fractions ξ mel and ξ sol , the dissipation function is introduced together with the inequality constraints derived from (42). The evolution of phase fractions is supposed to occur over a finite temperature range rather than jumping from zero to one at a specific temperature. Therefore, this development is considered dissipative. Altogether, two Biot-equations are carried out which govern the evolution of the independent volume fractions. In this case, refers to the driving force based on the free energy density Ψ .

Viscous strains
The evolution of viscous strains ε vis mel within the molten pool is determined via a visco-elastic and thus rate dependent primal dissipation function Therein, η vis mel denotes a viscosity related material constant. For this evolution equation, no inequality constraints have to be taken into account, therefore the right hand side in (55) equals zero when applied to the evolution of viscous strain contributions. In general, the driving force is calculated via (56), such that the standard driving force for the viscous strains is defined as Here, the driving force F vis equals the stress tensor weighted with ζ mel . However, this driving force is a quantity averaged over a domain due to the homogenisation framework applied, whereas the domain can contain multiple phases and changing mass fractions, as visualised in Fig. 4. This domain corresponds to an infinitesimal surrounding of a material point. As a consequence, a local driving force has to be defined that affects only the molten phase within the domain considered, see Fig. 4 and also Remark 1. This driving force is incorporated into the evolution equation, such thaṫ The viscous strains incorporate a fluid-like response, as a complete relaxation of the stresses is enabled.

Plastic strains
For the rate-independent evolution of the plastic strains ε pl sol as well as of the variable k hard sol related to the isotropic hardening, the following dual dissipation function is used. In the aforementioned equation, λ denotes the Lagrange multiplier and Φ refers to the yield function. Following the same argumentation as in the previous section, the yield function depends on the local driving forces which yields the associated evolution equation for the plastic strainṡ and for the hardening variablė The yield function is specifically chosen as where σ y sol refers to the yield stress of the solid phase and where • dev := • − 1 3 tr(•) I denotes the deviator of the quantity •. With this, the flow direction can then be specified as together with ∂ /∂κ loc = √ 2/3.

Remark 1
In general, a driving force F = −∂ V ψ is calculated for all internal variables V. However, the evolution of the internal state variablesV of the corresponding material only depends on the local driving force F loc , as discussed before. In contrast, the mechanical dissipation entry D mech = FV regards the complete driving force, as the thermodynamic consistent driving force F is already weighted within D mech due to the homogenisation approach applied for Ψ .

Heat effects
All standard heat transfer mechanisms, particularly heat expansion, conduction, radiation and convection, are considered within the model formulation and shall be briefly discussed in what follows. With respect to the heat expansion model and related thermal strains ε th sol , a standard linear relation In addition, a standard isotropic Fourier-type form is used for the heat conduction model, viz.
In this relation, represents the averaged heat conduction coefficient of the phase mixture according to a Reuss-type homogenisation with k • as the respective heat conductivity of each phase. Moreover, heat convection and radiation are regarded on the outer boundaries. In a standard manner, heat convection is considered by where h refers to the convective heat transfer coefficient, i.e. the reference film coefficient. Furthermore, θ sur and θ sink denote the surface and ambient (sink) temperature, respectively. In addition, radiation heat flow can be given by with σ boltz = 1.38064852 × 10 −23 m 2 kg s −2 K −1 corresponding to the Stefan-Boltzmann constant, emiss being the emissivity coefficient, θ Z defining the value of absolute zero on the specific temperature scale used and with θ amb referring to the ambient temperature of the surrounding media. For further insight into the heat transfer mechanism the interested reader is referred to, e.g., [6]. As introduced in (19) and (20), the balance equations depend on the effective specific heat capacity which is here given by It is worth noting that the specific heat capacity is dependent on the volume fraction ξ • rather than on the mass fraction ζ • . This is due to the direct dependence of the effective specific heat capacity with the averaged specific energy ψ. Besides, quantity χ stems from the thermoelastic coupling via the thermal strains in the solid phase and is defined as However, in the present framework, this contribution turns out to be negligible due to the dependence on the square of the heat expansion coefficient.

Implementation and algorithmic treatment
The framework at hand has been realised with help of the commercial FE-based software Abaqus. For the FE implementation, not only the time discretisation and the evolution of the volume fractions are necessary, but various user subroutines also need to be specified to implement the present modelling framework. For example, the constitutive model defined in Sect. 3 has to be incorporated into Abaqus via the subroutine UMAT, whereas the thermal material model is adapted within UMATHT. In general, the strong form of the energy balance is defined as in Abaqus, see [13], where r * symbolises all effective heat sources. Following the derivations in Sect. 2, some differences are visible. Overall, the heat sources in a thermomechanical problem can be additively decomposed into the sum of external heat supply r ext and the volumetric heat generation caused by mechanical working of the material This results in All terms contributing to r mech are calculated within the user subroutine UMAT as contributions to the variable RPL as defined in [13]. In a thermomechanical analysis, this contribution-if specified-is automatically added to the external heat sources. In doing so, it is possible to directly include the influence of mechanical working and dissipation in Abaqus. To complete the model, the moving heat flux representing the laser beam is added in the user subroutine DFLUX. In addition, an internal strategy of Abaqus (*Model Change) is used to model the layer construction during the SLM process. For further information, see [2] and the references cited therein.

Global finite-element-based algorithm
In this section, the discretised weak forms of (3) and (19) as introduced in Sect. 2 are defined, where the contributions from both mechanical equilibrium in the absence of mechanical body forces and heat equation for a thermomechanically fully coupled system are considered. The completely discretised weak form can be given by where A denotes the assembly operator, N u A , and N θ C symbolise the corresponding set of shape functions, ρ refers to the current density, see also Remark 2, surface tractions per unit area are represented by t = σ · n and the heat flux vector is given by q n = −q · n, with n describing the outward normal unit vector. The derivation of the discretised weak form is given in the Appendix A.

Remark 2
In analogy to (41), the current density is defined as ρ = ζ pow ρ pow + ζ mel ρ mel + ζ sol ρ sol . At the beginning of the simulation, the material completely consists of powder material, thus ρ := ρ pow . The effective density obviously changes during the process due to the evolution of the mass fractions. This can be implemented within the user subroutine UMATHT, where the density is artificially adapted, cf. the explanations in [2].

Local algorithm
In general, all time derivatives• are approximated by a forward difference quotienṫ wherein n • refers to the previous and n+1 • to the current time step of the quantity •. The time increment is defined as Δt = n+1 t − n t. This time discretisation is applied not only for the temperature and strain field rate,θ andε, respectively, but also for all rates of internal variablesV. The solution of the time-discretised versions of the evolution equations given by (66), (67), (72) and (78) with respect to n+1 V is discussed in the following. Evolution equation (72) can be explicitly solved for ε vis mel by using the backward Euler scheme, such that Consecutively, the viscous strains are replaced with this analytical expression throughout the algorithm. The constraint (45) can be further reformulated by inserting relation (60) and by applying the forward difference quotient so thaṫ With these conversion at hand, the distinct set of inequality constraints introduced in (62)-(65) is now defined by which has to be fulfilled by the algorithm. The first two residuals emanating from the evolution equations with respect to the volume fractions of the molten and solid phase are specifically given by where the dissipation function specified in (61) is analogously discretised by using (94), and where the explicit definition for ξ pow is not inserted because the terms are already rather lengthy. The compatibility conditions belonging to the constraints (97)-(100) shall be transformed into equality constraints by using the regularised Fischer-Burmeister nonlinear complementarity function, as applied in [2] and explained in the references cited therein. For a more general overview and further insight, cf. [5]. Thus, the residual entries can be implemented into the framework, where the perturbation parameter δ = 10 −10 is used. With this parameter, a continuously differentiable function exists. Furthermore, (77) and (78) are discretised in time by using a backward Euler scheme, resulting in To calculate the consistency parameter λ in (107) and (108), the complementarity conditions related to the evolution of plastic strains are incorporated via The final residual has to be solved with respect to where the tensorial quantities r 7 and ε pl sol are rewritten as vectors with the help of the Kelvin notation, resulting in vectorial entries r 7 and ε pl sol .
Using the Fischer-Burmeister nonlinear complementarity function has the advantage that a standard solver, such as the Newton-Raphson scheme, is applicable even though inequality conditions have to be fulfilled. This method is then used to solve (110) for (111), where the norm of the residuum is checked in every iteration with a tolerance = 10 −8 .

Numerical examples
In this section, different numerical examples are discussed. The material model is adapted to a Ti 6 Al 4 V titanium aluminium alloy with material parameters according to Table 1. The transformation strains are directly defined via the mass densities. The parameter η ξ incorporated in the dissipation potential D ξ is set to η ξ = 0.005 for the examples at hand. This parameter can be used to adjust the temperature range, in which the phase changes occur in comparison to experimental findings. Moreover, it can be used to potentially stabilise the global FE scheme in terms of numerical robustness. In addition, the viscous strain parameter is chosen as η vis mel = 70, so that a high relaxation within the molten (fluid) phase is feasible. The mechanical parameters, i.e. Poisson's ratio ν, Young's modulus E and yield limit σ y , are taken from the literature, see [32,36,46]. The hardening modulus H sol is, at this stage, chosen without particular literature reference. The thermal material parameters, i.e. the expansion coefficient α, the specific heat capacity c, the conductivity k and the latent heat L, as well as the initial temperature θ ini and the reference temperature θ ref are parameters based on [32,36]. The respective effective counterparts which are used to calculate the material response follow from the homogenisation approach introduced in the previous chapters in contrast to material models which directly incorporate temperature-dependent averaged material properties. Furthermore, it shall be remarked that the difference in the mass densities of the material's phases is rather large and, in consequence, the volume changes captured by the transformation strains. Therefore, one could argue that the used small strain approach may not be appropriate. However, we consider the small strain formulation acceptable, as little rotations are present (see Remark 3) and the rather large transformation strains are completely volumetric.

Remark 3
The deformation gradient F can be multiplicatively split into a rotation tensor R and a stretch tensor U, so that F = R · U. The general definition of the small strain tensor has been introduced in (1). In general, the strain measure is introduced via the deformation gradient, so that, e.g., ε :=

Proof of concept
In order to investigate the pure material response of the proposed framework, a temperature cycle is prescribed at material point level with zero Neumann boundary conditions as a proof of concept. At the beginning, the temperature is linearly increased in a time interval Δt = 1 ms from θ = 873 K to θ = 2873 K, then decreased back to θ = 873 K in the same time span, see Fig. 5a. The resulting strain evolution is shown in Fig. 5b, while Fig. 6 illustrates the corresponding evolution of mass and volume fractions. The significant changes in the strains occur simultaneously to the evolution of phase fractions. For constant phase fractions, the strains are therefore also constant, except for the re-solidified phase where the strains decrease while the material cools down due to heat expansion. Due to the zero Neumann boundary conditions, no stresses arise, so that no viscous or plastic strains are present in this conceptual proof. In Fig. 6, particularly the difference in the evolution of the mass and volume fractions for the three phases can be examined. The material initially purely consists of powder, thus ξ pow = ζ pow = 1. With increasing temperature, the first phase change towards the molten phase starts at approximately θ ≈ 1880 K and finishes at θ ≈ 2330 K. The solidification process begins at θ ≈ 1870 K and ends at θ ≈ 1520 K. As already indicated above, the parameter η ξ governs the time span in which the phase transition occurs and thereby significantly affects the temperature interval as well. In other words, it is possible to affect the mushy region between the solidus and liquidus temperature with this parameter, as it controls the rate dependent behaviour of the evolution equation (61). From this numerical example, it can be concluded that a volume change of approximately 36 % occurs during the process.

Simulation of basic SLM processes
In the following, the modelling and simulation of basic SLM processes are presented. The main geometrical model with all boundary conditions on ∂B is conceptually illustrated in Fig. 7 for the representative building chamber. For the example at hand, homogeneous Dirichlet boundary conditions are prescribed for the displacement, so that The temperature at the bottom surface is prescribed based on Additional heat flux terms due to convection q conv and radiation q rad are present at the top surface, i.e.
In addition, the laser beam heat input which is represented by the volume-distributed heat source r ext exists along the respective scanning path. Various models exist for the application of the heat flux r ext generated by the laser beam as accurately as possible, see e.g. [20,22,26,47]. For this model, a volumetric heat source as proposed in [19] for welding is adjusted for SLM, such that the external heat source is defined within the Abaqus subroutine DFLUX by where P defines the laser power, r 0 is the characteristic (focus) radius of the laser beam and x i refers to the moving coordinate system whose origin lies in the centre of the laser beam at the top surface of the current powder layer, see Fig. 7. Coordinates x 1 and x 2 are needed to define the laser movement in dependence of the laser velocity v lsr within the titanium alloy is set to emiss = 0.19, cf. [1]. A layer height h lyr = 50 µm is used, whereas the height of the base material is chosen to be larger (150 µm). Overall, a part made out of three layers n lyr = 3 will be simulated. To model the application of different layers, the element birth and death method implemented in Abaqus is used (*Model Change). More insight into the laser beam model and layer construction model can be found in [2]. Different models for the laser beam impact as well as models for the heat transfer mechanisms could be incorporated into the algorithm in a straight-forward manner.
The initial part only consists of powder material where the first layer of material is already activated. At first, the building chamber is homogeneously pre-heated to θ = 373 K subject to the respective boundary conditions. Afterwards, the laser beam is applied along the predefined scanning path and a cooling time is given before the next layer of powder is activated. This procedure is repeated until all three layers have been activated and until the laser beam has been applied. Finally, the work piece cools down to the building chamber temperature θ . For the simulation, the element type C3D8HT is chosen. Due to high temperature gradients and laser veloc- ity, a dense mesh is used close to the laser beam path. Within this region, a characteristic element length of l char ≈ 20 µm is applied, as visualised for the straight laser path in Fig. 8. For the L-shaped laser path, a constant mesh of l char ≈ 20 µm is used, as shown in Fig. 9. For both examples, three elements per layer are used in thickness direction. During the simulation, the automatic time incrementation included within Abaqus is used. The simulations are performed on the Linux HPC cluster (LiDO3) at TU Dortmund University, where one compute node using eight cores is taken for both simulations.

Straight laser path
For this example, the laser beam moves 1 mm along the x 1 -direction. The geometry considered, respectively half of the geometry, is specified by l = 3.3 mm, w = 1 mm and d = 0.3 mm for the boundary conditions introduced in (112)-(117). In Fig. 10, the temperature evolution θ(x, t) and the evolution of the mass fraction of the molten phase ζ mel is illustrated. The deformed mesh (with scale factor one) is plotted for all consecutive figures. The depth of the molten pool increases as can be seen in Fig. 10b, c compared to Fig. 10a due to the higher conductivity of the already resolidified layers. This allows the connection of newly added layers to previous ones so that there is one compound. The In addition, the corresponding mesh is visualised. (Color figure online) volume shrinkage of the material during the melting is indicated in Fig. 10. After the application of all three layers and the consecutive cooling, the final re-solidified part can be identified in terms of the volume fraction ξ sol , see Fig. 11a, b, where in the latter figure the solid part has been virtually extracted. These two Figures show the re-solidified material after applying the laser beam to three layers in contrast to Fig. 10c, where only the current state of the molten material during the third scanning is visualised. The maximum value of the volume fraction of the solid phase equals 0.633 (i.e. max ξ sol = ρ pow /ρ sol = 0.633) in contrast to the related mass fraction (max ζ sol = 1), compare Fig. 6c.
The incorporation of r mech in the subroutine UMAT as defined in (91) is quite important. Without the coupling term RPL, the maximum temperature within the molten pool would be overestimated by approximate 200 K. Within r mech , not only the dissipative effects due to the internal state variables V are included, but also the latent heat effects are thermomechanically consistent incorporated into the framework. For example, the latent heat terms in ψ pow and ψ sol change the heat conduction and cooling speed. In the literature, mostly the apparent heat capacity method as introduced by [8] is used to capture the effects of the latent heat in an artificial increase of the heat capacity, e.g. [11,15,18,42,45,48]. An in-depth study on different methods to include latent heat effects during phase changes can be found in [35]. In [42], the influence of the latent heat for the solely thermal problem has been specifically examined, where similar results were obtained compared to our studies. Overall, (91) mainly regulates the temperature evolution and thus the melt pool geometry rather than the absolute stress values. In general, the laser beam parameters influence the melt pool geometry which affects the hatching strategy and the layer height while simulating a complete part.
The residual stresses stemming from heat expansion and, more significantly, the volume changes due to the phase transitions, are illustrated in Fig. 12. These stresses are partic-ularly significant within the re-solidified part, as the powder material exhibits a small Young's modulus (Fig. 12d). In the contour plots, it is observable in Fig. 12a that particularly high tensile normal stresses are present in x 1 -direction of the third layer, in other words along the direction of the moving heat source. In addition, high compressive normal stresses exist in x 2 -direction which is perpendicular to the direction of the moving heat source, see Fig. 12b. Negligible stresses are present in x 3 -direction as presented in Fig. 12c. Not pictured are the shear stresses, as these stresses are negligible (for the coordinate system considered) compared to the highlighted normal stress contributions. To gain a better understanding of the stress evolution, a centred element in the first layer shall be considered in detail indicated by a red dot, see Fig. 11. In Fig. 13a, the related temperature and the molten mass fraction are illustrated. A steep temperature increase is found when the laser beam is applied. During the second scanning, the top material remelts, whereas during the scanning of the third layer, the first layer does not completely re-melt again. When examining the stress evolution in Fig. 13c, the relaxation of the stress σ 11 is noticeable during the scanning of the second layer. In contrast, stress σ 22 increases, as the expansion and contraction of the material in x 2 -direction is hindered.
The temporal evolution of various significant quantities is exemplarily presented for one element of the first layer in Fig. 13. The values of the nodes are averaged for the element. In Fig. 13c, the stress evolution is shown. High total strains are present for ε 22 and ε 33 in Fig. 13b. After cooling down for the third time, see Fig. 13a, a steady-state is reached for the strains and stresses, as visualised in Fig. 13b, c.
As mentioned before, various models use a phenomenological approach to define a macroscopic inelastic strain contribution that accumulates all irreversible processes. This is often denoted as the inherent strain method, see e.g. [24,31,40,41]. In a post processing step, the accumulated irreversible strains ε irr can be evaluated with the introduced model. These are then defined as The evolution of this averaged inherent strain of one element is exemplarily plotted in Fig. 13d. Due to the volumetric transformation strains, the evolution of inherent strains is (quasi) volumetric for the particular boundary value problem considered until the part solidifies. Then the accumulated inelastic strain evolution differs: larger strains are only found in x 2 -and x 3 -direction, while these strains almost coincide again during the second re-melting. After cooling, almost no further changes occur. Only one small peak is visible when the third layer is applied and molten. This is in accordance with Fig. 13a, as the element does not completely remelt during the third cycle. In the final state, high irreversible strains for the x 2 -and x 3 -direction are present.

L-shaped laser path
In addition to the aforementioned case, a more complicated laser beam path is simulated. For this example, the laser beam moves 0.75 mm along the x 1 -direction and 0.75 mm along the x 2 -direction. Similar to the example before, three layers are consecutively added. The boundary conditions defined in (112)-(117) are referred to a geometry with l = 1.6 mm, w = 1.6 mm and d = 0.3 mm, cf. Figs. 7 and 9. In theory, it is possible to model even more complex geometries. However, the computational time increases considerably with a larger number of elements and longer laser beam paths. When comparing both simulations, the straight laser path needs approximately 17 hours, where 27 693 elements are used. Overall, the model exhibits 151 821 degrees of freedom. The example of the L-shaped laser path, however, is applied using 108 800 elements and a total of 581 192 degrees of freedom. Altogether, the computational time increases to almost First of all, the normal and shear strain distribution is pictured in Fig. 14, where the deformed mesh (with scale factor one) is presented. Here, a couple of effects are striking: the dependence of the strains and the laser path orientation is obvious, compare Fig. 14a for ε 11 and ε 22 . While the beam moves along the x 1 -direction, small tensile strains ε 11 are present. The values of the strains change as soon as the laser beam turns and moves along the x 2 -direction. The strains ε 11 switch from tensile to compressive strains and small tensile strains ε 22 are present. In addition, especially high normal compressive strains ε 33 are found within the part, as it shrinks the most in depth. Due to the changing laser beam path, shear stresses are induced, as illustrated in Fig. 14b. Large negative shear strains ε 23 change to high positive shear strains ε 13 after the kink, as a symmetry is present with respect to a plane parallel to −x 1 = x 2 . Altogether, the strains are rather constant along the width of the re-solidified part and within one direction of each layer. In contrast-with an increasing number of layers-the magnitude of the strains increases.
Although the laser beam path, respectively the objective geometry, is symmetric with respect to the underlying symmetry plane and although the distribution of the strains is almost symmetric as well, the distribution of the resulting solid part is not completely symmetric. The widths of the two sides differ slightly. This can be explained by the rather coarse discretisation. In addition, the depth of the simulated workpiece increases at the turning point of the laser beam and the part is notably wider at the kink. Due to the turning of the laser, the heat influence takes place over a longer time period compared to the positions further afield and the straight lines considered in Sect. 5.2.1. This effect has to be borne in mind when manufacturing more complex components. Figure 15 shows the normal stresses which are significantly larger than the related shear stresses. In analogy to Fig. 14a, the stresses vary according to the laser beam movement. This becomes especially visible when comparing Fig. 15a, b. High tensile stresses σ 11 and σ 22 , respectively, can be found along the movement of the laser beam. High compressive strains exist perpendicular to the laser movement, especially in the lower portion of each layer. In x 3 -direction, less normal stresses σ 33 are computed.
Finally, the equivalent irreversible strain ε irr (in analogy to the equivalent von Mises stress) shall be analysed, as illustrated in Fig. 16. Overall, the results show a symmetric distribution of ε irr with respect to the loading path. This coincides with the inherent strain method, where the orientation of the laser beam path is of immense importance, cf. [24]. In summary, the values of strains and stresses themselves significantly change with respect to the loading path, but correlate to one another when taking into account the movement of the laser beam.

Conclusion
In this contribution we presented a three dimensional and thermomechanically consistent fully coupled framework for the simulations of Selective Laser Melting processes based on a phase transformation model. The different states of the material-powder, molten and re-solidified-are explicitly taken into account with the respective mass densities, where the mass fractions are determined via homogenisation and energy minimisation. The phase distribution is not coupled one-to-one to the temperature distribution, e.g. the melting point, but is a result of the calculated distribution of the mass fractions which is defined within the minimisation problem. Furthermore, we consider the volume change during the changes of the material state through physically sound motivated transformation strains. In addition, further inelastic strain contributions are considered, i.e. a viscous strain in the molten phase and a plastic strain in the re-solidified phase, which are derived via dissipation functionals. Moreover, a thermal strain is included in the re-solidified material. This constitutive model is then implemented into the commercial finite element program Abaqus to simulate the SLM process. To model the manufacturing process, further models are applied, e.g. for the laser beam heat source and the layer build-up. Due to the micromechanically properly motivated material model, it is possible to predict the process induced residual stresses and the temperature evolution in a reasonable manner. In a post-processing step, the accumu-lated inelastic strains present are calculated. This information can be incorporated in a phenomenological model by using an averaged inherent strain tensor.
In future work, the previously derived and computed irreversible strains shall be averaged to an inherent strain tensor and applied to more complex geometries. With the enhanced material model at hand, it is possible to establish an efficient simulation framework to predict residual stresses and the final geometry for larger process simulations. Thus, the idea of the inherent strain method, to model more complex geometries and laser beam paths, without the need for extremely time consuming simulations, can be applied based on the previously shown simulation results. In addition, a mesh optimisation and code optimisation could save further computational time. It shall be emphasised that the material model at hand can also be used for the simulation of other additive manufacturing processes based on metallic powder, for example laser cladding. Furthermore, the used phase transformation model allows the incorporation of multiple solid phases which could preferably form depending on the specific local cooling rate and stress state. In view of the titanium alloy considered here, the related αand βphase could also be taken into account. This is supposed to have a significant effect on the residual stress prediction and may therefore improve the results. In addition, the material model can be further enhanced, for example the necessity of a temperature-dependent yield limit for the plastic strains or of temperature-dependent material parameters needs to be investigated. In the long term, due to the highly changing mass densities, it may be necessary to enhance the constitutive framework by using a large strain framework. In view of the model verification, some parameters still need to be identified which requires comparisons to available experimental data. In this way, the material and process parameters can be calibrated accordingly to further improve the simulation results.
Acknowledgements The authors gratefully acknowledge the computing time provided on the Linux HPC cluster at TU Dortmund University (LiDO3), partially funded in the course of the Large-Scale Equipment Initiative by the German Research Foundation (DFG) as project 271512359.
Funding Open Access funding provided by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copy-right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Finite element framework
The discretised weak forms of (3) and (19) shall be derived in detail by considering the contributions from both mechanical equilibrium in the absence of mechanical body forces and heat equation for a thermomechanically fully coupled system as introduced in Sect. 2. The weak form is obtained by multiplying (3) and (19) by an admissible test function for the displacement and temperature field, δu and δθ, respectively. Subsequently, the resulting equations are integrated with respect to the volume of the body B. Finally, integration by parts and the divergence theorem are applied. This yields where ρ denotes the current density, and where n describes the outward normal unit vector. Using Cauchy's theorem, surface tractions per unit area are represented by t = σ · n and the heat flux vector q n = −q · n is analogously defined.
To discretise the problem in space, body B is decomposed into n el finite elements, thus for each element node n en , where N u A , N θ C denote the corresponding set of shape functions. For the problem at hand, the same type of shape functions are used for both fields. The contributions δu e A and δθ eC represent the respective values of the test function at the nodes. With these relations at hand, and by additionally defining the assembly operator A, the completely discretised weak form can finally be given, see (92) and (93).