Impact of Boundary Parameters Accuracy on Modeling of Directed Energy Deposition Thermal Field

Abstract

Within the large Additive Manufacturing (AM) process family, Directed Energy Deposition (DED) can be used to create low-cost prototypes and coatings, or to repair cracks. In the case of M4 HSS (High Speed Steel), a reliable computed temperature field during DED process allows the optimization of the substrate preheating temperature value and other process parameters. Such optimization is required to avoid failure during the process, as well as high residual stresses. If 3D DED simulations provide accurate thermal fields, they also induce huge computation time, which motivates simplifications. This article uses a 2D Finite Element (FE) model that decreases the computation cost through dividing the CPU time by around 100 in our studied case, but it needs some calibrations. As described, the identification of a correct data set solely based on local temperature measurements can lead to various sets of parameters with variations of up to 100%. In this study, the melt pool depth was used as an additional experimental measurement to identify the input data set, and a sensitivity analysis was conducted to estimate the impact of each identified parameter on the cooling rate and the melt pool dimension.

1. Introduction

Nowadays, several techniques exist to build metal pieces with high added value; for instance, it is possible to use subtraction or addition techniques. For high-cost material or topologically optimized pieces, it is preferable to use additive manufacturing [1]. In this field, many technologies exist, such as laser powder bed fusion (LPBF) or directed energy deposition (DED—see norms ASTM F3187-16 [2]). For DED, several techniques are available through the use of a wire feeding method, i.e., Wire Arc Additive Manufacturing (WAAM), which is well adapted to big parts at the expense of a low geometrical accuracy and surface finish [3], or a powder feeding method which provides a better surface roughness [4].

DED has high potential in cladding uneven substrates. It reaches high cooling rates and generates fine microstructures [5]. If optimized, the process can achieve little distortion and high-quality deposits [6]. The laser cladding process is advantageously used as a repair technology [7]. Production of functional prototypes is another possibility, alongside small series production [8,9]. However, some problems still need to be fixed, such as microstructure monitoring, porosity which can lead to bad mechanical behaviour [10], geometrical accuracy inducing the need of post-process rectification [11], and crack events during manufacturing, either at substrate-deposit interface or in the clad [12].

The melt pool size has been identified as a critical parameter for maintaining optimal building conditions [13]. In microstructure genesis, the melt pool size and the cooling rate determine the properties of the future clad as well [14]. However, to accurately model the melt pool size and shape, and to solve the issues linked to the fluid state in the melt pool such as balling effect, keyhole, and porosity formation, Computational Fluid Dynamic (CFD) analysis is required. For L-PBF process, Heeling et al. [15] present an advanced CFD model taking into account accurate surface tension, the Marangoni effect, recoil pressure due to melt ejection by vaporization, and buoyancy force. Their validations for different process parameters using cubic samples of stainless steel 316 L and nickel-based superalloy IN738LC are quite convincing. For 316L, with a similar type of CFD model, Yao et al. [16] also present very accurate results about track shape, melt pool morphology, and thermal field, showing that a quantitative agreement between predictions and experiments can be reached. However, due to their completeness, these types of simulations are quite greedy from a CPU point of view. They can only be applied to small samples and academic problems. Nevertheless, they provide input about the local thermal field for metallurgic models, such as Phase Field [17,18] or Cellular Automaton [19], that predict microstructure genesis which includes solidification, grain nucleation and growth, and texture formation.

At a macroscopic scale, numerous Finite Element simulations have been developed to predict the thermal history of the deposit within DED process [20], focusing on the solid phases only. Such an approach relies on the strong assumptions that, for instance, neglect all the internal phenomena at the melt pool scale. Therefore, some corrections are applied to improve the prediction accuracy. For instance, the experimental conductivity value above solidus temperature is often multiplied by a coefficient γ defined according to literature data or just trials and errors (γ = 2.5 in [21], γ = 5 in [22]). These different FE models allow for one to define the correct melt pool size, and give information on the heat-affected zone (HAZ) where the material undergoes cycles of temperature variations [23].

