A simple yet consistent constitutive law and mortar-based layer coupling schemes for thermomechanical macroscale simulations of metal additive manufacturing processes

This article proposes a coupled thermomechanical finite element model tailored to the macroscale simulation of metal additive manufacturing processes such as selective laser melting. A first focus lies on the derivation of a consistent constitutive law on basis of a Voigt-type spatial homogenization procedure across the relevant phases, powder, melt and solid. The proposed constitutive law accounts for the irreversibility of phase change and consistently represents thermally induced residual stresses. In particular, the incorporation of a reference strain term, formulated in rate form, allows to consistently enforce a stress-free configuration for newly solidifying material at melt temperature. Application to elementary test cases demonstrates the validity of the proposed constitutive law and allows for a comparison with analytical and reference solutions. Moreover, these elementary solidification scenarios give detailed insights and foster understanding of basic mechanisms of residual stress generation in melting and solidification problems with localized, moving heat sources. As a second methodological aspect, dual mortar mesh tying strategies are proposed for the coupling of successively applied powder layers. This approach allows for very flexible mesh generation for complex geometries. As compared to collocation-type coupling schemes, e.g., based on hanging nodes, these mortar methods enforce the coupling conditions between non-matching meshes in an $L^2$-optimal manner. The combination of the proposed constitutive law and mortar mesh tying approach is validated on realistic three-dimensional examples, representing a first step towards part-scale predictions.


Introduction
Metal additive manufacturing (AM) opens new opportunities in manufacturing technology [10,16]. Specifically, the family of powder bed fusion additive manufacturing (PBFAM) processes allows for high geometrical flexibility and the potential for a pointwise control of material properties. However, many of the underlying complex multiphysics phenomena are still insufficiently understood and a suboptimal choice of processing parameters might lead to poor part quality in terms of pores or high residual stresses, which in turn might induce cracks in the part [32]. Currently, extensive process tuning (via trial and error) is required, which severely hinders adoption by industry. This limitation could be overcome via adequate models that allow for sound and real physics-based simulation of the whole process in order to predict optimal processing parameters. Large efforts have been undertaken to establish such approaches, all of which still face a number of obstacles when it comes to manufacturing of realistic and complex parts.
Existing modeling approaches for PBFAM are often distinguished by the modeled length scales (see overview in [31,32]). Macroscale or part-scale models determine residual stresses or dimensional warping. Residual stresses are of special interest in the community and a general review and theoretical considerations can be found in [27,35]. Mesoscale models operate on a length scale from single powder particles up to one powder layer thickness in order to study mechanisms of defect generation arising from the melt pool [21,24,29,30,48,52,53] or the powder recoating process [15,33,34]. Microscale models investigate the formation of (typically anisotropic) metallurgical microstructures during solidification [11,28,39,45,49,55].
In this work we propose a coupled thermomechanical macroscale model based on the finite element method (FEM). The focus lies on two important aspects, namely consistent constitutive modeling in a homogenized macroscale sense and elaborate layer coupling schemes allowing for layerwise non-matching finite element (FE) meshes. Macroscale models commonly solve a thermomechanical (or sometimes a pure thermal) problem, where powder, melt and solid phase are all modeled as homogenized continua with specific temperature-and phase-dependent thermal and mechanical properties [6,7,8,17,18,22,47,51]. The heat of the incident laser beam is frequently modeled with the powder bed absorption model from Gusarov et al. [12,13,14]. Further approaches that are frequently followed for an efficient simulation of PBFAM processes include spatial mesh adaptivity [23,46,9,56], sometimes combined with code parallelization and load balancing techniques [37], reduction of the computational domain via equivalent thermal boundary conditions representing the powder phase [38], process layer agglomeration [54,18,56], as well as rather heuristic approaches such as inherent strain schemes [20,26,50].
In macroscale models, all three phases are typically modeled in a Lagrangian manner based on a quasi-solid constitutive law with artificial stiffness values in powder and melt phase that are orders of magnitude lower than the stiffness of the solid phase. This approximation is well-justified since the mechanical stresses arising in the powder and melt phase at the stress-free surface of each processed layer are typically very small compared to the stresses in the bulk solid material. The interfaces between these different phases are commonly modeled as diffuse interfaces, e.g., defined by the solidus and liquidus temperature of the alloy.
Critically, such modeling approaches have to ensure that comparatively large strains, as typically arising in the melt phase (with low stiffness), do not transfer to large stresses in the newly formed solid phase (with high stiffness), which is assumed to be (approximately) stress-free right after solidification. Some contributions [17,26] use a plastic constitutive law, which can achieve the desired effect of an (almost) stress-free configuration at solidification start [25]. Other works mention a reset or annealing of (plastic) strains and stresses if material (re)melts [3,8,41]. Unfortunately, most existing works do not go into detail how a material model, that was initially designed for a single solid phase, is applied to a three-phase mixture of powder, melt and solid, and how the aforementioned condition of a stress-free solidification start can be incorporated in such models. A notable exception are the recent contributions [2,40], where a constitutive model for the three-phase mixture powder-melt-solid has been consistently derived on basis of iso-stress homogenization and energy minimization. In this interesting contribution, the phase fractions for powder, melt and solid are treated as (independent) internal variables with associated evolution equations and compatibility (inequality) constraints.
The present work builds upon the purely thermal model developed in our previous contribution [43]. The stress state of the three-phase mixture powder-melt-solid is consistently described on basis of an iso-strain homogenization across these phases. As a main contribution we motivate and derive a constitutive model which ensures a stress-free state in newly solidified material by means of a reference strain term, formulated in rate form. Compared to existing approaches, the proposed scheme allows to consistently account for this stress-free solidification start while avoiding spatial and temporal discontinuities (jumps) in the stress-field, as typically resulting from (instantaneous) stress resetting procedures. Moreover, it results in a simple computational model that does not require any additional internal variables or constraint equations, and thus, can easily be integrated in existing material libraries of standard FEM codes.
The second aspect, that is addressed in a novel way in this contribution, is the problem of growing, complex geometries. In this context, various approaches have been utilized to model newly added powder layers, e.g. so-called quiet-element methods [9], element-birth methods [37] or combinations thereof [36], most of which use layerwise conforming FE meshes. Only very few approaches, e.g., based on the immersed boundary method [4], can be found that allow for more flexible FE discretizations that do not need to conform with the part shape In this contribution we propose to apply new powder layers with dual mortar mesh tying schemes, which enable efficient condensation of constraint equations from the global system of equations. This procedure offers superior flexibility with regard to spatial discretization by allowing for non-conforming meshes between subsequent layers. As compared to collocationtype coupling schemes, e.g., based on hanging nodes, mortar methods enforce the coupling conditions between non-matching meshes in an L 2 -optimal manner. Non-conforming interface discretizations between powder layers can significantly simplify mesh generation for complex part geometries, especially, when the cross section of the part changes rapidly. We demonstrate this aspect and the general capabilities and robustness of the mortar approach based on a larger three-dimensional example.
Finally, it should be noted that this work is mainly concerned with the numerical modeling of PBFAM, a family of processes with selective laser melting (SLM) as the most prominent member. However, the methodology presented herein may also be applied to the modeling of directed energy deposition or wire-feed AM processes.
The remainder of this article is structured as follows: Section 2 presents the underlying mathematical model of thermomechanics. More specifically, Section 2.1 reviews the modeling of temperature-and phase-dependent thermal material parameters and Section 2.2 presents the newly proposed solid material model. Its numerical treatment is detailed in Section 3 along with the general numerical solution procedure. The model is verified by elementary test cases allowing for analytic solutions in Section 4. Section 5 introduces the mortar mesh tying framework for adding new powder layers. Its capabilities are demonstrated in a variety of three-dimensional examples in Section 6. Finally, Section 7 concludes the article with a summary of the most important results and observations.

