Analytical Solution of the Non-Stationary Heat Conduction Problem in Thin-Walled Products during the Additive Manufacturing Process

The work is devoted to the development of a model for calculating transient quasiperiodic temperature fields arising in the direct deposition process of thin walls with various configurations. The model allows calculating the temperature field, thermal cycles, temperature gradients, and the cooling rate in the wall during the direct deposition process at any time. The temperature field in the deposited wall is determined based on the analytical solution of the non-stationary heat conduction equation for a moving heat source, taking into account heat transfer to the environment. Heat accumulation and temperature change are calculated based on the superposition principle of transient temperature fields resulting from the heat source action at each pass. The proposed method for calculating temperature fields describes the heat-transfer process and heat accumulation in the wall with satisfactory accuracy. This was confirmed by comparisons with experimental thermocouple data. It takes into account the size of the wall and the substrate, the change in power from layer to layer, the pause time between passes, and the heat-source trajectory. In addition, this calculation method is easy to adapt to various additive manufacturing processes that use both laser and arc heat sources.


Introduction
Additive technologies, in particular direct energy deposition (DED) technologies, are actively developing and are already used in modern production in the manufacture of metal parts and repair and restoration work [1][2][3][4][5]. Laser, electron beam, plasma, and arc are used as the main heat source, and the filler material can be used both in the form of powder and in the form of wire. A wide range of materials can be used: steel, titanium, aluminum, nickel alloys, and composites [6][7][8][9][10][11].
A characteristic feature of the direct deposition process is that the material undergoes multiple heating and cooling processes, including partial remelting of already-formed layers. Such a complex temperature field, which changes both in space and in time, significantly affects the local microstructure and its evolution, residual stresses, and deformations, as well as the distribution of defects, which can significantly affect the mechanical properties, and hence the entire product's reliability.
The problem of determining transient temperature fields can be investigated using experiments and mathematical modeling. Experimental determination is time-consuming and costly due to the large number of the process operating parameters. Typically, thermocouples and infrared thermography are used to determine temperature fields [12][13][14][15][16]. However, it is not possible to directly measure the object temperature using thermography. The accuracy of determining the surface temperature is limited by the unknown emissivity ε. In addition, emissivity can change during additive manufacturing (AM), as it depends on material, temperature, viewing angle, surface roughness, and presence of oxide films [17,18]. As a consequence, the instrument-calibration process is extremely difficult. The use of thermocouples makes it possible to determine the transient temperature much more accurately; however, thermocouples must be thin enough, and at the same time they only allow measuring the temperature locally, even if there are several of them.
Various methods of mathematical modeling are used to determine the three-dimensional temperature fields. Currently, one of the most common methods for calculating transient temperature fields in the AM process is the finite element method [19][20][21]. The objective of such calculations is usually to determine the temperature field evolution, temperature gradients, and their effect on residual stresses, which is related to a thermomechanical problem.
The main difficulty of modeling DED processes for large products is associated with the large temporal and spatial domain of the calculation. To determine the most accurate transient temperature fields, the time step should be microseconds, while the total deposition time is many hours. The melt pool has a characteristic size of the order of a few millimeters, while the model of the entire product is usually on a meter or sub-meter scale. As a consequence, such calculations require enormous computing power and time.
There are several approaches to solve such problems. The first approach assumes minimal time steps and a fine mesh of the part, which leads to fairly accurate temperature values, but the time spent on computations can be huge [22,23]. The second approach assumes a scheme according to which the material is added either in parts of a layer (a hatch-by-hatch), or in whole layers at once (layer-by-layer), or in several layers at once [21,24]. In this case, the deposited energy for a period of time corresponding to the trajectory traversed is distributed throughout the added material. Thus, the simulation time can be significantly reduced. As a result, only the history of the average temperature is recorded, but not the local thermal history.
In some works, the finite volume method is used to determine transient temperature fields, which also takes into account convective heat transfer. Due to the complexity of such calculations, the simulated samples are bodies with characteristic dimensions of the order of several tens of millimeters [25][26][27]. As already mentioned, real products have dimensions of the order of a meter, which ensures that the process conditions are different from the simulated ones, in particular, the temperature, cooling rates, and temperature gradients will differ.
In comparison with the above models, the peculiarity of this work is to provide a fast, simple, and universal, but at the same time reliable, analytical method for calculating transient temperature fields in the AM process. Maintaining stable temperatures and meltpool sizes is one of the key means of controlling the process stability. Then, as a result of using this calculation method, it becomes possible to theoretically determine the influence of the mode parameters on the formation of the deposited layers, as well as to select stable process modes. In addition, the application of this method will make it possible to study the degree of influence of local quasiperiodic temperature fields, which always accompany direct deposition processes, on structural phase transformations.

