Issue 
Int. J. Simul. Multidisci. Des. Optim.
Volume 12, 2021
Computation Challenges for engineering problems



Article Number  11  
Number of page(s)  10  
DOI  https://doi.org/10.1051/smdo/2021011  
Published online  22 July 2021 
Research Article
2D linear finite element simulation of laser metal heating for digital twins
^{1}
Laboratory of CAD CAM CAE, Universidad EAFIT, Medellín, Colombia
^{2}
Vicomtech Foundation, Basque Research and Technology Alliance (BRTA), San Sebastian, Spain
^{3}
Department of Mechanical Engineering, Universidad EAFIT, Medellín, Colombia
^{*} email: amoreno@vicomtech.org
Received:
7
January
2021
Accepted:
29
June
2021
In the context of laserbased 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 nonlinear analysis (considered as ground truth for this case). Future publications address nonlinear finite element simulations of the laser heating process (including convection and radiation and temperature dependent substrate properties).
Key words: Numerical simulation / heat transfer / finite element method
© D. MontoyaZapata et al., Published by EDP Sciences, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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 highvalued industrial pieces [2]. However, the characterization of the laserbased 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 nonlinear 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.
2 Literature review
2.1 Numerical studies of process parameters in additive manufacturing
In recent years, there has been an increasing amount of literature on the study of laserbased AM. Several studies have shown that the temperature history has a significant impact in the quality and mechanical properties of the parts manufactured via laserbased AM [3–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–10] analyzed the effects of the process speed, laser power and material feeding rate on the geometry of the meltpool (with, depth and burnin 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 thermalfluid [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 meltpool 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]).
2.2 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 cyberphysical 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 cyberphysical world (e.g. thermal behavior, meltpool 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 realtime 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].
2.3 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 nonlinear 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 designstages.
3 Methodology
3.1 Problem description
Laser Metal Deposition (LMD) is a manufacturing process in which a highpower 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 Ω for the thermal analysis is a 2D crosssection of the substrate with thickness Δz. The considered crosssection 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 Ω in Z (see Fig. 1).
The process studied in this paper can be applied to other laseraided processes, such as powderbed laser additive manufacturing.
Fig. 1 Simulation of the deposition of four parallel deposition tracks (Track 1, ..., Track 4). Graphical representation of the domain, reference frame and parameters involved. 
3.2 Governing equations
Let Ω denote our 2D domain of analysis (see Fig. 1). Let T = T(x, t) denote the temperature at position x ∈ Ω at time t. The temperature function satisfies the partial differential equation(1)where ρ 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(2) where κ is the thermal conductivity of the material. In general, the thermal conductivity κ is a secondorder tensor. However, we assume that the material is isotropic and, therefore, treat κ 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(3) 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 Ω as(4) (5) where and are known scalar functions, n(x) is the unitary outward normal to the boundary at x. In addition, ∂Ω_{T} ∩ ∂ Ω_{q} = ∅ and ∂Ω_{T} ∪ ∂ Ω_{q} = ∂ Ω, where ∂Ω denotes the boundary of Ω. The initial and boundary conditions we imposed are detailed in Section 3.7.
3.3 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]:
Find T^{h} ∈ ^{h} ⊂ ^{1}(Ω) such that for all w^{h} ∈ ^{h} ⊂ ^{1}(Ω):(6)
where(7) (8)and, if Δz is the thickness of Ω, the differential elements of volume and area in equation (6) become dV = ΔzdA when integrating over Ω and dA = ΔzdL when integrating over ∂Ω_{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.
3.4 Finite element discretization
The domain Ω is partitioned in finite elements Ω^{e} such that(9)where denotes the closure of Ω. After using an isoparametric formulation, we obtain the semidiscrete form of equation (6):(10)
where T(t) is the vector of the nodal temperatures at time t, is the vector of the nodal derivatives of the temperature w.r.t. time: , 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(11) (12)where and 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 3node 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(13)where is the component associated to the node a of the force in the element e and represents the nodes in the element e with Dirichlet boundary conditions. The symbol denotes the coordinates of the node b in element e.
The reader may notice that in equation (10), the vectors T and are continuous functions with respect to time. To complete the numerical scheme it is still necessary to perform the time discretization.
3.5 Time discretization
To execute the time discretization, the time interval [0, t_{max}] is divided into N subintervals of length Δt: [t_{0}, t_{1}], [t_{1}, t_{2}], …, [t_{N−1}, t_{N}], such that t_{N} − t_{N−1} = Δt, t_{0} = 0 and t_{N} = t_{max}. 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 for s = 0, 1, 2, … , N. We approximate the time derivative at s + 1 using the backward Euler method, as follows [21]:(14)
Assuming that M, K and f are timeindependent quantities, equation (14) is replaced into equation (10) to obtain(15)
From equation (15) can be obtained an expression for ^{s+1}T:(16)
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 Ω.
Assume the elements on top of the domain are e_{1}, e_{2}, … , e_{L}. Let 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 we measure the total power provided by the laser in a region which is the result of extruding the edge half of the thickness in +Z and −Z directions (see Fig. 2), as shown in the following equation:(17)where I(x, z) is the function that describes the laser intensity distribution and P^{ei} is the total power acting on the edge . Since the integration region is restricted to , only the power that acts on the 2D domain is considered.
Hence, to find the heat flux q^{ei} at due to the action of the laser, we divide P^{ei} by the area of the integration region, as shown below:(18)where 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. Figure 3 shows an example of the heat fluxes q^{ei} calculated for the elements on top of the domain Ω, considering a Gaussian intensity function for the time sequence t_{1}, … , t_{i}, t_{i+1}, … , t_{N}.
Fig. 2 Calculation of the inner heat flux provided by the laser at every finite element. Only the laser power that lies inside Ω (dotted lines) is considered. 
Fig. 3 Curves of the inner heat flux through the top side of Ω for the time sequence t_{1}, t_{2}, … , t_{i}, t_{i+1}, … , t_{N}. 
3.7 Finite element mesh and boundary conditions
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 nonzero 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.
Fig. 4 Surfaces that describe the laser power distributions for (a) Gaussian, (b) uniform circular and (c) uniform rectangular beams. Laser power P = 500 W and laser radius R = 2.5 mm. 
Fig. 5 Initial and boundary conditions for the simulations. 
Fig. 6 Finite element mesh used for the simulations. Mesh refined at neighborhood of laser spot. 
3.8 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.
Domain size and process parameters used for the numerical simulations.
4 Results
4.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.(19) (20) (21)
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 = 8Δt 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 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).(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.
Fig. 7 Substrate temperature. (i) Gaussian (left) vs. (ii) uniform circular (right) laser power distributions. t = 8Δt. 
Fig. 8 Substrate temperature. (i) Gaussian (left) vs. (ii) uniform rectangular (right) laser power distributions. t = 8Δt. 
Fig. 9 Comparison of the functions of laser power per unit of length I_{L} for the (i) Gaussian, (ii) uniform circular and (iii) uniform rectangular laser power distributions. Laser power (area under the curve) P = 500 W and laser radius R = 2.5 mm. 
Fig. 10 Temperature history at x_{0} (Fig. 1). (i) Gaussian, (ii) uniform circular and (iii) uniform rectangular laser power distributions. 
4.2 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 = 8Δt 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 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.
Fig. 11 Substrate temperature. Gaussian laser power distribution. (a) R = 2.5 mm (left), (b) R = 2.0 mm (right). t = 8Δt. 
Fig. 12 Substrate temperature. Gaussian laser power distribution. (a) R = 2.5 mm (left), (b) R = 3.0 mm (right). t = 8Δt. 
Fig. 13 Gaussian laser power distribution. Power per unit length I_{L}. (i) R = 2.0 mm, (ii) R = 2.5 mm and (iii) R = 3.0 mm. Laser power (area under the curve) P = 500 W. 
Fig. 14 Temperature history at x_{0} (Fig. 1). Gaussian laser power distribution. (i) R = 2.0 mm, (ii) R = 2.5 mm, (iii) R = 3.0 mm. 
4.3 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 = 8Δt and v = 10.4 mm/s (right) at time t = 10Δt. 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 = 8Δ vs. t = 10Δt).
Figure 16 compares the temperature distributions between the simulations with v = 13 mm/s (left) at time t = 8Δt and v = 15.6 mm/s (right) at time t = 7Δt. 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.
Fig. 15 Substrate temperature. Gaussian laser power distribution. Laser speed: (i) v = 13.0 mm/s (left), (ii) v = 10.4 mm/s (right). 
Fig. 16 Substrate temperature. Gaussian laser power distribution. Laser speed: (i) v = 13.0 mm/s (left), (ii) v = 15.6 mm/s (right). 
Fig. 17 Temperature history at x_{0} (Fig. 1). Gaussian laser power distribution. (i) v = 10.4 mm/s, (ii) v = 13.0 mm/s, (iii) v = 15.6 mm/s. 
4.4 Comparison with nonlinear simulations
In order to assess the suitability of the linear model in digital twin environments, we compared the simulations with a nonlinear model that included radiation (Eq. (24)) and natural convection (Eq. (23)) heat losses:(23) (24)where h_{c} is the natural convection coefficient, ϵ is the thermal emissivity, σ ≈ 5.67 × 10^{−8} W m^{−2} K^{−4} is the StefanBoltzmann 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 nonlinear 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 nonlinear 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 nonlinear model as ground truth.
5 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 nonlinear 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 laserbased 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 nonlinear FEA (considered as ground truth for this case). Future efforts must address nonlinear finite element simulations of the laser heating process.
Funding
This work has received funding from the Eusko Jaurlaritza/Basque Government under the grants KK2018/00115 (ADDISEND) and KK2018/00071 (LANGILEOK).
Glossary
Ω ⊂ ℝ^{2}: Studied domain with boundary ∂Ω
x ∈ Ω: Coordinates to represent the position of Ω [m]
T(x, t): Temperature at x ∈ Ω in the instant t [K]
q(x, t): Heat flux into or out of the medium at x ∈ Ω at time t [W / m^{2}]
s(x, t): Volumetric heat sources at x ∈ Ω in the instant t [W / m^{3}]
n(x): Outward unit normal to the boundary at x ∈ Ω
ρ: Density of the material [kg / m^{3}]
C_{p}: Specific heat capacity of the material [J / (kg K)]
κ: Thermal conductivity of the material [W / (m K)]
References
 U.M. Dilberoglu, B. Gharehpapagh, U. Yaman, M. Dolen, The role of additive manufacturing in the era of Industry 4.0, Proc. Manufactur. 11, 545–554 (2017) [Google Scholar]
 M. Leino, J. Pekkarinen, R. Soukka, The role of laser additive manufacturing methods of metals in repair, refurbishment and remanufacturing enabling circular economy, Phys. Proc. 83, 752–760 (2016) [Google Scholar]
 B. Cheng, S. Shrestha, K. Chou, Stress and deformation evaluations of scanning strategy effect in selective laser melting, Addit. Manufactur. 12, 240–251 (2016) [Google Scholar]
 L. Costa, T. Reti, A.M. Deus, R. Vilar, Simulation of layer overlap tempering kinetics in steel parts deposited by laser cladding, in Proceedings of International Conference on Metal Powder Deposition for Rapid Manufacturing. MPIF, Princeton, NJ, 2002, p. 172–176 [Google Scholar]
 J.C. Heigel, P. Michaleris, T.A. Palmer, In situ monitoring and characterization of distortion during laser cladding of Inconel^{®} 625, J. Mater. Process. Technol. 220, 135–145 (2015) [Google Scholar]
 G.A. Ravi, X.J. Hao, N. Wain, X. Wu, M.M. Attallah, Direct laser fabrication of three dimensional components using SC420 stainless steel, Mater. Des. 47, 731–736 (2013) [Google Scholar]
 F. Cordovilla, P. Álvarez, A. GarcíaBeltrán, M.A. Montealegre, J.L. Ocana, Nonlinear thermal model of the direct laser melting process considering the adhesion of the consolidated material to the substrate using a domain with discontinuous material properties, in Proceedings of Lasers in Manufacturing (2019) [Google Scholar]
 D.M. Goodarzi, J. Pekkarinen, A. Salminen, Effect of process parameters in laser cladding on substrate melted areas and the substrate melted shape, J. Laser Appl. 27, S29201 (2015) [Google Scholar]
 H. Tian, X. Chen, Z. Yan, X. Zhi, Q. Yang, Z. Yuan, Finiteelement simulation of melt pool geometry and dilution ratio during laser cladding, Appl. Phys. A 125 (2019) [Google Scholar]
 W. Ya, B. Pathiraj, S. Liu, 2D modelling of clad geometry and resulting thermal cycles during laser cladding, J. Mater. Process. Technol. 230, 217–232 (2016) [Google Scholar]
 D.J. Corbin, A.R. Nassar, E.W. Reutzel, A.M. Beese, N.A. Kistler, Effect of directed energy deposition processing parameters on laser deposited Inconel^{®} 718: External morphology, J. Laser Appl. 29, 022001 (2017) [Google Scholar]
 J.I. Arrizubieta, A. Lamikiz, F. Klocke, S. Martínez, K. Arntz, E. Ukar, Evaluation of the relevance of melt pool dynamics in laser material deposition process modeling, Int. J. Heat Mass Transfer 115, 80–91 (2017) [Google Scholar]
 H. Liu, X. Qin, S. Huang, Z. Hu, M. Ni, Geometry modeling of single track cladding deposited by high power diode laser with rectangular beam spot, Opt. Lasers Eng. 100, 38–46 (2018) [Google Scholar]
 J. Posada, C. Toro, I. Barandiaran, D. Oyarzun, D. Stricker, R. de Amicis, E.B. Pinto, P. Eisert, J. Dollner, I. Vallarino, Visual computing as a key enabling technology for Industrie 4.0 and industrial internet, IEEE Comput. Graph. Appl. 35, 2640 (2015) [Google Scholar]
 M. Garetti, P. Rosa, S. Terzi, Life cycle simulation for the design of productservice systems, Comput. Ind. 63, 361–369 (2012) [Google Scholar]
 T. DebRoy, W. Zhang, J. Turner, S.S. Babu, Building digital twins of 3D printing machines, Scr. Mater. 135, 119–124 (2017) [Google Scholar]
 A. Gaikwad, R. Yavari, M. Montazeri, K. Cole, L. Bian, P. Rao, Toward the digital twin of additive manufacturing: integrating thermal simulations, sensing, and analytics to detect process faults, IISE Trans. 52, 1204–1217 (2020) [Google Scholar]
 Z. Yang, D. Eddy, S. Krishnamurty, I. Grosse, P. Denno, Y. Lu, P. Witherell, Investigating greybox modeling for predictive analytics in smart manufacturing, in Volume 2B: 43rd Design Automation Conference. American Society of Mechanical Engineers (2017) [Google Scholar]
 G.L. Knapp, T. Mukherjee, J.S. Zuback, H.L. Wei, T.A. Palmer, A. De, T. DebRoy, Building blocks for a digital twin of additive manufacturing, Acta Mater. 135, 390–399 (2017) [Google Scholar]
 L. Zhang, X. Chen, W. Zhou, T. Cheng, L. Chen, Z. Guo, B. Han, L. Lu, Digital twins for additive manufacturing: a stateof theart review, Appl. Sci. 10 (2020) [Google Scholar]
 A. Ibrahimbegovic, Nonlinear solid mechanics: theoretical formulations and finite element solution methods (Springer Science Business Media, 2009), Vol. 160 [Google Scholar]
 D. Mejia, A. Moreno, A. Arbelaiz, J. Posada, O. RuizSalguero, R. Chopitea, Accelerated thermal simulation for threedimensional interactive optimization of computer numeric control sheet metal laser cutting, J. Manufactur. Sci. Eng. 140 (2017) [Google Scholar]
