Assuming a Lagrangian frame Ω and a material point located by r ( r ∈ Ω) as the reference, given thermal energy balance at time t, the governing equation can be formulated as follows: ρ ( T ) C P ( T ) ∂ T ( r , t ) ∂ t =   − ∇ ⋅ q ( r , t ) + Q ( r , t ) ,   r ∈ Ω (1) where ρ is the material density, C _p is the specific heat capacity, T is the temperature, t is the time, Q is the heat source, and q is the heat conduction flux vector, calculated as: q = − k ( T ) ∇ T ( r , t ) (2) Where k is the thermal conductivity of the material.

Table 1 presents the temperature dependent thermal properties for SS 316L (Mills and Mills, 2002). Linear interpolation is used to calculate the properties at any temperature.

The Single Ellipsoid (SE) heat input model (Goldak) can be used to describe an equivalent to the laser heat source (Goldak et al., 1984) as: Q SE =   6 3 AP abc π π exp ( − 3 ( x + v s t ) 2 a 2   − 3 y 2 b 2   − 3 z 2 c 2 ) (3)

The laser power is P and the laser absorption efficiency is A. The value for laser power P is based on measurement, as will be discussed in Experiment Set-Up Section. The value of A is calibrated using the method of reverse calibration by iteratively fitting the simulated temperature field at thermocouple location to match the experiment results described in Ref. (Denlinger et al., 2015). x, y and z are the local coordinates with the origin centred at the ellipsoid where the heat source reaches the maximum intensity and with a moving velocity v s . Parameters a, b and c represent the respective length of the longitudinal, transverse and through the depth semi-axes of the ellipsoid parallel to the local x, y and z axes. Generally, a is taken as the melt-pool length, b is taken as the half width of the deposition bead and c to the melt pool depth (Denlinger et al., 2015).

It has been shown previously that, using single ellipsoid heat source predicts a lower temperature gradient at the front and a higher temperature gradient at the trailing edge as compared to the experiments (Goldak et al., 1984). Therefore, it has been a common practise among the researchers to employ double ellipsoid heat source to be more computationally accurate. The double ellipsoid (DE) Heat Input model (Goldak) is used in the present work to describe the laser heat source (Goldak et al., 1984): Q DE =   6 3 A P f f , r a f , r bc π π exp ( − 3 ( x + v s t ) 2 a f , r 2   − 3 y 2 b 2   − 3 z 2 c 2 ) (4) f is a weighting fraction, that determines the energy partition among the front and rear ellipsoid. Typically, different values are employed in the front and rear of the heat source for the longitudinal axis length a (Goldak et al., 1984). The orientation of the double ellipsoid (DE) heat source and coordinate system are depicted in Figure 1A. For the single ellipsoid (SE) heat source the lengths of two ellipsoids are equal (a_f = a_r), hence depicting a single ellipsoid.

To accurately model the motion or movement of laser (input heat source), the simulation time increments should be small enough in such a way that the numerical heat source should move less than or equal to half of laser spot size i.e., spot size radius in one-time step (Irwin and Michaleris, 2016). Therefore, in the present work for all cases using Double Ellipsoid Heat Source, computation time-step does not exceed R / v s , where R is laser spot size radius and v s is laser source/travel speed. In LDED, track length is big and with the specified computation time increment, it can lead to high computation time. Therefore, to reduce the computation time, an elongated ellipsoid line heat input model can be used that averages the single ellipsoid heat source over its path (Irwin and Michaleris, 2016). Q AVERAGE =   1 Δ t ∫ t 0 t 0 +   Δ t Q SE d t (5) where t 0 is the time at beginning of the increment and Δ t is the duration of time increment over which SE model is averaged. But this averaging over its path will result in large thermal gradients at the edges of each segment especially when the linear segment is large ( v s Δ t ≫ a ) (Irwin and Michaleris, 2016). Therefore, to smooth out the discontinuities at the segment edges, elongated ellipsoid model is developed where peak value of the heat input is at the middle of the sub-track or segment. Power density of elongated ellipsoid at the beginning and end of each segment is half of its peak value that results in smooth distribution over successive linear segments. Thus, Elongated Ellipsoid (EE) line heat input model is also used in present work to describe the laser heat source (Irwin and Michaleris, 2016): Q EE =   6 3 AP a ^ bc π π exp ( − 3 ( x +   v s ( t   +   1 2 Δ t   ) ) 2 a ^ 2   − 3 y 2 b 2   − 3 z 2 c 2 ) (6) where the length of each elongated ellipsoid sub-track or segment a ^ is (Irwin and Michaleris, 2016): a ^ =   v s Δ t 2 3 log ⁡ 2 (7)

