2D linear finite element simulation of laser metal heating for digital twins

In the context of laser-based additive manufacturing, the thermal behavior of the substrate is relevant to define process parameters vis-à-vis piece quality. The existing literature focuses on two process variables: (a) lumped laser power and (b) process speed. However, this literature does not consider other variables, such as those related to the laser power distribution. To fill this vacuum, this manuscript includes the laser power spatial distributions (Gaussian, uniform circular and uniform rectangular) in addition to (a) and (b) above in 2D linear substrate heating simulations. The laser energy is modeled as a time dependent heat flux boundary condition on top of the domain. The total laser delivered power was identical for all spatial distributions. The results show that the laser intensity spatial distribution strongly affects the maximum temperature, and the depth and width of the heat affected zone. These 2D finite element simulations prove to be good options for digital twin based design environments, due to their simplicity and reasonable temperature error, compared to non-linear analysis (considered as ground truth for this case). Future publications address non-linear finite element simulations of the laser heating process (including convection and radiation and temperature dependent substrate properties).


Introduction
Metal Additive Manufacturing (AM) has enabled the fabrication of complex geometries that could not be build using traditional manufacturing techniques [1]. Laserbased AM has also grown up because of its applications in repair, reconditioning, coating and remanufacturing of high-valued industrial pieces [2]. However, the characterization of the laser-based AM is still a matter of research.
In this manuscript, we present an analysis of the influence of the laser intensity distribution, laser radius and process speed on the thermal behavior of the substrate. The analysis is carried out via numerical simulations of a 2D thermal model using the finite element method. The energy contributions of the laser into the substrate are modeled as time dependent heat flux (or Neumann) boundary conditions. We study the effects of three types of laser intensity distributions: Gaussian, uniform circular and uniform rectangular. In our simulations, we do not consider the addition of material, the phase change nor non-linear phenomena.
The remainder of this article is organized as follows: Section 2 provides a review of the relevant related work. Section 3 describes the governing equations, numerical scheme, and the materials used for the simulations. Section 4 presents and discusses the results obtained. Section 5 concludes the manuscript and makes suggestions for future work. In recent years, there has been an increasing amount of literature on the study of laser-based AM. Several studies have shown that the temperature history has a significant impact in the quality and mechanical properties of the parts manufactured via laser-based AM [3][4][5][6].
Physical experimental studies have been executed to assess the influence of the process parameters on the geometry of parts produced by laser metal deposition (LMD). Most of the research has focused on these three parameters: laser power, process speed and material feeding rate. References [7][8][9][10] analyzed the effects of the process speed, laser power and material feeding rate on the geometry of the melt-pool (with, depth and burn-in shape) and on the dilution ratio. References [10,11] studied the impact of these three process parameters in the clad geometry (width, height and angle of repose).
Numerical models have also been developed for this purpose. Pure thermal [7,10] and thermal-fluid [9,12] were implemented to model the influence of the laser power, process speed and material feeding rate on the geometry and thermal history of the melt-pool and the clad bead. In all of these works, the power intensity distribution of the laser beam was modeled using a Gaussian function. However, other beam spots shapes have been used in physical contexts (e.g. rectangular beam spot [13]).

Digital twins in additive manufacturing
In the context of Industry 4.0, the concept of digital twin is key in the manufacturing environment [14]. In few words, a digital twin is a virtual representation of a real system. This virtual representation resembles as much as possible the real system, allowing the knowledge transfer between the real and cyber-physical worlds [15]. Therefore, simulation over the virtual entity plays a major role, since the gained information can be fed into the real system.
Digital twins for AM are still at early stages of development. Some research has focused on the identification of the main features that should be included in the cyber-physical world (e.g. thermal behavior, melt-pool dynamics, distortions, geometry prediction) [16,17]. However, one of the major limitations of the current models is the large amount of computational resources associated to their use, which make them impractical for real-time or interactive applications [16].
Reference [18] presented an approach that integrated real data obtained by sensors and theoretical results in the context of smart manufacturing. Reference [17] used the prior approach to develop a model for AM that used numerical simulations and real data to detect manufacturing defects and deviations. Reference [19] presented a digital twin for AM that integrated some main features related to resemble the real process (e.g. obtained geometry, temperature history, cooling rates, microstructure). The main contribution was the assembly of those features, which are are commonly studied independently. A further review on the related works on digital twins for AM can be found in [20].