Mathematical problem statement
The considered thermomechanical problem consists of the dynamic heat equation coupled with the static balance of linear momentum: with the primary variables temperature T and displacement u. The magnitudes of displacements and rotations as arising from typical PBFAM processing conditions can be assumed as small, hence the problem is commonly modeled within the theory of geometrically linear continuum mechanics where the (engineering) strain tensor is defined by: The coupling between the two equations (1) and (2) is for now hidden in the structural material law σ = σ(ε(u), T ) which will be discussed in detail in Section 2.2. The heat flux q is specified by Fourier's law of heat conduction, The material parameters appearing in the thermal problem, namely volumetric heat capacity c and heat conductivity k, will in general depend on the temperature and phase. Their modeling is discussed in detail in the authors' publication [43] and is briefly reviewed in Section 2.1. The source termr is used to model the incident laser beam power based on [14].
Remark. Technically, the heat equation (1), which is derived from energy conservation, can contain coupled mechanical terms. The present, geometrically linear problem formulation without these coupling terms results in a one-way coupling, i.e., the temperature field influences the structural field, but not vice versa. This assumption is ubiquitous in the macroscale PBFAM simulation literature [3,8,17,18] and seems justified as strain-rate dependent heating effects are not relevant for a quasi-static solid mechanics problems. A model that considers the two-way coupled thermomechanical problem can be found in [40].
The initial boundary value problem is completed by initial conditions for the temperature field and Dirichlet and Neumann boundary conditions for the thermal and structural problem: in Ω for t = 0, q · n =q, on Γ q , u =û, on Γ u , where Ω is the problem domain, Γ T the Dirichlet and Γ q the Neumann boundary of the thermal problem, Γ u the Dirichlet and Γ σ the Neumann boundary of the solid problem and quantitieŝ (·) are prescribed values on the respective boundaries. No initial conditions are required for the quasi-static balance of linear momentum (2). In the present work an apparent capacity method accounts for the effects of latent heat. Essentially, this method modifies the heat capacity c, details on the derivation can again be found in [43].
Remark (Static vs. dynamic solid mechanics problem). The balance of linear momentum can either be treated as a static or dynamic problem. The PBFAM process can be assumed as quasi-static as no high accelerations takes place (in the solid phase) and inertia effects are thus negligible. Consequently, most of the literature focuses on the static structural problem [8,9,19,41]. Note that history-dependent, potentially irreversible phenomena, which play an important role for the considered class of phase change processes, are still captured by means of history information in the thermomechanical constitutive equations.

Temperature-and phase-dependent parameters
This section briefly summarizes the modeling of the three different phases powder, melt and solid. For a more detailed motivation and derivation, the reader is referred to [43]. The commonly used liquid fraction g is introduced as where T s and T l represent the solidus and liquidus temperature. The irreversibility of the powder-to-melt transition is captured via the consolidated fraction if r c (0) = 1 (i.e., initially consolidated) max t≤t g(T (t)), if r c (0) = 0 (i.e., initially powder) (11) The resulting fractions of powder (p), melt (m) and solid (s) are computed as These phase fractions can be used to interpolate arbitrary material parameters: where f interp is the interpolated parameter and f p , f s and f m are the single phase parameters. This technique is applied to the thermal conductivity k and the heat capacity c. For the mechanical material properties we refer to the next section.

