Reduced order modeling of selective laser melting: from calibration to parametric part distortion

Additive manufacturing is an appealing solution to produce geometrically complex parts, difficult to manufacture using traditional technologies. The extreme process conditions, in particular the high temperature, complex interactions and couplings, and rich metallurgical transformations that this process entails, are at the origin of numerous process defects. Therefore, the numerical simulation of the process is gaining the interest of both the scientific and the industrial communities. However, simulating that process demands impressive computational resources, limiting high resolution simulations to the microscopic and mesoscopic scales. This paper proposes a thermo-mechanical modeling framework at the process scale as well as its associated reduced order simulation counterpart, enabling the parametric evaluation of the part distortion. It deeply addresses the process calibration using a high-resolution computational procedure based on the use of an in-plane-out-of-plane separated representation at the heart of the so-called Proper Generalized Decomposition (PGD), as well as the analysis of the transient thermal effects, defining the conditions in which the thermal and mechanical analyses can be decoupled.


Introduction
In the past decade, additive manufacturing attracted interest because of its benefits and flexibility. It made possible the production of designs that were impossible to manufacture using traditional manufacturing techniques and opened the door to the advancement in topology optimization and design efficiency enhancement. Additive manufacturing is a wide technology with many variants. Among them, Selective Laser Melting-SLM-also known as Direct Metal Laser Sintering, that proceeds on a powder bed [1].
SLM consists of a concentrated laser that hits a metal powder layer and melts it according to the slicing of a CAD model. The powder is heated above the liquidus temperature of the metal hence producing parts with properties comparable to those of metal bulk, except in what concerns ductility.
Chady Ghnatios cghnatios@ndu.edu.lb Extended author information available on the last page of the article.
However, due to the extreme process conditions some defects occur, such as cracking [2] and part distortion [3] induced by installed residual stresses [4]. In order to account for these defects, scientists have tried to quantify them using either experimental methods as well as numerical simulation.
Experimentally, two types of measurements are used: with and without contact. In absence of contact, ultrasounds, thermal imaging such as infrared -IR-sensors [5] or X-ray tomography to quantify the part defects such as porosity and cracks [6], are widely considered, among many other techniques. This type of techniques can only inspect external exposed surfaces or its neighborhood inside the material depending on their penetration capability. Contactbased techniques use thermocouple joints to measure the temperature of the substrate (bottom surface or bed) [7] or other sensors enabling accessing to the installed stresses.
To circumvent the shortcomings of the experimental approach, researchers resorted to numerical techniques in an attempt to simulate the printing process. Among the numerical techniques, the finite elements analysis -FEMis mostly used to simulate the temperature field in additive manufacturing processes [8,9]. A thermal simulation is then coupled to a mechanical simulation in order to predict the distortion of 3D printed parts. Simulation concerns three different scales: Micro, Meso and Macro.
Microscale analysis focuses on the rich local physics including the interaction between the heat source and the substrate, the phase-change and the associated metallurgical transformations, as well as the localized thermal gradients. The fast cooling and the induced stresses along the printing path is usually described at the mesoscopic scale. These two approaches are usually coupled [10,11] for reproducing quite well experimental findings. However, the computational cost makes impossible envisaging its application at the part level (macroscopic scale).
In order to validate the numerical results, several approaches are used depending on the considered process. For SLM, the so-called "Twin Cantilever" test shown in Fig. 1 is usually considered. In this experiment a T-shaped part is printed on a support and the tip deflection after removing the support serves to calibrate and/or validate.
However in most additive manufacturing processes, distortion is difficult to predict because of the change of the mechanical properties during the process [12]. The high temperatures cause the microstructure to be altered, and the cyclic heating may cause for some materials a thermal fatigue leading to the material embrittlement [13]. These changes are difficult to model and take into account sincethey depend on the printing parameters such as the laser power and speed, the material itself, the powder characteristics, the heating spot size, ... affecting not only the microstructure but also the interlayer cohesion due to the varying porosity [6].
In these validation experiments, the printing parameters are varied and the distortion is studied. A quite good agreement between the experimental and numerical results are noticed when addressing the thermal modeling. Whereas in what concerns the distortion the error varies depending on the process parameters and materials employed, ranging from 10% to a total mismatch [14,15].
The present work aims at proposing a thermal simulation of SLM while studying the reheating effect when the heat source re-passes close to a solidified element. In addition, it will present a novel approach for simulating the distortion of a 3D printed part. This approach is used to generate a calibration model from which the stresses will be extracted and used to perform a macroscale simulation that gives the distortion as a function of the temperature reached by the molten metal and of the printing trajectories appropriately parametrized, in real time by using the non-intrusive Proper Generalized Decomposition approach [16,17], following the rationale initially introduced in [18] where the interested reader can find more references on the advanced simulation of those kind of processes.
The next section defines and implement the thermal model and analyses the effect of the thermal history in the transient temperature field. Section "Thermo-mechanical modelling" establishes a thermomechanical approach, enabling first the calculation of the installed residual stresses and then, its integration in a parametric solution of the manufactured part distortion.