Conclusions of the literature review
In the currently existing literature, metal laser heating is addressed by considering, in addition to the material properties, lumped laser power input at the material boundary and laser speed. The spatial distribution of the laser power follows a Gaussian profile. This manuscript addresses the role of the laser power spatial distribution on the temperature field at the metal substrate. This temperature distribution is located at the substrate cross section normal to the laser trajectory. Our linear 2D initiative is obviously less precise than the 2D non-linear counterparts. However, we contend that it has value for approximate simplified purposes, e.g. digital twin applications, which require a reasonable appraisal of the substrate temperatures, at early design-stages.

Problem description
Laser Metal Deposition (LMD) is a manufacturing process in which a high-power laser beam is used to melt a metallic material while it is being deposited on the surface of a metallic substrate.
In this work, we aim to analyze the thermal behavior of an LMD process consisting of several parallel linear deposition tracks (see Fig. 1). Along each track, we assume that the laser speed, laser power and material feeding rate are constant. Given these conditions, the process can be considered stable for points far enough of the start/end of the tracks. Therefore, we follow the approach presented in [10], in which the domain V for the thermal analysis is a 2D cross-section of the substrate with thickness Dz. The considered cross-section is perpendicular to the deposition tracks.
Despite the previous assumptions (constant speed, power and feeding rate), the process is still complex to model since it is affected by the properties of the laser beam, the properties of the deposited material and the substrate and the thermal conditions in which the process is executed. Therefore, we make the following simplifications: -The addition of material and the phase change (melting) are not taken into account. The analysis is limited to the thermal history of the substrate for temperatures below the melting point. -Heat loss is not considered. The following phenomena are not considered: convection between the air and the substrate, thermal radiation, conduction between the substrate and its supporting floor. -The properties of the substrate (specific heat, density and thermal conductivity) are assumed to be constant. -Heat phenomena in Z direction are neglected, including the conduction produced when the laser heats a neighborhood close to V in Z (see Fig. 1).
The process studied in this paper can be applied to other laser-aided processes, such as powder-bed laser additive manufacturing.

Governing equations
Let V denote our 2D domain of analysis (see Fig. 1). Let T = T(x, t) denote the temperature at position x ∈ V at time t. The temperature function satisfies the partial differential equation where r and C p are the density and specific heat capacity of the material, q = q(x, t) is the heat flux and s = s(x, t) is the body heat source. The heat flux q satisfies the following constitutive relation given by Fourier law where k is the thermal conductivity of the material. In general, the thermal conductivity k is a second-order tensor. However, we assume that the material is isotropic and, therefore, treat k as a scalar.
To complete the mathematical formulation of the thermal problem, it is still necessary to define the initial and the boundary conditions. The initial temperature field is given by where T 0 (x) = 300 K is the prescribed temperature at time t = 0. Temperature (Dirichlet) and flux (Neumann) conditions are imposed on the boundary of V as where T ðx; tÞ and qðx; tÞ are known scalar functions, n(x) is the unitary outward normal to the boundary at x. In addition, where ∂V denotes the boundary of V. The initial and boundary conditions we imposed are detailed in Section 3.7.

Galerkin weak form
Given the constitutive relation in equation (2), the boundary conditions in equations (4) and (5), and applying the Galerkin method on equation (1), the problem is stated as follows [21]: where and, if Dz is the thickness of V, the differential elements of volume and area in equation (6) become dV = DzdA when integrating over V and dA = DzdL when integrating over ∂V q . The function T h aims to approximate the exact solution T and w h is a weighting function. Notice that T h satisfies the Dirichlet boundary conditions and w h vanishes where Dirichlet boundary conditions are applied.