Mathematical formulation
An iso-strain homogenization (also known as Voigt-type homogenization) assumes that the strain in all phases is identical. Accordingly, the stress of the mixture is given by a weighted sum of the individual contributions, a procedure that is in fact similar to the interpolation scheme (15): Based on the iso-strain assumption, the total kinematic strain (3) is equal for all phases. For each single phase i ∈ p, m, s, it can be additively split according to although not all terms will be utilized for each phase. The first term on the right-hand side of (17) is the elastic strain ε σ,i which induces a stress σ i in each phase according to a linear hyper-elastic material where C i is the fourth-order constitutive tensor C i = λ i δ ab δ cd + µ i (δ ac δ bd + δ ad δ bc ), λ i = The artificial Young's modulus in powder, E p , and melt, E m , will be chosen orders of magnitude below the physically consistent value of the solid, E s . The Poisson's ratio is assumed to be the same in all phases. The remaining terms in (17) are inelastic contributions which are considered in more detail in the following. The plastic strains ε p,i , which are only relevant in the solid phase, could be calculated with standard approaches, e.g., an incremental problem formulation in combination with a return mapping algorithm. For simplicity, however, plastic strains will not be considered in the numerical examples in this work (ε p = 0). The strains due to thermal expansion ε T are assumed equal in all phases and read where α T is the (constant) coefficient of thermal expansion and T ref is a reference temperature.  sṙ sεref = 0 givenṙ s < 0. This case is necessary for a consistent notation but, as we will see later, can be circumvented in practice. Note, how the rate formulation in (21) causes a continuous change in the stress over the phase change interval [T s ; T l ], which is beneficial for a numerical solution, in contrast to existing approaches with an instantaneous reset of stresses at melting temperature.
For completeness, all introduced strain contributions can be inserted into (16), which after some rearrangement yields the following total stress of the phase mixture: The pre-factor of the first term in (22) is equivalent to an average of the single phase material parameters weighted with the phase fractions r i similar to (15).
Remark (Modeling assumptions). One of the main assumptions underlying the present and most existing thermo-mechanical PBFAM models is that mechanical stresses in the (open-surface) powder and melt phase domains are negligible. This behavior is approximated by applying a simple elastic constitutive law to these phases, with stiffness parameters that are considerably lower as compared to the solid phase, i.e., E p , E m Es . In practice, this approximation turns out to result in moderate, i.e., limited, strains, since the deformation of these powder and melt domains is mostly kinematically controlled by the motion of the significantly stiffer solid phase domains, thus yielding only small stress contributions as desired. Moreover, as compared to approaches exactly satisfying the zero-stress assumption in powder and melt, no additional means are required for tracking and discretization of sharp interfaces inside elements. Note, the assumption that thermal strains exist also in the powder and melt phase, and are equal to thermal strains in the solid phase, has been made for simplicity here. This assumption is neither necessary nor has it a significant influence on the resulting residual stresses due to the low stiffness of these phases and the definition of reference strains (21), which ensures that newly created solid material is stress-free. Further, we assume that from the beginning powder already has the volume and density it would have after consolidation, which has to be accounted for by defining correspondingly decreased layer thicknesses. Modeling solidification shrinkage, i.e., a density increase when powder consolidates, is deemed unnecessary due to the free surface at the top of the currently processed power layer, which allows for (approximately) stress-free consolidation and shrinkage in thickness direction when the powder melts.
The irreversibility of phase change and the reference strain (21) make the material behavior non-conservative. While the proposed approach is very general and can be combined with arbitrary (standard) solid material laws, e.g., elasto-plasticity, we purposefully restrict ourselves to purely elastic solid material behavior in the studied examples. In this scenario, the proposed reference strain term is the only non-conservative contribution to the overall material model, underlining that the creation of a stress-free state at solidification start is the most significant non-conservative aspect of the overall thermo-mechanical problem.
In the simple case of purely elastic material behavior (i.e, stresses do not exceed the yield stress) a heating to a maximum temperature below melting temperature and subsequent cooling to the initial temperature would not lead to any residual stress. Only when the melting temperature is exceeded, the reference strain term will cause an (approximately) stress-free configuration at solidification start, and thus, a residual stress will remain after cooling to the initial temperature. Therefore, we identify the reference strain contribution as the minimal necessary effect for residual stress prediction in such a simplified model. At time t + ∆t, the newly solidified fraction takes on the (now changed) current configuration (red) as reference configuration (c), such that the effective reference strain (which in fact is calculated in the proposed model) in the (total) solid phase can be interpreted as weighted average (purple) of the contributions from individual solid phase fraction increments (d).

Physical motivation and discussion
The reference strain (21) and stress (22) form a simple yet consistent solid constitutive law.
In the following, we want to discuss some properties of the material law and their physical motivation by means of analytically tractable cases. For simplicity, the plastic strain terms are neglected from now on. Note, that the reference strains ε ref in (21) only change when the solid phase fraction increases according toṙ s > 0, i.e., for temperatures T ∈ [T s ; T l ] in the phase change interval and negative temperature ratesṪ < 0. An elastic constitutive law with low stiffness values (i.e., E p , E m Es) as applied to powder and melt leads to small stresses yet considerable total strains in these phases. In this context, the reference strains according to (21) ensure that these strains do not translate into stresses during solidification. For the special case that kinematic ε and thermal strains ε T (as well as plastic strains ε p ) are constant during solidification, which approximately holds if the phase change interval T l − T s is sufficiently small, it can easily be verified from (17) and (21) that the elastic strain, and thus due to (18) the resulting stresses, in the evolving solid phase vanish. This fact corresponds to the physical intuition that newly formed solid should lose all history information and exhibit a new stress-free configuration when solidification starts.
From a stress homogenization point of view, the model assumptions underlying the reference strain formulation state that for each solidifying fraction of material, the reference strain contribution effectively creates its new, stress-free reference configuration, from which strains are calculated. This is illustrated in Figure 1.
Two more involved examples, which illustrate the evolution of the reference strain over multiple repeated melt and solidification cycles, can be found in Appendix A.

Numerical solution procedure
Spatial discretization of the partial differential equation is based on the FEM. The required weak form of the thermomechanical problem (1) and (2) is obtained via multiplication with test functions δu and δT and integration by parts, resulting in where the boundary conditions have already been inserted. (23) and (24) are equivalent to the strong forms (1) and (2)  T where x is the spatial position and N j (x) are the time-independent shape functions, which are the same for all solution and test functions. T j and δT j are the discrete nodal temperatures and their variations, and d j and δd j are the discrete nodal displacements and their variations, all of which depend on time. Although the static balance of linear momentum (2) is considered, i.e., inertia forces are neglected, the displacement field u(x, t) is depending on time via the coupling to the temperature field T (x, t), which is a solution to the dynamic heat equation (1). The thermal subproblem is discretized in time with a generalized trapezoidal rule. The discrete system of equations is consistently linearized and solved with Newton's method.

Discretization of the rate-based reference strain
Time integration of equation (21) is based on a backward Euler scheme and results in the following time-discrete form of the fraction-weighted reference strainε n+1 ref : which leads to the actual reference strains via relation (21) as ε n+1 Thus, the time-discrete total reference strain ε n+1 ref may also be obtained directly, viz.
Note, that in (28) the case of melting material (∆r n+1 s < 0) no longer needs to be treated separately due to the construction ofε ref as already discussed when it was first introduced.