Transient simulation of a moving heat source
This section discusses the finite element thermal modeling associated with a moving heat source hitting a powder bed and moving at a velocity V, as encountered in SLM.

Heat equation
The heat transfer mechanisms in SLM are conduction and natural convection as shown in Fig. 2. Note that radiation can be neglected since metallic materials have an excellent heat conduction ability, hence most of the heat will be diffused by conduction to the surrounding material [19]. In the reduced order modeling that we develops in the present work, the laser-matter interaction is not taken into account, only the effective heat transferred to the material is considered.
The transient heat equation that governs the temperature distribution T (x, y, z) reads where q latent is the latent heat, i.e. the heat released by the phase change, Cp refers to the heat capacity and k to the thermal conductivity.
The heat source that acts as a boundary condition, q(r), is considered having a gaussian distribution [19] and hence q source is expressed by in which P is the laser power (W ), r 0 is the laser spot radius and r is the distance from the spot center. The initial temperature is given by T (x, t = 0) = 293K and the free convection boundary condition reads q = h(T − T ∞ ), with h the exchange coefficient.
The implicit semi-discrete formulation reads (3) in which T n+1 is the temperature field at time step n + 1, i.e. t = (n + 1)Δt, with T 0 = 293K, and q = q latent /ρC p is assumed temperature dependent according to [20] where T l and T s are respectively the liquidus and the solidus temperatures, ΔH is the latent heat and f s the solid fraction. The value of ΔH is taken from the literature [20] assuming that the solidification in additive manufacturing is similar to that in casting, even if other expressions considered. The previous equation shows that there is no heat release when cooling in the liquid state or solid state, but only when there is a solid-liquid mixture.
Defining the solid fraction f s from [21] its temperature derivative involved in Eq. 4, df s dT , reads [21] df s dT In order to simulate the moving heat source, at each time step dt the laser moves a distance dx = Vdt. Hence the FEM discrete problem needs at each time step n, the solution of the linear problem expressed from where coefficient matrix at time t n , K n , depends on the temperature T n because of the nonlinearity.

Thermal properties
The thermal conductivity is a temperature dependent function and hence it must adapt to the temperature change at each iteration. The used material in the simulation is the 316L stainless steel. In [22] authors found correlations for the thermal conductivity and the temperature where T is the temperature in (K) and a, b, c and d are constants corresponding to iron, reported in Table 1.
The thermal conductivity is assumed to be the one of the bulk material expressed by Eq. 8 everywhere except in the powder layer in which it is reduced of one order of magnitude for taking into account the effects of porosity in the powder bed. The heat capacity can be expressed from the mixture rule