The main difference between Eq. 3 and Eq. 6 are the introduction of elongated length a ^ in place of a, that stretches the heat distribution in the local x direction and additional (1/2) Δ t which shifts the peak heat input from the end to the middle of the linear segment or sub-track (Irwin and Michaleris, 2016). Also, the double ellipsoid heat source moves continuously along the path replicating laser movement in experiments with computation time step less than R / v s , the elongated ellipsoid heat source also moves continuously along the path, but enables bigger time steps equal to Δ t that results in reduction of computation time. A comparison of the energy distribution with three different heat source models SE, DE and EE is shown in Figure 1B . A detailed analysis of the Elongated Ellipsoid Heat Source and its functioning is the presented in the reference work (Irwin and Michaleris, 2016).

Elongated Ellipsoid heat source model allows the simulation of an entire heat source scan in one time increment instead of hundreds of time increment when using Double Ellipsoid or Single Ellipsoid Heat Source. However, taking large time increments also leads to increase in errors as well (Irwin and Michaleris, 2016). Therefore, the track scan is divided into several linear segments (sub-tracks), for each of which, heat source is applied in one simulation time increment.

The thermal governing equations are supplemented by the initial condition and different types of boundary conditions: T ( r , 0 ) =   T 0 ( r , 0 )   r ∈ Ω (8)

Thermal radiation q r a d is accounted for using the Stefan–Boltzmann law: q rad =   ε   σ   ( T 4 − T amb 4 ) (9) where ε is the surface emissivity, σ is the Stefan–Boltzmann constant, T the surface temperature of the workpiece, and T a m b is the ambient temperature.

Newton’s law of cooling describes the heat loss due to convection q c o n v : q conv = h   ( T − T amb ) (10) where h is the convective heat transfer coefficient.

The effect of latent heat of fusion during the melting and solidification process is accounted by modifying the specific heat capacity C p * as in Eq. 1, where L f (J/kg) is latent heat of fusion, and T m is melting temperature (Yan et al., 2017). C p * ( T ) =   L f T m −   T amb +   C p ( T ) (11)