Linearization of constitutive law
The linearization of the total stress (22) with respect to the primary solution variables is required for the nonlinear solution procedure. As usual, the derivatives of the constitutive equation are written with respect to the kinematic strain and temperature, where H (x) is the Heaviside step function All phase fractions and their derivatives are evaluated at T n+1 . The first two terms in (30) stem from the phase interpolation in the first term of (22). The third term represents the contribution from a solidifying increment and is efficiently derived by inserting (28) into the total stress (22) and thereby canceling out the solid fraction r s . The last term represents the effect of a melting solid.

One-and two-dimensional numerical examples
This section focuses on the verification of the proposed material law with a series of elementary test cases. In all of them the temperature is effectively a prescribed time-dependent function that drives the structural simulation and no heat equation needs to be solved. The investigation starts out on a one-dimensional domain which is subject to different boundary conditions and temperature profiles with increasing complexity. These findings are reported and discussed in sections 4.1 and 4.2. Complexity is increased further by applying different temperature profiles to a two-dimensional domain in Section 4.3. The first series of examples examines a one-dimensional bar of length l = 1 mm, which is subject to a homogeneous temperature load, i.e., the temperature evolution is only a function of time but spatially constant, leading to melting and solidification of the material, possibly multiple times. Figure 2 illustrates the two considered types of boundary conditions and how the respective one-dimensional problem relates to a three-dimensional solid mechanics problem. The behavior of the bar is equivalent to the behavior of the depicted cube in x-direction. The bar is either constrained by a Dirichlet-type condition (Figure 2 left) prescribing the deformation or a Neumann-type condition (Figure 2 right) prescribing the traction. The prescribed values can be zero (homogeneous boundary conditions) or non-zero (inhomogeneous boundary conditions). The initial material phase is either solid or powder and all relevant material parameters are listed in Table 1.
Three representative examples are investigated: first, a full melt of the material, and second, a repeated partial melt followed by a full melt. Both scenarios use homogeneous boundary conditions. For completeness, the partial melt scenario is repeated with inhomogeneous boundary conditions as a third example. The numerical results for all examples are compared to an analytical reference solution, which assumes isothermal melting (T s = T l = T m ) and a zero stiffness in powder and melt phase.

Full melt
The material is completely molten by heating it to a peak temperatureT = 2200°C > T l . Figure 3 shows the temporal evolution of the (homogeneous) strain and stress state for both, Dirichlet and Neumann, scenarios as well as the evolution of the effective Young's modulus, E = r i E i , of the phase mixture. While the Neumann scenario is mostly shown for completeness, the Dirichlet scenario is more interesting regarding the difference between an initially powder and solid material: before melting takes place (t 2) the stress will stay close to a near-zero value in the case of initial powder (small stiffness) but take on significant values in the case of initial solid. When the material melts (t ≈ 2), the stress reduces to the same small value as in the powder case. For a finite but small phase change interval, as chosen for this problem, this change happens continuously. On the other hand, the analytical reference solution for the stress, assuming isothermal melting, exhibits a discontinuity at melting temperature. The proposed model can represent this discontinuity asymptotically correct when decreasing the phase change interval. After full melting, the two stress curves for initially powder and solid material are identical since in both cases all material is now completely (re)molten and has the same reference strain. At the end of the simulation the same final stress σ final = 2 MPa is obtained. This value can be calculated analytically, see Appendix B.1.
In the scenario with a homogeneous Neumann boundary condition no stress can occur. The  total strain is always equal to the thermal strain, ε = ε T , and thus the contributions to the reference strain in (21) will always be zero. Therefore, the strain is directly proportional to the change in temperature ∆T = T − T 0 , as shown on the right-hand side in Figure 3.

Repeated partial melt followed by full melt
In a second set of examples the material is partially molten with a peak temperatureT 1 inside the phase change interval, viz. T s <T 1 = 2000°C < T l , and then cooled down to T 0 . This cycle is repeated four times and, finally, the material is fully molten with a peak temperaturê T 2 = 2200°C > T l . The resulting strain and stress evolution is shown in Figure 4. Once again, attention is drawn to the stress evolution in the Dirichlet case which is equal to the results from the previous section until the first peak in the temperature profile is reached. This time, the stress does not reach a near-zero value for the initially solid case because the melting stops atT 1 < T l , i.e., there remains a phase fraction of solid material atT 1 , which exhibits thermal stresses due to the increased temperature level as compared to the stress-free configuration at T 0 . Subsequently, the stress rises to the same local peak value of 975 kPa for both initially solid and powder material after cooling to reference temperature. This can be explained as follows: the only relevant contribution from (22) is the last term containing the reference strain since the temperature is equal to the initial reference temperature and the strain is identical to zero. Although the solid fraction differs at this point (0.5 in powder case and 1 in solid case), multiplication with the reference strain (21) yields the same contributionε ref . Put differently: independent of the initial phase, the non-molten phase fraction yields a stress of zero, since T = T 0 ; the (re)molten phase fraction yields the same stress contribution, since the (re)molten fraction is equal for initially powder and solid material (same maximum temperature). The specific value of 975 kPa can be calculated analytically as demonstrated in Appendix B.2.  The same heating and cooling cycle is repeated three times. For an initial powder material the stress reaches a near-zero value atT 1 because the solid phase (from the previous partial melt) is completely remolten at this point. For an initially solid material, the situation was already discussed in Section 2.2.2: the material model contains no information about different solid phases but instead averages the reference contributions according to (28). As already outlined before, a repeated partial melt with the same peak value will lead to increasing reference strain and, therefore, also stresses (when cooled to reference temperature), a fact that is well illustrated by the simulation results in Figure 4. The analytical reference solution for isothermal melting of initially solid material is also included and again only differs for temperatures in the phase change interval (although not visible in Figure 4).
Finally, the material is fully molten with a peak temperatureT 2 = 2200°C > T l and cooled to T 0 . Both an initially solid and powder material yield the same final stress value of 2 MPa, which is the same value as in the previous example, where the material was fully molten directly in one thermal cycle. Once more, this confirms that all history information related to stresses is erased when all solid is (re)molten.

Inhomogeneous boundary conditions
The scenarios so far were only concerned with zero Dirichlet and Neumann boundary conditions. For completeness the partial melt case in Section 4.1.2 is repeated with an inhomogeneous Dirichlet-type constraint u(l) = 0.001 mm, which leads to a constant strain ε x = 0.001 for the given homogeneous temperature load. Figure 5 shows the resulting residual stresses, which, for an initially solid material, differs from Figure 4: the fixed strain induces an initial stress that decreases during the repeated partial melting cycles and, eventually, completely vanishes after the full melting, during which