Problem Statement
The following physical assumptions were made:

•
The physical properties of the substrate and the filler material (specific heat capacity c, density ρ, thermal conductivity λ, thermal diffusivity a) are temperature-independent.

•
The effect of convection of liquid metal is not considered.

•
Heat flux distribution of the heat source q h is presented as a surface normally distributed heat source.

•
Heat transfer occurs according to Newton's law.
In order to obtain the temperature field in the substrate and the deposited layers, it is necessary to solve the following linear non-stationary heat-conduction problem in a Cartesian coordinate system x, y, z: where λ-thermal conductivity of the material, c-specific heat capacity, and ρ-density.
The initial temperature of the substrate is equal to the ambient temperature: Boundary conditions on the front surface of the computational domain: where q h (x, y) is the heat flux density. The adiabatic boundary is set on other surfaces where there are no heat sources.
Heating of a product in the AM process is described as the action of a surface elliptical heat source with a power density q h (x, y). In the x0y plane, the power-density distribution is described by the Gaussian function: where Q h -heat-source power, η-heat efficiency, R H -effective radius of the heat source, and β-tilt angle.

Analytical Model of Non-Stationary Heat Transfer
The temperature increment at an arbitrary point with coordinates x, y, z (in a fixed coordinate system) at any time t from an elementary point source that acted at time t' on the surface of a semi-infinite body is known and is equal to dT [28]: where q-point heat-source power, v-heat-source moving speed (cladding speed), and a-thermal diffusivity, equal to a = λ/(c·ρ). When calculating the repeated heating and cooling process of thin-walled products, it is impossible to neglect heat transfer to the environment, since it leads to a noticeable error in determining the temperature. The smaller the deposited wall thickness, the greater the heat-transfer effect to the environment. Considering the above, it is necessary to obtain a non-stationary equation for the heat-propagation process, taking into account heat transfer to the air.
Let it be supposed that the heat transfer occurs according to Newton's law, but the temperature along the thickness is equalized instantly. Then the heat transfer from the wall side surfaces is taken into account by introducing the multiplier e −b(t−t ) into Equation (5). It means only a decrease in the average temperature in the section, but does not consider the temperature unevenness along the wall thickness. Thus, heat transfer is equivalent to a volumetric heat sink, while the condition of the adiabatic boundary is still satisfied. Then Equation (5) takes the form: where b = 2α cρh -coefficient of heat loss, α-coefficients of surface heat transfer, and h-wall thickness.
A moving heat source can be represented as elementary instantaneous sources acting sequentially and displaced relative to each other. Let us sum up the temperature increments from all elementary sources that acted in the general case during the time from t 1 to t 2 and make elementary transformations: where R(x, y, z) = (x − νt) 2 + y 2 + z 2 -distance from the heat source to the considered point of the body, t-considered the moment in time, t 1 -the start time of the source action, t 2 -the time when the source ends its action, and t > t 2 ≥ t 1 ≥0. A sufficient condition for determining the temperature field when the source has not yet stopped its action (when t = t 2 ) is the difference between t and t 2 by an infinitesimal value o(t).
To obtain the temperature field at any time ∆T (x, y, z, t, t 1 , t 2 ), it is necessary to calculate the integral in Equation (7) with the limits of integration t 1 and t 2 . For this, the integral in Equation (7) is represented as the difference of two integrals. Then the solution for a moving point source can be obtained using the substitution u 2 = 1/ √ t − t and the known integral 1.3.3.20 [29]: where Let us now consider the effect of the limited wall size under the assumption that its boundaries are adiabatic. This assumption allows the use of the method of images. To do this, it is necessary to mirror the actual heat source and each mirrored source from 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 (H = H s + H w ), and also from the side wall boundaries y = W/2, y = −W/2, where W is the wall width. As a result, we obtain a system of an infinite number of heat sources. A cylindrical wall (generally a closed wall) can be represented as a single wall by unwrapping the wall around one of its generatrices ( Figure 1). Figures 2 and 3 show the schematic of the reflection of sources along the x-axis for a single wall and a closed single wall, respectively. The red color denotes imaginary sources for which k = 1, and the blue color denotes that k = 1. The temperature field is calculated at an arbitrary point p.
is the wall height (H = Hs + Hw), and also from the side wall boundaries y = W/2, y = −W/2, where W is the wall width. As a result, we obtain a system of an infinite number of heat sources. A cylindrical wall (generally a closed wall) can be represented as a single wall by unwrapping the wall around one of its generatrices ( Figure 1). Figures 2 and 3 show the schematic of the reflection of sources along the x-axis for a single wall and a closed single wall, respectively. The red color denotes imaginary sources for which k = 1, and the blue color denotes that k = 1. The temperature field is calculated at an arbitrary point p.   Then, the temperature field in a bounded wall is determined by the following sum of the corresponding solutions for an unbounded semi-infinite body: -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 k and n considers the limited length; while for summation over j and p, over width and height, respectively.
Let us integrate Equation (9) from t1 to t2, repeating all the steps that have been done to obtain Equation (8) and take into account that the source can be distributed over the surface of the computational domain. The result is the equation:  Then, the temperature field in a bounded wall is determined by the following sum of the corresponding solutions for an unbounded semi-infinite body: -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 k and n considers the limited length; while for summation over j and p, over width and height, respectively.
Let us integrate Equation (9) from t1 to t2, repeating all the steps that have been done to obtain Equation (8) and take into account that the source can be distributed over the surface of the computational domain. The result is the equation: Then, the temperature field in a bounded wall is determined by the following sum of the corresponding solutions for an unbounded semi-infinite body: where X = k(x − 2nL * )-for the case of a single wall; X = k(x − nL * )-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πR w -for the case of a closed wall. Summation over k and n considers the limited length; while for summation over j and p, over width and height, respectively.
Let us integrate Equation (9) from t 1 to t 2 , repeating all the steps that have been done to obtain Equation (8) and take into account that the source can be distributed over the surface of the computational domain. The result is the equation: dT(x, y, z, xs, ys, t, t 1 , where R = (X − vt − xs) 2 + (y + jW − ys) 2 + (z + 2pH) 2 ; and xs, ys are the x, y coordinates, respectively, of the point source in the coordinate system associated with the source. By shifting the origin for each pass in the AM process, it is possible to set the times t 1 and t 2 in such a way that t 1 = 0 always, and . In this case, t 1 and t 2 are not arguments to the dT function. When using a distributed heat source, it is necessary to integrate Equation (10) over the source area (the area radius is equal to R b ). Then the heating temperature ∆T preh (x, y, z, t, t 1 , t 2 ) can be obtained as: dT(x, y, z, xs, ys, t, t 1 , t 2 ).
The function series in Equation (10) generally converge rapidly, so in practice, it is possible to limit the series to the first few terms. However, the longer the considered heating or cooling time t and/or thermal diffusivity a, the greater the number of series terms must be considered. In other words, the number of series terms is directly proportional to √ 4at. The criterion for choosing the number of series terms is the fact that the tangent of the slope of the tangent line to the temperature-distribution curve along the normal to the surfaces is equal to 0 at the adiabatic boundary, which corresponds to an angle of 0 • .
It should be noted that the practical application of Equation (10) is not convenient, since R can have a large value in absolute value, then in the product exp RvB , the multiplier exp RvB 2a tends toward a large value, which is taken as infinity (infinitely large value) when using mathematical and computational tends toward 0 (infinitesimal value).
Consequently, an indeterminate form of the type (0 × ∞) arises in Equation (10). 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 [30]: where k = 1 1+px , |ε(x)| ≤ 1.5 · 10 −7 , x = Using Equations (10) and (12) and performing elementary transformations, it is easy to get rid of the indeterminate form.

Influence of the Substrate on the Temperature Field
The next step is to consider the effect of the substrate and the wall height on the deposited wall temperature. The conditions for the formation of the layers differ as to their number increases. It is primarily due to the different heat-removal conditions caused by the heat accumulation in the wall and substrate during the initial deposition period and by an increase in the wall temperature. Upon reaching a steady state of the deposition process, the wall temperature stops rising.
It is convenient to use Green's function method to calculate the temperature field for bodies with a simple geometric shape. The cross-section of the wall deposited on the substrate has a T-shape, which complicates the use of this method.
In this regard, an equivalent computation scheme for calculating temperature fields was considered. Imagine a wall on a substrate as just a wall, removing the side parts of the substrate (Figure 4). To take into account the removed mass of the substrate, let us introduce heat sinks with certain energy equal to that accumulated in the substrate during deposition. These fixed sinks are activated as the real heat source moves, thereby the introduced sinks simulate the presence of a substrate. the substrate (Figure 4). To take into account the removed mass of the substrate, let us introduce heat sinks with certain energy equal to that accumulated in the substrate during deposition. These fixed sinks are activated as the real heat source moves, thereby the introduced sinks simulate the presence of a substrate. The temperature increment in the equilibrium state after the next pass is: where Vs-the substrate volume, and Vw-the wall volume.
To find the energy of the sinks, it is necessary 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. Thus, the total energy of the sinks in the nth passage is determined from the expression: where Vs'-truncated substrate volume, and Vw (n)-wall volume on the nth passage.
To find the power of each fixed sink, it is necessary to know the time of their action tsn and the total energy of all sinks En, which is already known. 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. So, the power of each sink is: where m-number of sinks. The total energy of the sinks decreases with an increase in the number of layers, and the action time increases, which indicates a decrease in the influence of the substrate as the wall height increases.
To take into account the inertia of the heat-propagation process, the sinks act with a delay n ts Δ relative to the real heat-source action (see Figure 4), while The equation describing the temperature field created by fixed-point sinks on the surface of a semi-infinite body can be obtained by setting v = 0 in Equation (8). Then, the The temperature increment in the equilibrium state after the next pass is: where Vs-the substrate volume, and V w -the wall volume.
To find the energy of the sinks, it is necessary 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. Thus, the total energy of the sinks in the nth passage is determined from the expression: where V s -truncated substrate volume, and V w (n)-wall volume on the nth passage.
To find the power of each fixed sink, it is necessary to know the time of their action ts n and the total energy of all sinks E n , which is already known. The action time of the sinks is proportional to the distribution time of uneven temperature; that is, ts n ∼ R n 2 4a , where R n is the characteristic size of the deposition wall after n passes. So, the power of each sink is: where m-number of sinks. The total energy of the sinks decreases with an increase in the number of layers, and the action time increases, which indicates a decrease in the influence of the substrate as the wall height increases.
To take into account the inertia of the heat-propagation process, the sinks act with a delay ∆ts n relative to the real heat-source action (see Figure 4), while ∆ts n ∼ H 2 n 4a . The equation describing the temperature field created by fixed-point sinks on the surface of a semi-infinite body can be obtained by setting v = 0 in Equation (8). Then, the equation describing heat propagation, considering the limited wall size, has the form: , i f t ≤ ts + ∆ts; ts + ∆ts, i f t > ts + ∆ts; In Equation (16), the summation over i is additionally introduced, which considers the action of each activated sink at the considered moment of time.
Taking into account the linearity of the thermal problem, the temperature field of heating on the nth layer T n is presented as the sum of the temperature fields as a result of the heat-source action and the heat-sink action at each pass.
In this work, the temperature field T n after the deposition of the nth number of layers, taking into account the linearity of the thermal problem, is represented as a temperature field as a result of the heat-source action depositing the nth layer, in front of which are n−1 heat sources, and also the action of sinks. These heat sources and sinks have equal or different power and operate at equal or different intervals of time, depending on the deposition strategy or scheme. Then, the heating temperature can be calculated using the following equation: where t pause -the pause time between passes, index i = 0 corresponds to the last pass, and index i = n − 1 corresponds to the first pass.

Results and Discussions
To verify the developed non-stationary heat-transfer model, two deposition processes were chosen and simulated. The first process was the DLD (direct laser deposition) process of a thin wall using Ti-6Al-4V powder [31], and the second process was the WAAM (wire arc additive manufacturing) process of a cylindrical thin wall using H08Mn2Si steel wire [32]. Type K thermocouples were used to measure the temperature. In the case of the DLD, the thermocouple was located on the substrate surface opposite the center of the wall at a distance of 0.5 mm. In the case of WAAM, the thermocouple was located on the substrate surface at a distance of 3 mm from the cylindrical wall. Both experiments were accompanied by calculations using numerical methods. The parameters of the modes are presented in Table 1. Details of the conditions for setting up experiments, placement of thermocouples, simulation parameters, and the values of thermophysical properties are described in detail in the corresponding works. It is worth mentioning that in the simulation in this work, the values of properties were adopted as constant, and averaged in the temperature interval between the ambient temperature and the melting point.  Figure 5 shows the compared results of the experimental data of temperature, as well as the data obtained using numerical methods with the calculations performed according to this work. Figure 6 shows the calculated temperature distribution in the deposited wall in its longitudinal section when the 20th layer was cladded. Figure 7 shows the comparing results of the experimental data of cylindrical wall temperature, as well as the data obtained using numerical methods with the calculations made according to this work.   Figure 5 shows the compared results of the experimental data of temperature, as well as the data obtained using numerical methods with the calculations performed according to this work. Figure 6 shows the calculated temperature distribution in the deposited wall in its longitudinal section when the 20th layer was cladded. Figure 7 shows the comparing results of the experimental data of cylindrical wall temperature, as well as the data obtained using numerical methods with the calculations made according to this work.     Figure 5 shows the compared results of the experimental data of temperature, as well as the data obtained using numerical methods with the calculations performed according to this work. Figure 6 shows the calculated temperature distribution in the deposited wall in its longitudinal section when the 20th layer was cladded. Figure 7 shows the comparing results of the experimental data of cylindrical wall temperature, as well as the data obtained using numerical methods with the calculations made according to this work.    The calculated transient temperature fields agreed satisfactorily with the experimental results both in the heating phase and in the free-cooling phase. In Figure 5, the calculated curve shows an increased temperature-change amplitude in comparison with the experimental data after 75 s. This can be explained by the fact that the experimental sample had a more complex cross-sectional shape, which allowed the heat flux to diverge in all directions at the junction of the wall and the substrate. Thus, it led to a decrease in the temperature change amplitude in the region below the plane of junction between the wall and the substrate. When measuring the wall temperature directly, this effect did not appear. It is worth mentioning that the calculations were performed under cardinally different conditions; namely, the heat sources, the process parameters, the movement trajectory, and the pause time were different. This fact suggests that the developed model made it possible to reproduce temperature fields in a wide range of technological process parameters and for various product configurations. Ignoring heat transfer to the environment greatly distorts the temperature, and given the specificity of additive manufacturing technologies in the form of multiple heating and cooling processes, it leads to a significant error. In addition, this method made it possible to calculate temperature fields in times in the order of several minutes using a personal computer. For this reason, this calculation method can be used to predict the temperature of large parts.

Conclusions
The three-dimensional, non-stationary heat-transfer model was developed for the direct deposition processes of thin-walled parts with various configurations. The calculated transient temperature fields obtained using the developed model agreed satisfactorily with the experimental results, both in the heating phase and in the free-cooling phase. At the same time, the presented calculations were performed with cardinally different input parameters of the model; namely, the types of heat sources, the process technological parameters, thermophysical properties of materials, movement trajectories, and pause times were different. These facts suggest that the model allowed reproducing three-dimensional temperature fields in a wide range of technological process parameters and for various product configurations.
To develop the model, a three-dimensional analytical solution of the non-stationary heat-conduction equation for a moving distributed heat source was obtained, taking into account heat transfer to the environment. The model made it possible to calculate all the temperature field characteristics; in particular, thermal cycles, temperature gradients, and the cooling rate in a path during the deposition process at any time. In this case, the wall The calculated transient temperature fields agreed satisfactorily with the experimental results both in the heating phase and in the free-cooling phase. In Figure 5, the calculated curve shows an increased temperature-change amplitude in comparison with the experimental data after 75 s. This can be explained by the fact that the experimental sample had a more complex cross-sectional shape, which allowed the heat flux to diverge in all directions at the junction of the wall and the substrate. Thus, it led to a decrease in the temperature change amplitude in the region below the plane of junction between the wall and the substrate. When measuring the wall temperature directly, this effect did not appear. It is worth mentioning that the calculations were performed under cardinally different conditions; namely, the heat sources, the process parameters, the movement trajectory, and the pause time were different. This fact suggests that the developed model made it possible to reproduce temperature fields in a wide range of technological process parameters and for various product configurations. Ignoring heat transfer to the environment greatly distorts the temperature, and given the specificity of additive manufacturing technologies in the form of multiple heating and cooling processes, it leads to a significant error. In addition, this method made it possible to calculate temperature fields in times in the order of several minutes using a personal computer. For this reason, this calculation method can be used to predict the temperature of large parts.

Conclusions
The three-dimensional, non-stationary heat-transfer model was developed for the direct deposition processes of thin-walled parts with various configurations. The calculated transient temperature fields obtained using the developed model agreed satisfactorily with the experimental results, both in the heating phase and in the free-cooling phase. At the same time, the presented calculations were performed with cardinally different input parameters of the model; namely, the types of heat sources, the process technological parameters, thermophysical properties of materials, movement trajectories, and pause times were different. These facts suggest that the model allowed reproducing threedimensional temperature fields in a wide range of technological process parameters and for various product configurations.
To develop the model, a three-dimensional analytical solution of the non-stationary heat-conduction equation for a moving distributed heat source was obtained, taking into account heat transfer to the environment. The model made it possible to calculate all the temperature field characteristics; in particular, thermal cycles, temperature gradients, and the cooling rate in a path during the deposition process at any time. In this case, the wall and the substrate sizes, the change in power from layer to layer, the pause time between passes, and the trajectory of the heat source were taken into account.
Thus, the developed model can be used to predict temperature fields, as well as to study the degree of influence of local quasiperiodic temperature fields, which always accompany direct deposition processes, on structural phase transformations.
The proposed calculation scheme applies to thin-walled products; however, it can be extended for the case of multi-pass thick walls. In addition, this model can be easily adapted to various direct deposition technologies that use a laser, electron beam, plasma, or arc as the main heat source.