In order to take the convective redistribution of heat in the melt pool due to the fluid flows into account, a higher value of thermal conductivity is considered, hence avoiding integrating complex analytic form in the model (Lampa et al., 1997). Thus, enhanced thermal conductivity factor k ^* is assumed in Eq. 2 (Ren et al., 2019): k ∗ ( T ) = { k ( T ) , T < T m 2.5 × k ( T ) , T ≥ T m (12)

In situ temperature is measured at two different locations on the bottom face of substrate, as shown in Figure 3A, using Ω GG-K-30 type K thermocouples of 250 µm diameter. The thermocouples have a measurement uncertainty of ± 0.75%. TC 1 and TC 2 are located on the bottom surface of the substrate along the deposition path. The thermocouple signals are read by National Instruments modules 9213. The module record data in Signal Express at a sampling rate of 200 Hz.

An Infra-Red imaging camera (NIT system) is used analyse the melt-pool stability and to estimate the melt-pool dimensions based on image analysis. Image acquisition was done at a sampling rate of 200 fps. Image analysis is performed with the Igor Pro software on these images and an average value of melt-pool width (W) and length (MP_L) is measured for each experiment case based on a specific contour level (L), as shown in Figure 3B. The laser penetration depth was also measured for each case by sectioning as recommended in (Goldak et al., 1984) and illustrated in the macrograph shown in Figure 3C. The melt-pool dimensions values given in the Table 3 are then used to calibrate the values of heat source as discussed in previous section.

The FEM analysis is performed using COMSOL Multiphysics based solver (PARDISO) with the implicit Backward Differentiation Formula (BDF) time stepping method. Adaptive time stepping method is employed rather than strict formulation with maximum time step of R / v s for DE and Δ t for EE heat sources. The solver is further adapted specifically to model additive manufacturing technologies in the present work. All simulation cases are performed on an Intel Xeon W-2275, 16 cores, with 128 GB RAM workstation.

Figure 4 displays the three-dimensional finite element mesh, generated in COMSOL Multiphysics, used for the thermal model. The same mesh is used for all cases and for different Heat Source models used in the present work. The mesh contains 33,852 Hex-8 elements and 49,790 nodes. Hex-8 elements were chosen because it has been proved that they give more accurate results as compared to tetrahedral elements for the plastic deformation (Benzley et al., 1995), therefore to have same mesh for thermal and mechanical analysis in future work, Hex-8 elements were used. The elements for the deposited material are allotted as 3 per laser spot size and 2 per deposition thickness, making the elements 0.75 × 0.75 × 0.19 mm in volume for experiment case D1, but varies for different experiment cases as track geometry changes. The mesh is coarsened at the substrate as it moves away from the wall builds. A mesh convergence study was done using three different mesh strategies to confirm the accuracy of thermal analysis.

The “quiet” element activation method is used to simulate the deposition of material during the DED process. The elements that represent the wall builds are pre-existing at the beginning of the analysis as shown in Figure 4. But their material properties are rescaled in such a way that it does not affect the computational analysis. In the present work, the thermal conductivity k and specific heat C _p are set to very low values to minimize the heat transfer from active to quiet elements, as follows: k quiet =   α k k (13) C P quiet =   α C P C P (14) where k quiet and C P quiet are the respective thermal conductivity and specific heat for the quiet elements; α k and α C P are the appropriate scaling factors that are usually assigned very small values, for instance, 10^–4 is employed in the present work. The properties of an element are switched from “quiet” to “active” when any Gauss point of the element is consumed by the heat source volume (Eq. 15). exp ( − 3 ( x + v s t ) 2 a f , r   2   − 3 y 2 b 2   − 3 z 2 c 2 )   ≥ 5 % (15)

This means that if the evaluated heat source value at any Gauss point is greater than 5% of the peak intensity, the element will be activated by switching material from quiet to active.

For the elongated ellipsoid line heat input model, material activation is done in the same way of heat source intensity as explained above, but in this case only the activation domain is larger over the deposition path. If the evaluated heat source value at any Gauss point is greater than 50% of the peak intensity, the element turned to active state. The activation criteria of 50% in the deposition direction is explained in the further sections. exp ( − 3 ( x +   v s ( t 0   +   1 2 Δ t   ) ) 2 a ^ 2   )   ≥ 50 % (16)

To simplify the model, deposited wall build is considered to be a flat rectangular shape with constant layer height and width. This is contrary to the experiments as layer height changes in the few first layers and then reaches a uniform layer thickness. As discussed previously, to develop an accurate thermal model, certain input parameters have to be calibrated against experiment results. The laser absorption efficiency ( A ) is taken as 0.4 with reverse calibration of fitting simulated temperature field to experiment results at thermocouple location iteratively as suggested in (Michaleris, 2014). Also, for all experiment cases, as laser spot size is 2.2 mm, laser spot size radius is taken as 1.1 mm. The double ellipsoid heat source dimensions parameters are dependent on experiment cases but follow the same rule i.e., front ellipsoid length a f =   W / 2   a n d   a r = 2 a f ,   b =   W / 2   a n d   c = H + M P D , that represent the double ellipsoid dimensions to be almost equal to experiment melt-pool dimensions for each experiment case. Emissivity ( ε ) is temperature dependent, but mean value of 0.6 is taken as widely reported in the literature (Biegler et al., 2018b). The convective heat transfer ( h ) coefficient is temperature dependent, but a mean value is taken at the substrate h = 5 (W.m ⁻² .K ⁻¹ ) and to consider the effect of forced convection due to the powder carrying argon gas at the wall builds, average convective heat transfer coefficient of h = 25   (W.m ⁻² .K ⁻¹ ) is taken as widely reported in the literature that it can be in the range of h = 20 − 60   (W.m ⁻² .K ⁻¹ ) (Heigel et al., 2015). At the clamped surfaces, the substrate is in contact with metallic fixtures clamps shown in Figure 1 (d). In the numerical model, to reduce the computation time, fixtures clamp is not included, therefore for these surfaces, a higher heat loss is modelled through conductive to convective equivalent heat loss expression similar to Newton’s law of cooling Eq. 10. The thermal conductance at the contact surfaces is thus h = 150 (W.m ⁻² .K ⁻¹ ) to model the conductive heat transfer from substrate to the metallic fixture clamps as suggested in (Lu et al., 2019).

The Quiet/Active material activation is well implemented in the model as can be seen in Figure 5. Indeed, during the deposition of the 12th layer, quiet elements of the layer and elements above are not activated and hence does not contribute to the heat transfer, as the material properties of quiet elements are assigned to a dummy material. The current deposition elements and previously deposited layer elements are correctly activated from quiet to active once the activation criteria are satisfied that depends upon the laser travel (Eq. 15). The input heat source absorptivity and heat losses parameters are kept same for all experiment cases as discussed in previous section Model Calibration. Double Ellipsoid heat source dimensions are dependent on experiment cases and is taken following the rule explained in the previous section Model Calibration. The thermal response of the workpiece is calculated by the Double Ellipsoid model and compared to the experimental measurements.

Figure 6 shows the comparison between experimental results, as measured by thermocouples 1 and 2, and the numerical results at the corresponding nodes in FEM analysis for experiment set 1 (effect of dwell time). As explained in the previous section, due to the experiment set-up, because thermocouple 1 and 2 are at different locations on the substrate but along the deposition line, they record almost same thermal histories with a time offset. Therefore, only Thermocouple 1 results are presented. The thermal response can be classified in 4 different zones depending upon the thermal history. In zone 1, peak temperature for successive build layers increases as layers are building up, in zone 2, peak temperature for successive layers stabilises and there is no further increase of temperature, and in zone 3, peak temperature for successive build layers starts to decrease due to the fact that heat source is moving vertically away from the thermocouple locations and fabricated material starts to increase that also conducts heat leading to the less heat transfer or temperature gradient. In zone 4, once the process is finished, temperature rapidly reduces with respect to time.

In the experiment set 1, analysis of effect of dwell time is studied, as shown in Figure 6. Peak temperature starts to decrease as dwell time increases from 0 to 30 s. Longer dwell times results in lower peak temperatures of 230°C with the 30 s dwell time (D4), while with 0 s dwell time (D1) exceeding 500°C. Therefore, dwell time has a strong influence on temperature evolution as well as the peak temperatures obtained during the deposition process.

For experiment set 1, analysis of effect of dwell time, it is observed that thermal response is changing drastically, not just in terms of peak temperatures, but also in terms of thermal evolution. Table 4 that counts the number of peaks (or layers) observed in the different zones (Z1 to Z4), the number of layers under Zone 1 are decreasing with an increase of dwell time, depicting peak temperature is stabilising at much earlier deposition stage (number of layers). But, the number of layers under Zone 2 and 3 are increasing with an increase of dwell time, depicting that stabilised peak temperature are maintained for longer duration of deposition. Numerical model with calibrated parameters captures the temperature evolution trend and peak temperatures correctly with change of dwell time.

In experiment set 2, analysis of effect of laser power is studied, as shown in Figure 7 . Peak temperature keeps on increasing as laser power increases from 800 to 1400 W. higher laser power results in higher peak temperatures of 480°C with 1400 W laser power (P4), and only 280°C with 800 W laser power (P1). Therefore, laser power has a strong influence on peak temperature, but does not influence the trend of temperature evolution obtained during the deposition process.

For experiment set 2 (analysis of effect of laser power), it is observed that thermal response is changing peak temperatures drastically, but does not influences the thermal gradient. As it can be seen in Table 5, number of layers under Zone 1, 2 and 3 remains same depicting that temperature trend is lifted upwards (increase of peak temperature). Higher peak temperatures are recorded because of the increase of laser power, but the thermal gradient remains the same. As can be seen in Table 5, peak temperature is stabilised at fifth layer, stabilises for next 4 layers and then starts to decrease as heat source is moving vertically away from thermocouple locations for all experiment cases (P1-P4).

The results of the transient thermal analyses are in close agreement with the experimental results as can be seen in Figure 6 and Figure 7. Errors between experiment and simulation results are calculated by comparing instances in time. %   Error =   100 ∑ i = 1 n | ( T exp ) i −   ( T sim ) i ( T exp ) i | n (17) where n is the total number of simulation time increments between the beginning and end of the deposition, i is the current time increment, T _sim is the simulated temperature, and T _exp is the measured temperature. The largest error at thermo-couple is found to be 3.91%.

Table 6 shows the computation time and percent error at both thermocouples TC1 and TC2 for all cases. For experiment set 1 (D1 to D4), it can be clearly seen that with the increase of dwell time, that leads to increase in number of time increments leading to increase in computation time. For experiment set 2 (P1 to P4), it can be observed that with an increase in laser power, that leads to an increase in simulation peak temperature at the melt-pool region i.e., higher thermal gradient, that leads to slight increase of computation time.

To reduce the computation time, elongated ellipsoid model is used to reduce the number of simulation time steps by dividing the complete track in number of linear sub-tracks, with each sub-track is solved in one simulation time step. Different track size (sub-track) is chosen for material activation presented in Table 7 and a comparison is done with the experiment results to make a compromise between computation speed and accuracy. K E =   v s Δ t a (18)

As shown in Figure 8A, 5 different simulation time increments are chosen, where K_E = 0.5 represents heat source movement of R with a time step of R / v s , hence the ellipsoid is not elongated and that is same as in the case of double ellipsoid heat source, K_E = 4 represents heat source averaging over a segment of 8R with a simulation time step of 8 R / v s and K_E = 27 represents heat source averaging over a segment of 54R with a simulation time step of 54 R / v s , that is averaging over the complete track, activating the complete track and performing heat transfer analysis in 1 simulation time increment.

The movement of the elongated ellipsoid heat source is such that there is a smooth distribution of power intensity over the successive scan segments. The power intensity at the start and end of each scan segment is half of its peak value, resulting in a smooth distribution over the successive segments of a track as shown in Figure 8B. And this is the reason why activation criteria for elongated ellipsoid (EE) heat source model is 50% that accounts for the half of its peak value in the deposition direction. Of course, increasing the factor K_E helps in reducing the total computation time, but it also leads to the increase of computation error. Thus, effect of sub-track size activation on temperature evolution at thermocouple location is shown in Figure 9. K_E of 0.5 (∆t = R / v s ) gives the same accurate results as we observed with Double Ellipsoid (DE) heat source with a simulation time step of R / v s . As explained in the section further, computation error increase with K_E, while computation time drops, as shown in Figure 10.

For all experiment cases, simulation is performed to find out the effect of K_E on computation time and accuracy. Figure 10A shows the results for D1-D4 (effect of dwell time). For case D1, with increase in K_E from 0.5 to 27, computation time reduces drastically from 90 to 35 min, but leads to an increase of error from 2 to 22%. For case D2, with respect to increase of K_E, computation time reduction is more drastic from 130 to 36 min and computation error increases from 3 to 21%. For case D3 and case D4 as well, computation time reduction is more as dwell time is increasing but computation error is also increasing.

Figure 10B shows the results for all cases P1-P4 (effect of laser power). For case P1-P4, with increase in K_E from 0.5 to 27, computation time reduces from 175 to 25 min, but leads to an increase of error from 2 to 21%. Computation time reduction is exponential (3–4 times reduction) when K_E is increased from 0.5 to 4. But then the trend in computation time becomes linear when K_E is increased from 4 to 27 (not even half). Computation error increases linearly with an increase in K_E. So, an intelligent compromise needs to be done to reduce the computation time but also keeping in mind the computation error as well. In this work, the objective is that computational error should be less than 10%, that is well accepted in the scientific and industrial community. Keeping this in mind, K_E = 4 seems to satisfy both the objectives of reducing the computation time drastically but also keeping computation error less than 10%. Further increase of K_E results in computation time reduction, but computation error exceeds 20% when K_E = 27. However, with the introduction of elongated ellipsoid (EE) heat source, it is noticed that there is a drop-down of the temperature at melt-pool scale (local scale) as well as part scale (global scale) e.g., thermocouple locations that is far away from the deposition region as shown in Figure 9. This modification of the thermal behaviour (especially, melt-pool scale) will surely impact the local mechanical response of the sample. Therefore, same analysis needs to be performed to verify if the optimised value of K_E found in thermal analysis is also valid for mechanical response e.g., distortion and residual stresses.

Due to the increase in K_E (sub-track length), computation time reduces exponentially, but it also leads to increase in computation errors. To compensate those errors, a variable source correction factor K_Q is introduced that leads to an increase in source intensity in elongated ellipsoid model, as shown in the equation below. Q EE =   6 3 APK Q a ^ bc π π exp ( − 3 ( x +   v s ( t 0   +   1 2 Δ t   ) ) 2 a ^ 2   − 3 y 2 b 2   − 3 z 2 c 2 ) (19) With the introduction of K_Q, distribution and peak intensity over the ellipsoid is increased artificially. Thus, the correct value of K_Q should be dependent upon the value of K_E, so calibration of source correction factor is done with respect to different sub-track size for all experiment cases. Then for the correct value of correction factor that should be dependent upon the size of sub-track ( K E   ), Calibration of correction factor is done with respect to sub-track size ( K E ) for all cases. K_E = 0.5 does not elongate the ellipsoid. So, to analyse the effect of correction factor, values of K_E starting from 1 to 27 is used, because K_E = 0.5 using Elongated Ellipsoid does averaging over a sub-track of R, i.e., SE or DE with a simulation time step of R / v s i.e., conventional simulation increment.

As it can be seen in Figure 11 with K_E values (1, 4, 8 and 27) shown in black vertical lines, for experiment cases D1-D4 (effect of dwell time), with an increase in K_E, computation error increases as we have seen in the previous section, source correction factor K_Q value is also increasing to compensate for this increase in computation error with an increase in K_E. Therefore, there is a direct relation between K_E and K_Q, that is helpful in achieving the second objective to reduce the computation error. With the introduction of K_Q, K_E = 8, 16 can also be utilised as computation error falls below 10% which was not the case previously. As it can be seen in Figure 12, for experiment cases P1-P4 (effect of laser power), same relation is established between K_E and K_Q. In this set of experiments also, K_E = 8, 16 can now be used in heat transfer analysis as computation error falls below 10%. So, with the introduction of K_Q, computation error is reduced along with reduction of computation time as higher values of K_E yields desired accurate results. With the introduction of correction factor, there is no significant impact on the computation time of analysis.

As it can be seen in Figure 13, for experiment case D1 with K_E = 4, with the introduction of correction factor K_Q using EE heat source model, temperature evolution at thermocouple location shifts upwardly due to the increase of heat intensity. K_Q values of 1.1 and 1.2 seems to be well compensating the effect of K_E that leads to underprediction of temperature evolution, but with further increase of K_Q values of 1.3, 1.4 and 1.5, it can be seen that it over-compensates the thermal error generated due to K_E, and starts to yield higher computation error. Therefore, a correlation between K_E and K_Q should be established for different process parameters that yields minimum computation error for elongated ellipsoid (EE) heat source.

For all experiment cases D1-D4 and P1-P4, it can be seen in Figure 14 that for most experiment cases, with an increase of K_E from 0.5 to 8, there is also an increase in correction factor from K_Q = 1 to 1.25, that leads to minimum computation error for temperature history at thermocouple location. But with further increase of K_E from 8 to 27, further increase in values of K_Q is not required. Therefore, a correlation between K_E and K_Q is recommended in this work that yield to minimum computation error.