One-dimensional domain: inhomogeneous temperature load
The next scenario focuses on a more complex, transient and inhomogeneous temperature profile and a homogeneous Dirichlet boundary condition u(l) = 0 (see left-hand side in Figure 2). A moving temperature peak sketched in Figure 6, that starts outside the initially powder domain (length l = 1 m) and moves through it to the other side, mimics the effect of a moving laser beam in PBFAM. For completeness, the mathematical form of such a temperature profile is stated as where the parameters for the present example are chosen as T 0 = 0, T max = 2200°C, w = 0.1 m, x 0 = −w and v = 1.0 m/ sec. With these parameters every point in the domain melts and solidifies. Figure 7 shows the temperature profile and the resulting displacement and strain field at select time steps. The stress, which is constant over the domain, is plotted over the whole simulation time. In the beginning, the low stiffness of the remaining powder at the right side of the domain enables an almost free thermal expansion and the stress is close to zero. It increases to significant values only when the whole bar has a solid fraction r s > 0, i.e., when the temperature peak reaches the right end of the bar, where the displacement is fixed to zero.   The current example is also used to judge the influence of the stiffness ratio Es Ep,m between the solid and powder/melt Young's moduli by varying the (artificial) powder/melt stiffness E p,m . Figure 8 shows how the final stress in the bar (after cooldown to the reference temperature) converges with a refined spatial discretization and for different stiffness ratios. For a high stiffness ratio the limit value lies at around 0.1 MPa, which one would obtain as the analytical value assuming the theoretical value of zero powder and melt stiffness, i.e. a stiffness ratio of infinity. A proof is sketched in Appendix B.3.

Two-dimensional domain: interaction of boundary condition and temperature load
The final verification example investigates a two-dimensional setting under the assumption of plane-stress 1 . As in the previous section, the geometry consists of a slender, bar-like structure of powder material which is now constrained with a homogeneous Dirichlet boundary condition on the bottom edge (z = 0), while the left, right and top edges are free (homogeneous Neumann boundary condition), see Figure 9. The motivation for this setup is to mimic the processing of a powder layer which is applied atop an already solid domain. For this example we use the realistic material parameters from Table 2.
As a baseline the domain is heated uniformly across the entire domain from initial temperature T 0 = 303 K to a maximum temperature T max = 2000 K > T l and then cooled down to the initial temperature. This case is compared to the scenario of a moving temperature peak profile as in (32)      z-direction, which is a rough approximation to the temperature profile encountered in a PBFAM application. Figure 10 compares the resulting residual stresses σ xx in the different scenarios: a uniform temperature profile yields the highest tensile stresses and almost no variation of the stress σ xx in z-direction (the "powder layer height") while a moving temperature peak profile leads to a significant drop-off of σ xx over the height, meaning that stresses are lower at the non-constrained upper surface. This drop-off becomes more pronounced with a decreased peak width w and can be more clearly observed in Figure 11 (left). These observations can be explained by the size of the domain which has a temperature above T 0 and is cooling down: if this region is small (for a small w), the deformation occuring during cooldown is almost unconstrained since these solid material regions with elevated temperatures are very close to the (approximately) stress-free liquid domain, thus resulting in a comparatively small tensile stress.
If the cooling region is wider, large portions of the cooling solid material are far away from the stress-free solid-liquid boundary and thus are strongly constrained in their ability to deform, which consequently results in a higher tensile stress. Next, the temperature profile shall vary in z-direction as well. Given a time and x-dependent function T z0 (x, t) for the temperature profile at z = 0 the following function allows to control the variation of the temperature over the height h = 0.05 mm. The scaling factor s describes the ratio of the temperature at the upper edge z = h compared to the constrained bottom z = 0. For better comparison, it is calculated from the width w in (32) as such that the partial derivative of (33) with respect to z is comparable to the slope of the moving temperature peak of width w in (32). First, the temperature only varies in z-direction but is constant in x-direction, which is achieved by (33) and the purely time-dependent function The scaling factor s is computed from w according to (34). The results in Figure 11 (right) indicate that a variation only in z-direction does not lead to significantly different stresses (compared to the uniform temperature field), even for extreme values such as w = 0.01 mm, i.e., a scaling factor of s > 7 between upper and lower edge temperatures. All material points with the same z-coordinate solidify at the same time.
The situation becomes more interesting when the moving peak is combined with a variation in z-direction, which is closest to the real temperature profiles expected in PBFAM processes. This is easily achieved by inserting (32) as T z0 into (33) to obtain another kind of temperature profile, see Figure 12 (right). The width (in x-direction) of the peak again varies between w ∈ {0.4, 0.2, 0.1} K. Figure 12 (left) shows how the stress σ xx varies over the height z. Especially for small widths, the additional variation in z-direction has an influence on the resulting stress compared to the corresponding cases without z-variation.
We do not claim that the results in this section are directly transferable to a PBFAM process simulation. Still, they suggest that the actual temperature profile might play a more important role than often assumed. In the authors' opinion, the significant differences in residual stress from a uniform and a moving temperature field raise the question how well-suited common simplifications in PBFAM simulation, such as layer-agglomeration, inherent-strain or flashing complete tracks at once, are for an accurate prediction of residual stress distributions.