The layer height is higher in the DED process (close or larger than 1 mm) than in the L-PBF one (usually a few tens of µm [24,25]) and smaller than that in the case of WAAM (above 1 mm [26]). Nevertheless, for any real part manufactured by a LPBF, DED, or WAAM process, the length of the laser path is of a magnitude larger than the equivalent diameter of the track (for this study, the typical transversal dimensions of a track were about 0.75 mm, while the total path was about 40,000 mm long). So, if the FE model simulates the total laser path, the computation time can reach days, and even months, and becomes unaffordable. In order to reduce computation time, many numerical techniques exist, such as using an artificially elongated heat source (part of a track or even a few layers) to increase the zone that is printed at once, with the major drawback being that temperature peaks are lost [27], or data driven models (mostly deep neural network enriched with physical laws) which need a large amount of data for their training [28,29], or an equivalent mathematical data driven technique such as the Proper Orthogonal Decomposition (POD). This latter technique can be applied on parts of the studied volume with low thermal gradients in order to reduce the total number of degrees of freedom (replacing FE computed thermal field by thermal modes). However, the limitation of such a method lies in the need to update these modes to avoid increasing error [30]. For these last two numerical models, the huge drawback resides in the necessity to have previous FE results in order to exploit them. Another direction followed in the present article is to decrease CPU time through the replacement of the 3D thermal field computation with a 2D assumption for thermal FE simulations. The 2D approach can be used for a 3D bulk sample, considering only the thermal field of its middle section. Extensive tests on such 2D simulations [31] confirmed the substrate temperature for parallelepiped bulk samples, as well as the melt pool size predictions, but also confirmed the different type of carbides precipitation associated to the thermal histories in different points of interest along the height of the sample were validated. Ref. [31] has also investigated the impact of thermophysical properties (conductivity, capacity) and boundary conditions (convection and emissivity coefficient) on the FE 2D predicted thermal history to conclude on the highest impact of the emissivity on the substrate temperature.

The main goal of the present study was to estimate the impact of the identified numerical parameters during the 2D FE simulations, such as the absorption coefficient of the laser power $\beta$, the convection, and radiation coefficients on the melt pool dimensions. The selected words “numerical parameters” underline the fact that in order to have compensated the thermal flow in the transversal direction, which is not taken into account in the 2D FE simulations, these parameters have lost their physical meaning and need careful identification to recover the experimental observations, including substrate temperature at thermocouples and melt pool size at the last layer. Note that further validations of the 2D FE computed temperature fields, such as the heterogeneity of the carbide type along the height of the bulk sample studied, as well as of the Vickers hardness profile and nano hardness information, were provided respectively within reference [5]. The present research highlights the fact that different sets of numerical parameters and associated results could be identified only if the thermocouple temperature measurements are exploited, or if poor accuracy is accepted within the identification methodology. Therefore, it is important to consider various experimental measurements to reach an accurate and robust fitting of the numerical 2D FE model.

The intensive research field around accurate thermal fields in additive manufacturing can be understood, as it is the main loading field in the process. As pointed out by Bayat et al. [32] regarding DED, if thermal simulations are validated, the next steps can focus on thermomechanical simulations. Let us be reminded that these calculations are not often totally coupled. Usually, the thermal field is computed for a chosen time increment, and then the thermal field is imposed to the mechanical FE model. However, a quicker way has been often applied, which is called the inherent strain method. With the aim of calculating the residual distortions and stresses at the end of the process, this method solves a limited number of mechanical equilibrium problems in which the loading field, known as the inherent strain, is used to compute the stress field, generally in an incremental way. Developed initially for the welding process, this method is extensively used in LPBF, and its application in DED can also be found (see, for instance, ref. [33], or recent PhD of J. Keumo-Tematio [34]).

To conclude this introduction, it should be emphasized that several important modelling results, such as the final microstructure or the field of residual stresses, need a well-defined thermal field, which leads to the potential of 3D simulations with high computing time. Therefore, 2D simulations that reduce the time-consuming issue are being investigated. This paper is dedicated to the potentiality of receiving erroneous results in the case of an incorrect calibration of 2D simulations.

2. Experimental Aspect