Results and discussion
Numerical simulations are performed by using our own finite element method implementation. For that purpose, the mesh consists of 14553 hexahedral elements (with 66 000 degrees of freedom). The problem is solved using an incremental implicit time integration scheme with a timestep of 1.25 ms. The converged solution (with resperct to the mesh size) using the parameters reported in Table 1, is visualized in the Paraview open source software.
Using these values, the simulation is performed in a thin 3D plate of length 15 mm, width 5 mm and thickness 1 mm, where a temperature snapshot is depicted in Fig. 3.
The objective is to obtain the characteristic temperature time evolution in a representative node.

Heating
When studying the maximum temperature at each point along the laser path, it results the profile depicted in Fig. 4. This result clearly proves that except at the beginning, all the points on the laser trajectory will reach almost the same Temperature field at a given time t, with the source moving from the left boundary towards the right one

Reheating
Printing implies considering close trajectories as depicted in Fig. 6. When the heat source comes close to a certain point during the next printing path, this point will be reheated, effect known as neighboring effect.
As shown in Fig. 5, the characteristic cooling time in order to reach a temperature verifying the condition is of about 44 ms, time that is referred as critical time t cr . Thus, if the time elapsed before coming back to the vicinity of a certain point is higher than that critical time t cr , this latter has cooled down to the almost ambient temperature T amb and it will experience the same reheating as any point in the neighbor path-line, as the marked point in Fig. 7.  In these circumstances the temperature evolution of the referred point during its reheating is representative of the reheating of all the points located in the contiguous path-line, and consequently all them experience the same maximum temperature (obviously lower than the one reached when the laser applied directly on them), except at the beginning of the heating path-lines, as shown in Fig. 8, with the associated representative cooling profile given in Fig. 9.
Thus, if the elapsed time between two passes is higher than t cr , the previous analysis proves that all the points start their reheating from an almost ambient temperature and thus will reach almost the same maximum temperature when reheated, and will follow a similar cooling.
However, when the elapsed time between two passes is lower than t cr , reheating applies on a still hot initial state. This situation applies when a printing path starts close to the end of the previous one as illustrated Fig. 10 or when the  The simulation is done on the printing path shown in Fig. 10 in which x = 0 represents the start of the second path. The maximum temperature produced by the reheating exhibits values much higher than the ones obtained when reheating proceeds from the almost ambient temperature. Figure 11 proves that this deviation only applied in the vicinity of x = 0 where the end of one path coincided with the beginning of the next one.
If we assume that the printing occurs as shown in Fig. 10, t cr /2 = 22 ms corresponds about 2 mm for V = 0.1m/s, then points x > 2 mm will be reheated from the ambient temperature, in agreement with Fig. 11 that clearly states that 2 mm suffice for reaching the temperature characteristic of a reheating from the ambient temperature.  The previous analysis allows neglecting the thermal heterogeneities induced by elapsed times lower than t cr as soon as the printed part are large enough.

Parametric distortion
In what follows we assume a thermal load ΔT = T max − T amb , with T max discussed in the previous section, and the deposition trajectories parametrized as proposed in [18] where the trajectories coincide with the curves Φ = cte, with Φ(x, y) solution of the problem ΔΦ(x, y) = μ 5 (11) with (x, y) ∈ ω = [0, L] × [0, H ], and where the printed part Ω ⊂ ω. Equation 11 is solved with the boundary condition expressed from In the previous expressions μ i , i = 1, . . . , 5, represent the model parameters, the first four parametrizing the Dirichlet boundary conditions, and μ 5 representing the source term in problem (11).
In [18] authors assume a local residual tension aligned with the deposition trajectory. This means that at each point, if one considers a frame aligned with the deposition trajectory, the one consisting of x 0 and y 0 depicted in Fig. 12, the resulting tensor σ 0 | (x 0 ,y 0 ) will have a single non-zero component.
After expressing the residual stress at each point by using the laboratory frame, and parametrically with respect to μ i , i = 1, . . . , 5, the distortion results from the virtual work principle or by using the elastic behavior with C the fourth order elasticity tensor, and the residual tension scaling with ΔT as discussed in [18].

