Skip to main content

On phase change and latent heat models in metal additive manufacturing process simulation

Abstract

This work proposes an extension of phase change and latent heat models for the simulation of metal powder bed fusion additive manufacturing processes on the macroscale and compares different models with respect to accuracy and numerical efficiency. Specifically, a systematic formulation of phase fraction variables is proposed relying either on temperature- or enthalpy-based interpolation schemes. Moreover, two well-known schemes for the numerical treatment of latent heat, namely the apparent capacity and the so-called heat integration scheme, are critically reviewed and compared with respect to numerical efficiency and overall accuracy. Eventually, a novel variant of the heat integration scheme is proposed that allows to directly control efficiency and accuracy by means of a user-defined tolerance. Depending on the chosen tolerance, it is shown that this novel approach offers increased numerical efficiency for a given level of accuracy or improved accuracy for a given level of numerical efficiency as compared to the apparent capacity and the original heat integration scheme. The investigation and comparison of all considered schemes is based on a series of numerical test cases that are representative for application scenarios in metal powder bed fusion additive manufacturing.

Introduction

Additive manufacturing (AM) is widely considered to be a key technology for future advances in engineering. AM offers highest flexibility in part design while still achieving the mechanical properties required for functional parts [1]. In metal powder bed fusion additive manufacturing (PBFAM) multiple layers of metal powder are successively molten at selected positions, which eventually form the cross-sections of the final part after solidification. Energy is commonly deposited by a laser or electron beam giving rise to the respective names selective laser melting (SLM) and electron beam melting (EBM). These processes come with very challenging thermophysical phenomena on multiple length and time scales [2]. Accordingly, existing modeling approaches can be classified with respect to the resolved length scales: Macroscale approaches commonly aim at determining spatial distributions of physical fields such as temperature, residual stresses or dimensional warping on the scale of the design part. Mesoscale approaches resolve the length scale of individual powder particles on domains smaller than one powder layer to either study the melt pool thermo-fluid dynamics during the melting process [3,4,5,6,7,8] or the cohesive powder flow during the previous powder recoating process [9,10,11]. Last, microscale approaches predict the evolution of the metallurgical microstructure during solidification [12,13,14,15,16]. A broad overview of state-of-the-art modeling approaches on these different length scales can be found in [17]. The present article focuses on the development of a thermal macroscale model for metal PBFAM processes.

Macroscale PBFAM models typically treat the powder phase as a homogenized continuum described via spatially averaged thermal and mechanical properties, without resolving individual powder particles. Also, the complex fluid dynamics within the melt pool is typically not explicitly resolved. Instead, a pure thermo-(solid-)mechanical problem is solved, usually based on a Lagrangian description and a spatial finite element discretization, with specific temperature- and phase-dependent thermal and mechanical constitutive parameters for the (homogenized) powder phase, the melt phase and the solidified phase. On the one hand, from a modeling point of view, such a procedure considerably simplifies the coupling of the different phase domains. On the other hand, this approach seems to be well-justified for certain simulations and quantities of interest since the mechanical forces transferred from powder and molten phase onto the solid phase are often negligible in good approximation.

In their works [18,19,20], Gusarov et al. proposed a model for powder bed laser absorption, which has been incorporated in many existing macroscale modeling approaches. For example, by using this absorption model, [21, 22] proposed a thermo-mechanical finite element (FE) model accounting for temperature- and phase-dependent thermal and mechanical constitutive behavior. Further developments in this field consider e.g. the accuracy of the physical model by adding additional physical effects such as residual stress relaxation [23], improved models for temperature and phase-dependent thermal conductivity of the powder [24, 25] and melt [26] phase, anisotropic conductivity [27], phase-dependent laser absorptivity [28], thermodynamically consistent constitutive modeling based on phase energies [29], or by explicitly modeling the melt pool fluid dynamics [30], an approach inspired by similar schemes in the context of laser and electron beam welding [31,32,33]. Another important aspect for macroscale PBFAM models is computational efficiency, which has been addressed, among others, by applying dynamic mesh adaptivity schemes [34,35,36,37], code parallelization and load balancing techniques [38] as well as process layer agglomeration approaches [22, 37, 39].

The present work addresses two important aspects of macroscale PBFAM models, namely the modeling of phase change and latent heat effects. Concerning the first aspect, we propose phase fraction variables which allow to formulate temperature- and phase-dependent material parameters in phase transition regions by consistent interpolation of the single phase parameters. While the definition of phase fraction variables is often somehow hidden in existing works, the present contribution defines these phase fraction variables in a transparent and systematic manner. Moreover, as basis of these phase fraction variables, different interpolation strategies, e.g. temperature-based or enthalpy-based interpolation, are discussed in detail. An alternative formulation of phase fractions based on energy minimization can e.g. be found in [29].

Concerning the modeling of latent heat effects, the two schemes that are most widely used in PBFAM models, namely the simple apparent (heat) capacity approach [35, 40,41,42] and the more involved heat integration method [21, 43, 44], will be investigated in the present work. It is known in the context of PBFAM process simulation that the specific choice of the latent heat model might considerably influence the overall efficiency of the numerical model [21]. In the present work, the two aforementioned schemes, namely the apparent capacity approach and the heat integration method, are critically reviewed and compared with respect to numerical efficiency and overall accuracy. Eventually, a novel variant of the heat integration scheme is proposed that allows to directly control efficiency and accuracy by means of a user-defined tolerance. Depending on the chosen tolerance, it is shown that this novel approach offers increased numerical efficiency for a given level of accuracy or improved accuracy for a given level of numerical efficiency as compared to the apparent capacity and the original heat integration scheme. One example where high accuracy, e.g. in the prediction of melt pool shape and thermal gradients, is essential is the modeling of microstructure evolution on the basis of temperature solutions provided by macroscale PBFAM models. While also the prediction of residual stresses can be considered as one of the major objectives of macroscale PBFAM process simulation, the present study intentionally focuses on purely thermal problems to isolate the effects of primary interest, namely the formulation and comparison of different phase change and latent heat models. For the same reason, local effects in the melt pool such as evaporation and fluid flow have been purposely neglected.

This article is structured as follows: “General thermal model” section presents the mathematical problem statement and summarizes the main model constituents in space- and time-continuous form. Specifically, “Modeling of the laser heat source” section introduces the heat source modeling of a laser beam. Next, in “Modeling of phase change and latent heat” section the phase change problem is introduced. It is shown how both mentioned methods for modeling latent heat, namely the apparent capacity and the heat integration scheme, can be derived from a Lagrange multiplier potential. In “Modeling of temperature- and phase-dependent parameters” section temperature-history dependent material parameters are interpolated on the basis of properly defined phase fraction variables that allow to distinguish powder, solid and melt phase. The spatial and temporal discretization schemes as well as the general numerical solution procedure are outlined in “General numerical solution procedure” section. In “Discretization and algorithm of heat integration scheme” section the fully discretized version as well as algorithmic details of the heat integration method are discussed and eventually the novel tolerance-based variant of the heat integration scheme is proposed. Numerical experiments are presented in “Numerical results” section with a focus on accuracy and efficiency of the considered methods. Finally, a summary of the present contribution and a brief outlook on future research work is given in “Conclusion” section.

General thermal model

For the purpose of this study it is sufficient to focus on the following transient purely thermal problem described by the heat equation for the temperature T and appropriate boundary conditions. The problem statement in strong form is:

$$\begin{aligned} \begin{aligned} c(T)\,{\dot{T}}+\nabla \cdot \varvec{q}&= {\hat{r}},&\ \text {in}\ \Omega ,\\ T&= {\hat{T}},&\ \text {on}\ \Gamma _T, \\ \varvec{q}\cdot \varvec{n}&= {\hat{q}},&\ \text {on}\ \Gamma _{\varvec{q}}. \end{aligned} \end{aligned}$$
(1)

Temperatures T are prescribed on the boundary part \(\Gamma _T\) and heat fluxes on \(\Gamma _{\varvec{q}}\). The heat flux \(\varvec{q}\) is specified by Fourier’s law for isotropic material,

$$\begin{aligned} \varvec{q} = -k(T)\ \nabla T. \end{aligned}$$
(2)

Material properties, namely (volumetric) heat capacity c and conductivity k, may in general depend on the temperature T but also on the phase r (see “Modeling of temperature- and phase-dependent parameters” section for the definition of phase fractions and the interpolation of material parameters). In this contribution, spatial discretization will be based on the finite element method (see “General numerical solution procedure” section), i.e. (1) has to be transferred into its weak form via multiplication with a test function \(\delta T\) and integration by parts, viz.

$$\begin{aligned} \int _\Omega \delta T\,c(T){\dot{T}}\, \mathrm {d}\Omega - \int _\Omega \nabla \delta T\cdot \varvec{q}\,\mathrm {d}\Omega + \int _{\Gamma _{\varvec{q}}}\delta T\,{\hat{q}}\,\mathrm {d}\Gamma -\int _\Omega \delta T\, {\hat{r}}\,\mathrm {d}\Omega = 0, \end{aligned}$$
(3)

where the boundary conditions have already been inserted. By introducing the trial space \({\mathcal {V}} \!=\! \{T \, | \, T \!\in \! {\mathcal {H}}^1, T\!=\!{\hat{T}} \, \text {on} \, \Gamma _T \}\) as well as the weighting space \({\mathcal {W}} \!=\! \{\delta T \, | \, \delta T \!\in \! {\mathcal {H}}^1, \delta T\!=\!0 \, \text {on} \, \Gamma _T \}\), where \({\mathcal {H}}^1\) denotes the Sobolev space of functions with square-integrable first derivatives, (3) is equivalent to the strong form (1).

Modeling of the laser heat source

The source term \({\hat{r}}\) represents the energy deposited by the laser beam as a volumetric heat source based on the model for radiative and conductive heat transfer in powder beds by Gusarov et al. [20]. In the following, a summary of the main model constituents is given. Let a powder layer be distributed in the xy-plane, where the powder material extends in positive z-direction from the powder layer surface at \(z=0\) up to the layer thickness L at \(z=L\). A laser beam of nominal power W and size R is applied normal to this plane. The laser beam as well as the local coordinate system move in x-direction with a velocity v. The source term is then given in this local coordinate system relative to the laser beam center by