Multilayer mesh tying strategy
From the PBFAM process perspective, it is natural to slice the complex part geometries into layers, which may either represent a single powder layer or an upscaled process layer, i.e., an accumulation of several powder layers. It seems that many existing approaches, even in commercial codes, suffer from the need to generate matching meshes between layers for real application scenarios and complex geometries. In this work, we want to propose a flexible, layerwise spatial discretization strategy without the need of matching meshes across the interface, thus allowing for easy mesh generation without distorted elements or large size differences between elements. Dual mortar methods allow to couple such non-matching FEM meshes in a manner that results in optimal discretization errors (measured in the L 2 -norm) and that leads to an efficient condensation of the constraint equations from the global system of equations.
This section summarizes the most important aspects of the approach for the application to PBFAM, while a general overview on mortar methods for contact and mesh tying problems Figure 13: Mesh tying concept for multilayer PBFAM simulations.
is found in [42,44]. Figure 13 shows an exemplary part geometry which is sliced into layers Ω ( ) with thickness equal to the powder thickness h p . Each layer can be meshed separately and the finite element nodes may be non-matching at the interface between two layers. The continuity of the primary solution variables across the interface Γ mt, is enforced in a weak sense by additional (mesh tying) constraint equations. The -th mesh tying interface Γ mt, consists of the boundaries Γ These conditions are enforced by introducing a Lagrange multiplier potential for each constraint: The total variation of the two potentials is found as  (40) and (41) represents the original mesh tying constraint and the second term contributes to the weak forms (23) and (24). From dimensional analysis it becomes apparent that the unknown structural Lagrange multiplier can be interpreted as a force vector acting on the interface, λ u = t mt . The Lagrange multiplier fields are discretized with special dual shape functions that allow for an efficient condensation of the Lagrange multiplier DOFs from the global system of equations [44]. After spatial discretization of all primary variable fields, the consistent link between discrete nodal DOFs on both sides of the interface is established via constant mortar matrices D and M [42], which can be combined into constant projector matrices P = D −1 M . The projector matrices map discrete solution increments (from the nonlinear solution procedure), where ∆d and ∆T are increments of the discrete solution vectors of discrete displacements d and discrete temperatures T . Given that the mesh tying constraints (36) and (37) hold initially, these incremental mappings will enforce the constraints also in later configurations. Typically this is easy to ensure in problem settings which utilize a tied interface from the beginning. However, in PBFAM applications, new powder layers are added over time, thus introducing new mesh tying constraints. In general, a new powder layer's initial temperature and displacement will not conform to the, already processed, last layer. A simple approach to initialize temperature and displacements on side Γ ( ) mt, of a newly added layer is to apply mappings (42) and (43) to the solution vectors on Γ In essence, the nodal solution at the bottom surface of a newly added layer is set to the (consistently mapped) nodal solution at the top surface of the last processed layer. After this initialization the meshes are correctly tied and the primary solution variables are continuous across the interface. The remaining DOFs in the newly added domain Ω ( ) \Γ ( ) mt, are set as follows: temperature DOFs are set to the initial temperature T 0 , while displacement DOFs are initialized to zero. The first time step after a new mesh tying interface is activated will solve the (static) mechanical equilibrium equation (2) for a consistent displacement state.
Remark. Technically, the direct mapping (44) violates conservation of energy when adding a new powder layer: temperature and displacement on Γ ( ) mt, are modified without an associated external force. An alternative initialization of the temperature field for element-birth methods is discussed in [36] and could be transferred to the proposed mesh tying approach. However, if the cooldown phase after processing of one layer is long enough for the top surface to reach the initial temperature, the energy error from the temperature modification vanishes. The error introduced by modification of the displacements of a newly added layer is generally small due to the low stiffness of the newly added powder, i.e., possible artificial strains in the new layer only cause a small stress response. Furthermore, the material formulation from Section 2.2 ensures that these artificial stresses vanish after melting.

Three-dimensional numerical examples
This section presents larger three-dimensional numerical examples that are simulated with the research code BACI [1], jointly developed at the Institute for Computational Mechanics.