Cite this article as: Diego MontoyaZapata, Juan M. Rodríguez, Aitor Moreno, Jorge Posada, Oscar RuizSalguero, 2D linear finite element simulation of laser metal heating for digital twins, Int. J. Simul. Multidisci. Des. Optim. 12, 11 (2021)
All Tables
All Figures
Fig. 1 Simulation of the deposition of four parallel deposition tracks (Track 1, ..., Track 4). Graphical representation of the domain, reference frame and parameters involved. 

In the text 
Fig. 2 Calculation of the inner heat flux provided by the laser at every finite element. Only the laser power that lies inside Ω (dotted lines) is considered. 

In the text 
Fig. 3 Curves of the inner heat flux through the top side of Ω for the time sequence t_{1}, t_{2}, … , t_{i}, t_{i+1}, … , t_{N}. 

In the text 
Fig. 4 Surfaces that describe the laser power distributions for (a) Gaussian, (b) uniform circular and (c) uniform rectangular beams. Laser power P = 500 W and laser radius R = 2.5 mm. 

In the text 
Fig. 5 Initial and boundary conditions for the simulations. 

In the text 
Fig. 6 Finite element mesh used for the simulations. Mesh refined at neighborhood of laser spot. 

In the text 
Fig. 7 Substrate temperature. (i) Gaussian (left) vs. (ii) uniform circular (right) laser power distributions. t = 8Δt. 

