1. Introduction
The development of additive manufacturing (AM), particularly direct energy deposition (DED) technologies, creates more and more opportunities with the lowest economic and time costs for enterprises of different industries in the manufacture of large-sized, non-standard metal products or parts with a complex shape, as well as in repair and restoration work [
1,
2,
3,
4]. In the DED process, a heat source is used to form a melt pool on the surface of a solid metal substrate. The filler material is fed into the melt pool, which causes it to melt. Thus, there is an increase in the volume of the molten material and the formation of the deposited bead. Then, the process is repeated layer-by-layer. There are various direct deposition technologies, differing in the type of heat source used and the filler material used. Technologies using laser and metal powder are well-developed and widespread. The powder is fed into the treatment area by means of protective gas (such as argon), most often either coaxially or multi-jet [
5]. A laser, electron beam, and arc can be used when using filler wire. When using a laser or electron beam, the wire is usually fed from the side, which imposes corresponding restrictions on the product shape [
6]. The process of using an electric arc as a heat source is similar to fusion welding [
7].
The formation of a product layer-by-layer means that the deposited material undergoes quasi-periodic heating and cooling processes, including partial remelting of already formed layers. The resulting non-stationary thermal conditions significantly affect the local microstructure, residual stresses, and deformations, as well as the distribution of defects [
8]. Gushchina et al. [
9] investigated the effect of the interpass temperature on the Ti-6Al-4V structure in the direct laser deposition (DLD) process. The temperature values varied from approximately 250 to 550 °C. A high interpass temperature corresponds to a lower cooling rate. It leads to the formation of an equilibrium structure, to the elongation values increased by two times, but the microhardness values are underestimated. In [
10,
11], reduced flow stress was found for material obtained by AM (selective laser melting (SLM), electron beam melting (EBM), and DED) compared to wrought material in compression tests at high temperatures. However, it was found in [
12] that the Young’s modulus of the DLD-processed Ti-6Al-4V alloy is approximately 46% higher at 500 °C than that of the wrought alloy. In the study [
13], deposited Ti-6Al-4V samples showed low values of relative elongation and uneven grain size across the sample width. On the other hand, when a wobbling or oscillating laser beam is used, an even distribution of grain sizes, a uniform distribution of hardness, and a high value of mechanical properties are already observed. This allows the use of such products without additional heat treatment. There are also attempts to develop maps of the microstructure evolution based on experimental results, which showed different areas of microstructures depending on the cooling rate and repeated thermal cycles [
14,
15]. In the AM process, defects are also encountered, such as non-fusion, the main cause of which is insufficient heat source energy while simultaneously feeding an excess amount of filler material into the melt pool [
16]. The incorrect choice of the mode can also lead to the formation of cracking and delamination [
8]. At the same time, the above effects can be combated by optimizing the process, which includes strict control of operating parameters, the presence or absence of preheating, control of the cooling rate, and the choice of a scanning strategy. Thus, the problem of determining transient temperature fields is not new, but it is still relevant. Experimental determination of temperature fields, cooling rates, as well as related regularities is time-consuming and costly due to the large number of technological parameters in the process. Generally, temperature measurement is carried out by a contact method based on thermocouples and a non-contact method based on infrared thermography [
17,
18,
19]. However, infrared thermography does not allow direct measurements of the object’s surface temperature. The accuracy of determining the surface temperature is limited by the unknown emissivity ε. Moreover, emissivity is not constant but depends on many factors, such as the material, surface temperature, surface roughness, viewing angle, and presence of oxide films [
20,
21]. Thermocouples can significantly improve the accuracy of measurements. For this reason, they are used both for measuring temperature and for calibrating the thermal camera.
Another way to determine three-dimensional temperature fields is to use different methods of mathematical modeling. One of the most common methods for calculating transient temperature fields in the AM process is the finite element method (FEM) [
22,
23]. The material deposition in the DED process is often modeled using quiet or inactive elements, which are activated as the added filler material solidifies [
24,
25]. In the quiet approach, the elements are present in the calculations but are assigned special properties. For this reason, such elements do not affect the result. In the inactive element approach, elements are not included in the calculations until the corresponding material has been added. Michaleris [
26] proposed a numerical model using the hybrid quiet/inactive element method. This approach assumes that the elements are initially inactive and then switch to quiet mode. As a result, equivalent results can be obtained with less computation time. Peyre et al. [
27] developed a three-stage DED model, taking into account the powder temperature and the deposited bead geometry. Before solving the heat transfer problem using FEM, the layer height is calculated based on the law of conservation of mass and energy, and the layer width is calculated by the iterative method. Heigel et al. [
28] compared the results of modeling temperature fields under various convection conditions. According to the results, the forced convection model gives the best agreement with the experimental data as compared to the free convection model. Chiumenti et al. [
29] showed that heat loss occurs mainly through thermal radiation. This effect cannot be neglected but it can be replaced by a temperature-dependent convective heat transfer coefficient that will exceed the free convection coefficient. A great deal of research is aimed at reducing the calculation time. A widespread approach is that the material is added either in parts of a layer, in whole layers at once, or in several layers at once [
23,
25]. Denlinger [
30] introduced layer-by-layer coarsening of the mesh at a certain distance from the heat source. Patil, Pal et al. [
31] used a fine mesh only locally near the heat source. With these strategies, it is possible to reduce the number of elements without losing the high accuracy of the solution.
There are also numerical approaches that take into account convective heat transfer and mixing in the melt pool [
32]. Using the control volume method, Manvatkar et al. [
33] demonstrated that neglecting the liquid metal flow leads to an overestimation of the temperature field—in particular, peak values of the temperatures. This result was confirmed by Gan et al. [
34], who also indicated that the Marangoni convection causes mixing in the molten pool, as a result of which the distribution of the composition evens out rapidly to a uniform value. This must be taken into account when solving the crystallization problem.
There is also an analytical approach to determining temperature fields [
35,
36,
37]. However, such works have several characteristic features. First, heat transfer to the environment is not taken into account, which, as is known, has a significant effect on temperature [
29]. In addition, the adiabatic boundary of surfaces (or convective heat transfer) is not specified in the calculations, or the adiabatic boundary is specified only on some surfaces. The physical interpretation of the adiabatic boundary is that the heat transferred to the body remains inside the body. The case in the absence of an adiabatic boundary corresponds to a body that is infinite in all directions; the case with adiabatic boundaries on some surfaces corresponds to a body with an infinite size in the direction where there is no adiabatic boundary. From a physical point of view, it corresponds to a different problem statement. In addition, no works have been found in the literature that would correctly apply the analytical solution to the determination of transient temperature fields over long periods.
The objective of this work is to develop a simple and fast but, at the same time, reliable method for calculating transient three-dimensional temperature fields in the AM process of multi-track thick-walled products. This work is a continuation of a previous work [
38], where the heat conduction problem was solved for single-pass thin-walled products. The developed calculation scheme assumes its further inclusion in the optimization system for selecting the technological parameters. It presupposes the processes of multiple direct calculation (direct problem). As a consequence, the issue of computation time is fundamental. Analytical methods are inferior in accuracy to numerical methods, but they make it possible to directly determine the relationship between individual parameters of a mode or a group of parameters. In such conditions, high accuracy of the calculation is not a fundamental point, since the process is extremely complex and some assumptions will have to be made in any case. In addition, due to the lack of accurate values of the thermophysical properties of materials in high-temperature regions, the calculation accuracy using any method is, in principle, limited. Thus, in decision support systems for technologists or engineers, the issue of quickly determining relationships comes to the fore. As a result of using the extended solution, it becomes possible to study the degree of influence of the process parameters and different heat source trajectories on the temperature distribution and cooling rate in the product. From a practical point of view, it is applicable because maintaining stable temperatures and melt pool sizes is one of the key means of controlling the process stability.
2. Methods and Model Description
2.1. Problem Statement
The following physical assumptions were made in the heat transfer model:
The thermophysical properties of the filler material and substrate are constant and temperature-independent;
Heat flux distribution of the heat source qh is presented as a surface normally distributed heat source;
Heat transfer occurs according to Newton’s law;
Phase transition is not considered;
The effect of convection of liquid metal is neglected.
Transient heat transfer throughout the model is determined by a three-dimensional linear heat conduction equation in a Cartesian coordinate system with appropriate boundary conditions. It is worth considering that heat is released into the environment from the wall side surfaces under conditions of convective heat transfer. Then, convective heat transfer can be represented as a heat sink
, where α—surface heat transfer coefficient, and
W—wall width. This approach considers only the decrease in the average temperature in the section, and does not consider the temperature unevenness along the wall width. The heat conduction equation is expressed as:
where
λ—thermal conductivity,
ρ—density, and
c—specific heat capacity.
The initial temperature distribution of the substrate and the wall at time t = 0 is equal to the ambient temperature T0.
Boundary conditions on the front surface of the computational domain are as follows:
where
qh(
x, y) is the heat flux density.
The adiabatic boundary is set on other surfaces of the wall and substrate where there are no heat sources.
2.2. Analytical Model of Transient Heat Transfer
The solution to Equation (1) can be obtained using the Green’s function method for the heat equation. Then, the temperature increment at an arbitrary point with coordinates
x,
y,
z at any time
t from an elementary point source (moving with speed
v) that acted at time
t′ on the surface of a semi-infinite body in a stationary coordinate system is equal to
dT:
where
q—point heat source power,
a—thermal diffusivity,
vx—the projection of the heat source speed
v onto the
x-axis (
Figure 1),
vy—the projection of the heat source speed
v onto the
y-axis, and
—coefficient of heat loss.
The solution for a moving point source can be obtained by taking the integral over the source action time
t′ from
t1 to
t2 (the source started acting at time
t1 and finished acting at time
t2):
where
,
—distance from the heat source to the considered point,
,
, and
t >
t2 ≥
t1 ≥ 0.
Let us pass from the solution for a semi-infinite body to the solution for a bounded body. In order to satisfy homogeneous boundary conditions of the second kind in the form of the absence of heat flow through the wall boundaries, it is necessary to continue the function, defining the heat source as an even function relative to the body boundaries. Then, the resulting even function is continued periodically on the entire numeric axis. The physical interpretation of this action consists of the introduction of imaginary heat sources symmetrically relative to the body’s boundaries. In this case, the planes of symmetry are the planes
x = 0,
x =
L*, where
L* is the wall length, from the planes
z = 0 and
z = H, where H is the wall height (sum of substrate height
Hs and wall height
Hw), and also from the side wall boundaries
y = 0,
y =
W.
Figure 2 and
Figure 3 show the schemes of introducing imaginary sources along the
x-axis and
y-axis for a thick wall and a closed wall. The scheme for a closed wall (in particular, a cylindrical wall) can be obtained by unwrapping the wall around one of its generatrices (
Figure 3b). The temperature field is calculated at an arbitrary point
P.
Then, the temperature field in the deposited wall is defined as an infinite series of corresponding solutions for an unbounded semi-infinite body (in this case, the source is usually distributed over the computational domain surface):
where
,
xs,
ys—the coordinates of the point source in the coordinate system associated with the source,
—for the case of a single wall,
—for the case of a closed wall,
,
k = −1, 1—for the case of a single wall,
k = 1—for the case of a closed wall,
L* = L—for the case of a single wall, and
L* = 2πRw—for the case of a closed wall. Summation over
n considers the limited length, while for summation over
j and
p—over width and height, respectively.
By shifting the origin for each pass in the AM process, it is possible to set the times t1 and t2 in such a way that t1 = 0 always, and . In this case, t1 and t2 are not arguments to the dT function.
If the source acts over an area with radius
Rb, then it is necessary to integrate Equation (5) over this area to obtain the temperature field. Then, the temperature Δ
Tpreh(
x,y,z,t,t1,
t2) can be obtained as:
The number of series terms in Equation (5), which must be taken into account in the calculations, depends on two aspects: the value of the considered heating or cooling time t in the DLD process, and the value of the thermal diffusivity of the material a. Thus, the number of series terms is directly proportional to .
It is worth noting that the variable R in Equation (5) can have a large value in modulus. This leads to the fact that in the product the multiplier tends to be a large value, which is taken equal to infinity (infinitely large value). This result occurs when using mathematical and computational software in the calculation process. The multiplier tends toward 0 (infinitesimal value). Consequently, an indeterminate form of the type (0 × ∞) arises in Equation (5).
For the evaluation of the indeterminate form, let us use, for example, the known approximation of the error function erf (
x) using elementary functions 7.1.26 [
39]:
where
,
,
,
p = 0.3275911,
a1 = 0.254829592,
a2 = −0.284496736,
a3 = 1.421413741,
a4 = −1.453152027, and
a5 = 1.061405429.
Then, it is easy to eliminate the arising indeterminate form using this replacement and performing elementary transformations.
2.3. Influence of the Substrate on the Temperature Field
At the initial stage of the deposition process, the conditions for the formation of layers differ as their number increases. The heat accumulation in the wall during the initial deposition period causes different heat removal conditions. A steady state of the process is achieved when the wall temperature stops rising, and the conditions for heat removal become constant. At the same time, the conditions for the formation of each subsequent layer, all other factors being equal, do not differ from the conditions for the formation of the previous layer.
It is convenient to calculate temperature fields using the Green’s function method for bodies of simple geometric shape. The shape of the wall deposited onto the substrate is not simple but rather T-shaped. In this regard, several approximate computational schemes for calculating temperature fields were considered. According to the first scheme, the side parts of a substrate that are outside the wall construction area are removed (
Figure 4). To take into account the removed mass of the substrate, heat sinks of the corresponding energy are introduced. One of the variants of this approach is described in [
38].
The second approach consists in replacing the real horizontal substrate with a vertical substrate (rib), the section of which coincides with the wall section. It should be borne in mind that the masses of these substrates must be equivalent to each other. However, this approach is applicable only when the wall width exceeds the thickness of the substrate by two times or more. This is due to the fact that the heat removal conditions change with a smaller wall width.
The third approach is a combination of the first two. In this case, the vertical substrate height is half the horizontal substrate width. The remaining removed part of the substrate mass is taken into account by introducing heat sinks.
The temperature increment in the equilibrium state after the next pass is
where
L′—the trajectory length of the heat source on the current pass,
Vs—the substrate volume, and
Vw—the wall volume.
Based on the condition that the increment in the wall temperature in the absence of a substrate corresponds to the increment in the wall temperature in the presence of a substrate, it is possible to determine the total energy of the sinks at each pass:
where
Vs′—truncated substrate volume, and
Vw(
n)—wall volume on the
n-th passage.
The power of each fixed sink is determined based on the time of their action
tsn and the total energy of all sinks
En. Let us take into account that the power of each sink is not constant, but decreases exponentially in time according to the following law:
The action time of the sinks is proportional to the distribution time of uneven temperature; that is,
, where
Rn is the characteristic size of the deposition wall after
n passes. The maximum power of each sink can be determined using the expression
where
m—number of sinks.
The inertia of the heat transfer process in the wall is taken into account due to the introduction of the delay time
of the action of the sinks in relation to the actual heat source action (see
Figure 4), while
.
To obtain the equation of the temperature field as a result of the action of sinks, suppose
v = 0 in Equation (3):
where
,
Another way to take into account the removed mass of the substrate is to introduce a volumetric sink (
m = 1), the temperature field of which can be determined by the following equation:
The total temperature field
Tn after the deposition of the
n-th number of layers can be represented as a superposition of temperature fields from the action of heat sources and the action of sinks at each pass. In this case, the process parameters (power, cladding speed, trajectory, pause time) can differ at each pass. Then, the heating temperature can be calculated using the following equation:
where
tpause—the pause time between passes, index
i = 1 corresponds to the last pass, and the index
i =
n corresponds to the first pass.
It should be noted that the above approaches to determining the temperature field are applicable since the difference in heat removal condition is insignificant. As a result, it can be ignored. Let us demonstrate this by the example of the repeated heating of a steel substrate with a size of 200 × 100 × 10 mm (length × width × height) in its middle part. The calculation results are compared with the results of using approximate calculation schemes. Case 0 corresponds to the heating of the horizontal substrate in its middle part; case 1 is a vertical substrate (rib) with a wall width
w0 =
h0; case 2 is a truncated substrate with a wall width
w0 =
h0; case 3 is a vertical substrate with a wall width
w0 = 2
h0; case 4 is a vertical substrate with a wall width
w0 =
h0/2 and a height equal to the half-width of the substrate (
Figure 5).
The difference between the curves for cases 2, 3, and 4 compared to the heating of the real substrate (case 0) is insignificant. As mentioned earlier, case 1 has reduced heat removal. The temperature deviation is observed upward at a short time interval for cases 2 and 4 after the passage of the temperature peak. The smaller the wall width, the greater the difference between the actual temperature and the calculated one. However, all curves converge to the same value. Moreover, this effect is observed only at the initial deposition stage, when the wall height is insignificant in comparison with the substrate thickness. The curves completely coincide for the wall width w ≥ 2h. The real wall will increasingly approach the approximate calculation scheme with an increase in the number of layers. Thus, in this work, we will use a calculation scheme with the replacement of the horizontal substrate by a vertical one at w ≥ 2h and a calculation scheme with the replacement of the horizontal substrate by a vertical one with a height equal to the half-width of the substrate and heat sinks at w < 2h.
Let us show that the proposed calculation method also satisfies the law of conservation of energy. Consider several arbitrary modes (10 passes) of heating a body with different powers (pure powers), heat source movement speeds, and pause times between passes (
Table 1). However, the total input amount of heat energy remains the same for all cases. The heating body is a steel substrate and a truncated substrate (similar to case 0 and case 2). The truncated substrate width is equal to the half-width of the actual substrate. In this case, heat transfer to the environment is absent, so we will consider the completely adiabatic boundaries of all substrate surfaces. The choice of arbitrary modes will make it possible to draw conclusions about the independence of the result from the specific parameters of the mode. The calculated results are compared with the steady-state temperature (
Figure 6) obtained from the equation
.
Despite the different heating modes (the total invested thermal energy is the same), the body temperature in the calculations has always tended to its equilibrium state at a given energy value. Thus, the coincidence of the temperatures of all curves (in the absence of heat exchange with the environment) indicates the consistency and adequacy of the developed method for calculating non-stationary temperature fields.
4. Results and Discussion
The thermocouple data were used to verify the developed transient heat transfer model. For each sample, two thermocouples were installed at different distances from the substrate. The position of the thermocouples is shown in
Figure 7. For mode 1—
H1 = 12 mm (15 layers),
H2 = 78 mm (96 layers); for mode 2—
H1 = 10 mm (12 layers),
H2 = 73 mm (91 layers); for mode 3—
H1 = 12 mm (15 layers),
H2 = 73 mm (91 layers). The modeling parameters were taken according to
Table 2. The calculations used the absorption coefficient of laser radiation equal to 0.35, the convective heat transfer coefficient for the mode without a pause was 70 W·K
−1·m
−2, and that for the modes with a pause was 35 W·K
−1·m
−2. The values of the thermophysical properties were taken from [
40] and were adopted as constant, and averaged in the temperature interval between the interlayer temperature of the previous layer and the melting point.
Figure 8 shows thermal cycles obtained using thermocouples (point
H1 and
H2), as well as calculation cycles for modes with a pause time between passes of 0, 10, and 25 s. However, it is worth clarifying that a pause of 0 s (as with 10 and 25 s) is a pause that is set in the control program of the robot. The actual pause time was more than 0 (more than 10 and 25 s) since the robot needs a certain time to return from the point of the end of the trajectory to the point of the beginning of the trajectory.
Thermocouples were installed at the deposition of 15th layer (point
H1) and 96th layer (point
H2) for the mode without a pause. For point
H1, the thermocouple readings were considered from the moment the deposition of the 16th layer begins (
Figure 8a), which corresponds to 106 passes. For point
H2, the thermocouple readings were considered from the moment the deposition of the 97th layer begins (
Figure 8b), which corresponds to 673 passes. Similarly, data were obtained in modes with a pause of 10 and 25 s. However, with a pause of 10 s, thermocouples were installed during deposition of 12 layers and 91 layers, and thermocouples were installed during deposition of 15 layers and 91 layers with a pause of 25 s.
The decrease in the temperature of the fixed points of the wall can be observed over time with a pause time between passes of 0 s (
Figure 8a,b). A less intense decrease in temperature is also observed with a pause of 10 s (
Figure 8c,d). However, the temperature of the fixed points does not decrease with a pause of 25 s but remains constant (
Figure 8e,f), which is completely consistent with the experimental results. A decrease in temperature at a fixed point, as in the case of a mode without a pause, indicates high temperature gradients along with the wall height. On the contrary, in the case of a pause of 25 s, when the temperature hardly changes, the temperature gradient along the height is almost absent. Knowing the change in temperature over time, heating and cooling rates can be found.
Figure 9 shows the calculated cooling rates in the near-surface layer of the wall (point
H1) during the deposition of 106th pass (16th layer) at different pause times. The highest cooling rate is observed in the mode with the longest pause. Thus, the longer the pause time, the faster the heating and cooling rate. This is directly related to the value of the interpass or interlayer temperature. At low temperatures, the heat removal process proceeds faster. Cooling rates may differ by 30% according to calculations.
Let us consider the change in the interlayer temperature depending on the pause time between passes. The temperature of the first five layers after attaching the thermocouple was used to fix the interlayer temperature since the interlayer temperature is actually the wall surface temperature. Moreover, at a shallow depth from the surface, the temperature has time to equalize. The temperature values that were used to determine the interlayer temperature are highlighted with red dots on thermal cycles (see
Figure 8a–f). In this case, the temperature values were used after the sixth pass, since, after the seventh pass, outliers are observed. Therefore, these points are ignored.
Figure 10 shows the calculated curves of the interlayer temperature depending on the number of layers and the pause time between passes with plotted experimental points.
The temperature is significantly higher when depositing without a pause, compared to modes with a pause. The highest rate of increase in the interlayer temperature is observed when the first layers are deposited near the substrate. The presence of a pause leads to the fact that the temperature rises rapidly to around 10 layers, after which it changes insignificantly. However, the rapid increase in temperature continues after the 10th layer and almost stops after the 25th layer in the non-pause mode. The relative difference in the interlayer temperature in comparison with experimental points in all modes and at various distances from the substrate does not exceed 15% in absolute value. At the same time, deviations are observed both in the direction of overestimated values and in the direction of underestimated values.
Figure 11 shows the calculated distribution of the wall surface temperature during the depositing of the second, fourth, and sixth passes of the 16th layer. It can be seen that the closer the melt pool is located to the side surface of the wall, the more the isotherms elongate. More elongated isotherms, other factors being equal, correspond to a lower cooling rate. This may cause the appearance of heterogeneity of grain sizes and properties along the wall width.
It should be noted once again that since the specificity of additive manufacturing technologies presupposes repetitive heating and cooling processes, ignoring convective heat transfer with the environment leads to a significant error. In addition, the time spent on the above calculations is on the order of several minutes or several tens of minutes using a personal computer. However, functional–analytical methods make it possible to obtain the value of a function at a certain moment in time; without calculating the values of the function in previous moments, the calculations are seconds.
5. Conclusions
The three-dimensional transient heat transfer model was developed for the DED processes of multi-track walls. The proposed model makes it possible to reproduce temperature fields at various values of the technological process parameters. This model is based on the obtained three-dimensional solution of the non-stationary heat conduction equation for a moving distributed heat source, considering convective heat transfer to the environment. Since the model made it possible to calculate the temperature change both in space and in time, the model was allowed to calculate all the temperature field characteristics (heat cycles, temperature gradients, and cooling rates). In this case, the wall and the substrate sizes, the change in power from layer to layer, the change in the cladding speed, the interpass dwell time, and the heat source trajectory were taken into account.
The relative difference in the interlayer temperature fluctuated within 15% when comparing the experimental and calculated data. The temperature change over time at fixed points of the samples was completely consistent with the experimental data. A decrease in temperature at a fixed point, as in the case of a mode without a pause, indicates the presence of a temperature gradient along the height of the wall, while it is practically absent with a pause of 25 s. According to calculations, the highest cooling rate is observed in the mode with the longest pause. Thus, the heating and cooling rate is proportional to the value of the pause time between passes and, consequently, inversely proportional to the values of interpass temperature.
Considering the features of beam or non-beam heat sources, it is possible to adapt the model to various direct deposition technologies. In addition, based on the proposed calculation scheme, it is possible to develop a method that takes into account any curved trajectories of the heat source.