Single tracks per layer
The first three-dimensional example investigates ten processed powder layers with a single track per layer each on top of a solid base plate. A schematic of the geometry and the boundary conditions is shown in Figure 14. The computational domain consists of a 0.1 mm high base  plate that has a prescribed temperature T 0 = 303 K and zero displacements u = 0 at the bottom (indicated in blue in Figure 14). This temperature T 0 is also used as the initial temperature throughout the domain and in all newly activated powder layers. The powder layers of height h p are connected via the mesh tying algorithm introduced in the last section as soon as they are activated. The tracks are scanned atop each other with a unidirectional or serpentine pattern, i.e., the laser tracks are either oriented in the same direction in all layers (unidirectional) or the orientation alternates between successive layers (serpentine). To save computational resources, a symmetry plane along the laser track center line is used to half the size of the computational domain. All lateral faces (indicated in red) are assumed to be thermally insulating (q = 0) and subject to homogeneous Dirichlet boundary conditions.. The built part will not attach to these boundaries for the chosen laser beam path and diameter and, since the remaining powder close to the boundary is modeled with low stiffness and conductivity, the zero displacement boundary conditions will not influence the results in the final part geometry significantly. Every scanned track (pure scanning time of 0.005 s) is followed by a cooldown time of 1.0 s (which is long enough to reach a homogeneous temperature state close to the initial value T 0 in good approximation) during which significant residual stresses form. After all ten layers have been scanned, the process of cutting the part from the base plate can be simulated as follows: Figure 15: Unidirectional scanning pattern: distribution of different stress components in symmetry plane at selected time steps. Blue color represents tension, red color represents compression. Liquidus isotherm (solid black line) approximates melt pool shape.
The boundary conditions on the red surfaces are removed everywhere except on the leftmost surfaces (with normal vector (−1, 0, 0) T ) to avoid rigid body modes. The mesh tying constraint between the base plate and the first layer is removed. A structural analysis is run to obtain a static equilibrium solution for the stress and deformation state after detaching the part from the base plate.
Material parameters for the thermal and mechanical problem as well as other processing parameters are listed in Table 2. A time step of ∆t = 5 × 10 −6 s is used for the scanning phase and the initial cooling phase of 0.002 s, while the remaining 0.998 s of cooling use a time step of ∆t = 1 × 10 −3 s. The whole domain is discretized with 61 440 linear, hexahedral 8-node finite elements with edge size h ele = hp 4 = 12.5 µm. Figure 15 shows an overview of the six components of the stress tensor at selected points in time: the first and third column show stresses during scanning of layer 6 and 10 with a comparable relative laser position. The second and fourth column show stresses after these layers have been cooled down. The fifth column shows the residual stress after the part is detached from the base plate. The laser paths of the individual tracks remain clearly visible in the plots of the σ xx component (first row). In agreement with the literature this stress in scanning direction is a large tensile stress [3,5]. The highest stresses can be observed in the heavily constrained base plate. Specifically, it can be observed how geometrical compatibility during cool down lead to high tensile normal stresses σ xx in the first track (and remolten portion of the base plate) and high compressive normal stresses σ xx in the (non-molten) lower part of the build plate. In the last row of Figure 15 the shear stresses σ xz can be observed, which are -according to the mechanical equilibrium equation (2) -necessary to balance gradients of the normal stress σ xx at the beginning and end of a track, especially close to the base plate. These σ xz shear stresses, Figure 16: Serpentine scanning pattern: distribution of different stress components in symmetry plane at selected time steps. Blue color represents tension, red color represents compression. Liquidus isotherm (solid black line) approximates melt pool shape.
in turn, also induce σ zz normal stresses. Note, that the different stress components at many locations exceed the yield stress of typical materials used in PBFAM, hence the inclusion of an elasto-plastic material model, as already sketched in Section 2.2.1, is desirable for future work. The inclusion of such material models, or even of more elaborate microstructure-informed constitutive laws, will allow for a detailed quantitative validation based on experimental results.
All stresses decrease in their magnitude when the part is cut from the base plate. Figure 16 shows the resulting stress distribution for an alternating, serpentine scanning pattern. The absolute values are comparable to Figure 15. To better highlight the influence of the different scanning patterns, Figure 17 visualizes the normal stresses σ xx for the unidirectional and serpentine scanning pattern with a rescaled color map.
Selected time steps are plotted in Figure 18 for a more detailed view of the stress evolution of σ xx during scanning of track 7 with a unidirectional pattern. The following considerations start from a fully cooled layer 6, which exhibits a tensile stress throughout the consolidated material. The subsequent snapshots (from left to right) in the first and second row of Figure 18 show the processing of track 7. To visualize the melt pool size, the temperature isoline corresponding to the liquidus temperature T l is depicted in these snapshots (solid black line). While no stresses occur in the powder material of layer 7 in front of the laser, the thermally induced volume expansion during heating leads to negative (compressive; red color) stresses in the solid material of the previously processed track 6, mostly pronounced in the direct vicinity of the T l -isoline. As desired, these stresses rapidly drop to zero in the narrow temperature interval T ∈ [T s ; T l ] such that no visible stresses remain in the melt pool domain. This strong gradient between vanishing stresses in the melt pool and high compressive stresses in the solid material beneath remains after solidification and is superimposed by additional tensile stress contributions due to  the thermally induced volume shrinkage during cooling. After cooling down (see snapshot at bottom right of Figure 18), this process results in high tensile stresses in the upper, re-molten part of track 6, and stresses close to zero in its lower part. The same characteristics are observed in the previously processed tracks below. Even though the base plate has the same stiffness as the solidified tracks above, this characteristic band structure of the normal stresses is much more pronounced for the first track, i.e., the highest overall tensile stresses occur in the first track, accompanied by compressive stresses of comparable magnitude in the base plate below. This observation can be explained by the fact that the base plate is initially stress-free, while solidified melt tracks are subject to tensile stresses after cooling down (e.g., snapshot at bottom right of Figure 18), which partly compensate the compressive stresses arising from thermal expansion when processing the subsequent track above. The tensile stresses in layer 1 are transferred to the base plate as concentrated tensile stresses in the vicinity of the part edges, which may lead to notch effects in practice. It should be mentioned once more that the displacement Dirichlet boundary conditions on the lateral faces of the base plate are not responsible for the high stresses and results are almost identical if these conditions are replaced with zero Neumann boundary conditions.
Finally, we investigate a slightly modified scenario, where the base plate is scanned first with the same heat source as the powder layers. Figure 19 shows the evolution of stress σ xx . Large   Figure 19). The subsequent processing of layer 1, using a unidirectional pattern, leads to remelting in the uppermost region of the base plate, which causes a reduction of the tensile stress after the next cooldown phase. The strong gradient between tensile and compressive stresses, which was observed close to the interface between base plate and layer 1 in Figure 18, shifts down in Figure 19 while the magnitude stays roughly the same. For comparison, Figure 20 shows a detailed view of the evolution of stress σ xx when the first layer is scanned atop a stress-free (not scanned) base plate, i.e., the scenario that was initially investigated in this section. The stress distribution in layer 1 after cooldown for this case (right-most plot in Figure 20) is equal to the stress distribution in layer 1 after cooldown with a pre-processed base plate (bottom right in Figure 19), which demonstrates that the chosen pre-processing of the base plate does not influence the stress in layer 1.

Multiple tracks per layer
In the second three-dimensional example, a closer look is taken at the classic in-plane serpentine scanning pattern. The setup is shown in Figure 21: two powder layers are deposited on top of a base plate and five tracks following a serpentine pattern are scanned in each layer. Afterwards the domain is allowed to cool down for 1.0 s. The boundary conditions are essentially the same as in the last example and the material parameters from Table 2 are reused. As before, a time step of ∆t = 5 × 10 −6 s is used for the scanning phase and the initial cooling phase of 0.002 s, while the remaining 0.998 s of cooling use a time step of ∆t = 1 × 10 −3 s. The domain is discretized with 153 600 linear, hexahedral 8-node finite elements with edge size h ele = hp 4 = 12.5 µm.
In this example, we look at the stress distribution after each layer was processed. Figure 22 visualizes the stresses σ xx and σ yy in the consolidated material, which represents the built part. One quarter of the consolidated volume is cut out for visualization of the stress distribution in z-direction. In the first layer, the laser travels in x-direction, which is clearly reflected in the stress distributions shown in the first column of Figure 22. When the second layer is applied and scanned in y-direction, the melt pool penetrates into the first layer. Therefore, after the second layer was scanned, the stress distribution in the first layer reflects the laser beam movement in y-direction, see the second column in Figure 22). At the same time, stresses in the second layer also align with the y direction, see the third column in Figure 22. Overall, only tensile stresses remain after the final cool down.