Thermo-mechanical coupled model
In this section we will present a new macroscale approach to simulate the distortion of a part processed by additive manufacturing. In this approach the cooling of the deposited metal will be divided into two stages: T > T melting and T < T melting .
In the first stage, the metal is in a liquid phase and no residual stresses are induced because the liquid metal adapts to the substrate deformation that will expand when heated from the ambient temperature. Then, in the second stage, the metal solidifies and hence sticks to the substrate Fig. 13 Simulation flowchart for the calculation of both displacement components d 1 and d 2 contribution to the total one d surface. The whole system (deposited later and substrate) continues the cooling from the temperature existing when the deposited layer started its solidification, until reaching the final ambient temperature. The flowchart is presented in Fig. 13. Predictions were validated experimentally as described in the next section.

Experimental validation
The performed experiment consists of evaluating the distortion of a cantilevered plate with one deposit line on its top surface. The geometry is shown in Fig. 14.
The experiment is performed using Laser Metal Deposition -LMD-, with the deposit dimensions being: width 0.8 mm, thickness 175 μm and the substrate thickness 2 mm with its left boundary clamped (all the other remaining free). In addition the material used is stainless steel 316L. The maximum plate deflection is measured after the deposition process finished.
For the thermal simulation, the transient heat transfer equation is solved with a prescribed initial temperature in the deposit layer T (t = 0) = 2000K and T (t = 0) = 293K elsewhere. As discussed in the previous section, the thermal simulation provides the temperature field in the substrate when the deposit-substrate interface reaches the melting temperature T s . The substrate heating produces the substrate plate deflection depicted in Fig. 15.
When the metal begins to solidify, both thermal strains, in the substrate and the deposited layer, apply, and result in the deformed part depicted in Fig. 16, that added to the pre-deformed configuration (Fig. 15) results in the final configuration depicted in Fig. 17.
The quantity of interest is the maximum deflection at the tip of the plate. When comparing the simulated results with the experimental measurements of the maximum tip deflection, an excellent agreement was observed. Note that the modulus of elasticity of the deposit is taken to be slightly lower than that of the substrate due to the porosities and phase change during solidification (E deposit = 175 GPa and E substrate = 200 GPa) [13,23].

Calibration model for the macroscale simulation
In order to perform fast simulations on large scale models, a high-resolution calibration model is proposed from which the residual stresses to be used as local stresses in the whole part (to proceed as discussed in Section "Parametric distortion") are extracted.
For addressing the solution of the associated thermomechanical problems in the extremely degenerated domains    with thicknesses orders of magnitude much smaller that the other in-plane characteristic dimensions, the in-planeout-of-plane separated representation is considered. The last expresses the displacement field as a finite sum of terms, each involving a function of the in-plane coordinates and the second depending on the thickness coordinate. Thus, extremely rich out-of-plane discretizations can be considered to capture all the 3D effects event in such very degenerated domains, without impacting the computational efficiency. That separated representation expresses the displacement field according to where • refers to the Hadamard product and P i and T i are vectors involving respectively the in-plane and out-of-plane coordinates.  The calibration model geometry consists of three deposited layers with three printed lines per layer as shown in Fig. 18, all of them along direction y. In this model we assume that the deposit of interest, highlighted in red in Fig. 18 of 0.8 mm width and 200 μm thickness, is a repre sentative one accounting for the effects of the neighboring deposits and hence its stress state can be used at any point in a simulated part.
The simulation is performed by adding a deposit after the other to accumulate the distortions and the installed stresses. In the deformed model shown in Fig. 19 we can see that the highest distortion is observed between the tracks and between the layers and mainly on the edges. In addition, in the middle of the representative domain the distortions reached an almost uniform value, meaning that the deposit length is sufficient to extract the installed residual stresses.
From the computed strains we can then calculate the stresses and also the Von Mises effective stress. Figure 20 proves that there is a stress concentration on the interface between the layers which is a problem widely encountered in additive manufacturing [24]. Figure 21 extracts the stresses in the deposit of interest. The plot shows a region in the middle of the domain in Fig. 26 Resulting concentric circular trajectories which stresses are almost uniform and all in-plane stresses vanish. These stresses will be used as local stresses to enrich the parametric model proposed in [18] whose assumption (residual tension aligned with the printing path-line) seems too simplistic when compared with the stresses reported in Fig. 21.
To better appreciate the high-resolution reached by the in-plane-out-of-plane separated representation, Figs. 22 and 23 zoom-in respectively the deformation and the associated equivalent mesh resolution. It is important to note that few minutes in a standard laptop allowed for a resolution equivalent to 2.5 millions of 3D finite elements.

