proell@lnm.mw.tum.de Institute for Computational This article proposes a coupled thermomechanical ﬁnite element model tailored to the Mechanics, Technical University macroscale simulation of metal additive manufacturing processes such as selective of Munich, Boltzmannstraße 15, 85748 Garching b. München, laser melting. A ﬁrst focus lies on the derivation of a consistent constitutive law on basis Germany 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 conﬁguration 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 solidiﬁcation scenarios give detailed insights and foster understanding of basic mechanisms of residual stress generation in melting and solidiﬁcation 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 ﬂexible 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 -optimal manner. The combination of the proposed constitutive law and mortar mesh tying approach is validated on realistic three-dimensional examples, representing a ﬁrst step towards part-scale predictions. Keywords: Constitutive model, Thermomechanical coupling, Metal additive manufacturing, Mortar method, Finite element method Introduction Metal additive manufacturing (AM) opens new opportunities in manufacturing technol- ogy [1,2]. Speciﬁcally, the family of powder bed fusion additive manufacturing (PBFAM) © The Author(s) 2021. 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 copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 0123456789().,–: volV Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 2 of 37 processes allows for high geometrical ﬂexibility and the potential for a pointwise control of material properties. However, many of the underlying complex multiphysics phenom- ena are still insuﬃciently 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 [3]. 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 eﬀorts 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 [3,4]). 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 [5,6]. 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 [7–13]or the powder recoating process [14–16]. Microscale models investigate the formation of (typically anisotropic) metallurgical microstructures during solidiﬁcation [17–22]. In this work we propose a coupled thermomechanical macroscale model based on the ﬁnite 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 ﬁnite 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 spe- ciﬁc temperature- and phase-dependent thermal and mechanical properties [23–30]. The heat of the incident laser beam is frequently modeled with the powder bed absorption model from Gusarov et al. [31–33]. Further approaches that are frequently followed for an eﬃcient simulation of PBFAM processes include spatial mesh adaptivity [34–37], some- times combined with code parallelization and load balancing techniques [38], reduction of the computational domain via equivalent thermal boundary conditions representing the powder phase [39], process layer agglomeration [27,37,40], as well as rather heuristic approaches such as inherent strain schemes [41–43]. In macroscale models, all three phases are typically modeled in a Lagrangian manner based on a quasi-solid constitutive law with artiﬁcial stiﬀness values in powder and melt phase that are orders of magnitude lower than the stiﬀness of the solid phase. This approx- imation is well-justiﬁed 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 diﬀerent phases are com- monly modeled as diﬀuse interfaces, e.g., deﬁned 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 stiﬀness), do not transfer to large stresses in the newly formed solid phase (with high stiﬀness), which is assumed to be (approximately) stress-free right after solidiﬁcation. Some contributions [26,42] use a plastic constitutive law, which can achieve the desired eﬀect of an (almost) stress-free conﬁguration at solidiﬁ- cation start [44]. Other works mention a reset or annealing of (plastic) strains and stresses Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 3 of 37 if material (re)melts [25,45,46]. 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 solidiﬁcation start can be incorporated in such models. A notable excep- tion are the recent contributions [47,48], where a constitutive model for the three-phase mixture powder-melt-solid has been consistently derived on basis of iso-stress homog- enization 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 [49]. The stress state of the three-phase mixture powder-melt-solid is con- sistently 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 solidiﬁed 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 solidiﬁcation start while avoiding spatial and temporal discontinuities (jumps) in the stress-ﬁeld, as typically resulting from (instantaneous) stress resetting pro- cedures. 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 [36], element- birth methods [38] or combinations thereof [50], most of which use layerwise conforming FE meshes. Only very few approaches, e.g., based on the immersed boundary method [51], can be found that allow for more ﬂexible 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 eﬃcient condensation of constraint equations from the global system of equations. This procedure oﬀers superior ﬂexibility with regard to spatial dis- cretization by allowing for non-conforming meshes between subsequent layers. As com- pared to collocation-type coupling schemes, e.g., based on hanging nodes, mortar methods enforce the coupling conditions between non-matching meshes in an L -optimal man- ner. Non-conforming interface discretizations between powder layers can signiﬁcantly 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 mod- eling 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 “Mathematical prob- lem statement” presents the underlying mathematical model of thermomechanics. More speciﬁcally, section “Temperature- and phase-dependent parameters” reviews the mod- eling of temperature- and phase-dependent thermal material parameters and section “Mechanical constitutive law” presents the newly proposed solid material model. Its Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 4 of 37 numerical treatment is detailed in section “Numerical solution procedure” along with the general numerical solution procedure. The model is veriﬁed by elementary test cases allowing for analytic solutions in section “One- and two-dimensional numerical exam- ples”. section “Multilayer mesh tying strategy” introduces the mortar mesh tying frame- work for adding new powder layers. Its capabilities are demonstrated in a variety of three- dimensional examples in section “Three-dimensional numerical examples”. Finally, sec- tion “Conclusion” 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 cou- pled with the static balance of linear momentum: c(T) T +∇ · q = rˆ (1) ∇· σ = 0 (2) with the primary variables temperature T and displacement u. The magnitudes of dis- placements and rotations as arising from typical PBFAM processing conditions can be assumed as small, hence the problem is commonly modeled within the theory of geo- metrically linear continuum mechanics where the (engineering) strain tensor is deﬁned by: ε = ∇u + (∇u) (3) The coupling between the two Eqs. (1)and (2) is for now hidden in the structural material law σ = σ(ε(u),T) which will be discussed in detail in section “Mechanical constitutive law”. The heat ﬂux q is speciﬁed by Fourier’s law of heat conduction, q =−k(T)∇T. (4) The material parameters appearing in the thermal problem, namely volumetric heat capac- ity c and heat conductivity k, will in general depend on the temperature and phase. Their modeling is discussed in detail in the authors’ publication [49] and is brieﬂy reviewed in section “Temperature- and phase-dependent parameters”. The source term rˆ is used to model the incident laser beam power based on [33]. Remark Technically, the heat equation (1), which is derived from energy conservation, can contain coupled mechanical terms. The present, geometrically linear problem formu- lation without these coupling terms results in a one-way coupling, i.e., the temperature ﬁeld inﬂuences the structural ﬁeld, but not vice versa. This assumption is ubiquitous in the macroscale PBFAM simulation literature [25–27,45] and seems justiﬁed as strain-rate dependent heating eﬀects are not relevant for a quasi-static solid mechanics problems. A model that considers the two-way coupled thermomechanical problem can be found in [48]. The initial boundary value problem is completed by initial conditions for the temperature ﬁeld and Dirichlet and Neumann boundary conditions for the thermal and structural Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 5 of 37 problem: T = T , in for t = 0, (5) T = T, on , (6) q · n = q, ˆ on , (7) u = uˆ, on , (8) σ · n = t, on (9) where is the problem domain, the Dirichlet and the Neumann boundary of the T q thermal problem, the Dirichlet and the Neumann boundary of the solid problem u σ and quantities (·) 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 eﬀects of latent heat. Essentially, this method modiﬁes the heat capacity c, details on the derivation can again be found in [49]. 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 eﬀects are thus negligible. Consequently, most of the literature focuses on the static structural problem [25,36,46,52]. 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 brieﬂy summarizes the modeling of the three diﬀerent phases powder, melt and solid. For a more detailed motivation and derivation, the reader is referred to [49]. The commonly used liquid fraction g is introduced as 0,T < T T −T g(T) = (10) ,T ≤ T ≤ T T −T ⎪ l 1,T > T where T and T represent the solidus and liquidus temperature. The irreversibility of the s l powder-to-melt transition is captured via the consolidated fraction 1, if r (0) = 1 (i.e., initially consolidated) r (t) = (11) max g(T(t)), if r (0) = 0 (i.e., initially powder) t≤t The resulting fractions of powder (p), melt (m)and solid(s) are computed as r = 1 − r , (12) p c r = g, (13) r = r − g. (14) s c Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 6 of 37 These phase fractions can be used to interpolate arbitrary material parameters: f = r (T)f (T) + r (T)f (T) + r (T)f (T), (15) interp p p m m s s where f is the interpolated parameter and f , f and f are the single phase parameters. interp p s m 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. Mechanical constitutive law 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): σ = r σ with i ∈{p, m, s}. (16) i i 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 ε = ε = ε + ε + ε + ε , (17) i σ ,i p,i T,i ref,i although not all terms will be utilized for each phase. The ﬁrst term on the right-hand side of (17) is the elastic strain ε which induces a stress σ in each phase according to a σ ,i i linear hyper-elastic material σ = C : ε , (18) i i σ ,i where C is the fourth-order constitutive tensor E ν E i i C = λ δ δ + μ (δ δ + δ δ ), λ = , μ = . i i ab cd i ac bd ad bc i i (1 + ν)(1 − 2ν) 2(1 + ν) (19) The artiﬁcial Young’s modulus in powder, E ,and melt, E , will be chosen orders of p m magnitude below the physically consistent value of the solid, E . 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 ε , which are only relevant in the solid phase, p,i 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 (ε = 0). The strains due to thermal expansion ε are assumed equal in all phases and read ε = ε = I α dT = α (T − T )I, (20) T,i T T T ref ref where α is the (constant) coeﬃcient of thermal expansion and T is a reference tem- T ref perature. Finally, the following reference strain ε =: ε , which is only relevant for the ref,s ref solid phase, i.e., ε = ε = 0 , is proposed in rate form: ref,p ref,m Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 7 of 37 (ε − ε − ε )˙r , if r˙ > 0 ⎪ p T s s r˙ ˙ s ε = εˆ , with εˆ = , εˆ (0) = 0, (21) εˆ , if r˙ < 0 ref ref ref ref ref s r ⎪ 0, otherwise where ε represents an accumulation of reference strain contributions weighted by solid ref fractions, which is used as an intermediate variable. The ﬁrst case in (21) refers to a solidifying material point (r˙ > 0) and is motivated by physics as discussed in the next section, while the second case, for a melting material point (r˙ < 0), ensures that the 1 1 reference strain ε does not change during melting, i.e., ε˙ = εˆ − r˙ εˆ = 0 ref ref ref 2 ref s r given r˙ < 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 ; T ], which is beneﬁcial 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: σ = (r C + r C + r C ):(ε − α (T − T )I) − r C : ε (22) p p m m s s T ref s s ref The pre-factor of the ﬁrst term in (22) is equivalent to an average of the single phase material parameters weighted with the phase fractions r 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 approxi- mated by applying a simple elastic constitutive law to these phases, with stiﬀness param- eters that are considerably lower as compared to the solid phase, i.e., E ,E Es .In p m 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 signiﬁcantly stiﬀer 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 signiﬁcant inﬂuence on the resulting residual stresses due to the low stiﬀness of these phases and the deﬁnition 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 deﬁning correspondingly decreased layer thicknesses. Modeling solidiﬁcation 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 com- bined with arbitrary (standard) solid material laws, e.g., elasto-plasticity, we purposefully Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 8 of 37 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 solidiﬁca- tion start is the most signiﬁcant 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 conﬁguration at solidiﬁcation 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 eﬀect for residual stress prediction in such a simpliﬁed model. 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 ε in (21) only change when the solid phase fraction ref increases according to r˙ > 0, i.e., for temperatures T ∈ [T ; T ] in the phase change s s interval and negative temperature rates T < 0. An elastic constitutive law with low stiﬀness values (i.e., E ,E Es) as applied to powder and melt leads to small stresses yet p m 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 solidiﬁcation. For the special case that kinematic ε and thermal strains ε (as well as plastic strains ε )are T p constant during solidiﬁcation, which approximately holds if the phase change interval T − T is suﬃciently small, it can easily be veriﬁed from (17)and (21) that the elastic l s 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 conﬁguration when solidiﬁcation starts. From a stress homogenization point of view, the model assumptions underlying the ref- erence strain formulation state that for each solidifying fraction of material, the reference strain contribution eﬀectively creates its new, stress-free reference conﬁguration, from which strains are calculated. This is illustrated in Fig. 1. Two more involved examples, which illustrate the evolution of the reference strain over multiple repeated melt and solidiﬁcation cycles, can be found in Appendix A. Numerical solution procedure Spatial discretization of the partial diﬀerential 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 δTc(T)T d − ∇δT · q d + δT qˆ d − δT rˆ d = 0, (23) δε : σ d − δu · t d =0(24) σ Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 9 of 37 (a) (b) (c) (d) Fig. 1 Eﬀect of the reference strain contribution during solidiﬁcation in a one-dimensional setting. Diﬀerent phase fractions are distinguished in vertical direction together with their respective reference length, visualized in horizontal direction. At time t, starting from a fully molten phase (a), a fraction of material solidiﬁes (b) and due to the calculation of the reference strain eﬀectively takes on a new reference conﬁguration equal to the current conﬁguration (blue). At time t + t, the newly solidiﬁed fraction takes on the (now changed) current conﬁguration (red) as reference conﬁguration (c), such that the eﬀective 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) where the boundary conditions have already been inserted. Equations (23)and (24)are equivalent to the strong forms (1)and (2) if the solution functions are chosen from the trial 1 1 spaces V ={T ∈ H (): T = T on } and V ={u ∈ H (): u = uˆ on }, and the test T T u u functions are chosen from the weighting spaces W ={δT ∈ H (): δT =0on } and δT T 1 1 W ={δu ∈ H (): δu = 0 on }. Here, H () denotes the Sobolev space of functions δu u with square-integrable ﬁrst derivatives. The solution and test functions are discretized with a Bubnov–Galerkin ansatz in space. T(x,t) = N (x)T (t), δT(x,t) = N (x)δT (t) (25) j j j j j j u(x,t) = N (x)d (t), δu(x,t) = N (x)δd (t) (26) j j j j j j where x is the spatial position and N (x) are the time-independent shape functions, which are the same for all solution and test functions. T and δT are the discrete nodal temper- j j atures and their variations, and d and δd are the discrete nodal displacements and their j j 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 ﬁeld u(x,t) is depend- ing on time via the coupling to the temperature ﬁeld T(x,t), which is a solution to the dynamic heat Eq. (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 Eq. (21) is based on a backward Euler scheme and results in the n+1 following time-discrete form of the fraction-weighted reference strain εˆ : ref ⎪ n n+1 n+1 n+1 n+1 εˆ + r (ε − ε )if r > 0, ⎪ s s ref T n+1 n+1 r n+1 n n+1 ˆ ˆ ˆ ε = ε + ε if r < 0, (27) n+1 s ref ref ref ⎪ s εˆ otherwise, ref Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 10 of 37 n+1 n+1 which leads to the actual reference strains via relation (21)as ε = ε . Thus, the n+1 ref ref n+1 time-discrete total reference strain ε may also be obtained directly, viz. ref 1 n n+1 ⎨ n n+1 n+1 n+1 ε r + r (ε − ε ) , if r > 0 n+1 s s T s n+1 ref ε = (28) ref ⎩ n ε , otherwise ref n+1 Note, that in (28) the case of melting material ( r < 0) no longer needs to be treated separately due to the construction of ε as already discussed when it was ﬁrst introduced. ref 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, ∂σ n+1 n+1 = (r C + r C + r C ) − r C H r (29) p p m m s s s s s n+1 ∂ε ∂r ∂σ ∂r ∂r p m s n+1 n+1 = C + C + C : ε − α (T − T )I p m s T 0 n+1 ∂T ∂T ∂T ∂T − (r C + r C + r C ): α I p p m m s s T ∂r n+1 n+1 n+1 n+1 − C : ε − ε − r α I H r s T T s s ∂T ∂r n n+1 − C : ε H − r (30) ref ∂T where H x is the Heaviside step function ( ) 1, if x > 0 H (x) = (31) 0, otherwise. n+1 All phase fractions and their derivatives are evaluated at T . The ﬁrst two terms in (30) stem from the phase interpolation in the ﬁrst term of (22). The third term represents the contribution from a solidifying increment and is eﬃciently derived by inserting (28) into the total stress (22) and thereby canceling out the solid fraction r . The last term represents the eﬀect of a melting solid. One- and two-dimensional numerical examples This section focuses on the veriﬁcation of the proposed material law with a series of elementary test cases. In all of them the temperature is eﬀectively 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 diﬀerent boundary conditions and temperature proﬁles with increasing complexity. These ﬁndings are reported and discussed in sections “One-dimensional domain: homo- geneous temperature load” and “One-dimensional domain: inhomogeneous temperature load”. Complexity is increased further by applying diﬀerent temperature proﬁles to a two-dimensional domain in section “Two-dimensional domain: interaction of boundary condition and temperature load”. Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 11 of 37 Fig. 2 Relation of investigated one-dimensional veriﬁcation examples to their respective three-dimensional solid mechanics problem Table 1 Material parameters for one-dimensional simulations Parameter Description Value Unit T Solidus temperature 1900 C T Liquidus temperature 2100 C T Reference temperature 0 C E Young’s modulus in solid 1 GPa E Young’s modulus in powder 10 MPa E Young’s modulus in melt 10 MPa −6 −1 α Coeﬃcient of thermal expansion 1 × 10 K Thermal properties (e.g. thermal conductivity) are not necessary for these simulations One-dimensional domain: homogeneous temperature load The ﬁrst 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 func- tion of time but spatially constant, leading to melting and solidiﬁcation of the material, possibly multiple times. Figure 2 illustrates the two considered types of boundary con- ditions 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 condi- tion (Fig. 2 left) prescribing the deformation or a Neumann-type condition (Fig. 2 right) prescribing the traction. The prescribed values can be zero (homogeneous boundary con- ditions) 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: ﬁrst, 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 inho- mogeneous boundary conditions as a third example. The numerical results for all exam- ples are compared to an analytical reference solution, which assumes isothermal melting (T = T = T ) and a zero stiﬀness in powder and melt phase. s l m Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 12 of 37 Dirichlet: displacement u(l) = 0 Neumann: traction t(l) = 0 -3 initial phase powder solid analytical -2 0.5 time t [s] time t [s] Fig. 3 Stress and strain over time for the Dirichlet and Neumann scenario with initially powder or solid material induced by a full melt. Simulation is driven by a homogeneous temperature load with a single peak T > T . Phase change interval in gray peak l Full melt The material is completely molten by heating it to a peak temperature T = 2200 C > T . 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 eﬀective Young’s modulus, E = r E , of the phase mixture. While the Neumann scenario is mostly i i shown for completeness, the Dirichlet scenario is more interesting regarding the diﬀerence 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 stiﬀness) but take on signiﬁcant 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 ﬁnite 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 ﬁnal stress σ = 2 MPa is obtained. This value can ﬁnal 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, ε = ε , 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 , as shown on the right-hand side in Fig. 3. eff. Young's m. E [MPa] stress [MPa] strain [-] temperature T [°C] Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 13 of 37 Dirichlet: displacement u(l) = 0 Neumann: traction t(l) = 0 -3 initial phase powder solid analytical -2 0.5 0 5 10 15 20 0 5 10 15 20 time t [s] time t [s] Fig. 4 Stress and strain over time for the Dirichlet and Neumann scenario with initially powder or solid ˆ ˆ material induced by repeated partial melting with T < T < T andaﬁnalfullmeltwith T > T . Phase s 1 l 2 l change interval in gray Repeated partial melt followed by full melt In a second set of examples the material is partially molten with a peak temperature T inside the phase change interval, viz. T < T = 2000 C < T , and then cooled down to s 1 l T . This cycle is repeated four times and, ﬁnally, the material is fully molten with a peak temperature T = 2200 C > T . The resulting strain and stress evolution is shown in 2 l Fig. 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 ﬁrst peak in the temperature proﬁle is reached. This time, the stress does not reach a near-zero value for the initially solid case because the melting stops at T < T , i.e., there remains a phase fraction of solid 1 l material at T , which exhibits thermal stresses due to the increased temperature level as compared to the stress-free conﬁguration at T . 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 diﬀers at this point (0.5 in powder case and 1 in solid case), multiplication with the reference strain (21) yields the same contribution εˆ . Put diﬀerently: independent of ref the initial phase, the non-molten phase fraction yields a stress of zero, since T = T ;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 speciﬁc 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 mate- rial the stress reaches a near-zero value at T because the solid phase (from the previous partial melt) is completely remolten at this point. For an initially solid material, the situ- eff. Young's m. E [MPa] temperature T [°C] stress [MPa] strain [-] Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 14 of 37 Dirichlet: displacement u(l) = 0.001 mm initial phase powder solid analytical -1 0 5 10 15 20 time t [s] Fig. 5 Stress over time induced by a repeated partial melt and a ﬁnal full melt. Inhomogeneous Dirichlet boundary conditions. Phase change interval in gray ation was already discussed in section “Physical motivation and discussion”: the material model contains no information about diﬀerent solid phases but instead averages the ref- erence 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 simu- lation results in Fig. 4. The analytical reference solution for isothermal melting of initially solid material is also included and again only diﬀers for temperatures in the phase change interval (although not visible in Fig. 4). Finally, the material is fully molten with a peak temperature T = 2200 C > T and 2 l cooled to T . Both an initially solid and powder material yield the same ﬁnal 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 conﬁrms 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 bound- ary conditions. For completeness the partial melt case in section “Repeated partial melt followed by full melt” is repeated with an inhomogeneous Dirichlet-type constraint u(l) = 0.001 mm, which leads to a constant strain ε = 0.001 for the given homoge- neous temperature load. Figure 5 shows the resulting residual stresses, which, for an initially solid material, dif- fers from Fig. 4: the ﬁxed strain induces an initial stress that decreases during the repeated partial melting cycles and, eventually, completely vanishes after the full melting, during which all pre-existing residual stresses are relaxed. Again, the behavior of the solid after full melting is identical to the examples in “Full melt” and sections “Repeated partial melt followed by full melt”. For an initially powder material, the results are almost unaf- fected by the inhomogeneous Dirichlet value, which is consistent with the basic modeling assumption that mechanical stresses resulting from a deformation of the powder phase are negligible. stress [MPa] temperature T [°C] Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 15 of 37 Fig. 6 Schematic of temperature peak function, which will move with speed v in x-direction One-dimensional domain: inhomogeneous temperature load The next scenario focuses on a more complex, transient and inhomogeneous temperature proﬁle and a homogeneous Dirichlet boundary condition u(l) = 0 (see left-hand side in Fig. 2). A moving temperature peak sketched in Fig. 6, that starts outside the initially powder domain (length l = 1 m) and moves through it to the other side, mimics the eﬀect of a moving laser beam in PBFAM. For completeness, the mathematical form of such a temperature proﬁle is stated as w − x, ˆ if 0 ≤ xˆ < w T − T max 0 T(x, t) = T + × , with xˆ = x − vt − x 0 w + x, ˆ if − w < xˆ < 0 0 0, otherwise (32) where the parameters for the present example are chosen as T = 0, T = 2200 C, 0 max w = 0.1m, x =−w and v = 1.0 m/ sec. With these parameters every point in the domain melts and solidiﬁes. Figure 7 shows the temperature proﬁle and the resulting displacement and strain ﬁeld at select time steps. The stress, which is constant over the domain, is plotted over the whole simulation time. In the beginning, the low stiﬀness 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 signiﬁcant values only when the whole bar has a solid fraction r > 0, i.e., when the temperature peak reaches the right end of the bar, where the displacement is ﬁxed to zero. The current example is also used to judge the inﬂuence of the stiﬀness ratio between p,m the solid and powder/melt Young’s moduli by varying the (artiﬁcial) powder/melt stiﬀness E .Figure 8 shows how the ﬁnal stress in the bar (after cooldown to the reference p,m temperature) converges with a reﬁned spatial discretization and for diﬀerent stiﬀness ratios. For a high stiﬀness 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 stiﬀness, i.e. a stiﬀness ratio of inﬁnity. A proof is sketched in Appendix B.3. Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 16 of 37 t=0.20 t=0.60 t=1.10 t=1.25 0.2 -0.2 -3 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1 x [m] x [m] x [m] x [m] 0.15 0.1 0.05 0 0.2 0.4 0.6 0.8 1 1.2 time t [s] Fig. 7 Temperature, displacement and strain at selected time steps (ﬁrst three rows) and stress evolution over time for the inhomogeneous temperature load example. Dashed black lines in the time-stress diagram indicate the selected time steps for the upper diagrams 0.5 0.5 E /E 0.4 elements 0.4 s p,m 20 100 0.3 0.3 50 1000 100 10000 200 100000 0.2 0.2 0.1 0.1 2 4 6 10 10 10 10 20 50 100 200 E /E elements s p,m Fig. 8 Convergence of the ﬁnal stress in a constrained bar subject to a moving peak temperature load. Convergence is inﬂuenced by the number of (ﬁnite) elements and the stiﬀness ratio between solid and powder/melt. Both diagrams contain the same information in diﬀerent representations Two-dimensional domain: interaction of boundary condition and temperature load The ﬁnal veriﬁcation example investigates a two-dimensional setting under the assump- tion of plane-stress. 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 Fig. 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 = 303 K to a maximum temperature T = 2000 K > T and then 0 max l cooled down to the initial temperature. This case is compared to the scenario of a moving temperature peak proﬁle as in (32) with diﬀerent peak widths w ∈{0.4, 0.2, 0.1}Kand −1 T = 2000 K. The peak still moves only in x-direction with speed v = 1.0mm s , i.e., max A plane-strain setting, which was also simulated by the authors, would change the absolute values but not the qualitative trends presented in this section. stress [MPa] strain [-] disp. u [mm] temp. T [°C] stress [MPa] stress [MPa] Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 17 of 37 Fig. 9 Geometry and boundary conditions of two-dimensional domain: bottom edge ﬁxed; left, right and top edge free. Dimensions in mm (not to scale) Table 2 Material and simulation parameters for three-dimensional examples Quantity Description Value Unit E Young’s modulus in solid 200 GPa E ,E (Artiﬁcial) Young’s modulus in powder and melt 2 GPa p m ν Poisson ratio 0.3 – −6 −1 α Coeﬃcient of thermal expansion 15 × 10 K −3 −1 c Volumetric speciﬁc heat, solid 4.25 M J m K −3 −1 c Volumetric speciﬁc heat, powder 2.98 M J m K −3 −1 c Volumetric speciﬁc heat, melt 5.95 M J m K −1 −1 k ,k Conductivity in solid and melt 20 W m K s m −1 −1 k Conductivity in powder 0.2@200, 0.3@1600 W m K @K −3 h Latent heat of fusion 2.18 GJ m T Initial/reference temperature 303 K T Solidus temperature 1500 K T Liquidus temperature 1900 K ρ Hemispherical reﬂectivity 0.7 – −1 β Extinction coeﬃcient 60 mm h Powder layer thickness 50 µm W Eﬀective laser power 30 W R Eﬀective laser radius 0.08 mm v Laser scan speed 100 mm/ sec scan the temperature is considered constant in z-direction, which is a rough approximation to the temperature proﬁle encountered in a PBFAM application. Figure 10 compares the resulting residual stresses σ in the diﬀerent scenarios: a uniform temperature proﬁle xx yields the highest tensile stresses and almost no variation of the stress σ in z-direction xx (the “powder layer height”) while a moving temperature peak proﬁle leads to a signiﬁcant drop-oﬀ of σ over the height, meaning that stresses are lower at the non-constrained xx upper surface. This drop-oﬀ becomes more pronounced with a decreased peak width w and can be more clearly observed in Fig. 11 (left). These observations can be explained by the size of the domain which has a temperature above T 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. Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 18 of 37 Fig. 10 Two-dimensional example: residual stress distribution σ after subjecting domain to a uniform (ﬁrst xx row) or moving (second to fourth row) peak temperature proﬁle uniform heating uniform heating z-var w=0.4 moving w=0.4 4000 z-var w=0.2 moving w=0.2 z-var w=0.1 moving w=0.1 z-var w=0.01 1000 3800 00.02 0.04 0 0.02 0.04 z [mm] z [mm] Fig. 11 Two-dimensional example: ﬁnal stress σ over height z (zero is constrained bottom) at x = 0.5mm xx (see dashed line in Fig. 10): induced by a temperature peak moving in x-direction, constant in z-direction (left); induced by a temperature proﬁle varying in z-direction, constant in x-direction (right). Note the diﬀerent ranges of the ordinates Next, the temperature proﬁle shall vary in z-direction as well. Given a time and x- dependent function T (x, t) for the temperature proﬁle at z = 0 the following function z0 z(s − 1) + h T(x, z, t) = T (x, t), (33) z0 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 T − T h max 0 s = 1 + , (34) T − T w m 0 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 2t, if t ≤ 0.5s T (t) = T × . (35) z0 max 2 − 2t, if t > 0.5s Stress [MPa] xx Stress [MPa] xx Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 19 of 37 Fig. 12 Two-dimensional example. Left: ﬁnal stress σ over height z (zero is constrained bottom) at xx x = 0.5 mm depending on width of temperature peak and (possibly) z-variation. Solid lines represent moving temperature proﬁles, dashed lines represent moving and z-varying temperature proﬁles. Right: snapshots of temperature proﬁles a–f which induce the associated ﬁnal stress response on the left. Solidus and liquidus isotherms indicated in red The scaling factor s is computed from w accordingto(34). The results in Fig. 11 (right) indicate that a variation only in z-direction does not lead to signiﬁcantly diﬀerent stresses (compared to the uniform temperature ﬁeld), 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 proﬁles expected in PBFAM processes. This is easily achieved by inserting (32)as T into (33)toobtain z0 another kind of temperature proﬁle, see Fig. 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 inﬂuence 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 proﬁle might play a more important role than often assumed. In the authors’ opinion, the signiﬁcant diﬀerences in residual stress from a uniform and a moving temperature ﬁeld raise the question how well-suited common simpliﬁcations in PBFAM simulation, such as layer-agglomeration, inherent-strain or ﬂashing 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, suﬀer from the need to generate matching meshes between layers for real application scenarios and complex geometries. In this work, we want to Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 20 of 37 Fig. 13 Mesh tying concept for multilayer PBFAM simulations propose a ﬂexible, 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 diﬀerences 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 -norm) and that leads to an eﬃcient 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 is found in [53,54]. Figure 13 shows an exemplary part geometry which is sliced () into layers with thickness equal to the powder thickness h . Each layer can be meshed separately and the ﬁnite element nodes may be non-matching at the interface between two layers. The continuity of the primary solution variables across the interface is mt, enforced in a weak sense by additional (mesh tying) constraint equations. The -th mesh (−1) (−1) tying interface consists of the boundaries (top surface of domain )and mt, mt, () () (bottom surface of domain ). In this framework, applying a new powder layer mt, boils down to activating a new mesh tying constraint. In order to have a continuous part, degrees of freedom (DOFs) associated with nodes () on the bottom surface of powder domain and on the top surface of powder domain (−1) need to be constrained. Displacements and temperatures are the primary solution variables of the thermomechanical problem (1)and (2) and thus the necessary constraint equations read, () (−1) u = u , on , (36) mt, () (−1) T = T , on . (37) mt, These conditions are enforced by introducing a Lagrange multiplier potential for each constraint: () (−1) = λ · u − u d , (38) u,mt u mt, () (−1) = λ · T − T d . (39) T,mt T mt, The total variation of the two potentials is found as () (−1) () (−1) δ = δλ · u − u d + λ · δu − δu d (40) u,mt u mt, u mt, Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 21 of 37 () (−1)) () (−1) δ = δλ · T − T d + λ · δT − δT d (41) T,mt T mt, T mt, The space and time dependency of u, T, their variations δu and δT, and the Lagrange multipliers λ and λ is not written out in Eqs. (36)–(41) for brevity. The Lagrange u T multipliers λ and their variations δλ are chosen from M , the dual space of the trace −1/2 space of V , i.e., M = H ( ), where ∈{u,T }. The ﬁrst term in (40)and mt (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 () (−1) the interface, λ = t =−t , and the thermal Lagrange multiplier as a heat ﬂux mt mt () (−1) across the interface, λ = q =−q . The Lagrange multiplier ﬁelds are discretized mt mt with special dual shape functions that allow for an eﬃcient condensation of the Lagrange multiplier DOFs from the global system of equations [54]. After spatial discretization of all primary variable ﬁelds, the consistent link between discrete nodal DOFs on both sides of the interface is established via constant mortar matrices D and M [53], which can −1 be combined into constant projector matrices P = D M . The projector matrices map discrete solution increments (from the nonlinear solution procedure), () (−1) d = P d , on (42) u mt, () (−1) T = P T , on , (43) T mt, where d and T are increments of the discrete solution vectors of discrete displace- ments 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 conﬁgurations. 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 of mt, (−1) a newly added layer is to apply mappings (42)and (43) to the solution vectors on , mt, viz. () (−1) () (−1) d ← P d , T ← P T on . (44) u T mt, 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 () () \ are set as follows: temperature DOFs are set to the initial temperature T , mt, while displacement DOFs are initialized to zero. The ﬁrst time step after a new mesh tying interface is activated will solve the (static) mechanical equilibrium Eq. (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 are modiﬁed without an mt, associated external force. An alternative initialization of the temperature ﬁeld for element- birth methods is discussed in [50] and could be transferred to the proposed mesh tying Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 22 of 37 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 modiﬁcation vanishes. The error introduced by modiﬁcation of the displacements of a newly added layer is generally small due to the low stiﬀness of the newly added powder, i.e., possible artiﬁcial strains in the new layer only cause a small stress response. Furthermore, the material formulation from section “Mechanical constitutive law” ensures that these artiﬁcial stresses vanish after melting. Three-dimensional numerical examples This section presents larger three-dimensional numerical examples that are simulated with the research code BACI [55], jointly developed at the Institute for Computational Mechanics. Single tracks per layer The ﬁrst 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 Fig. 14. The computational domain consists of a 0.1 mm high base plate that has a prescribed temperature T = 303 K and zero displacements u = 0 at the bottom (indicated in blue in Fig. 14). This temperature T is also used as the initial temperature throughout the domain and in all newly activated powder layers. The powder layers of height h 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 stiﬀness and conductivity, the zero displacement boundary conditions will not inﬂuence the results in the ﬁnal part geometry signiﬁcantly. 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 in good approximation) during which signiﬁcant 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: The boundary conditions on the red surfaces are removed everywhere except on the leftmost surfaces (with normal vector (−1, 0, 0) ) to avoid rigid body modes. The mesh tying constraint between the base plate and the ﬁrst 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 −6 parameters are listed in Table 2. A time step of t = 5 × 10 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 Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 23 of 37 Fig. 14 Geometry, boundary conditions and simulation setup of ten layer example. Dimensions in mm −3 atimestepof t = 1 × 10 s. The whole domain is discretized with 61 440 linear, ele hexahedral 8-node ﬁnite elements with edge size h = = 12.5 µm. Figure 15 shows an overview of the six components of the stress tensor at selected points in time: the ﬁrst 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 ﬁfth 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 σ component (ﬁrst row). In agreement with the literature xx this stress in scanning direction is a large tensile stress [45,56]. The highest stresses can be observed in the heavily constrained base plate. Speciﬁcally, it can be observed how geometrical compatibility during cool down lead to high tensile normal stresses σ in xx the ﬁrst track (and remolten portion of the base plate) and high compressive normal stresses σ in the (non-molten) lower part of the build plate. In the last row of Fig. 15 the xx shear stresses σ can be observed, which are—according to the mechanical equilibrium xz Eq. (2)— necessary to balance gradients of the normal stress σ at the beginning and xx end of a track, especially close to the base plate. These σ shear stresses, in turn, also xz induce σ normal stresses. Note, that the diﬀerent stress components at many locations zz exceed the yield stress of typical materials used in PBFAM, hence the inclusion of an elasto-plastic material model, as already sketched in section “Mathematical formulation”, 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. Fig- ure 16 shows the resulting stress distribution for an alternating, serpentine scanning pat- tern. The absolute values are comparable to Fig. 15. To better highlight the inﬂuence of the diﬀerent scanning patterns, Fig. 17 visualizes the normal stresses σ for the unidirectional xx and serpentine scanning pattern with a rescaled color map. Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 24 of 37 Fig. 15 Unidirectional scanning pattern: distribution of diﬀerent 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 Fig. 16 Serpentine scanning pattern: distribution of diﬀerent 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 Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 25 of 37 Fig. 17 Stress distribution of σ in symmetry plane at selected time steps for unidirectional (ﬁrst row) and xx serpentine (second row) scanning pattern. Blue color represents tension, red color represents compression. Liquidus isotherm (solid black line) approximates melt pool shape. For improved visibility of the relevant stresses the color bar is rescaled Selected time steps are plotted in Fig. 18 for a more detailed view of the stress evolution of σ during scanning of track 7 with a unidirectional pattern. The following considerations xx 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 ﬁrst and second row of Fig. 18 show the processing of track 7. To visualize the melt pool size, the temperature isoline corresponding to the liquidus temperature T 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 -isoline. As desired, these stresses rapidly drop to zero in the narrow temperature interval T ∈ [T ; T ] such that no visible stresses remain in the s l melt pool domain. This strong gradient between vanishing stresses in the melt pool and high compressive stresses in the solid material beneath remains after solidiﬁcation 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 Fig. 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 stiﬀness as the solidiﬁed tracks above, this characteristic band structure of the normal stresses is much more pronounced for the ﬁrst track, i.e., the highest overall tensile stresses occur in the ﬁrst 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 solidiﬁed melt tracks are subject to tensile stresses after cooling down (e.g., snapshot at bottom right of Fig. 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 eﬀects 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 modiﬁed scenario, where the base plate is scanned ﬁrst with the same heat source as the powder layers. Figure 19 shows the evolution of stress σ . xx Large tensile stresses are visible in the scanned regions of the base plate after cooldown Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 26 of 37 Fig. 18 Unidirectional scanning pattern: detailed view of the evolution of the longitudinal stress σ in the xx symmetry plane during scanning of layer 7. Blue color represents tension, red color represents compression. Liquidus isotherm (solid black line) approximates melt pool shape. Time progresses from left to right and top to bottom. For improved visibility of the relevant stresses the color bar is rescaled Fig. 19 Modiﬁed scenario: scanning base plate and layer 1. Detailed view of the evolution of the longitudinal stress σ in the symmetry plane. Blue color represents tension, red color represents compression. Liquidus xx isotherm (solid black line) approximates melt pool shape. Time progresses from left to right and top to bottom. For improved visibility of the relevant stresses the color bar is rescaled Fig. 20 Uniform scanning pattern: detailed view of the evolution of the longitudinal stress σ in the xx symmetry plane during scanning of layer 1. Blue color represents tension, red color represents compression. Liquidus isotherm (solid black line) approximates melt pool shape. Time progresses from left to right. For improved visibility of the relevant stresses the color bar is rescaled (top right in Fig. 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 Fig. 18, shifts down in Fig. 19 while the magnitude stays roughly the same. For comparison, Fig. 20 shows a detailed view of the evolution of stress σ when the ﬁrst layer xx 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 Fig. 20) is equal to the stress distribution in layer 1 after cooldown with a pre-processed base plate (bottom right in Fig. 19), which demonstrates that the chosen pre-processing of the base plate does not inﬂuence the stress in layer 1. Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 27 of 37 Fig. 21 Geometry, boundary conditions and simulation setup for two layer example with multiple tracks per layer (indicated in green for layer 1 and blue for layer 2). Dimensions in mm 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 Fig. 21: two powder layers are deposited on top of a base plate and ﬁve 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 −6 are reused. As before, a time step of t = 5 × 10 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 −3 of t = 1 × 10 s. The domain is discretized with 15 3600 linear, hexahedral 8-node ele ﬁnite elements with edge size h = = 12.5 µm. In this example, we look at the stress distribution after each layer was processed. Fig- ure 22 visualizes the stresses σ and σ in the consolidated material, which represents xx yy the built part. One quarter of the consolidated volume is cut out for visualization of the stress distribution in z-direction. In the ﬁrst layer, the laser travels in x-direction, which is clearly reﬂected in the stress distributions shown in the ﬁrst column of Fig. 22. When the second layer is applied and scanned in y-direction, the melt pool penetrates into the ﬁrst layer. Therefore, after the second layer was scanned, the stress distribution in the ﬁrst layer reﬂects the laser beam movement in y-direction, see the second column in Fig. 22). At the same time, stresses in the second layer also align with the y direction, see the third column in Fig. 22. Overall, only tensile stresses remain after the ﬁnal cool down. Two pyramids example The ﬁnal numerical example investigates a more complex geometry. Yet, the geometry is deﬁned simple enough so that it can easily be adopted by other researchers. The powder material that is not supposed to contribute to the ﬁnal part is not included, an assumption which is justiﬁed by the low thermal conductivity and stiﬀness of the powder. The ﬁnal part’s cross section changes drastically over the layers which demonstrates the ﬂexibility of the mortar mesh tying approach. The geometry, see Fig. 23, consists of two pyramid Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 28 of 37 stress in layer 1 stress in layer 1 stress in layer 2 after scan of layer 1 after scan of layer 2 after scan of layer 2 xx yy 0 MPa 3100 MPa Fig. 22 Multiple tracks per layer: stress distribution on consolidated part geometry after scanning layer 1 (ﬁrst column) and layer 2 (second and third column). A quarter of the geometry is cut out for visualization purposes bodies on a build plate: the ﬁrst pyramid has curved surfaces and consists of bulk material, while the second pyramid is hollow. These geometries are chosen to illustrate diﬀerent behavior in bulk geometries versus slender geometries. In each layer, ﬁrst 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 ﬁrst 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 , in analogy to the tracks depicted in Fig. 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 ﬁxed to the initial temperature T .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 signiﬁcant convective surface heat transfer would occur in the relevant time scales of this problem. −6 Again, a time step of t = 5×10 s is used for the scanning phase and the initial cooling −3 phase of 0.002 s, while the remaining 0.998 s of cooling use a time step of t = 1 × 10 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 “Multilayer mesh tying strategy”. Figure 24 depicts the normal stress components as well as the von Mises equivalent stress 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 stiﬀ hollow pyramid compared to the bulk pyramid. The curvature follows a circular arc with radius R = 33 2/80 mm, which is tangent to the z-axis at the pyramid top. The von Mises equivalent stress is calculated as 2 2 2 2 2 2 σ = σ + σ + σ − σ σ − σ σ − σ σ + 3(σ + σ + σ ). vM xx yy zz xx yy xx zz yy zz xy xz yz Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 29 of 37 Fig. 23 Geometry of pyramids example. Dimensions in mm The stress distribution σ in vertical direction exposes a zone of compressive stresses zz inside the bulk pyramid and tensile stresses in the outer region, which is in agreement with the results in section “Single tracks per layer” and results in literature [27,42]. Note, that parts of the material close to the edges in the ﬁrst layers were only partially molten, hence the artifacts in the stress distributions. The observations for stresses are complemented by Fig. 25: here, the displacement magnitude is visualized on a deformed mesh (with deformation scaled by a factor of 5 for improved visibility). In the stiﬀer 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. Of course, this example could also be scaled up to more realistic part dimensions when using a layer-agglomeration approach [27] or eﬀective heat track- or layerwise heat sources [37]. Conclusion A thermomechanically coupled computational model for the macroscale simulation of PBFAM processes was presented. A ﬁrst emphasis of the present work was the consis- Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 30 of 37 (a) (b) (c) (d) Fig. 24 Two pyramids example: stress distribution in pyramids after ﬁnal cooldown. Normal stresses a σ , b xx σ , c σ ,and d von Mises equivalent stress σ . Part of the geometry is cut out above ﬁrst powder layer and yy zz vM at y = 0.6 mm for visualization of stress inside the pyramid Fig. 25 Two pyramids example: displacement magnitude after ﬁnal cooldown displayed on warped geometry. The distortion is scaled up by a factor of 5 tent 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 veriﬁed with academic test cases. Crucially, this inelastic reference strain contribution, formulated in rate form, was introduced to allow for a stress-free state when solidiﬁcation of molten material starts. Compared to existing (instantaneous) stress-reset procedures, which might result in jumps in the stress-ﬁeld, the formulation of the reference strain term in rate form guarantees a smooth stress ﬁeld. Compared to more involved material models, e.g. [45,48], this simple approach can easily be integrated in standard material libraries Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 31 of 37 and FEM codes while still consistently capturing the most essential aspects of residual stress generation in PBFAM. This capability has been veriﬁed 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 ﬂexibility 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 eﬀectiveness of the proposed mortar mesh tying strategy. Acknowledgements This work was supported by the German Research Foundation (DFG) and the Technical University of Munich (TUM) in the framework of the Open Access Publishing Program. Authors’ contributions All authors contributed to the derivation of model equations, the discussion of results, and writing the manuscript. In addition, SP performed the speciﬁc code implementation and the shown numerical simulations. CM and WAW worked out the general conception of the proposed modeling approaches. All authors read and approved the ﬁnal manuscript. Funding Open Access funding enabled and organized by Projekt DEAL. This work was supported by funding of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within project 437616465 and project 414180263. Availability of data and materials The research code, numerical results and digital data obtained in this project are held on deployed servers that are backed up daily. Fundamental work results are also stored on a NAS server managed by the Leibniz Rechenzentrum (LRZ) in Garching. The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. Declarations Competing interests The authors declare that they have no competing interests. Appendix A: Behavior of the rate-based reference strain formulation For further insight into the behavior of the reference strain term, the proposed rate-based reference strain formulation (21) is integrated over the phase change interval to obtain rem rem ε (r ) = r ε + (ε − ε )H(˙r )dr , (45) ref s T s ref r rem s r rem with a remaining solid fraction r at solidiﬁcation start, which has a reference strain rem ε from a previous solidiﬁcation process, as initial value for the integration. Here, H (x) ref is the Heaviside step function (31). A.1 Example scenario 1 Table 3 lists a temperature history and the corresponding reference strain evolution com- puted with the integral form (45) under the assumption that the integrand ε − ε stays constant during solidiﬁcation, which is a good approximation for suﬃciently small phase change intervals T − T .At t the material consists solely of powder phase and the refer- l s ence strain is zero. After heating to 0.5(T +T ) (50% of material molten) at t , the material l s Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 32 of 37 Table 3 Exemplary evolution of reference strain for repeated cooling and heating assuming constant strains during solidiﬁcation Time Temperature Reference strain r r r p s m 0 0 t T < T ε = 0 1.0 0.0 0.0 ref 1 1 t T = 0.5(T + T ) ε = 0 0.5 0.0 0.5 l s ref 2 2 1 1 1 t T < T ε = 0.5(ε − ε ) 0.5 0.5 0.0 ref 0.5 T 3 3 1 2 2 t T = 0.25T + 0.75T ε = 0.25ε = ε 0.5 0.25 0.25 l s ref 0.25 ref ref 4 4 1 3 3 3 1 3 1 3 3 t T < T ε = (0.25ε + 0.25(ε − ε )) = ε + (ε − ε ) 0.5 0.5 0.0 ref 0.5 ref T 2 ref 2 T 5 5 r 4 4 t T = 0.75T + 0.25T ε = lim ε = ε 0.25 0.0 0.75 s r →0 l s ref r ref ref 6 6 1 5 5 5 5 5 t T < T ε = (0 · ε + 0.75(ε − ε )) = ε − ε 0.25 0.75 0.0 ref 0.75 ref T T Table 4 Exemplary evolution of reference strain for a cyclic, partial remelting of solid phase assuming ε − ε is constant in all phase change intervals time temperature reference strain r r s m 0 0 t T < T ε =01.00.0 ref 1 1 t T = 0.5(T + T ) ε =00.50.5 l s ref 2 2 1 1 1 t T < T ε = (0.5ε + 0.5(ε − ε )) = (ε − ε)1.0 0.0 s T T ref 1 ref 2 3 3 1 2 2 t T = 0.5(T + T ) ε = 0.5ε = ε 0.5 0.5 l s ref 0.5 ref ref 4 4 1 3 3 t T < T ε = (0.5ε + 0.5(ε − ε )) = (ε − ε)1.0 0.0 s T T ref 1 ref 4 5 5 1 4 4 t T = 0.5(T + T ) ε = 0.5ε = ε 0.5 0.5 l s ref 0.5 ref ref 6 6 1 5 7 t T < T ε = (0.5ε + 0.5(ε − ε )) = (ε − ε)1.0 0.0 s T T ref 1 ref 8 2 2 cools to a temperature below T at t . A reference strain ε evolves during this solidiﬁca- ref 1 1 tion process based on the strains ε − ε . Another heating up to 0.25T + 0.75T (25% of material molten) does not change the reference strain. However, note that the last term in the stress (22) changes continuously during melting due to the multiplication with the decreasing solid fraction (which melts before powder can melt). Not all previously created solid melts and all three phases are present at t . The next cooldown to a temperature 4 4 below solidus temperature at t leads to ε , computed as a weighted average of the old ref 3 3 3 reference strain ε in the remaining solid and the current strains ε − ε in the melt. As ref T expected, the powder material does not contribute to the reference strain. A ﬁnal heating to 0.75T + 0.25T (75% of material molten) at t melts all of the previously created solid l s phase and, additionally, some more of the powder phase. Thus, the subsequent cooling at 6 6 5 5 t leads to a reference strain ε , which only depends on the current strains ε − ε ,as ref T no solid fraction (with a reference strain contribution) remained at t . A.2 Example scenario 2 The ﬁnal special case looks at cyclic, partial remelting of an initially solid material point as outlined in Table 4. This case is also investigated as a numerical example in section “Repeated partial melt followed by full melt”. The material is repeatedly heated up to a T +T s l maximum temperature T = (i.e., 50% of the material melts) and afterwards cooled to a temperature below T . In the ﬁrst cycle, half of the initial solid will melt (at t )and obtain a new reference strain after cooldown (at t ), where the strain ε − ε is assumed (approximately) constant during the complete solidiﬁcation 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 ε = 0.5(ε − ε ). In the next cycle, again 50% of the solid phase melts ref 3 4 (at t ). After cooldown below T (at t ) the new reference strain is computed according to (45) as a weighted average of the reference strain of the non-molten solid (at t ) and the Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 33 of 37 4 3 new strain contribution ε − ε for the remolten fraction, ε = 0.5ε + 0.5(ε − ε ) = T T ref ref 5 6 0.75(ε − ε ). The same logic applies to the next heating at t and cooling t 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 ε − ε . If this process is repeated many times, this has the (at ﬁrst glance) paradoxical result that, although never fully molten, the complete initial solid will asymptotically convert to newly solidiﬁed 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 suﬃcient partial melting cycles each solid material fraction has molten at least once such that the total solid phase eﬀectively exhibits a new stress-free reference state. Appendix B: Analytical reference solutions for one-dimensional problem This section presents analytical solutions for veriﬁcation of the one-dimensional scenario described in section “One- and two-dimensional numerical examples”. B.1 Final stress after full melt and solidiﬁcation The homogeneous Dirichlet boundary condition and homogeneous temperature load used in sections “Full melt” and “Repeated partial melt followed by full melt” imply ε ≡ 0. Thus, the ﬁnal 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: 1 T T − T σ =−E ε = E α (T − T )dr =−E α dT ﬁnal s ref s T 0 s s T T − T 0 T l s T + T = E α − T (46) s T 0 where the average of solidus and liquidus temperature can be interpreted as the melting T +T s l temperature T := . Result (46) is equivalent to the expected solution for isothermal phase change at melting temperature T and the exact value of the (artiﬁcial) powder and melt stiﬀness is irrelevant for the ﬁnal stress value in this speciﬁc case. B.2 Stress after partial melt and solidiﬁcation The computation in (46) can be generalized to ﬁnd the ﬁnal stress after a partial melt with a peak temperature T, lying in the phase change interval [T ; T ], and subsequent cooldown: s ˆ ˆ T − T T − T T + T 0 s s σ =−E ε =−E α dT = − T . (47) ﬁnal s ref s T 0 ˆ T − T T − T 2 T l s l s =g(T) B.3 Final stress for moving temperature peak The following derivations consider the example in section “One-dimensional domain: inhomogeneous temperature load”, where the space and time-dependent, peak-shaped temperature proﬁle T(x, t) given in (32) is prescribed on the whole domain. We assume zero stiﬀness for powder and melt phase and a ﬁnite value for the solid phase. Let t be the Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 34 of 37 instant of time, when, for the ﬁrst time, there is a solid fraction r > 0 everywhere in the bar. No stress occurs before and at this time (due to zero stiﬀness in powder and melt) and consequently σ (t) = 0. With these simpliﬁcations we ﬁnd from (22) σ (t) = r (x, t)E (ε(x, t) − α (T(x, t) − T ) − ε (x, t)) , ∀t ≥ t (48) s s T 0 ref where the spatial dependency on x is left out for the stress. Due to geometric compatibility the displacement on both ends of the domain must vanish, and therefore 0 = ε(x, t)dx l l l dx = σ (t) + α T(x, t) − T dx + ε (x, t)dx. (49) ( ) T 0 ref r (x, t)E 0 s s 0 0 ˆ ˆ Since σ (t) = 0 as stated above, evaluating (49)at t yields l l ˆ ˆ 0 = α T(x, t) − T dx + ε (x, t)dx, (50) T 0 ref 0 0 and subtracting (50)from(49) gives l l dx 0 = σ (t) + α T(x, t) − T(x, t) dx r (x, t)E s s 0 0 + ε (x, t) − ε (x, t)dx. (51) ref ref Let t be the ﬁnal time, where all material is solid, r (x, t ) = 1, and the temperature is f f equal to the initial temperature, T(x, t ) = T . Then the ﬁrst term in (51)istrivial to integrate and, after some rearrangement, one ﬁnds for the ﬁnal stress: l l E E s s ˆ ˆ σ = σ t = α T(x, t) − T dx + ε (x, t) − ε (x, t )dx . (52) ﬁnal f T 0 ref ref f l l 0 0 I I 1 2 The ﬁrst term depends on the integrated diﬀerence in temperature when solidiﬁcation starts (at time t) compared to the ﬁnal temperature proﬁle (at time t ). For the temperature proﬁle (32), it can be computed in analogy to (46)and (47)as: α (T − T ) T l 0 I = w . (53) 2 T − T max 0 Result (53) demonstrates that the stress after cooldown depends on the shape of the temperature proﬁle (here described by w and T ) and not just the temperature diﬀerence max between liquidus temperature T and ﬁnal temperature T . l f To compute the second term in (52), the deﬁnition of the reference strain (21) is inserted: I = ε (x, t) −ε (x, t )dx ref ref f =0, since for t≤t: ε=α T l t ˜ ˜ ˜ ˜ =− H(˙r ) ε(x, t) − α T(x, t) r˙ (t)dt dx (54) s T s r (x, t ) ˆ 0 f t Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 35 of 37 ˜ ˆ Assume that the strain stays constant during solidiﬁcation, i.e., ε(x, t) = ε(x, t) = α T(x, t), which is a good approximation for a small phase change interval [T ; T ]. T s l The ﬁnal solid fraction is again given as r (x, t ) = 1 and its rate during solidiﬁcation as s f r˙ = T /(T − T ). Thus, (54) simpliﬁes to: s s l t ˆ ˜ ˜ I =− H(˙r )α (T(x, t) − T(x, t)) dt dx (55) 2 s T T − T 0 t s l ˙ ˜ ˜ ˆ After a change of variables, H(˙r )T dt → H T(x, t) − T dT, the inner integral can be s s computed as: l T s ˜ T − T(x, t) I =−α H T(x, t) − T dT dx 2 T s ˆ T − T 0 T(x,t) l s l 2 α (T − T(x, t)) T s =− H T(x, t) − T dx (56) 2 T − T 0 l s Finally, the speciﬁc temperature proﬁle (32) is inserted into (56) which allows to compute the ﬁnal integral, yielding α (T − T ) T s I =− w , (57) 6 T − T max 0 which is a signiﬁcantly smaller value compared to I for typical values encountered in PBFAM application, e.g. for the parameters used in section “One-dimensional domain: inhomogeneous temperature load” a ratio of |I /I |≈ 0.3% is obtained. 2 1 Received: 30 July 2021 Accepted: 28 September 2021 References 1. Gibson I, Rosen D, Stucker B, Khorasani M. Additive manufacturing technologies, vol. 17. New York: Springer; 2014. 2. Herzog D, Seyda V, Wycisk E, Emmelmann C. Additive manufacturing of metals. Acta Mater. 2016;117:371–92. 3. Meier C, Penny RW, Zou Y, Gibbs JS, Hart AJ. Thermophysical phenomena in metal additive manufacturing by selective laser melting: fundamentals, modeling, simulation, and experimentation. Annu Rev Heat Transf. 2017;20:241–316. 4. Meier C, Fuchs SL, Much N, Nitzler J, Penny RW, Praegla PM, Proell SD, Sun Y, Weissbach R, Schreter M, Hodge NE, John Hart A, Wall WA. Physics-based modeling and predictive simulation of powder bed fusion additive manufacturing across length scales. Surv Appl Math Mech GAMM-Mitteilungen. 2021. https://doi.org/10.1002/gamm.202100014. 5. Li C, Liu Z, Fang X, Guo Y. Residual stress in metal additive manufacturing. Procedia Cirp. 2018;71:348–53. 6. Mercelis P, Kruth J-P. Residual stresses in selective laser sintering and selective laser melting. Rapid Prototyping J. 2006;12(5):254–65. 7. Khairallah SA, Anderson A. Mesoscopic simulation model of selective laser melting of stainless steel powder. J Mater Process Technol. 2014;214(11):2627–36. 8. Körner C, Attar E, Heinl P. Mesoscopic simulation of selective beam melting processes. J Mater Process Technol. 2011;211(6):978–87. 9. Markl M, Ammer R, Rüde U, Körner C. Numerical investigations on hatching process strategies for powder-bed-based additive manufacturing using an electron beam. Int J Adv Manuf Technol. 2015;78(1–4):239–47. 10. Meier C, Fuchs SL, Hart AJ, Wall WA. A novel smoothed particle hydrodynamics formulation for thermo-capillary phase change problems with focus on metal additive manufacturing melt pool modeling. Comput Methods Appl Mech Eng. 2021;381: 113812. 11. Russell M, Souto-Iglesias A, Zohdi T. Numerical simulation of laser fusion additive manufacturing processes using the sph method. Comput Methods Appl Mech Eng. 2018;341:163–87. 12. Wessels H, Weißenfels C, Wriggers P. Metal particle fusion analysis for additive manufacturing using the stabilized optimal transportation meshfree method. Comput Methods Appl Mech Eng. 2018;339:91–114. 13. Yan W, Qian Y, Ge W, Lin S, Liu WK, Lin F, Wagner GJ. Meso-scale modeling of multiple-layer fabrication process in selective electron beam melting: inter-layer/track voids formation. Mater Des. 2018;141:210–9. 14. Herbold E, Walton O, Homel M. Simulation of powder layer deposition in additive manufacturing processes using the discrete element method. Technical report, Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States); 15. Meier C, Weissbach R, Weinberg J, Wall WA, Hart AJ. Critical inﬂuences of particle size and adhesion on the powder layer uniformity in metal additive manufacturing. J Mater Process Technol. 2019;266:484–501. Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 36 of 37 16. Meier C, Weissbach R, Weinberg J, Wall WA, Hart AJ. Modeling and characterization of cohesion in ﬁne metal powders with a focus on additive manufacturing process simulations. Powder Technol. 2019;343:855–66. 17. Gong X, Chou K. Phase-ﬁeld modeling of microstructure evolution in electron beam additive manufacturing. JOM. 2015;67(5):1176–82. 18. Lindgren L-E, Lundbäck A, Fisk M, Pederson R, Andersson J. Simulation of additive manufacturing using coupled constitutive and microstructure models. Addit Manuf. 2016;12:144–58. 19. Nitzler J, Meier C, Müller KW, Wall WA, Hodge NE. A novel physics-based and data-supported microstructure model for part-scale simulation of laser powder bed fusion of Ti-6Al-4V. Adv Model Simul Eng Sci. 2021;8(16):1–39. 20. Rai A, Markl M, Körner C. A coupled cellular automaton-lattice Boltzmann model for grain structure simulation during additive manufacturing. Comput Mater Sci. 2016;124:37–48. 21. Salsi E, Chiumenti M, Cervera M. Modeling of microstructure evolution of Ti6Al4V for additive manufacturing. Metals. 2018;8(8):633. 22. Zhang J, Liou F, Seufzer W, Newkirk J, Fan Z, Liu H, Sparks TE. Probabilistic simulation of solidiﬁcation microstructure evolution during laser-based metal deposition. In: Proceedings of 2013 annual international solid freeform fabrication symposium–an additive manufacturing conference; 2013. 23. Cervera GBM, Lombera G. Numerical prediction of temperature and density distributions in selective laser sintering processes. Rapid Prototyping J. 1999;5:12–26. 24. Childs T, Hauser C, Badrossamay M. Selective laser sintering (melting) of stainless and tool steel powders: experiments and modelling. Proc Inst Mech Eng Part B J Eng Manuf. 2005;219(4):339–57. 25. Denlinger ER, Heigel JC, Michaleris P. Residual stress and distortion modeling of electron beam direct manufacturing Ti-6Al-4V. Proc Inst Mech Eng Part B J Eng Manuf. 2015;229(10):1803–13. 26. Hodge N, Ferencz R, Solberg J. Implementation of a thermomechanical model for the simulation of selective laser melting. Comput Mech. 2014;54(1):33–51. 27. Hodge N, Ferencz R, Vignes R. Experimental comparison of residual stresses for a thermomechanical model for the simulation of selective laser melting. Addit Manuf. 2016;12:159–68. 28. Kollmannsberger S, Carraturo M, Reali A, Auricchio F. Accurate prediction of melt pool shapes in laser powder bed fusion by the non-linear temperature equation including phase changes. Integr Mater Manuf Innov. 2019;8(2):167–77. 29. Roy S, Juha M, Shephard MS, Maniatty AM. Heat transfer model and ﬁnite element formulation for simulation of selective laser melting. Comput Mech. 2018;62(3):273–84. 30. Shen N, Chou K. Numerical thermal analysis in electron beam additive manufacturing with preheating eﬀects. In: Proceedings of the 23rd solid freeform fabrication symposium, Austin, TX, p. 774–84; 2012. 31. Gusarov A. Homogenization of radiation transfer in two-phase media with irregular phase boundaries. Phys Rev B. 2008;77(14): 144201. 32. Gusarov A, Kruth J-P. Modelling of radiation transfer in metallic powders at laser treatment. Int J Heat Mass Transf. 2005;48(16):3423–34. 33. Gusarov A, Yadroitsev I, Bertrand P, Smurov I. Model of radiation and heat transfer in laser-powder interaction zone at selective laser melting. J Heat Transf. 2009;131(7):10. 34. Kollmannsberger S, Özcan A, Carraturo M, Egger J, Schröder A, Rank E. A multi-level model for the simulation of AM processes. In: Simulation for additive manufacturing; 2017. 35. Riedlbauer D, Scharowsky T, Singer RF, Steinmann P, Körner C, Mergheim J. Macroscopic simulation and experimental measurement of melt pool characteristics in selective electron beam melting of Ti-6Al-4V. Int J Adv Manuf Technol. 2017;88(5–8):1309–17. 36. Denlinger ER, Irwin J, Michaleris P. Thermomechanical modeling of additive manufacturing large parts. J Manuf Sci Eng. 2014;136(6):8. 37. Zhang Y, Guillemot G, Bernacki M, Bellet M. Macroscopic thermal ﬁnite element modeling of additive metal manu- facturing by selective laser melting process. Comput Methods Appl Mech Eng. 2018;331:514–35. 38. Neiva E, Badia S, Martín AF, Chiumenti M. A scalable parallel ﬁnite element framework for growing geometries. application to metal additive manufacturing. Int J Numer Methods Eng. 2019;119(11):1098–125. 39. Neiva E, Chiumenti M, Cervera M, Salsi E, Piscopo G, Badia S, Martín AF, Chen Z, Lee C, Davies C. Numerical modelling of heat transfer and experimental validation in powder-bed fusion with the virtual domain approximation. Finite Elem Anal Des. 2020;168: 103343. 40. Zaeh MF, Branner G. Investigations on residual stresses and deformations in selective laser melting. Prod Eng. 2010;4(1):35–45. 41. Keller N, Ploshikhin V. New method for fast predictions of residual stress and distortion of am parts. In: Solid Freeform Fabrication Symposium, vol. 25; 2014. 42. Li C, Liu Z, Fang X, Guo Y. On the simulation scalability of predicting residual stress and distortion in selective laser melting. J Manuf Sci Eng. 2018;140(4):10. 43. Setien I, Chiumenti M, van der Veen S, San Sebastian M, Garciandía F, Echeverría A. Empirical methodology to determine inherent strains in additive manufacturing. Comput Math Appl. 2019;78(7):2282–95. 44. Kruth J-P, Deckers J, Yasa E, Wauthlé R. Assessing and comparing inﬂuencing factors of residual stresses in selective laser melting using a novel analysis method. Proc Inst Mech Engi Part B J Eng Manuf. 2012;226(6):980–91. 45. Bruna-Rosso C, Mergheim J, Previtali B. Finite element modeling of residual stress and geometrical error formations in selective laser melting of metals. Proc Inst Mech Eng Part C J Mech Eng Sci. 2020;235:2022–38. 46. Parry L, Ashcroft I, Wildman RD. Understanding the eﬀect of laser scan strategy on residual stress in selective laser melting through thermo-mechanical simulation. Addit Manuf. 2016;12:1–15. 47. Bartel T, Guschke I, Menzel A. Towards the simulation of selective laser melting processes via phase transformation models. Comput Math Appl. 2019;78(7):2267–81. 48. Noll I, Bartel T, Menzel A. A computational phase transformation model for selective laser melting processes. Comput Mech. 2020;66(6):1321–42. 49. Proell SD, Wall WA, Meier C. On phase change and latent heat models in metal additive manufacturing process simulation. Adv Model Simul Eng Sci. 2020;7(24):1–32. Proell et al. Adv. Model. and Simul. in Eng. Sci. (2021) 8:24 Page 37 of 37 50. Michaleris P. Modeling metal deposition in heat transfer analyses of additive manufacturing processes. Finite Elem Anal Des. 2014;86:51–60. 51. Carraturo M, Jomo J, Kollmannsberger S, Reali A, Auricchio F, Rank E. Modeling and experimental validation of an immersed thermo-mechanical part-scale analysis for laser powder bed fusion processes. Addit Manuf. 2020;36: 52. Hussein A, Hao L, Yan C, Everson R. Finite element simulation of the temperature and stress ﬁelds in single layers built without-support in selective laser melting. Mater Des. 2013;1980–2015(52):638–47. 53. Popp A, Gitterle M, Gee MW, Wall WA. A dual mortar approach for 3D ﬁnite deformation contact with consistent linearization. Int J Numer Methods Eng. 2010;83(11):1428–65. 54. Puso MA, Laursen TA. Mesh tying on curved interfaces in 3d. Eng Comput. 2003;20:305–19. 55. BACI: A comprehensive multi-physics simulation framework. https://baci.pages.gitlab.lrz.de/website. Accessed 28 July 2021. 56. Carraturo M, Kollmannsberger S, Reali A, Auricchio F, Rank E. An immersed boundary approach for residual stress evaluation in selective laser melting processes. Addit Manuf. 2021;46:102077. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional aﬃliations.
"Advanced Modeling and Simulation in Engineering Sciences" – Springer Journals
Published: Oct 19, 2021
Keywords: Constitutive model; Thermomechanical coupling; Metal additive manufacturing; Mortar method; Finite element method
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.