Finite element discretization
The domain V is partitioned in finite elements V e such that where V denotes the closure of V. After using an isoparametric formulation, we obtain the semi-discrete form of equation (6): where T(t) is the vector of the nodal temperatures at time t, _ TðtÞ is the vector of the nodal derivatives of the temperature w.r.t. time: _ T a ¼ ∂T a =∂t, M and K denote the global mass and conductivity matrices, and f is the global force vector. The components of the corresponding element (local) mass and conductivity matrices are where M e ab and K e ab are the components of the mass and conductivity matrices that relate nodes a and b in the finite element e, and {N a } are the shape functions. We use 3-node linear triangular finite elements, therefore the functions {N a } , a = 1, 2, 3 are linear.
The components of the element force vector are given by where f e a is the component associated to the node a of the force in the element e and B T ¼ {b : x e b ∈ ∂V e T } represents the nodes in the element e with Dirichlet boundary conditions. The symbol x e b denotes the coordinates of the node b in element e.
The reader may notice that in equation (10), the vectors T and _ T are continuous functions with respect to time. To complete the numerical scheme it is still necessary to perform the time discretization.

Time discretization
To execute the time discretization, the time interval The goal in this section is to obtain the solution at time t s+1 given the solution at time t s .
Let s T = T(t s ) and s _ T ¼ _ Tðt s Þ for s = 0, 1, 2, … , N. We approximate the time derivative _ T at s + 1 using the backward Euler method, as follows [21]: Assuming that M, K and f are time-independent quantities, equation (14) is replaced into equation (10) to obtain From equation (15) can be obtained an expression for s+1 T: 3.6 Modeling of the heat provided by the laser The energy provided by the laser is modeled as a time dependent heat flux boundary condition. At every time step, we calculate the total influx through the boundary of each element on top (i.e. maximum Y direction) of the domain V. Assume the elements on top of the domain are e 1 , e 2 , … , e L . Let ∂V e i q be the edge of the element e i , i = 1, 2, … , L that lies on top of the domain. To calculate the influx energy on each of the ∂V e i q we measure the total power provided by the laser in a region R e i q which is the result of extruding the edge ∂V e i q half of the thickness in +Z and ÀZ directions (see Fig. 2), as shown in the following equation: where I(x, z) is the function that describes the laser intensity distribution and P e i is the total power acting on the edge ∂V e i q . Since the integration region is restricted to R e i q , only the power that acts on the 2D domain is considered.
Hence, to find the heat flux q e i at ∂V e i q due to the action of the laser, we divide P e i by the area of the integration region, as shown below: where j∂V e i q j denotes the length of the segment. Figure 2 shows a scheme of the process and the entities involved in it. In this figure, a Gaussian intensity function is depicted, but other intensity functions can be considered.
Since the laser is moving, the function I(x(t), z(t)) is also a function of time. Therefore, the previous procedure must be repeated at every time step of the simulation.   shows an example of the heat fluxes q e i calculated for the elements on top of the domain V, considering a Gaussian intensity function for the time sequence t 1 , … , t i , t i+1 , … , t N . Figure 5 presents the initial and boundary conditions imposed on our 2D domain for the thermal analysis. The initial temperature was 300 K over all the domain. We set constant ambient temperature of 300 K on the left and right hand sides of our domain. The bottom boundary was subjected to an adiabatic boundary condition, i.e. the heat flux was zero. Regarding the top boundary, as mentioned in Section 3.6, the laser energy input was model as a heat flux (Neumann) boundary condition. Therefore, those elements that interacted with the laser were subjected to a non-zero flux boundary condition. On the other hand, the elements that did not interact with the laser were under adiabatic boundary conditions (zero flux). Figure 6 shows the finite element mesh used for the simulations. The mesh was formed by 3671 linear triangular elements and 1922 nodes. The mesh was refined at the center of the top boundary, since it was the zone that interacted directly with the laser beam. The changes in the temperature obtained with a finer mesh or a smaller time step were negligible for the analysis performed in this work.

Material properties and process parameters for the numerical simulation
For the numerical simulations, we used the AISI 4140 steel. The material properties (thermal conductivity, density, specific heat and melting point) are given in Table 1 [10].
We executed seven numerical simulations in which we studied the influence of different laser intensity distributions, laser radii and process speeds. The simulations aimed to represent the deposition of four parallel tracks (as shown in Fig. 1). The domain configuration and process parameters used for the simulations are listed in Table 2. For a graphical representation of these parameters, we refer the reader to Figure 1.