In the text 
Fig. 8 Substrate temperature. (i) Gaussian (left) vs. (ii) uniform rectangular (right) laser power distributions. t = 8Δt. 

In the text 
Fig. 9 Comparison of the functions of laser power per unit of length I_{L} for the (i) Gaussian, (ii) uniform circular and (iii) uniform rectangular laser power distributions. Laser power (area under the curve) P = 500 W and laser radius R = 2.5 mm. 

In the text 
Fig. 10 Temperature history at x_{0} (Fig. 1). (i) Gaussian, (ii) uniform circular and (iii) uniform rectangular laser power distributions. 

In the text 
Fig. 11 Substrate temperature. Gaussian laser power distribution. (a) R = 2.5 mm (left), (b) R = 2.0 mm (right). t = 8Δt. 

In the text 
Fig. 12 Substrate temperature. Gaussian laser power distribution. (a) R = 2.5 mm (left), (b) R = 3.0 mm (right). t = 8Δt. 

In the text 
Fig. 13 Gaussian laser power distribution. Power per unit length I_{L}. (i) R = 2.0 mm, (ii) R = 2.5 mm and (iii) R = 3.0 mm. Laser power (area under the curve) P = 500 W. 

In the text 
Fig. 14 Temperature history at x_{0} (Fig. 1). Gaussian laser power distribution. (i) R = 2.0 mm, (ii) R = 2.5 mm, (iii) R = 3.0 mm. 

In the text 
Fig. 15 Substrate temperature. Gaussian laser power distribution. Laser speed: (i) v = 13.0 mm/s (left), (ii) v = 10.4 mm/s (right). 

In the text 
Fig. 16 Substrate temperature. Gaussian laser power distribution. Laser speed: (i) v = 13.0 mm/s (left), (ii) v = 15.6 mm/s (right). 

In the text 
Fig. 17 Temperature history at x_{0} (Fig. 1). Gaussian laser power distribution. (i) v = 10.4 mm/s, (ii) v = 13.0 mm/s, (iii) v = 15.6 mm/s. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.