Two pyramids example
The final numerical example investigates a more complex geometry. Yet, the geometry is defined simple enough so that it can easily be adopted by other researchers. The powder material that is not supposed to contribute to the final part is not included, an assumption which is justified by the low thermal conductivity and stiffness of the powder. The final part's cross section changes drastically over the layers which demonstrates the flexibility of the mortar mesh tying approach. The geometry, see Figure 23, consists of two pyramid bodies on a build plate: the first pyramid has curved surfaces 2 and consists of bulk material, while the second pyramid is hollow. These geometries are chosen to illustrate different behavior in bulk geometries versus slender geometries. In each layer, first the bulk pyramid is scanned with a serpentine pattern with a maximal hatching space (space between two neighboring tracks' center lines) of 0.12 mm. Then, the hollow pyramid is scanned with a closed quadrilateral track, starting from the corner with minimal y-and z-coordinate and first moving into y-direction. After each processed layer, the domain is allowed to cool down for 1.0 s. In the next layer, the serpentine scanning pattern in the bulk pyramid is rotated by 90 degrees, in analogy to the tracks depicted in Figure 21. However, the quadrilateral track along the hollow pyramid wall always follows the same qualitative shape and only increases in size.
In terms of boundary conditions, the displacements on the bottom of the base plate are constrained to zero and the temperature is fixed to the initial temperature T 0 . All other faces are unconstrained, i.e.,thermally insulating and not loaded by external forces. Therefore, the only cooling mechanism is heat transfer over the base plate and all energy supplied through the heat source must eventually be dissipated this way. The example geometry is still small enough that this approach is feasible and no significant convective surface heat transfer would occur in the relevant time scales of this problem.
Again, a time step of ∆t = 5 × 10 −6 s is used for the scanning phase and the initial cooling phase of 0.002 s, while the remaining 0.998 s of cooling use a time step of ∆t = 1 × 10 −3 s. Figure 23 gives a qualitative view of the spatial discretization. Elements are relatively undistorted and roughly equally sized. The use of a non-matching discretization between layers and between base plate and the bulk pyramid is clearly visible, and is enabled by the mortar mesh tying approach presented in Section 5. Figure 24 depicts the normal stress components as well as the von Mises equivalent stress 3 after the full part was processed. Again, the highest stresses appear close to the strongly constrained base plate. With increasing distance to the base plate the stresses relax: this happens faster in the less stiff hollow pyramid compared to the bulk pyramid. The stress distribution σ zz in vertical direction exposes a zone of compressive stresses inside the bulk pyramid and tensile stresses in the outer region, which is in agreement with the results in Section 6.1 and results in literature [18,26]. Note, that parts of the material close to the edges in the first layers were only partially molten, hence the artifacts in the stress distributions. The observations for stresses are complemented by Figure 25: here, the displacement magnitude is visualized on a deformed mesh (with deformation scaled by a factor of 5 for improved visibility). In the stiffer bulk pyramid the overall deformation is small, while the hollow pyramid can deform much easier and an undesired indentation of the side walls occurs as consequence of residual stresses.

Conclusion
A thermomechanically coupled computational model for the macroscale simulation of PBFAM processes was presented. A first emphasis of the present work was the consistent modeling of material behavior in a macroscale sense. To this end, a thermo-elastic material model with an additional inelastic reference strain contribution was derived and verified with academic test cases. Crucially, this inelastic reference strain contribution, formulated in rate form, was introduced to allow for a stress-free state when solidification of molten material starts. Compared to existing (instantaneous) stress-reset procedures, which might result in jumps in the stress-field, the formulation of the reference strain term in rate form guarantees a smooth stress field. Compared to more involved material models, e.g. [3,40], this simple approach can easily be integrated in standard material libraries and FEM codes while still consistently capturing the most essential aspects of residual stress generation in PBFAM. This capability has been verified in a series of elementary test cases and via reference solutions stated for the sharp-interface limit. The material model did not consider plasticity, although its inclusion is possible in a straight-forward manner as already demonstrated in the stated model equations.
For the coupling of successively processed powder layers we suggested a dual mortar mesh tying approach. Compared to other popular methods, such as element-birth or quiet-element methods, it provides additional flexibility with respect to the spatial discretization, as it allows for layerwise non-matching meshes. The applicability of both novelties to larger three-dimensional problems was demonstrated. Especially the rapid cross section changes of a rather complex geometry, involving a hollow pyramid as well as a solid pyramid with curved surfaces, illustrated the effectiveness of the proposed mortar mesh tying strategy. Table 3: Exemplary evolution of reference strain for repeated cooling and heating assuming constant strains during solidification. Table 4: Exemplary evolution of reference strain for a cyclic, partial remelting of solid phase assuming ε − ε T is constant in all phase change intervals. strain ε−ε T is assumed (approximately) constant during the complete solidification process. The reference strain of the non-molten solid remains at its initial value, i.e., zero. The total reference strain after cooldown is thus ε 2 ref = 0.5(ε − ε T ). In the next cycle, again 50% of the solid phase melts (at t 3 ). After cooldown below T s (at t 4 ) the new reference strain is computed according to (45) as a weighted average of the reference strain of the non-molten solid (at t 3 ) and the new strain contribution ε − ε T for the remolten fraction, ε 4 ref = 0.5ε 3 ref + 0.5(ε − ε T ) = 0.75(ε − ε T ). The same logic applies to the next heating at t 5 and cooling t 6 and further heating-cooling cycles. Apparently, the evolution of the pre-factor in the reference strain over these cycles follows a geometric series that converges to ε − ε T . If this process is repeated many times, this has the (at first glance) paradoxical result that, although never fully molten, the complete initial solid will asymptotically convert to newly solidified solid. From a microscale perspective this observation can be interpreted as follows: within a representative volume, the choice of the solid material fractions (i.e., ensembles of molecules) that will actually melt is random and will change during the repeated partial melting cycles. Thus, after sufficient partial melting cycles each solid material fraction has molten at least once such that the total solid phase effectively exhibits a new stress-free reference state.

B Analytical reference solutions for one-dimensional problem
This section presents analytical solutions for verification of the one-dimensional scenario described in Section 4.

B.1 Final stress after full melt and solidification
The homogeneous Dirichlet boundary condition and homogeneous temperature load used in Sections 4.1.1 and 4.1.2 imply ε ≡ 0. Thus, the final stress after a full melt and cooldown to initial temperature is found from (22) by integrating (21) after a change of variables from solid fraction to temperature: where the average of solidus and liquidus temperature can be interpreted as the melting temperature T m := Ts+T l 2 . Result (46) is equivalent to the expected solution for isothermal phase change at melting temperature T m and the exact value of the (artificial) powder and melt stiffness is irrelevant for the final stress value in this specific case.