The 5-axis Irepa Laser Cladding system, with a Nd-YAG laser of maximum power capacity of 2000 W from Sirris Research Centre (Seraing, Belgium), is schematically described in Figure 1 [35]. The laser has a wavelength of 1064 µm and operates continuously. The metal powder is injected with an angle of 38 to 45 degrees, while the angle of the nozzle is 90°. The laser power has a top-hat energy distribution with a diameter of 1400 µm (it is considered a constant in the simulations [36]). A home-made protection cap filled with argon avoids contamination of the melt pool through interstitial elements such as oxygen or nitrogen [37]. The oxygen concentration within the protective atmosphere is below 10 ppm. The substrate of 50 mm height and 100 mm diameter is made of 42CrMo4 steel and preheated to 300 °C to prevent cracks in the deposit. The substrate top is surrounded by a 10 mm aluminium plate while, at the bottom, steel and aluminium plates are present, helping to diffuse heat.

The composition of the HSS M4 commercial powder with particle size ranging from 50 µm to 150 µm is provided in

Table 1. Chemical composition of HSS M4 (wt. %).

C	Mn	Cr	Mo	V	W	Ni	Si	Fe
1.35	0.34	4.30	4.64	4.10	5.60	0.9	0.33	Balance

, while

Table 2. Process parameters for the samples manufactured by laser cladding.

Laser beam speed [mm/s]	6.67
Laser power [W]	1100
Preheating [°C]	300

provides the different process parameters for the additive manufacturing.

Figure 2a shows the 40 × 40 × 27.5 mm³ (length × width × height) built sample made with 36 layers. The targeted height per layer was 0.76 mm. Figure 2c describes the position of the four thermocouples used: 5 mm deep in the substrate and at 20 mm from the edge of the substrate. For the calibration of this study, the thermocouple t3 was used as a reference (see Figure 3). Each layer consisted of 27 tracks, and the deposit strategy is shown in Figure 2b (note that only two layers and seven tracks per layer are shown for simplification purpose).

In Figure 3, the 36 large waves correspond to successive layers of the deposit. For each layer, in the zoom figure, one can identify the 14 local maxima associated with the moments when the laser is at the closest position to the thermocouple t3 for each track, and the 13 local minima that represent moments when the laser is at the opposite side of the clad.

To match experimental data with the 2D modelling assumption, a time shrinking factor had to be used. Indeed, the time to build a layer corresponds to the time to build a track in the 2D FE simulation, so one has to post-process the data (a Butterworth filter of order 5 was used in this study). As the simulation models the behaviour at the middle of the clad, it is necessary to use criteria to neglect oscillations due to tracks out of the middle plane. It is not possible for the model to reproduce these oscillations. According to the selected criterion, slightly different experimental targeted temperature curves are defined even if the whole trend is conserved. It is possible to refer to the extreme temperatures (minimum and maximum values) of the central track [38], or to refer to the extreme temperatures of the layer (by considering that these extremes happen when the laser is at the closest point to the thermocouple, or when the layer is finished (maximal distance)). To get a smooth curve, a spline was used to pass through the local minima and maxima of temperature values, see Figure 4.

During the process, the metal that is close in region to the laser melts, which represents the melt pool. Its dimensions will be largely influenced by the process parameters. Due to the remelting of the previous clad zone (Figure 5), it is impossible to measure the melt pool height, except for the last layer [39]. For this 36th layer, an average value of the melt pool depth measured on different samples manufactured with process parameters of Table 2 was 2.3 mm.

3. Numerical Methods

3.1. Implementation

In thermal simulations, the governing equation is the heat conservation equation

$$\rho C_p\frac{dT}{dt}=-\frac{\mathsf{\partial }\left(q_i\left(x_j,t\right)\right)}{\mathsf{\partial }{x}_i}+Q(x_j,t)$$

(1)

where T is the temperature, $\rho$ the volumic mass, $C_p$ the specific heat, $x_j$ the position vector, $q_i$ the heat flux vector, Q the body heat source, and t the time.

The heat conduction is described by the Fourier law

$$q_{cond}=-k(T,x)\frac{\mathsf{\partial }T}{\mathsf{\partial }{x}_i}$$

(2)

where k is the conduction coefficient, which depends on the temperature [40] with a neglected Marangoni effect within the melt pool [41,42]. However, it is possible to take it into account in a global way by multiplying the value of k by a factor [43] to artificially recover this effect.

As for boundary conditions, two different contributions are taken into account: heat loss due to convection, which is described by Newton’s law of cooling