$$\begin{aligned} {\hat{r}}(r_h,z) = -\beta _h Q_0 \frac{\partial q}{\partial \xi '}, \end{aligned}$$
(4)

where \(r_h\) is the distance in the xy-plane from the laser beam center and \(\beta _h\) is the extinction coefficient. The nominal power density \(Q_0\) is radially distributed around the laser beam center as

$$\begin{aligned} Q_0 = {\left\{ \begin{array}{ll} \frac{3W_e}{\pi R^2}\left( 1-\frac{r_h}{R}\right) ^2 \left( 1+\frac{r_h}{R}\right) ^2,\ &{}0< r_h < R \\ 0, &{}\text {otherwise} \end{array}\right. }. \end{aligned}$$
(5)

The nominal laser power W has been averaged and reduced to an effective power \(W_e\) to account for various losses. Thus, (5) describes the spatial distribution in x- and y-direction. The normalized power density q is given in terms of the dimensionless coordinate \(\xi ' = \beta _h z\) as

$$\begin{aligned} q&= \frac{\rho _h a}{(4\rho _h-3)D}\nonumber \\&\quad \left\{ e^{-\lambda }(1-\rho _h^2)\left[ e^{-2a\xi '}(1-a) + e^{2a\xi '}(a + 1)\right] \right. \nonumber \\&\quad -\left[ e^{2a({\xi '-\lambda })}(1-a - \rho _h(a + 1))\right. \nonumber \\&\quad +\left. \left. e^{2a(\lambda - \xi ')}(a + \rho _h(a - 1) + 1)\right] (\rho _h e^{-2\lambda } + 3)\right\} \nonumber \\&\quad -\frac{3(1-\rho _h)(e^{-\xi '} - \rho _h e^{\xi '-2\lambda } )}{4\rho _h - 3}, \end{aligned}$$
(6)

with hemispherical reflectivity \(\rho _h\), constant \(a=\sqrt{1-\rho _h}\), optical thickness \(\lambda = \beta _h L\) and the constant

$$\begin{aligned} D&= (1-a)\left[ 1-a-\rho _h(1+a)\right] e^{-2a\lambda }\nonumber \\&-(1-a)\left[ 1+a-\rho _h(1-a)\right] e^{2a\lambda }. \end{aligned}$$
(7)

Finally, the derivative of (6) with respect to \(\xi '\) as required in (4) reads

$$\begin{aligned} \frac{\partial q}{\partial \xi '}&= \frac{(3-3\rho _h)(e^{-\xi '} + \rho _h e^{\xi '-2\lambda } )}{4\rho _h - 3}\nonumber \\&\quad +\frac{2a^2\rho _h}{D(4\rho _h - 3)}\left\{ e^{-\lambda }(1-\rho _h^2)\left[ e^{-2a\xi '}(a - 1) + e^{2a\xi '}(a + 1)\right] \right. \nonumber \\&\quad +(\rho _h e^{-2\lambda } + 3)\nonumber \\&\quad \times \left. \left[ e^{-2a(\lambda - \xi ')}(a + \rho _h(a + 1) - 1) + e^{2a(\lambda - \xi ')}(a + \rho _h(a - 1) + 1)\right] \right\} . \end{aligned}$$
(8)

Apart from the optical properties, (8) only depends on the z-coordinate \(z = \xi ' / \beta _h \).

Modeling of phase change and latent heat

So far we have not addressed the (crucial) phase change problem (i.e. melting and solidification), first from powder to molten and eventually from molten to solid phase. Consider a domain containing both solid and liquid phase separated by an interface \(\Gamma _m\), which is defined by the isotherm \(T = T_m\). Based on an energy balance at the interface, the following so-called Stefan-Neumann equation has to hold:

$$\begin{aligned} \underbrace{\varvec{n}_{sl}\cdot \left( k_s\left. \frac{\partial T}{\partial \varvec{x}}\right| _s-k_l\left. \frac{\partial T}{\partial \varvec{x}}\right| _l\right) }_{\Delta q_m} = h_m\varvec{n}_{sl}\cdot \varvec{v}_{sl}\quad \text {on}\ \Gamma _m, \end{aligned}$$
(9)

where \(\varvec{n}_{sl}\) is the interface normal vector, \(h_m\) the (volume-specific) latent heat of melting and the subscripts \((\cdot )_s\) and \((\cdot )_l\) denote quantities evaluated in the solid or liquid phase, respectively. The absorbed or released heat flux \(\Delta q_m\) is proportional to the velocity \(\varvec{v}_{sl}\) of the evolving interface, in general leading to discontinuous heat fluxes across the phase interface. The free boundary condition (9) is especially suitable for sharp interface models in Eulerian description when combined with explicit interface tracking schemes, e.g. via level set functions [37], whose temporal evolution is defined by \(\varvec{v}_{sl}\). In this work, as typical for macroscale PBFAM models, the thermal problem is described in a Lagrangian manner in combination with a diffuse interface model, i.e. it is assumed that phase change takes place across an extended interface volume of finite thickness. Within this volume, the temperature is constrained until the melting enthalpy \(h_m\) has been absorbed or released:

$$\begin{aligned} T-T_m&= 0 \quad \text {if } h_0 \le h \le h_0+h_m \end{aligned}$$
(10a)
$$\begin{aligned} {\dot{T}}&= 0 \quad \text {if } h_0 \le h \le h_0+h_m, \end{aligned}$$
(10b)

where \(h_0\) represents the enthalpy level at which melting starts. The left part of Fig. 1 (solid line) illustrates this concept. When deriving the weak form (3) in a variational manner, constraint Eq. (10a) [or alternatively its rate form (10b)] can be enforced via the following Lagrange mutliplier potential

$$\begin{aligned} \Pi _m = \int _\Omega \lambda (T-T_m)\,\mathrm {d}\Omega , \end{aligned}$$
(11)
Fig. 1
figure 1

Left: enthalpy-temperature diagram for isothermal (solid line) and mushy (dotted line) phase change. Right: capacity-temperature diagram shows a singularity for isothermal phase change (solid line) at \(T_m\). For the AC method the modified capacity \(c^*\) includes the effects of latent heat \(h_m\) within a regularized phase change interval \(\left[ T_s,T_l\right] \) similar to the mushy phase change type

where \(\lambda \) represents the Lagrange multiplier enforcing (10a). Its total variation yields

$$\begin{aligned} \delta \Pi _m = \int _\Omega \delta \lambda (T-T_m)\,\mathrm {d}\Omega + \int _\Omega \delta T\, \lambda \,\mathrm {d}\Omega . \end{aligned}$$
(12)

The first term in (12) represents the original constraint Eq. (10a), while the second term yields an additional contribution to the weak form (3):

$$\begin{aligned} \int _\Omega \delta T\left( c(T){\dot{T}} +\lambda \right) \mathrm {d}\Omega - \int _\Omega \nabla \delta T\cdot \varvec{q}\, \mathrm {d}\Omega + \int _{\Gamma _{\varvec{q}}}\delta T\,{\hat{q}}\, \mathrm {d}\Gamma -\int _\Omega \delta T\, {\hat{r}}\,\mathrm {d}\Omega = 0. \end{aligned}$$
(13)

Since \({\dot{T}} = 0\) during phase change, the Lagrange multiplier equals the enthalpy rate

$$\begin{aligned} {\dot{h}}(t) = \lambda (t) \quad \text {if} \quad h_0 \le h \le h_0+h_m, \end{aligned}$$
(14)

which leads to the following expression for the enthalpy during phase change

$$\begin{aligned} {h}(t) = h_0 + \int _{t(h_0)}^t \lambda ({\tilde{t}})\, \mathrm {d}{\tilde{t}} \quad \text {if} \quad h_0 \le h \le h_0+h_m, \end{aligned}$$
(15)

and eventually to an integral limiting condition for the Lagrange multiplier:

$$\begin{aligned} 0 \le \int _{t(h_0)}^{t} \lambda ({\tilde{t}})\, \mathrm {d}{\tilde{t}} \le h_m. \end{aligned}$$
(16)

Hu and Argyropoulos [45] review various methods to account for the latent heat associated with phase change. In the following two sections, the basics of two methods especially popular in the field of PBFAM modeling are briefly presented. In particular, we propose that these two schemes can be interpreted as different realizations of the constrained weak form (13).

Remark

So far we have considered pure materials, i.e. isothermal phase change at a fixed melting temperature \(T_m\). For alloys, phase change typically happens gradually, i.e. the latent heat is absorbed or released within a rather narrow temperature interval between solidus temperature \(T_s\) and liquidus temperature \(T_l\), in the following denoted as mushy phase change.

Remark

Throughout this work, the phase interface is implicitly defined by isotherms, a common choice in the context of PBFAM process simulation. Depending on the type of phase change (isothermal or mushy), one might take the isotherm at melting temperature (supplemented by a proper tolerance) or the isotherms at solidus and liquidus temperature to represent the interface.

Apparent capacity method

One of the simplest approaches to capture the effects of latent heat is the so-called apparent capacity (AC) method. The basic idea is to regularize the constraint (10) such that phase change takes place within a finite temperature interval of width 2d given by \(T \in [T-d; T+d]\). Throughout this work, the temperature bounds confining this regularized phase change interval for isothermal phase change will be represented by the same variables \(T_s=T-d\) and \(T_l=T+d\) as the solidus and liquidus temperature for mushy phase change. This choice seems reasonable as it will turn out that the algorithmic treatment of both cases is identical. Considering now the weak from (13), the apparent capacity method results from setting

$$\begin{aligned} \lambda = c_m(T) {\dot{T}}, \end{aligned}$$
(17)

where the factor \(c_m(T) \ge 0\) penalizes non-zero temperature rates \({\dot{T}} \ne 0\), i.e. violation of constraint (10b). Thus, according to (17) the Lagrange multiplier is not considered as an independent primary variable and (10b) is not satisfied exactly anymore. Inserting (17) into (13) allows to identify a modified capacity,

$$\begin{aligned} c^*(T) = c_m(T) +c(T), \end{aligned}$$
(18)

justifying the name apparent capacity method. The choice of \(c_m(T)\) is only restricted by the limiting condition (16), which after a change of variable eventually reads

$$\begin{aligned} \int _{t(T_s)}^{t(T_l)} c_m(T) {\dot{T}} \, \mathrm {d}t = \int _{T_s}^{T_l} c_m(T)\,\mathrm {d}T \, \dot{=} \, h_m. \end{aligned}$$
(19)

According to (19), small values of \(c_m(T)\) require a large regularized phase change interval typically resulting in a more good-natured numerical algorithm at the cost of lower solution accuracy and vice versa. The original phase change constraint (10) (solid line) as well as the regularized phase change constraint based on the apparent capacity (18) (dashed line) are illustrated in Fig. 1. In this contribution a smoothed triangular distribution is chosen for \(c_m\), as illustrated in the right part of Fig. 1.

As simple as the AC method may be, it suffers from one major drawback: An accurate representation of the phase change constraint, i.e. the choice of a small phase change interval 2d yields large values and steep gradients of the function \(c_m(T)\), which are not only challenging for numerical solution schemes (e.g. nonlinear solvers, see “General numerical solution procedure” section) but also prone to large time integration errors. In this case, the absorbed or released enthalpy during phase change (19), and thus the overall energy balance of the phase change problem, is only captured with low accuracy.

Heat integration method

Another popular and more advanced procedure is what Hu and Argyropoulos [45] call the heat integration (HI) method. It has first been applied to a FE setting by Rolph and Bathe [43] and is still used in more recent contributions [21, 44]. Here, we present the basic concept by proposing an alternative interpretation of the heat integration method as an augmented Lagrange constraint enforcement scheme [46]. The full algorithmic details of the method in the spatially and temporally discretized problem setting will be presented in “Discretization and algorithm of heat integration scheme” section.

In a first step, the Lagrange multiplier in the weak form (3) is replaced by an augmented Lagrange formulation for the constraint (10a) of the form

$$\begin{aligned} \lambda ={\tilde{\lambda }} + \epsilon (T-T_m) \end{aligned}$$
(20)

where \(\epsilon >0\) represents a penalty parameter. Comparable to an augmented Lagrange version based on the so-called Uzawa algorithm [46], the new Lagrange multiplier \({\tilde{\lambda }}\) is not considered as an independent primary variable but rather as a history variable within an iterative solution procedure of the form \({\tilde{\lambda }}^i=\lambda ^{i-1}\) leading to

$$\begin{aligned} \lambda ^i=\lambda ^{i-1} + \epsilon (T^{i}-T_m) \quad \text {for} \quad i \ge 1 \quad \text {with} \quad \lambda ^{0}:=0, \end{aligned}$$
(21)

where i represents an iteration counter to be defined in “General numerical solution procedure” section. In the final algorithm, the enthalpy rate \(\lambda ^i\) at iteration i has to fulfill the constraint Eq. (10a) and the enthalpy inequality (16) up to a user-defined tolerance, which will eventually define a stopping criterion for the iterative procedure (see “Tolerance-based heat integration scheme” section).

Remark

The HI scheme, which has originally been proposed in the time-discrete problem setting [43], can be recovered from (21) by choosing the penalty parameter according to \(\epsilon = c(T) / \Delta t\), where \(\Delta t\) represents the time step size to be introduced in “General numerical solution procedure” section. Now, to motivate this specific choice and to illustrate the working principle of the HI scheme, consider the (theoretical) case of a material with constant capacity \(c(T)=\text {const.}\) and vanishing conductivity \(k(T)=0\). Then, the time-discrete version of heat Eq. (1) without any additional constraint, i.e. \(i=0\) and \(\lambda ^0=0\) in (21), yields the solution \(c\Delta T^1 = {\hat{r}} \Delta t\) with \(\Delta T^1 = {\hat{r}} \Delta t / c \ne 0\). If the considered material point is undergoing phase change in the current (and previous) time step, the solution can be expressed as \(\Delta T^1 = T^1-T_m = {\hat{r}} \Delta t / c \ne 0\), which violates constraint (10a). Thus, an additional iteration \(i=1\) with \(\lambda ^1=\epsilon (T^1-T_m)\) has to be conducted. For this example, the specific choice \(\epsilon = c / \Delta t\) allows to find the correct temperature solution that is consistent with (10a) already in the first iteration \(i=1\):

$$\begin{aligned} \Delta T^2 = {\hat{r}}{\Delta t}/{c} - \lambda ^1{\Delta t}/{c} = {\hat{r}}{\Delta t}/{c} - \underbrace{(T^1-T_m)}_{{\hat{r}} \Delta t / c} = 0 \end{aligned}$$

Even though general scenarios with \(c(T) \ne \text {const.}\) and \(k(T) \ne 0\) will typically require a higher number of iterations, the choice \(\epsilon = c / \Delta t\) still seems reasonable and is well-established also in this case.

Remark

The original definition of the Uzawa algorithm defines an additional fixed-point iteration wrapped around the iteration loop of the nonlinear solver (e.g. a Newton–Raphson scheme as discussed in “General numerical solution procedure” section) to perform the iterative update (21). In contrast, the HI algorithm presented in “Discretization and algorithm of heat integration scheme” section will perform these updates directly during the iterations of the nonlinear solver without any additional “outer” iteration loop.

Modeling of temperature- and phase-dependent parameters

As common in macroscale PBFAM models, all three phases are modeled by varying material parameters of a solid material law. On the one hand, this procedure considerably simplifies the numerical schemes for coupling the different phases. On the other hand, this approach seems also to be justified for the thermo-mechanical problem since the mechanical forces transferred from powder and molten phase (with vanishing stiffness) onto the solid phase are negligible in good approximation for many questions. It is also important to note that the change from powder phase to molten phase is irreversible. In the following, we will consider the different types of phase changes in more detail.

Fig. 2
figure 2

Temperature history diagram illustrating the two-dimensional nature of material parameter interpolation between powder (p), melt (m) and solid (s). Arrows indicate possible movement within the diagram

Figure 2 gives an illustration of the possible states of a material point. The current temperature T is on the ordinate, while the abscissa shows the highest temperature ever reached, \(T_\text {max}\). The different areas correspond to different mixtures of powder, melt and solid. By definition there cannot be any possible state for \(T>T_\text {max}\) and this area is blanked out. If \(T_\text {max}<T_s\), all material must still be in powder form (p). If \(T_\text {max}>T_l\), there is no more powder, all material is consolidated and thus must be solid (s) or molten (m). The exact constitution of the two then depends on the current temperature. The same reasoning applies to the current temperature. If \(T>T_l\), all material must be molten. If \(T<T_s\), material is a mixture of powder and solid, where the exact constitution is depending on \(T_\text {max}\). Perhaps the most interesting scenario is obtained when \(T_s<T<T_\text {max}<T_l\). In this case some powder is still left and the consolidated phase consists of melt and solid. The arrows in Fig. 2 indicate the possible evolution of phases. Due to the definition of \(T_\text {max}\) the only way to increase its value is an irreversible movement along the black diagonal line in Fig. 2. For all other locations in the diagram a reversible horizontal movement is possible. With these considerations a temperature-based interpolation procedure for any material parameter can be derived.

Temperature-based interpolation

First, the focus is on the transition from (non-powder) solid to melt. The liquid fraction g is introduced as

$$\begin{aligned} g(T) = {\left\{ \begin{array}{ll} 0, &{} T < T_s\\ \frac{T-T_s}{T_l-T_s}, &{}T_s \le T \le T_l\\ 1, &{} T > T_l. \end{array}\right. } \end{aligned}$$
(22)

If only solid and molten phase were present, any material parameter f could be interpolated from the solid and melt values \(f_s\) and \(f_m\). The history-dependent material behavior is captured by the fraction of consolidated material \(r_c\) defined as

$$\begin{aligned} r_c(t) = {\left\{ \begin{array}{ll} 1, \quad &{} \text {if } r_c(0) = 1 \text { (i.e. initially consolidated)}\\ \max \limits _{{\tilde{t}}\le t}\, g(T({\tilde{t}})), \quad &{}\text {if } r_c(0) = 0 \text { (i.e. initially powder)}.\\ \end{array}\right. } \end{aligned}$$
(23)

The time argument is explicitly stated to emphasize the history-dependence of \(r_c(t)\). The consolidated fraction is initialized with a proper start value \(r_c(0)\) (zero/one for locations initially covered with powder/consolidated material). For example, in a region that has initially been covered with powder, definition (23) equals the all-time maximum of the liquid fraction g at this location which, according to (22), carries the same information as the maximum temperature \(T_\text {max}\). In a region that has initially been covered with solid material, definition (23) equals one for all times since solid material can never transform back to powder. The monotonously increasing fraction of consolidated material \(r_c\) together with the liquid fraction g allow to define volume fractions for powder, melt and solid phases according to:

$$\begin{aligned} r_p&= 1-r_c, \end{aligned}$$
(24)
$$\begin{aligned} r_m&= g, \end{aligned}$$
(25)
$$\begin{aligned} r_s&= r_c - g. \end{aligned}$$
(26)

Their physical motivation is as follows: The powder fraction \(r_p\) given in (24) is by definition the complement of the consolidated fraction \(r_c\). The molten fraction \(r_m\) in (25) is independent of the history and is always determined by (22). The solid fraction \(r_s\) defined in (26) is the part of the consolidated fraction which is not molten. Note that definitions (24), (25) and (26) satisfy partition of unity and are thus suitable for interpolation.

Fig. 3
figure 3

Evolution of a material parameter depending on temperature history. Left: material is completely molten and consolidated. Right: temperature history stays below liquidus temperature and thus a fraction of material stays powder

Any material parameter f can now be interpolated from the single phase values \({f_p, f_m, f_s}\) weighted by the corresponding fractions

$$\begin{aligned} f(T)&= r_p(T)f_p + r_m(T)f_m + r_s(T)f_s. \end{aligned}$$
(27)

For the special choice \(f_p = f_s\), a two phase interpolation without history-dependent behavior would be recovered. The dependence on temperature is explicitly stated in (27) since this requires a consistent linearization to achieve robust convergence of the nonlinear solver. Figure 3 shows the evolution of an exemplary material parameter f over the temperature history for the chosen liquid fraction definition (22). In the left diagram powder melts completely and consequently the parameter f takes on the value of a solid after cooling down below \(T_s\). The right diagram shows partial melting: After cooling down below \(T_s\) the parameter is a weighted average of the powder and solid value based on the consolidated fraction which is now smaller than 1. Figure 4 shows the same interpolation in a three-dimensional representation which makes the history-dependence on \(T_\text {max}\) more explicit and shows that the parameter interpolation is also continuous over the history. This is important in the modeling of PBFAM processes since it ensures a continuous transition of material parameters between regions of molten or solidified material and regions that are still covered with unmolten powder. Note that the definition of the liquid fraction (22) determines the exact shape of the interpolating curve (27). The kinks at \(T_s\) and \(T_l\) could also be smoothed out if this seems necessary for an improved numerical behavior (e.g. robust convergence of the nonlinear solver, see “General numerical solution procedure” section).

Remark

The definitions (24), (25) and (26) imply that in a material mixture containing solid and powder the solid fraction always melts first, i.e. when increasing the temperature at a scenario \(T_s<T<T_\text {max}<T_l\) the powder fraction does not increase before the previous maximum temperature \(T_\text {max}\) is exceeded. With more history variables, also a different behavior could be realized (e.g. powder should melt first) but this is not deemed necessary in the context of macroscale PBFAM.

Enthalpy-based interpolation

In the case of isothermal phase change (if not regularized by an AC scheme), the definition of the liquid fraction based on temperature as in (22) is not directly applicable. A temperature-based parameter interpolation would still be possible based on a definition of solidus and liquidus temperature \(T_s=T_m-d\) and \(T_l=T_m+d\) as numerical regularization values as done in the AC scheme. However, a more consistent alternative approach utilizes the accumulated melt enthalpy h(t) in (15) to define

$$\begin{aligned} g(t) = {\left\{ \begin{array}{ll} 0, &{} h = h_0\\ \frac{h(t) - h_0}{h_m}, &{} h_0< h < h_0 + h_m\\ 1, &{} h = h_0 + h_m. \end{array}\right. } \end{aligned}$$
(28)

Essentially, the liquid fraction is defined as the ratio of already absorbed melt enthalpy to the latent heat required for phase change. All other parts of the interpolation (27) remain unchanged. This liquid fraction will prove useful in combination with the HI method as further discussed in “Discretization and algorithm of heat integration scheme” section.

Fig. 4
figure 4

Three-dimensional visualization of parameter interpolation. Material parameters are continuous in T and \(T_\text {max}\) direction. The blue curve shows the evolution of the parameter when \(T_\text {max}<T_s\)

General numerical solution procedure

The weak form of the heat Eq. (3) is discretized in space by the finite element method (FEM) following a Bubnov–Galerkin approach:

$$\begin{aligned} T({\varvec{x}},t) = \sum _{j} N_j({\varvec{x}}) T_j(t), \quad \delta T({\varvec{x}}) =\sum _j N_j ({\varvec{x}}) \delta T_j, \end{aligned}$$
(29)

where \({\varvec{x}}\) is the spatial position, and \(T_j(t)\) as well as \(\delta T_j\) are the nodal temperatures and temperature variations. For time discretization a one-step theta time integration scheme is employed, which, for the model equation \({\dot{\phi }}=f(\phi ,t)\), is defined as

$$\begin{aligned} \phi _{n+1} = \phi _n + \theta \Delta t f(\phi _{n+1},t_{n+1}) + (1-\theta ) \Delta t f(\phi _n,t_n). \end{aligned}$$
(30)

Here, \(\Delta t:= t_{n+1} - t_{n}\) is the time step size and the subscript \((.)_n\) indicates a quantity that is evaluated at the discrete time step \(t_{n}\). Together, the spatial and temporal discretization result in a fully discrete, nonlinear system of equations \({\varvec{R}}({\varvec{T}}_{n+1}) = 0\), where \({\varvec{R}}\) is the global residual vector and \({\varvec{T}}_{n+1}\) the global vector of nodal temperatures at \(t_{n+1}\). The equations are nonlinear due to the temperature-dependence of the heat capacity and thermal conductivity and due to the underlying phase change subproblem represented by (the non-smooth) constraint Eq. (10a). The system of nonlinear equations is linearized and solved iteratively with a Newton–Raphson scheme, which yields the following iterative update procedure

$$\begin{aligned} {\varvec{T}}_{n+1}^{i+1} = {\varvec{T}}_{n+1}^{i} + \Delta {\varvec{T}}_{n+1}^{i+1} \quad \text {with} \quad \left. \frac{\partial {\varvec{R}}}{\partial {\varvec{T}}}\right| _{{\varvec{T}}_{n+1}^i} \Delta {\varvec{T}}_{n+1}^{i+1} = -{\varvec{R}}({\varvec{T}}_{n+1}^i). \end{aligned}$$
(31)

The following two convergence criteria \(\vert \vert {\varvec{R}}({T}_{n+1}^{i+1})\vert \vert _2 < \varepsilon _R\) and \(\vert \vert \Delta {\varvec{T}}_{n+1}^{i+1}\vert \vert _2 < \varepsilon _T\) are considered, i.e., for convergence both the norm of the residual and the iterative solution vector increment have to fall below given tolerances. All implementations and simulations presented in the next sections have been performed in the in-house research code BACI [47], a parallel, multi-physics finite element framework.

Discretization and algorithm of heat integration scheme

Requirements and objective: In “Numerical results” section, it will turn out that the HI method typically describes the phase change problem more accurately [e.g. in terms of constraint (10a)] as the AC scheme. Unfortunately, the HI method in its original form is known to typically result in a slow Newton–Raphson convergence [21, 45], which can be traced back to residual manipulations in subsequent iteration steps i on the basis of (21), which are not accompanied by an associated consistent linearization contribution. In order to satisfy the constraint Eq. (10a), the original HI scheme typically leads to such manipulations, i.e. to changes of the Lagrange multiplier \(\lambda ^i\) in subsequent iterations, throughout the entire Newton–Raphson loop, which slows down convergence considerably. While the basics of this original HI scheme are presented in “Original heat integration scheme” section, we propose a novel realization of the HI scheme in “Tolerance-based heat integration scheme” section. This variant allows to control the accuracy of constraint enforcement in (10a) via a user-defined tolerance. For practically relevant choices of this tolerance, \(\lambda ^i\) in (21) typically takes on a constant value after only a few iterations, thus allowing for quadratic convergence of the Newton–Raphson scheme afterwards.

The basic idea of the HI method has been introduced in an abstract manner in “Heat integration method” section for the case of isothermal phase change. Essentially, the method performs the update (21) for the Lagrange multiplier \(\lambda ^i\) occurring in the discrete version of the weak form (13) for each Newton–Raphson iteration i at which constraint Eq. (10a) is violated. It is important to note that the HI method in the spatially discretized setting enforces this constraint at element nodes (and not at integration points). For this purpose, the nodal volume \(V_k\) is defined for each node k as

$$\begin{aligned} V_k = \int \limits _\Omega N_k\, \mathrm {d}\Omega , \end{aligned}$$
(32)

Taking advantage of the partition of unity \(\sum _k N_k=1\), (32) implies that the sum of all the nodal volumes \(V_k\) associated with an element indeed equals the volume occupied by this element. Considering the Lagrange multiplier term in (13), the integral over \(\Omega \) can then be approximated by a nodal evaluation of the integrand \(\lambda \) times the nodal volume. Thus, the contribution to the residual entry \(R_k\) yields:

$$\begin{aligned} R_k = \int _\Omega N_k \lambda \, \mathrm {d}\Omega \approx \lambda ({\varvec{x}}_k) \int _\Omega N_k \, \mathrm {d}\Omega = \lambda _k V_k =: {\dot{H}}_{k} \end{aligned}$$
(33)

Since \(\lambda _k\) is a volume-specific enthalpy rate, \({\dot{H}}_{k}\) can be interpreted as an absolute enthalpy rate. Similarly, the latent heat of melting associated with node k is

$$\begin{aligned} {H}_{m,k} = h_m V_k, \end{aligned}$$
(34)

Employing (21) to express the Lagrange multiplier in (33) at \(t_{n+1}\) yields:

$$\begin{aligned} \!\!\!\!\!\!\!\! {\dot{H}}_{k,n+1}^i = V_k \lambda _{k,n+1}^i = V_k [\lambda _{k,n+1}^{i-1} + \epsilon (T_{k,n+1}^{i}-T_m)] = {\dot{H}}_{k,n+1}^{i-1} + \frac{\Delta H_{k,n+1}^i}{\Delta t} ,\!\!\!\! \end{aligned}$$
(35)

where the penalty parameter was chosen as \(\epsilon = {c}/{\Delta t}\) (see “Heat integration method” section) and hence

$$\begin{aligned} \Delta H_{k,n+1}^i := c(T_{k,n+1}^{i}-T_m)V_k. \end{aligned}$$
(36)

The iterative update rule (35) together with (36) is the discrete counter part to (21). Employing a backward Euler scheme for time integration of \({\dot{H}}_{k,n+1}^i\) in (35) and assuming (without loss of generality) \(h_0=0\) in (15), the nodal enthalpy \({H}_{k,n+1}^i\) is

$$\begin{aligned} \begin{aligned} {H}_{k,n+1}^i&= {H}_{k,n} + \Delta t {\dot{H}}_{k,n+1}^i = \underbrace{{H}_{k,n} +\Delta t {\dot{H}}_{k,n+1}^{i-1}}_{{H}_{k,n+1}^{i-1}} + \Delta H_{k,n+1}^i \\&= {H}_{k,n+1}^{i-1} + \Delta H_{k,n+1}^i, \end{aligned} \end{aligned}$$
(37)

where \({H}_{k,n}\) represents the nodal enthalpy at the converged configuration of the last time step \(t_n\). Thus, the iterative update rule (37) for the nodal enthalpy \({H}_{k,n+1}^i\) has the same form as the update rule (35) for the nodal enthalpy rate \({\dot{H}}_{k,n+1}^i\) but with the different initial values \({H}_{k,n+1}^0={H}_{k,n}\) and \({\dot{H}}_{k,n+1}^0=0\) required in the first iteration \(i=1\) of a time step. Also the enthalpy limiting condition (16) during phase change can be transferred to the discrete problem setting, which reads (for \(h_0=0\)):

$$\begin{aligned} 0 \le {H}_{k,n+1}^i \le H_{m,k}. \end{aligned}$$
(38)

Essentially, the working principle of the heat integration scheme is to add residuum contributions in (33) defined by the iterative update scheme (35) as long as the discrete limiting condition (38) is satisfied. The detailed algorithm of the original HI scheme as well as the required adaptions for the proposed tolerance-based HI scheme will be presented in the next two sections.

Remark

As noted by the authors of the original work [43] and in accordance with (33), the HI method requires the consistent use of the nodal lumped capacity. Thus, when using this algorithm, the capacity matrix must enter the residual and Jacobian in the Newton–Raphson algorithm (31) in lumped form to obtain a robust scheme. Note also that no linearization contribution associated with the residual term in (33) is considered in the context of the Newton–Raphson method.

Original heat integration scheme

We will first introduce our understanding of the original HI method [43]. The specific notation was inspired by the formulation in [44]. So far, the considerations have been restricted to isothermal phase change but they can be extended to the treatment of mushy phase change. The only difference is the fact that during melting the temperature is not constrained to the constant value \(T_m\) but rather to a gradually evolving intermediate temperature \(T'\) between \(\left[ T_s,T_l\right] \). Thus, in the following, we will present the more general algorithm for mushy phase change and point out the relevant differences to the isothermal scenario.

Each node k stores the nodal enthalpy \({H}_{k,n+1}^{i-1}\) as a history-variable. At the beginning of the simulation, it is initialized to zero for solid material. At the first Newton–Raphson iteration of a time step it is initialized with the converged value from the last time step, i.e. \({H}_{k,n+1}^{0}={H}_{k,n}\). Now, after each Newton–Raphson iteration i in time step \(n+1\) the following calculations are performed:

  1. 1.

    Skip node k if

    $$\begin{aligned}&\left[ T_{k,n+1}^i< T_s\ \text {and}\ T_{k,n} < T_s\right] \ \mathbf{or} \nonumber \\&\left[ T_{k,n+1}^i> T_l\ \text {and}\ T_{k,n} > T_l\right] \end{aligned}$$
    (39)

    which means it is not undergoing phase change. In this case the increment \(\Delta H_{k,n+1}^{i}\) in (35) is set equal to zero. Here, \(T_{k,n}\) denotes the converged temperature value of node k at the last time step n.

  2. 2.

    Else, for each node k which is undergoing phase change compute the increment

    $$\begin{aligned} \Delta H_{k,n+1}^{i} = c'\left( T_{k,n+1}^{i} -T'\right) V_k. \end{aligned}$$
    (40)

    Here, \(T'\) is an intermediate temperature given by

    $$\begin{aligned} T' = T_s + \left| \frac{{H}_{k,n+1}^{i-1}}{H_{m,k}}\right| \, (T_l-T_s), \end{aligned}$$
    (41)

    calculated from the amount of latent heat already absorbed (released) during melting (freezing). The modified capacity is computed as

    $$\begin{aligned} c' = \left( {\frac{T_l-T_s}{h_m} + \frac{2}{c_s+c_l}}\right) ^{-1}, \end{aligned}$$
    (42)

    where \(c_s\) and \(c_l\) are the values of heat capacity at \(T_s\) and \(T_l\), respectively. In case of isothermal phase change, the intermediate temperature simplifies to the melting temperature, i.e. \(T' = T_m = T_s = T_l\), and the modified capacity simplifies to the averaged capacity at the melting point, i.e. \(c' = \frac{1}{2}(c_s+c_l)\) such that (40) becomes identical to (36).

  3. 3.

    Limit each increment \(\Delta H_{k,n+1}^{i}\) such that the following condition is fulfilled:

    $$\begin{aligned} 0 \le H_{k,n+1}^{i-1} + \Delta H_{k,n+1}^{i} \le H_{m,k}. \end{aligned}$$
    (43)

    If condition (43) is violated, the increment \(\Delta H_{k,n+1}^{i}\) is limited such that the respective bound in (43) is exactly met (e.g. \(\Delta H_{k,n+1}^{i} = H_{m,k} - H_{k,n+1}^{i-1}\) if the right bound in (43) was exceeded). Afterwards update the nodal enthalpy:

    $$\begin{aligned} {H}_{k,n+1}^i = {H}_{k,n+1}^{i-1} + \Delta H_{k,n+1}^i \quad \text {with} \quad {H}_{k,n+1}^{0}={H}_{k,n}. \end{aligned}$$
    (44)
  4. 4.

    For each node k where \(\vert \Delta H_{k,n+1}^i\vert > 0\) reset the temperature

    $$\begin{aligned} T_{k,n+1}^i = T'. \end{aligned}$$
    (45)
  5. 5.

    Calculate the updated enthalpy rate according to (35),

    $$\begin{aligned} {\dot{H}}_{k,n+1}^i = {\dot{H}}_{k,n+1}^{i-1} + \frac{1}{\Delta t} \Delta H_{k,n+1}^i \quad \text {with} \quad {\dot{H}}_{k,n+1}^{0}=0, \end{aligned}$$
    (46)

    and add the residual contribution (33). The start value \({\dot{H}}_{k,n+1}^0 = 0\) in the first iteration \(i=1\) of each time step is the equivalent of \(\lambda ^0\) in (21).

It is emphasized that limiting condition (43) allows to use the same algorithm for both melting and solidification as discussed in the remark at the end of this section.

Tolerance-based heat integration scheme

As stated at the beginning of “Discretization and algorithm of heat integration scheme” section, the original HI method is known to result in a slow Newton–Raphson convergence due to the residual manipulations (46), which are not accompanied by an associated consistent linearization contribution and which would typically occur throughout the entire Newton–Raphson loop if no additional, tolerance-controlled abort criterion is utilized. To improve the algorithm, we propose to introduce a tolerance \(\varepsilon _\text {tol}\) to control the accuracy of the phase change representation. Hence, we stop adding increments \(\Delta H_{k,n+1}^i\) in the algorithm above when these are small compared to the latent heat of melting, i.e.,

$$\begin{aligned} \left| \frac{ \Delta H_{k,n+1}^i}{{H}_{m,k}}\right| < \varepsilon _\text {tol}. \end{aligned}$$
(47)

where \(\varepsilon _\text {tol}<1\) is a relative tolerance describing the relative amount of latent heat which is not absorbed (released) after melting (solidification) is complete. Inserting (34) and (40) into (47) yields an alternative for step 1 in the algorithm above:

1\(^*\):

Skip node k if

$$\begin{aligned} \left| T_{k,n+1}^i - T' \right| < \varepsilon _\text {tol}\frac{ h_m}{ c'}. \end{aligned}$$
(48)

The new criterion provides a way to stop the HI algorithm when constraint (10a) is satisfied with a certain accuracy. At first glance, the new criterion (48) seems to significantly change the outcome of the algorithm since (40) will be evaluated for nodes that are far away from a phase change. However, it can easily be verified that the signed limiting condition (43) will automatically skip all nodes that are not meant to undergo phase change on the basis of their current temperature value (see also the discussion in the remark below).

Remark

Consider a node that is heated up and the material is undergoing a melting process. Initially, the node is in solid state, i.e. \(H_{k,n+1}^{i-1}=0\), with a temperature below the solidus temperature, i.e. \(T_{k,n+1}^i < T_s\). The calculated increment (40) would be negative. Since there is no phase change yet, no increment should be added. Indeed, limiting according to (43) will not allow a negative increment with a zero history and the increment is set to zero. When the temperature rises above the solidus temperature \(T_s\), however, the increments become positive and will be considered according to (43). Positive increments are added to the the nodal latent heat \(H_{k}\) until it reaches the allowed maximum \(H_{m,k}\). Then the material is fully molten.

Next, consider the inverse process, i.e. cooling of a node with material initially in the molten state. If material is molten, then \(H_{k,n+1}^{i-1} = H_{m,k}\). Looking at the limiting condition (43) shows that only negative increments are allowed, when material is molten. As expected, negative increments (40) are obtained, and the solidification process is initiated, as soon as temperature drops below the liquidus temperature \(T_l\). The negative increments are added to \(H_{k}\) until it reaches a value of zero. Then the material has returned to the solid state. The phase change as it was just described is fully reversible.

Remark

Now that the details of the HI scheme have been introduced, the enthalpy-based liquid fraction (28) used for parameter interpolation in “Modeling of temperature- and phase-dependent parameters” section can be further specified. For node k it reads

$$\begin{aligned} g_k = \frac{H_{k,n+1}^i}{H_{m,k}}. \end{aligned}$$
(49)

In the numerical examples of the following section, this liquid fraction based on latent heat will be employed in combination with the isothermal HI scheme. All latent heat schemes employing a finite phase change interval \([T_s;T_l]\), will use the temperature-based parameter interpolation with liquid fraction (22).

Numerical results

Solidification front in a 1D slab

To validate the implementation, first a series of numerical experiments are conducted on a one dimensional domain for which an analytic solution is available [45]. This example has already been used to show the validity of methods for capturing latent heat [40]. A pseudo one-dimensional slab (material properties of ice/water) with length \(L=1\,{\hbox {m}}\) is subject to a fixed temperature \({\hat{T}}=253\,{\hbox {K}}\) on its left edge at \(x=0\). The initial temperature in the whole slab is \(T_0=283\,{\hbox {K}}\). The left part of Fig. 5 illustrates this scenario. The interface separating frozen and liquid water will slowly travel from left to right. Material parameters for water are given in Table 1, they are taken directly from [40]. The problem is discretized with 25, 50 or 100 linear finite elements in space and three different fixed step sizes \(\Delta t \in \{200\,{\hbox {s}},400\,{\hbox {s}},800\,{\hbox {s}}\}\) in time. Total simulation time is \(t_f = 72\cdot 10^3\,{\hbox {s}}\). At this point in time the temperature on the right edge is still at the initial level and the analytic solution (which is calculated on a semi-infinite domain) remains valid.

Fig. 5
figure 5

Geometry, thermal loads, boundary and initial conditions for solidification front example (left) and melting volume example (right)

Table 1 Material parameters of water [40] for “Solidification front in a 1D slab” and “Melting volume in a 1D slab” sections

The introduced HI method will be used in four variants by distinguishing (a) isothermal and mushy phase change as well as (b) the original criterion (39) and the novel tolerance-based criterion (48). In this example, phase change of water is isothermal and thus the isothermal HI methods with either original or tolerance-based criterion can be applied directly. For the AC method an artificial phase change interval is chosen with \(T_s=270\,{\hbox {K}}\) and \(T_l=276\,{\hbox {K}}\), i.e., \(d=3\,{\hbox {K}}\). The same interval is used to apply the original and tolerance-based mushy HI methods. Additionally, for tolerance-based HI the tolerance is chosen as \(\varepsilon _\text {tol}=0.001\), i.e., up to 0.1% of latent heat will be neglected.

Fig. 6
figure 6

Melting front example: temperature profiles at \(t_f\) of analytic solution compared to numerical results with AC and HI method. 100 elements, \(\Delta t = 200{\hbox {s}}\)

Fig. 7
figure 7

Melting front example: maximum error in temperature \(\left\| \frac{T-T_\text {ref}}{T_0-{\hat{T}}}\right\| _\infty \) of different methods for latent heat depending on spatial and temporal discretization. Missing data points for the AC scheme indicate that no convergence was achieved for these parameter choices

All approaches yield results that agree very well with the analytic solution. Figure 6 shows the solutions obtained with AC and original isothermal HI method on the finest mesh as an example. The maximum errors in the numerical solutions provided by the different methods are shown for the three investigated meshes and time step sizes in Fig. 7. The maximum errors lie below 4% for the coarsest mesh and around 2% for the finer meshes, which is deemed accurate enough for the intended use case. Within the considered scope, the time step size seems to have only minor effect on the accuracy. All versions of HI schemes produce errors that are slightly higher compared to the error from the AC scheme. However, for the large time step sizes \(\Delta t = 400\,{\hbox {s}}\) and \(\Delta t = 800\,{\hbox {s}}\) the AC method does not always converge.

Fig. 8
figure 8

Melting front example: average number of iterations of different methods for latent heat depending on spatial and temporal discretization. Missing data points for the AC scheme indicate that no convergence was achieved for these parameter choices

Fig. 9
figure 9

Melting volume example: temperature profiles at \(t_f\) obtained with five different methods for latent heat. 100 elements, \(\Delta t = 200{\hbox {s}}\)

The real difference between the methods comes to light when numerical efficiency is investigated in terms of Newton iterations needed per time step. We choose to analyze Newton iterations as a measure for the efficiency of the proposed methods instead of computational time, which typically leads to a stronger dependency on the specific code implementation. Still, the resulting computation time, which is usually of practical interest, scales directly with the number of Newton iterations. Figure 8 shows a strong dependence of the original HI method [43] on the spatial discretization with heavily increased iterations per time step in case of finer spatial resolution. On top of that a larger time step leads to a further increase in iterations. Both isothermal and mushy version of the original method suffer from this effect. Hodge et al. [21] mention that small time steps had to be used because of the original HI method. However, our proposed tolerance-based method does not only require significantly less iterations per step in every scenario but is also less sensitive to spatial and temporal discretization. This seems beneficial when moving to simulation of PBFAM processes on a part-scale.

Melting volume in a 1D slab

We investigate the same variants on a slightly modified second example (Fig. 5, right). The Dirichlet condition on the left boundary is dropped and all faces are assumed to be insulating. Instead a spatially varying source term \({\hat{r}} = 20,000(1-x)\,\nicefrac {\hbox {W}}{\hbox {m}^3}\) is applied to the whole slab of water which is initially frozen at \(T_0=263\,{\hbox {K}}\). Material parameters for solid and liquid water are again given in Table 1. In contrast to the previous example, melting will not take place at a single node representing the phase interface. Instead a whole volume can be in phase transition (i.e. melting). The same spatial and temporal discretizations from before are used, total simulation time is \(t_f = 20\cdot 10^3\,{\hbox {s}}\).

Again, the AC method and the four variants of HI are applied with the settings from above. No analytic solution is available for this scenario. Figure 9 shows the obtained temperature profiles along the one dimensional slab for the fine discretization with 100 elements and a step size of \(200\,{\hbox {s}}\). Obviously, AC and mushy HI methods will not keep temperatures fixed at the melting temperature of 273 K during phase change. When one is concerned about exact representation of isothermal phase change, only the original and tolerance-based version of isothermal HI are accurate enough, although some oscillation around the melting temperature is observed. Looking at the final temperatures on the left edge, when melting is already finished, reveals that the latent heat of melting still is captured with good accuracy by all methods, and the predicted temperatures agree well. Given the large temperature range prevalent in the targeted application (PBFAM simulation on part-scale), a highly accurate representation of temperature profiles around the melting point is not of highest practical importance.

Fig. 10
figure 10

Melting volume example: temperature profiles at \(0.5\,{t_f}\) obtained from the original and the proposed tolerance-based HI scheme (100 elements, time step varying)

A short look on the accuracy of isothermal HI schemes is taken in this paragraph. Solutions with the HI scheme can become quite inaccurate around the melting temperature especially when larger time steps are used. Figure 10 shows such solutions of the original scheme at \(t=0.5t_f\) for different step sizes and compares it to our proposed tolerance-based method. Both methods cannot exactly enforce isothermal phase change which would be characterized by a horizontal plateau region at \(T_m\). The original method’s criterion (39) for determining nodes undergoing phase change proves to be ill-suited for this scenario. Larger time steps lead to larger undershoots in temperature. The fluctuations around melting temperature obtained with the proposed tolerance-based HI scheme on the contrary are independent of step size and only controlled by the tolerance \(\varepsilon _\text {tol}\). The temperature profiles resulting from three different tolerances (0.01, 0.001 and 0.0001) can be seen in Fig. 11.

Next we turn to the efficiency of all methods and examine the average number of Newton iterations per time step shown in Fig. 12. Again, the original variants of the HI scheme need the most Newton iterations and are especially sensitive to temporal and spatial discretization. They are no longer considered in the remaining examinations. Figure 13 only compares iterations of the AC and tolerance-based HI methods. The iteration count increases with increasing time step size for all three methods but stays more or less constant over all spatial discretizations. In the rightmost graph showing the largest time step, the AC for the first time requires more iterations than the proposed tolerance-based HI methods.

Fig. 11
figure 11

Melting volume example: temperature profiles at \(0.5\,{t_f}\) obtained from the proposed tolerance-based HI scheme for different tolerances (100 elements, \(\Delta t=200{\hbox {s}}\))

Fig. 12
figure 12

Melting volume example: average number of iterations of all investigated methods for latent heat depending on spatial and temporal discretization

Fig. 13
figure 13

Melting volume example: average number of iterations of best performing methods for latent heat depending on spatial and temporal discretization

Fig. 14
figure 14

Melting volume example: temperature profiles at \(t_f\) obtained with AC methods for three time step sizes and with different phase change interval widths d. 100 elements

Fig. 15
figure 15

Geometry and moving laser heat source for single track scan example (all dimensions given in mm) [21]

This melting volume example reveals an already mentioned problem of AC methods, namely the possibility to neglect much of latent heat by stepping over or passing through the phase change interval too fast. The AC method is used with three widths \(d\in \{1{\hbox {K}},2{\hbox {K}},3{\hbox {K}}\}\) to compute artificial solidus and liquidus temperatures \(T_s = T_m -d\) and \(T_l=T_m+d\). The final temperature profiles are graphed in Fig. 14 for three time step sizes and the finest mesh with 100 finite elements. Obviously, the profiles differ in the respective phase change intervals which would not be the main concern in PBFAM. Instead focus lies on the maximum temperature predicted on the left edge. For the smallest time step \(\Delta t = 200\,{\hbox {s}}\) all widths reach almost the same maximum temperature. Increasing the step size leads to larger discrepancies in the maximum temperature. A smaller width d correlates with higher maximum temperatures which in turn implies that not all latent heat has been absorbed.

Summary of preliminary simulations: Taking into account all experience gathered from the two numerical examples, the authors recommend the use of the proposed tolerance-based HI method or an AC method. We found that the new criterion (48) to determine nodes that undergo phase change is superior to the one originally introduced by [43]. This is due to two reasons. First, the accuracy is user-controllable by setting a respective tolerance. Second, the new stopping criterion (48) typically leads to a significant reduction of nonlinear Newton iterations except for scenarios with very strict tolerances, which are, however, not expected to be necessary in the targeted application PBFAM.

Table 2 Material parameters for single track scan example [21]

Naturally, it is hard to predict a-priori which method will lead to less Newton iterations since this depends on the specific problem, tolerances and the (artificial) phase change interval width. Therefore, in the following, an actual example in the context of PBFAM will be investigated.

Single track scan

The following example simulates the scanning of a single track in a PBFAM process and was introduced in [20] and has also been simulated elsewhere, see [21]. A schematic sketch of the setup is shown in Fig. 15. A volumetric heat source described by Eq. (4) with effective power \(W=30\,{\hbox {W}}\) and size \(R=0.06\,{\hbox {mm}}\) is applied to the powder domain.

The geometry consists of a cuboid of \(0.6\times 0.2\times 0.2\, {\hbox {mm}^3}\). The top layer of 0.05 mm in z-direction is in powder form and rests on top of a solid substrate domain. The material is a 316L type steel with parameters summarized in Table 2. The whole domain is initialized at a fixed temperature \(T_0=303\,{\hbox {K}}\). All surfaces are insulating, only \(x=0.6\,{\hbox {mm}}\) is subject to an essential boundary condition \({\hat{T}}=T_0=303\,{\hbox {K}}\).

The laser beam center moves from an initial position at \(x=-0.06\,{\hbox {mm}}\) (one laser beam radius outside the domain) with a scanning speed of \(v=\nicefrac 120\,{\hbox {mms}}\) in x-direction along the symmetry plane \(y=0\). The powder layer is discretized by a regular hexahedral mesh with element size \(h_0= 0.0025\,{\hbox {mm}}\), i.e. \(n_\text {ele}^z=20\) elements across the powder layer height. In the substrate domain, a mesh with double element height in z-direction is applied. Moreover, an adaptive time stepping scheme is applied that halves the time step size if no convergence is achieved by the employed Newton–Raphson scheme (within a prescribed maximal number of iterations), and that doubles the step size again after four convergent time steps on the smaller step size level. As initial step size a value of \(\Delta t^{(0)} = 1\,{\upmu \hbox { s}}\) has been chosen, which has been verified to yield a sufficiently small time discretization error.

In prior simulations of this example, latent heat effects have been taken into account via an enthalpy method [20] and an isothermal HI method [21]. Here, we will use an AC method and subsequently the newly proposed tolerance-based isothermal HI scheme to simulate the process. An artificial melting interval is introduced for the AC method. As a baseline we chose \(T_s = 1600\,{\hbox {K}}\) and \(T_l = 1800\,{\hbox {K}}\), i.e., \(d=100\,{\hbox {K}}\). Isothermal HI is applied with a tolerance of \(\varepsilon _\text {tol}=0.001\). The isothermal HI scheme is used in combination with enthalpy-based parameter interpolation, while the AC method is used with temperature-based interpolation. In a first step, qualitative characteristics of the solution shall be discussed. After a short period of time the melt pool shape reaches a steady-state. Its geometry can be visualized by the isotherm \(T=T_m\). Figure 16 compares the results obtained with the AC and tolerance-based HI method to the results reported in [20, 21]. Both the AC and HI solution show good agreement with the reference. The melt pool dimensions and peak temperatures are compared quantitatively to the reference solutions in Table 3. All compared quantities show good agreement.

Fig. 16
figure 16

Surface temperature profile and melt pool shape in steady-state. Current results from HI scheme (top, left) and AC scheme (top, right) as well as reference results from Gusarov et al. [20] (bottom, left) and Hodge et al. [21] (bottom,right)

Table 3 Comparison of maximum temperature and melt pool dimensions resulting from the different latent heat models
Fig. 17
figure 17

Surface temperature profiles for AC and HI along laser path (i.e. \(y=0.0\,{\hbox {mm}}\)) for different spatial discretizations. Zoomed in to melting front of melt pool

Since no more quantitative data is provided by the reference papers, we compare AC and HI method with each other. In the preliminary examples in “Solidification front in a 1D slab” and “Melting volume in a 1D slab” sections a strong dependency on spatial discretization was recognized. Three additional, coarser spatial discretizations with elements of size \(2h_0\), \(4h_0\) and \(\frac{20}{3}h_0\) (which results in \(n_\text {ele}^z \in \{10,5,3\}\) elements over the powder layer height) are introduced to investigate this phenomenon for the single track scan. The accuracy of both methods can be judged by looking at characteristic temperature profiles in the steady-state. All results presented in the following are shown for (approximately) \(t=4.6\,{\upmu \hbox { s}}\), which is the same point in time for which the melt pool shape has been illustrated in Fig. 16. First, Fig. 17 shows the surface temperature profiles for all discretizations plotted along the laser path (i.e. \(y=0\)) in the vicinity of the melting front. This front is characterized by high temperature gradients and rates. With increasing element size, larger nonphysical oscillations in the temperature profile are observed for the AC scheme. These oscillations can be traced back to the limitation of the employed first-order finite elements in representing strong gradients and material nonlinearities, here mainly caused by the extreme gradients of the thermal conductivity across the phase interface. Employing finite elements based on higher order shape functions can remedy this numerical issue [48, 49].

No such oscillations occur with the HI scheme, which enforces the temperature in the phase transition region to lie within a temperature interval (implicitly) prescribed through the tolerance \(\varepsilon _\text {tol}\). Although the shape functions used with the HI scheme are still first-order, the reset of temperature to a consistent melt temperature as described in (45) seems to prohibit temperature oscillations in the phase interface region as observed for AC, even though the HI scheme performs phase change and parameter interpolation within a considerably smaller temperature interval. It is important to note that so far the oscillations have not been observed to cause stability issues (e.g. a significant energy increase in the discrete system) and they remain small compared to the overall temperature range.

A second phase change happens along the laser path (i.e. \(y=0.0\,{\hbox {mm}}\)) when material cools down again at the tail of the melt pool. Figure 18 gives a detailed view of the temperature profile in this region. The HI method produces a kink in the temperature profile at \(T_m\), which is to be expected for a phase interface. The AC method produces a mushy phase change region and no kink is observed. However, further away both temperature profiles are in good agreement again.

Fig. 18
figure 18

Surface temperature profiles for AC and HI along laser path (i.e. \(y=0.0\,{\hbox {mm}}\)) for different spatial discretizations. Zoomed in to solidification region in melt pool tail

Fig. 19
figure 19

Surface temperature profiles for AC and HI transverse to laser path at \(x=0.4\,{\hbox {mm}}\) showing transition from melt pool to powder

Fig. 20
figure 20

Surface temperature profiles for AC and HI transverse to laser path at \(x=0.2\,{\hbox {mm}}\) showing transition from solid to powder

Another aspect to investigate is the resulting temperature profile transverse to the laser scanning direction. Change from melt to powder can be investigated with a cut in y-direction at \(x=0.4\,{\hbox {mm}}\) as shown in Fig. 19. Again oscillations appear in the AC solution for coarser meshes but not in the HI solution. A similar picture emerges for a cut in y-direction through solid and powder, e.g. at \(x=0.2\,{\hbox {mm}}\) as illustrated in Fig. 20. In both transverse cuts temperature profiles differ for AC and HI method in proximity to \(T_m\) and agree well in some distance to the phase interface.

As suspected, in PBFAM applications the chosen method for latent heat only has a very local influence on the resulting temperature field, but does not significantly affect the global temperature characteristics from a rather macroscopic point of view. Another aspect of importance for PBFAM process simulations is numerical efficiency, here assessed in terms of nonlinear solver performance. Figure 21 depicts the total number of Newton iterations accumulated over the whole simulation time for the different spatial discretizations. The results of the preliminary numerical examples in the previous sections seem to transfer to a larger example: The HI method shows a strong dependency on the spatial resolution. The difference in the number of Newton iterations between AC and HI method is less pronounced in the practically relevant range of discretizations (e.g. three elements across the powder layer thickness).

Fig. 21
figure 21

Numerical efficiency of AC and HI in terms of total Newton iterations required

Fig. 22
figure 22

Surface temperature profiles for AC with different phase change intervals along laser path (i.e. \(y=0.0\,{\hbox {mm}}\)). Zoomed in to melting front

Fig. 23
figure 23

Surface temperature profiles for AC with different phase change intervals transverse to laser path at \(x=0.4\,{\hbox {mm}}\) showing transition from melt pool to powder

Fig. 24
figure 24

Surface temperature profiles for AC with different phase change intervals transverse to laser path at \(x=0.2\,{\hbox {mm}}\) showing transition from solid to powder

Fig. 25
figure 25

Total number of Newton iterations required for different initial time step sizes \(\Delta t^0\) in the time step halving scheme

In a next step, it shall be investigated how the accuracy and numerical efficiency of the two considered phase change schemes can be influenced by the respective numerical parameters d (phase change interval of the AC scheme) and \(\varepsilon _\text {tol}\) (tolerance of the HI scheme). In a first step, the phase change interval for the AC method is varied for the two coarsest spatial discretizations to see how it affects the temperature oscillation. It has to be noted that this also changes the temperature interval for temperature-based parameter interpolation. The intervals are given by \(T_s = T_m - d\) and \(T_l = T_m +d\), where \(d \in \{100,250,500\} {[K]}\). Moreover, two additional versions of HI with relaxed tolerances of \(\varepsilon _\text {tol}=0.01\) and \(\varepsilon _\text {tol}=0.1\) are simulated. The resulting surface temperatures are plotted in Figs. 22, 23 and 24 in the regions that showed oscillations in the previous plots. While the increased phase change interval for the AC decreases the amplitude of these oscillations by a certain amount, the overall accuracy of the temperature profiles decreases as well. Thus, it has to be concluded that for a given spatial discretization the width d of the phase change interval is not a suitable parameter to control the accuracy of AC schemes. The solutions obtained with different tolerances for the HI method show no oscillation. The solution from HI with tolerance \(\varepsilon _\text {tol}=0.01\) is indistinguishable from the one with stricter tolerance \(\varepsilon _\text {tol}=0.001\). The solution from HI with tolerance \(\varepsilon _\text {tol}=0.1\) deviates slightly from the ones with stricter tolerances in the solid-powder transition region, see Fig. 24. However, it still seems to be considerably more accurate than the solution from AC with \(d=500\,{\hbox {K}}\).

Finally, the initial time step size \(\Delta t^0\) for the three considered HI and AC variants is increased up to 64 times of the original value. The total number of Newton iterations that is now required is depicted in Fig. 25. It is emphasized that the time step halving scheme is still employed. It was optimized individually for each method to yield low iteration counts. Moreover, it has been checked that the time discretization error is still sufficiently small, i.e. there is no visible difference in the results for base time step sizes up to 16 times \(\Delta t^0\) for all variants. Higher step sizes lead to visible albeit small deviations in the temperature profiles, which are smaller, however, than the deviations resulting from the different HI and AC variants. According to Fig. 25, an increased phase change interval d in AC allows for larger step sizes and thus less iterations. Especially on the coarsest mesh, HI with a high tolerance (\(\varepsilon _\text {tol}=0.1\)) yields a comparable number of accumulated Newton iterations, but at a considerably increased accuracy, as compared to the AC method with the (unphysically) large phase change interval \(d=500\,{\hbox {K}}\). All in all, for a given spatial discretization, the proposed HI scheme allows to directly control numerical efficiency and accuracy by means of a user-defined tolerance. Within the considered scope of numerical examples and practically relevant spatial and temporal discretizations, this property made the novel tolerance-based HI scheme preferable as compared to the original HI scheme and the investigated AC method.

Conclusion

In this work, an extension of phase change and latent heat models for the simulation of metal powder bed fusion additive manufacturing processes on the macroscale has been proposed and different models have been compared with respect to accuracy and numerical efficiency. In this context, a systematic formulation of phase fraction variables has been proposed relying either on temperature- or enthalpy-based interpolation schemes. Moreover, latent heat has been considered either by means of an apparent capacity (AC) or heat integration (HI) method. For the latter, a novel phase change criterion has been proposed, which combines superior accuracy and numerical efficiency (in terms of an improved nonlinear solver performance allowing for larger time step sizes and fewer iterations per time step) as compared to the original HI scheme. Compared to the AC approach, the numerical efficiency of the proposed tolerance-based HI scheme is comparable while offering an increased level of accuracy. Numerical results from the literature have been reproduced, which shows the validity of the proposed scheme. In summary, both the AC and the proposed tolerance-based HI scheme perform well when considering the accuracy requirements as well as practically relevant spatial and temporal discretization resolutions for PBFAM process simulation. Specifically, global temperature characteristics, such as the peak temperature, can be accurately captured with both methods. Still, the authors believe that the new tolerance-based HI method is advantageous over AC schemes due to the user-controllable tolerance, which allows to directly control numerical efficiency and accuracy of the scheme, and which can directly be interpreted as the error in latent heat made during a phase change process.

For part-scale PBFAM models thermo-mechanical interaction is one of the primary interests. Structural parameters such as Young’s modulus or the thermal expansion coefficient may depend upon the temperature history in a similar fashion as proposed in the current work for the thermal parameters. Future research work will focus on the question of how the proposed methods for latent heat and parameter interpolation will behave in the thermo-mechanically coupled scenario.

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 analysed during the current study are available from the corresponding author on reasonable request.

References

  1. Gibson I, Rosen D, Stucker B. Additive manufacturing technologies. New York: Springer; 2015. https://doi.org/10.1007/978-1-4939-2113-3.

    Google Scholar 

  2. Herzog D, Seyda V, Wycisk E, Emmelmann C. Additive manufacturing of metals. Acta Mater. 2016;117:371–92. https://doi.org/10.1016/j.actamat.2016.07.019.

    Google Scholar 

  3. Khairallah SA, Anderson A. Mesoscopic simulation model of selective laser melting of stainless steel powder. J Mater Process Technol. 2014;214(11):2627–36.

    Google Scholar 

  4. Körner C, Attar E, Heinl P. Mesoscopic simulation of selective beam melting processes. J Mater Process Technol. 2011;211(6):978–87.

    Google Scholar 

  5. 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.

    MathSciNet  MATH  Google Scholar 

  6. 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.

    Google Scholar 

  7. Russell MA, Souto-Iglesias A, Zohdi TI. Numerical simulation of Laser Fusion Additive Manufacturing processes using the SPH method. Comput Methods Appl Mech Eng. 2018;341:163–87. https://doi.org/10.1016/j.cma.2018.06.033.

    MathSciNet  MATH  Google Scholar 

  8. Yan W, Qian Y, Ge W, Lin S, Liu WK, Lin F, et al. Meso-scale modeling of multiple-layer fabrication process in selective electron beam melting: inter-layer/track voids formation. Mater Des. 2018;141:210–9.

    Google Scholar 

  9. Herbold EB, Walton O, Homel MA. Simulation of powder layer deposition in additive manufacturing processes using the discrete element method. Livermore: Lawrence Livermore National Laboratory (LLNL); 2015. http://www.osti.gov/servlets/purl/1239200/.

  10. Meier C, Weissbach R, Weinberg J, Wall WA, John Hart A. Modeling and characterization of cohesion in fine metal powders with a focus on additive manufacturing process simulations. Powder Technol. 2019;343:855–66.

    Google Scholar 

  11. Meier C, Weissbach R, Weinberg J, Wall WA, Hart AJ. Critical influences of particle size and adhesion on the powder layer uniformity in metal additive manufacturing. J Mater Process Technol. 2019;266:484–501.

    Google Scholar 

  12. Zhang J, Liou F, Seufzer W, Newkirk J, Fan Z, Liu H, et al. Probabilistic simulation of solidification microstructure evolution during laser-based metal deposition. In: Proceedings of the solid freeform fabrication symposium; 2013. p. 739–48. https://sffsymposium.engr.utexas.edu/Manuscripts/2013/2013-59-Zhang.pdf.

  13. Gong X, Chou K. Phase-field modeling of microstructure evolution in electron beam additive manufacturing. JOM. 2015;67(5):1176–82. https://doi.org/10.1007/s11837-015-1352-5.

    Google Scholar 

  14. 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.

    Google Scholar 

  15. Salsi E, Chiumenti M, Cervera M. Modeling of microstructure evolution of Ti6Al4V for additive manufacturing. Metals. 2018;8(8):633. http://www.mdpi.com/2075-4701/8/8/633.

    Google Scholar 

  16. Lindgren LE, 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.

    Google Scholar 

  17. 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. Ann Rev Heat Transf. 2017;20(1):241–316. http://www.dl.begellhouse.com/references/5756967540dd1b03,562e7b3835dec96e,0860d4f32b9f248d.html.

    Google Scholar 

  18. Gusarov AV, Kruth JP. Modelling of radiation transfer in metallic powders at laser treatment. Int J Heat Mass Transf. 2005;48(16):3423–34.

    MATH  Google Scholar 

  19. Gusarov AV. Homogenization of radiation transfer in two-phase media with irregular phase boundaries. Phys Rev B Condens Matter Mater Phys. 2008;77(14):1–14.

    Google Scholar 

  20. Gusarov AV, 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):072101.

    Google Scholar 

  21. Hodge NE, Ferencz RM, Solberg JM. Implementation of a thermomechanical model for the simulation of selective laser melting. Comput Mech. 2014;54(1):33–51. https://doi.org/10.1007/s00466-014-1024-2.

    MathSciNet  Google Scholar 

  22. Hodge NE, Ferencz RM, Vignes RM. Experimental comparison of residual stresses for a thermomechanical model for the simulation of selective laser melting. Addit Manuf. 2016;12:159–68.

    Google Scholar 

  23. 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. https://doi.org/10.1177/0954405414539494.

    Google Scholar 

  24. Cervera GBM, Lombera G. Numerical prediction of temperature and density distributions in selective laser sintering processes. Rapid Prototyp J. 1999;5(1):21–6. https://doi.org/10.1108/13552549910251846.

    Google Scholar 

  25. Childs THG, Hauser G, 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.

    Google Scholar 

  26. Shen N, Chou K. Numerical thermal analysis in electron beam additive manufacturing with preheating effects. In: Proceedings of the 23rd solid freeform fabrication symposium. 2012. p. 774–84.

  27. 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. https://doi.org/10.1007/s40192-019-00132-9.

    Google Scholar 

  28. Roy S, Juha M, Shephard MS, Maniatty AM. Heat transfer model and finite element formulation for simulation of selective laser melting. Comput Mech. 2018;62(3):273–84.

    MathSciNet  MATH  Google Scholar 

  29. 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.

    MathSciNet  Google Scholar 

  30. Jamshidinia M, Kong F, Kovacevic R. Temperature distribution and fluid flow modeling of electron beam melting® (EBM). In: Volume 7: lrts A, B, C, and D. ASME; 2013. p. 3089. http://proceedings.asmedigitalcollection.asme.org/proceeding.aspx?doi=10.1115/IMECE2012-89440.

  31. Chang B, Allen C, Blackburn J, Hilton P, Du D. Fluid flow characteristics and porosity behavior in full penetration laser welding of a titanium alloy. Metall Mater Trans B. 2015;46(2):906–18. https://doi.org/10.1007/s11663-014-0242-5.

    Google Scholar 

  32. Geiger M, Leitz KH, Koch H, Otto A. A 3D transient model of keyhole and melt pool dynamics in laser beam welding applied to the joining of zinc coated sheets. Prod Eng. 2009;3(2):127–36.

    Google Scholar 

  33. Rai R, Burgardt P, Milewski JO, Lienert TJ, Debroy T. Heat transfer and fluid flow during electron beam welding of 21Cr-6Ni-9Mn steel and Ti-6Al-4V alloy. J Phys D Appl Phys. 2009;42(2):025503.

    Google Scholar 

  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.

    Google Scholar 

  36. Denlinger ER, Irwin J, Michaleris P. Thermomechanical modeling of additive manufacturing large parts. J Manuf Sci Eng. 2014;136(6):061007. https://doi.org/10.1115/1.4028669.

    Google Scholar 

  37. Zhang Y, Guillemot G, Bernacki M, Bellet M. Macroscopic thermal finite element modeling of additive metal manufacturing by selective laser melting process. Comput Methods Appl Mech Eng. 2018;331:514–35.

    MathSciNet  MATH  Google Scholar 

  38. Neiva E, Badia S, Martín AF, Chiumenti M. A scalable parallel finite element framework for growing geometries. Application to metal additive manufacturing. Int J Numer Methods Eng. 2019;. https://doi.org/10.1002/nme.6085.

    MathSciNet  Google Scholar 

  39. Zaeh MF, Branner G. Investigations on residual stresses and deformations in selective laser melting. Prod Eng. 2010;4(1):35–45. https://doi.org/10.1007/s11740-009-0192-y.

    Google Scholar 

  40. Comini G, Del Guidice S, Lewis RW, Zienkiewicz OC. Finite element solution of non-linear heat conduction problems with special reference to phase change. Int J Numer Methods Eng. 1974;8(3):613–24. https://doi.org/10.1002/nme.1620080314.

    MATH  Google Scholar 

  41. de Moraes D, Czekanski A. Parametric thermal FE analysis on the laser power input and powder effective thermal conductivity during selective laser melting of SS304L. J Manuf Mater Process. 2018;2(3):47.

    Google Scholar 

  42. Morgan K, Lewis RW, Zienkiewicz OC. An improved algorithm for heat conduction problems with phase change. Int J Numer Methods Eng. 1978;12(7):1191–5. https://doi.org/10.1002/nme.1620120710.

    Google Scholar 

  43. Rolph WD, Bathe KJ. An efficient algorithm for analysis of nonlinear heat transfer with phase changes. Int J Numer Methods Eng. 1982;18(1):119–34. https://doi.org/10.1002/nme.1620180111.

    MATH  Google Scholar 

  44. Oliveira B, Afonso JC, Zlotnik S. A Lagrangian–Eulerian finite element algorithm for advection-diffusion-reaction problems with phase change. Comput Methods Appl Mech Eng. 2016;300:375–401.

    MathSciNet  MATH  Google Scholar 

  45. Hu H, Argyropoulos SA. Mathematical modelling of solidification and melting: a review. Model Simul Mater Sci Eng. 1996;4(4):371–96.

    Google Scholar 

  46. Laursen TA. Computational contact and impact mechanics. Berlin: Springer; 2002.

    MATH  Google Scholar 

  47. Wall WA, Popp A, Kronbichler M, Mayr M, Meier C, Vuong AT, et al. BACI: A multiphysics simulation environment. Institute for Computational Mechanics, Technische Universität München. 2019.

  48. O’Hara P, Duarte CA, Eason T. Transient analysis of sharp thermal gradients using coarse finite element meshes. Comput Methods Appl Mech Eng. 2011;200(5–8):812–29. https://doi.org/10.1016/j.cma.2010.10.005.

    MathSciNet  MATH  Google Scholar 

  49. Kollmannsberger S, Özcan A, Carraturo M, Zander N, Rank E. A hierarchical computational model for moving thermal loads and phase changes with applications to selective laser melting. Comput Math Appl. 2018;75(5):1483–97. https://doi.org/10.1016/j.camwa.2017.11.014.

    MathSciNet  MATH  Google Scholar 

Download references

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.

Funding

Not applicable.

Author information

Authors and Affiliations

Authors

Contributions

All authors contributed to the derivation of model equations, the discussion of results, and writing the manuscript. In addition, SP conducted the specific code implementation and the shown numerical simulations. CM and WAW worked out the general conception of the proposed modeling approach. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Christoph Meier.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Proell, S.D., Wall, W.A. & Meier, C. On phase change and latent heat models in metal additive manufacturing process simulation. Adv. Model. and Simul. in Eng. Sci. 7, 24 (2020). https://doi.org/10.1186/s40323-020-00158-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40323-020-00158-1

Keywords