Influence of the laser intensity function
We executed three simulations to study the influence of the laser intensity distribution function on the thermal behavior of the substrate. The intensity distributions used   were (1) Gaussian in equation (19), (2) uniform circular in equation (20) and (3) uniform rectangular in equation (21). Figure 4 shows the corresponding laser intensity distributions to a laser power P = 500 W and a laser radius R = 2.5 mm.  Figure 7 shows a comparison of the temperature distribution of the Gaussian intensity distribution (on the left) vs. the uniform circular intensity distribution (on the right). The figure presents the temperature at time t = 8Dt for the region of the substrate most affected by the laser. We can see that the temperature values and the shapes of the heat affected zone (HAZ) are similar. However, the maximum temperature in the case of the Gaussian distribution is 1210 K while in the case of the uniform circular distribution is 1141 K. This temperature difference is caused by the peak that we can observe in the Gaussian distribution (Fig. 4) which causes an energy concentration on the center of the HAZ. Figure 8 presents the comparison of the temperature distributions between the Gaussian (on the left) and the uniform rectangular (on the right) intensity distributions. It is noticeable that the HAZ of the uniform rectangular Table 1. Material properties of the AISI 4140 steel used in the numerical simulations [10].

Property Value
Thermal conductivity (k) 45 W/(m K) Density (r) 7800 kg/m 3 Specific heat (C p ) 500 J/(kg K) Melting point 1689 K Table 2. Domain size and process parameters used for the numerical simulations.