$$q_{cond}=-k(T,x)\frac{\mathsf{\partial }T}{\mathsf{\partial }{x}_i}$$

(3)

and heat loss due to radiation. The latter is described by the Stephan-Boltzmann’s law

$$q_{rad}=ϵ\sigma (T_{surf}^4-T_{\infty }^4)$$

(4)

where $T_{surf}$ is the temperature of the surface, $T_{\infty }$ the temperature of the environment (supposed constant), $\sigma$ the Stephan-Boltzmann constant (5.67 × 10⁻⁸ W/m²K⁴). h is the convection coefficient, and $ϵ$ is the emissivity.

Some analytical models neglect these boundary contributions in AM processes, such as the well-known simple Rosenthal approach [44], even if their effect is important on the temperature of the clad. The difficulty of measuring the convection coefficient that evolves with the geometry explains why this coefficient is often assumed constant, or just with piecewise constant values [45]. The latter approach leads to a good trade-off, and will be applied for the present study as in [46]. A different value of the coefficient was used for the substrate and the clad. The same assumption was applied for the emissivity.

These different equations were implemented in an updated Lagrangian thermomechanical-metallurgical FE software named Lagamine (http://www.lagamine.uliege.be/dokuwiki/doku.php, accessed on 26 January 2024), developed in 1986 at the University of Liège. It was first developed to model continuous casting process, and this FE software has been applied in several different processes such as 3D printing, centrifugal casting of rolling mill, deep drawing, and incremental sheet forming, among others. It is also exploited to model material science phenomena in steels, Aluminium, Titanium, and Copper alloys like creep-fatigue damage at hot temperature, scale effect, and phase transformation. Its extensive library of constitutive laws is one of its advantages, as well as its flexibility. Within a High Performance Computing (HPC) environment, it is used to train Deep Learning models, or to identify macroscopic behaviour based on Representative Volume Element (RVE). The development of this software always relies on different experiments to offer serious validations, and to comparisons with other FE software packages.

The apparition of new added material is made using the birth element technique, which consists of neglecting a part of the domain during the FE computation, i.e., the elements that are still not activated. Once the laser passes over one element, it is activated and taken into account for the rest of the simulation, representing an additional printed volume of material. One exception concerns the boundary condition elements, which needs to be turned off once it is no longer on the boundary.

3.2. Identification Method

Based on the temperature measured by the thermocouple t3 (see Figure 3), an optimization was made with the iterative algorithm of Levenberg-Marquardt (LM) [47,48]. This method allows the identification of the parameters by minimizing the sum of the squares of the errors between the experimental and numerical curves through a sequence of updates to the set of parameter values. This objective function is often nonlinear and depends on several variables. This algorithm is based on two methods: the gradient descent method and the Gauss-Newton method. It has been implemented into a module specifically designed to allow the exploitation of multiple experiments [49], see Figure 6.

4. Numerical Results

4.1. Identification Results

The use of this LM optimization module permits one to identify three different sets of numerical parameters that could be considered as optimal (see

Table 3. Sets of parameters obtained with the Levenberg-Marquardt algorithm.

	$\mathit{\beta }$ [/]	${\mathit{h}}_{\mathit{c}\mathit{l}\mathit{a}\mathit{d}}$ $[\mathbf{W}/{\mathbf{m}}^{\mathbf{2}}\mathbf{K}$]	${\mathit{\epsilon }}_{\mathit{c}\mathit{l}\mathit{a}\mathit{d}}$ [/]	${\mathit{h}}_{\mathit{s}\mathit{u}\mathit{b}\mathit{s}\mathit{t}\mathit{r}\mathit{a}\mathit{t}\mathit{e}}$ $[\mathbf{W}\mathbf{/}{\mathbf{m}}^{\mathbf{2}}\mathbf{K}$]	${\mathit{\epsilon }}_{\mathbf{s}\mathbf{u}\mathbf{b}\mathbf{s}\mathbf{t}\mathbf{r}\mathbf{a}\mathbf{t}\mathbf{e}}$ [/]	Mean Error [K]
Set 1	0.03955	479.93	1.83	72.23	0.95	4.66
Set 2	0.04097	823.47	0.8	82.49	0.8	4.27
Set 3	0.04172	707.6	0.9	120	0.9	4.39

, where grey cells represent parameters that were chosen as a fixed value for the current identification), demonstrating the non-unicity of the solution for this 2D FE model. The optimal data sets generate quite a close temperature evolution at the thermocouple position (see Figure 7). The optimized parameters are $\beta$, h, and $ϵ$, which are the absorption, the convection, and the emissivity coefficients, respectively. For Set 1, the five parameters were mathematically tuned without any constraint. For Set 2, emissivity coefficients were fixed to physical values [50], and for Set 3, both emissivity and substrate convection coefficients were imposed. The mean error column corresponds to the average of the absolute difference between experimental data and simulation temperature.

As expected, the 2D FE model cannot accurately predict the oscillations, but it can simulate the thermal trends in a quantitative way. In Figure 7, the predicted oscillations are larger for the first layers than for the next layers, and disappear at the end of the simulation. Experimentally, the corresponding large oscillations are related to successive layers (see Figure 3). The larger oscillations at the start of the process are due to the small distance between the laser and the thermocouple, which reduces the effect of the thermal diffusion. This effect also explains why the model does not reproduce oscillations after a certain number of layers. The post treatment of the experimental curve itself, discussed at the end of Section 2, justifies why one should not focus on the exact match of the temperature at the substrate level. A more important validation is the melt pool depth that confirms the accuracy of the local thermal field.

4.2. Validation of the Predicted Results

The melt pool depth was computed over the entire simulation with the three different sets, as shown in Figure 8. Its value has been estimated as the middle of the track to minimize edge effects.

A first observation can be made on the predicted temperature at the thermocouple compared to the computed melt pool of the different sets, even if the predicted temperature is the lowest during the first 80 seconds for Set 2 (see Figure 7), the smallest computed melt pool belongs to the Set 1 (see Figure 8).This can be related to the fact that the different parameters do not have the same effect. An increase in the values of the substrate convection and radiation parameters will decrease the temperature of the substrate more efficiently (close to the thermocouple), but will change melt pool dimensions less. This also explains why those sets provide close melt pool depths until the fifth layer. Indeed, the clad temperature and the clad height at this moment are not high enough to generate any difference between the three cases. Conversely, the increased $\beta$ value of the Set 2 compared to Set 1 brings more heat, close to the melt pool, which tends to increase its size. For the Set 3, the same approach can be applied, where an increase of the heat loss flux is compensated by an increase of $\beta$, leading to a similar temperature at the position of the thermocouple. However, closer to the melt pool, the quantity of heat received is increased, which conducts to a deeper melt pool.

4.3. Sensitivity of the Simulation Results to the Boundary Conditions

To evaluate the impact of these boundary parameters, some perturbed sets have been defined using the set 2 as a validated reference. The value of the perturbation is about 0.1% of the identified parameter, which gives the new sets of

Table 4. Sets selected used to evaluate the sensitivity of simulations results.

Perturbed Parameter Sets	$\mathit{\beta }$ [/]	${\mathit{h}}_{\mathit{c}\mathit{l}\mathit{a}\mathit{d}}$ $[\mathbf{W}\mathbf{/}{\mathbf{m}}^{\mathbf{2}}\mathbf{K}$]	${\mathit{\epsilon }}_{\mathit{c}\mathit{l}\mathit{a}\mathit{d}}$ [/]	${\mathit{h}}_{\mathit{s}\mathit{u}\mathit{b}\mathit{s}\mathit{t}\mathit{r}\mathit{a}\mathit{t}\mathit{e}}$ $[\mathbf{W}\mathbf{/}{\mathbf{m}}^{\mathbf{2}}\mathbf{K}$]	${\mathit{\epsilon }}_{\mathbf{s}\mathbf{u}\mathbf{b}\mathbf{s}\mathbf{t}\mathbf{r}\mathbf{a}\mathbf{t}\mathbf{e}}$ [/]
Set 2	0.04097	823.5	0.8	82.49	0.8
Perturbed $\beta$	0.04101	823.5	0.8	82.49	0.8
Perturbed $h_{clad}$	0.04097	824.3	0.8	82.49	0.8
Perturbed ${\epsilon }_{clad}$	0.04097	823.5	0.801	82.49	0.8
Perturbed $h_{substr}$	0.04097	823.5	0.8	82.57	0.8
Perturbed ${\epsilon }_{\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{s}\mathrm{t}\mathrm{r}}$	0.04097	823.5	0.8	82.49	0.801

, the perturbed parameter is represented with a grey cell.

From these perturbed sets of parameters, the sensitivity of the depth of the melt pool and the cooling rates were estimated at three different point of interest (POI) positions (Figure 9). Their positions are defined by taking the upper surface as a reference. The height of these POI is sufficiently different enough to represent the material heterogeneity of the DED processed part proved by their hardness profile (see [31]). A measure of the result sensitivity to the perturbed parameter is defined in Equation (5):

$$S_f=\left(\frac{\frac{f\left({BC(x}^{*})\right)-f\left(BC\left(x\right)\right)}{f\left(BC\left(x\right)\right)}}{\frac{x^{*}-x}{x}}\right)$$

(5)

$S_f$ is the sensitivity to the parameter x, which belongs to a set of values that represents the boundary conditions BC, x the identified value, $x^{*}$ is the perturbed condition, and f defines the FE prediction values. This definition of the sensitivity is quite similar to the one used in the LM algorithm to define the Jacobian matrix.

As shown in Figure 10, increasing the power of the laser leads to an increase in the melt pool depth. Substrate parameters have nearly no impact on the melt pool dimension, as they do not influence the local thermal field in a region close to the laser. Conversely, an increase in the clad parameters (convection and radiation) tends to decrease the melt pool depth.

Concerning the cooling rates, two different values can be estimated based on the temperature history (see Figure 11). The average cooling rate between peak max and min temperature (corresponding to the red crosses in Figure 11) and the average cooling rate in the liquidus-solidus zone (zone between red dashed lines), if liquidus or solidus is not reached, the maxima are taken to get the average slope (green crosses in Figure 11). For POI3, there is no third peak (process is over), so the minimum is taken after a certain time that is equal to the time between the maximum and the minimum of the first peak. This assumption is made for consistency.

In Figure 12 and Figure 13, the sensitivity is evaluated on the four first peaks of temperature, except for POI3 which undergoes only 2 peaks. It is important to consider melt pool dimensions (see Figure 8), which shows that the three first peaks happen when the POI is in the melt pool (each layer of the clad is modelled with a height of 0.75 mm which corresponds to the experimental mean height), which explains why the sensitivity of the fourth peak does not always follow the same trend as others. For instance, in Figure 12 for POI1 with the set Perturbed $\beta$, for the three first peaks, increasing $\beta$ will lead to an increase in the cooling rate (higher local temperature which leads to higher gradients and consequently to higher cooling rates). While, for the fourth peak, the cooling rate will decrease due to an increase of the amount of heat provided to the clad.

5. Conclusions

Within the large amount of models applied to DED, this article is focused on how a 2D FE simulation can be reliable, as saving CPU time must not prevent accuracy. Using a bulk sample (with a size of 40 × 40 × 27 mm³) in M4 high speed steel, it is demonstrated that a 2D FE model requires the use of several virtual parameters for the boundary conditions, indeed, avoiding the recognition that the thermal flow in the out of plane directions decreases the physical meaning of the parameters. The correct prediction of the temperature measured at a substrate thermocouple is considered to be global information, which is not sufficient enough to reach an accurate local temperature field. The article demonstrates the need of an improved fitting procedure including the melt pool size, which represents a local measure.

A sensitivity study was made to determine the impact of the accuracy of the fitted boundary conditions on the melt pool depth and the cooling rate, which are two quantities determining the clad microstructure.

The main conclusion of this study is that the laser power and the clad boundary parameters (convection and radiation) are critical in defining the melt pool dimensions. However, the sensitivity values highly depend on the reference peak of temperature and the positions that are considered.

This study was based on the 2D assumption, which is a rough approximation as it considers that transversal heat fluxes are negligible. This hypothesis can be applied on a bulk clad with low surface-volume ratio but, to our knowledge, it has not yet been investigated on thin parts.

Some future work will be done about the influence of the thickness of the deposit on the temperature field for 3D cases compared to their 2D model. These interactions will be analyzed, thanks to artificial intelligence, as it has already been done previously in [51] to reproduce the different results with a surrogate model that can be rapidly trained thanks to the 2D assumption.