Parametric macroscopic modeling of part distortion
In this section we will present the application of the proposed methodology for the simulation of the distortion of a barrel.

Simulated geometry and parametrized trajectories
The simulated part geometry is shown in Fig. 24. The part is made of Titanium Grade 5 and printed on a Titanium substrate. The cross sectional area being circular we calculate the solution of the Poisson equation (11) in a slightly bigger rectangular domain. Then, the solution will be projected from the rectangular embedding domain onto the part.
In order to obtain circular trajectories, Eq. 11 is solved with μ 5 = 1 enforcing Φ = 0 outside the part domain, by using a penalty formulation. Figure 25 shows the penalized domain and the Poisson equation solution.
Hence, after calculating the solution, curves Φ = cte define the trajectories shown in Fig. 26.
The technique proposed in [18] is extended by incorporating the richer stress state that is parametrized by the deposition temperature T . Thus, a parametric solution representing the part distortion u(μ 1 , μ 2 , μ 3 , μ 4 , μ 5 , T ) is computed. The deposition temperature can be easily correlated with the process parameters: deposition velocity and laser power. For computing such a parametric solution several runs (about one thousand) for different values of the parameters were performed and the computed results combined for defining the solution for any other choice by interpolation [25].
Moreover, to avoid axial distortion occurring when all trajectories start at the same location, the starting positions will be randomly chosen in each layer. Figure 27 shows the difference between the two cases and it reveals the benefits of randomly choosing the starting point in each layer. Figure 28 shows the distortion depending on the temperature of the deposited layer in the case of concentric trajectories with random starting position at each layer. It can be noticed that as the temperature increases from near liquidus temperature the distortion increases up to a certain temperature after which the distortion starts to decrease again. This behavior can be explained by resorting the consideration discussed when addressing the calibration model, since as the temperature increases the gap between the liquidus temperature and the maximum temperature increases, hence the effect of the pre-heating induced deformation can compensate a cooling induced distortion. Thus, increasing the deposit temperature could be beneficial, however there are other many issues, out of the scope of the present work, pointing in the opposite direction. Figure 29 shows the distortion in the case of horizontal printing (μ 1 = 1 and μ 2 = 1 all the others vanishing). In this case the distortion is maximal on the right and left parts of the domain. Figure 30 shows the distortion in the case of vertical printing (μ 1 = 1 and μ 4 = 1, all the other vanishing). In this case the distortion is maximal on the top and bottom of the domain.
Using this parametric solution one can modify the trajectory or process parameters and evaluate their impact in almost real-time.
The methodology performs in larger parts without major difficulties, as the one illustrated in Fig. 31 consisting of a cylinder of inner diameter 35 mm, thickness 10 mm and height 65 mm.

Conclusion
This paper started by discussing transient thermal simulations of SLM. The reheating effect when the heat source re-passes next to an already printed point of a previous pass was quantified. The analysis proved that if the metal had enough time to cool down to reach the ambient temperature, the heat source will reheat it to an almost representative temperature ensuring a certain uniformity of the mechanical state all along the part. Then a thermo-mechanical coupled approach for calculating the distortion in a 3D printed part was presented. This model was used to propose a calibration procedure in order to extract the characteristic stresses to be used for calculating the distortion at the part level.
Finally, several simulations were done by varying the trajectory and the temperature of the deposited material to define the parametric part distortion with respect to the process parameters.