Parameter Value
Width (size in X) of the domain (W in Fig. 1) 100 mm Height (size in Z) of the domain (H in Fig. 1) 3 0 m m Thickness of the domain (Dz in Fig. 1) 5 m m Length of the tracks (L T in Fig. 1) 100 mm Separation distance between tracks (d T in Fig. 1   laser has a larger width and a smaller depth. Likewise, the maximum temperature in the case of the uniform rectangular intensity distribution is 1014 K, almost 200 K less than for the Gaussian intensity distribution. To obtain a measure of the energy that goes into the domain along the points of the X axis, we calculated the power per unit of length (I L ) corresponding to each laser intensity function, as given by equation (22).
The resulting functions are shown in Figure 9. It is noticeable the relation between the temperature distributions in Figures 7-9. The peak of the Gaussian distribution for I L is reflected in the large value for the maximum temperature. Likewise, the larger width in the HAZ corresponding to the rectangular laser beam is explained by the larger width of its corresponding I L . It is important to remark that the areas under the three curves in Figure 9 are equal to the total power of the laser (P = 500 W). Figure 10 shows the thermal history of the point at the top of the domain along the second track (see Fig. 1) for the three simulations executed with different laser intensity distributions. In this figure, we can observe that the maximum temperatures are reached when the laser describes the trajectory of the second track. We can see that the maximum temperature corresponds to the simulation performed with the Gaussian intensity distribution. However, before and after the temperature peak, the resulting temperatures for the three laser intensity distributions are very similar. This can be explained by the fact that during these phases, the thermal behavior is determined by the material properties.

Influence of the laser radius
To study the influence of the laser radius, we performed three simulations using a Gaussian intensity distribution with laser radii R = 2.0 mm, R=2.5 mm and R = 3.0 mm. Figure 11 compares the temperature distributions between the simulations with R = 2.5 mm (left) and R = 2.0 mm (right). On the other hand, Figure 12 compares the temperature distributions between the simulations with R = 2.5 mm (left) and R = 3.0 mm (right). The two figures present the temperature at time t = 8Dt for the region of the substrate most affected by the laser. In both figures, we can observe that the simulation with the smallest laser radius produces the largest temperature and depth of the HAZ. Since the delivered power of the laser is kept constant, a     12, 11 (2021) smaller radius means that the power is more concentrated at the center of the beam. This power concentration provokes the temperature differences.
We also calculated the amount of power per unit of length (I L ) along the points on the X axis (Eq. (22)). The resulting functions are depicted in Figure 13. It is noticeable how the peak of the function at the center of the laser beam (X = 0) increases while the laser radius decreases. Figure 14 shows the thermal history of the point at the top of the domain along the second track (see Fig. 1) for the simulations performed with different laser radii. The behavior is similar to the observed in our previous analysis of the impact of the laser intensity distribution: when the laser beam is close to the studied point, temperature differences are noticeable. However, when the laser beam is not close, the thermal behavior is determined by the material properties and the temperatures obtained for the three simulations become very similar.

Influence of the process speed
To study the influence of the process speed, we executed three simulations using a Gaussian intensity distribution with process speed v = 13.0 mm/s, v = 10.4 mm/s and R = 15.6 mm/s. Figure 15 presents the temperature distribution of the simulations with v = 13 mm/s (left) at time t = 8Dt and v = 10.4 mm/s (right) at time t = 10Dt. The figures show the configurations with maximum temperature on the first track. First, notice that, since the initial point of the laser was the same (see Tab. 2), the maximum temperature for the simulation with v = 13 mm/s was reached before than the one in the case of v = 10.4 mm/s (t = 8D vs. t = 10Dt). Figure 16 compares the temperature distributions between the simulations with v = 13 mm/s (left) at time t = 8Dt and v = 15.6 mm/s (right) at time t = 7Dt. As in the previous case, those images correspond to the configuration of maximum temperature for the first track.
In Figures 15 and 16, the maximum temperatures of the simulations with v = 10.4 mm/s, v = 13 mm/s and v = 15.6 mm/s were 1276 K, 1210 K and 1160 K, respectively. When the process speed is lower, the interaction time between the laser beam and the domain is larger. Therefore, the energy delivered by the laser (heat influx through the top boundary) is larger for lower process speeds. Consequently, larger temperatures and depths of the HAZs are obtained for lower process speeds. Figure 17 shows the thermal history of the point at the top of the domain along the second track (see Fig. 1) for the simulations performed with different process speeds.    The observed peaks of maximum temperature coincide with the previous analysis: lower process speed generates larger maximum temperatures. On the other hand, since the starting point of the laser is the same in the three simulations, we observe that the peaks of maximum temperature occur at different time. It occurs because the point of analysis is reached by the laser at different times, as a consequence of the process speed. The reader may notice that the larger the process speed, the faster the appearance of the temperature peak.

Comparison with non-linear simulations
In order to assess the suitability of the linear model in digital twin environments, we compared the simulations with a non-linear model that included radiation (Eq. (24)) and natural convection (Eq. (23)) heat losses: where h c is the natural convection coefficient, e is the thermal emissivity, s ≈ 5.67 Â 10 À8 W m À2 K À4 is the Stefan-Boltzmann constant and T ∞ is the ambient temperature. These two processes were included as Neumann boundary conditions on the top boundary (red boundary in Fig. 5). The details of the implementation of the non-linear model are out of the scope of this manuscript. We executed the same simulations that in Section 4.1 with the Gaussian, circular and rectangular laser intensity distributions, but including convection and radiation. The value of the parameters for the simulations are reported in Tables 1-3. Since the non-linear model includes heat losses, the temperature calculated are consistently lower than for the linear model. The temperature difference at every time step was below 5% in all cases, using the non-linear model as ground truth.

Conclusions and future work
This manuscript presents an analysis of the influence of the laser intensity distribution, the laser spot radius and the process speed on the thermal history of a substrate that is heated by the action of the laser. For the analysis, we implemented a 2D linear transient thermal model using the finite element method. The energy provided by the laser was represented as a time/space dependent heat influx (Neumann) boundary condition.
We executed simulations with three types of intensity distribution functions, namely Gaussian, uniform circular and uniform rectangular. Likewise, we performed simulations with laser radius R = 2.0 mm, R = 2.5 mm and R = 3.0 mm. In all the simulations, the total laser delivered power was identical and equal to P = 500 W. Results showed that these two parameters (the type of the intensity distribution function and the laser radius) strongly affect the shape (width and depth) of the HAZ and the maximum temperature.
The comparison of the linear and non-linear models proved that the linear simulations are good options for digital twin based design environments, due to their simplicity and reasonable temperature error. Apart of convection and radiation, further work is required to include important aspects to resemble the real process such as phase change and temperature dependent material properties.

Implications and influences
In the context of laser-based additive manufacturing, the thermal behavior of the substrate is relevant to define process parameters vis-à-vis piece quality. This manuscript studies the laser power spatial distributions (Gaussian, uniform circular and uniform rectangular) in addition to the laser power and process speed in 2D linear substrate heating simulations. The laser energy is modeled as a time dependent heat flux boundary condition on top of the domain. The results show that the laser intensity spatial distribution strongly affects the maximum temperature and the depth and width of the heat affected zone. These 2D finite element simulations prove to be good options for digital twin based design environments, due to their simplicity and reasonable temperature error, compared to