A Semi-Analytical Solution for Temperature Distribution in Friction Stir Welding


 In this paper, a semi-analytical thermal model of friction stir welding processes is developed in which the heat-generating regions are divided into many elements as point heat sources. The heat generation in each element involves the friction and plastic deformation, and the temperature rise caused by each element is calculated by solving the heat conduction equation of a moving heat source in a solid body. The heat loss through the top and bottom surfaces are considered in the model as heat sinks. The asymmetric distribution of the temperature is calculated through the whole process and over the whole volume of the workpiece by integrating the effects of all heat sources and sinks. The temperature-dependent material properties are updated by a numerical routine. The comparison between the calculated results and the experimental data clearly approved the validity of the proposed method for some aluminum and steel alloys.


Introduction
One of the solid-state joining processes for producing joints with acceptable properties in aluminum alloys and steels is the friction stir welding (FSW), which was developed at The Welding Institute (TWI), UK, in 1991 [1]. The FSW process has lower heat inputs than the conventional fusion welding reducing metallurgical changes in the heat affected zone (HAZ), does not produce gases or fumes, and can join materials that are difficult to fusion weld.
Currently, this welding technique is used in various industries, such as aviation industry, automotive, railway manufacturing, and shipbuilding [2].
In the FSW process, as shown in Fig. 1(a), a non-consumable cylindrical tool with a profiled probe (pin) rotates and slowly plunges between two clamped plates until the tool shoulder touches the surface of the plates. The localized heating due to friction between the rotating shouldered tool and workpiece, and plastic deformation of the material softens the materials near the tool surfaces without reaching the melting point (a solid-state process). As the tool traverses along the weld line, the plasticized material transfers from the front to the back of the tool pin. Since the local relative motion of the tool and the workpiece is a combination of the translational and rotational motions, the heat generation rate varies from one side to the other side of the workpiece, which are referred to as the advancing and retreating sides ( Fig.1(b)). The temperature distribution directly affects the mechanical properties, the residual stress, the distortions, and the microstructure of the joints; therefore, it is important to predict the distribution of temperature during the FSW process. Up to now, numerous models have been developed to simulate the heat transfer and materials flow during the FSW process. The nonlinear thermal equation can be solved either analytically or numerically for finding temperature variations during the FSW process. Vilaca et al. [3] developed an analytical thermal model to simulate the asymmetric temperature field The previous semi-analytical methods had some limitations, such as not considering all the components of heat generation, the temperature-dependent material properties, or the heat loss due to convection and radiation. In addition, the numerical methods are time-consuming, or are not applicable for later mechanical analysis; therefore, in the present paper, a semianalytical thermal model is developed which considers all components of heat generation and heat loss due to convection and radiation, and can quickly and accurately predict the temperature distribution in the FSW process.

Governing equation
The governing heat conduction equation for the transient heat transfer analysis during the FSW process is given as follows: where ρ, Cp , and k are the density, the specific heat capacity and the thermal conductivity of the material, respectively; the term P Q is the volumetric heat generation caused by plastic dissipation in the FSW process. By removing P Q from Eq. (1) and considering it as a surface flux boundary condition at the tool/workpiece interface, Eq. (1) can be simplified and then expanded as follows: (2) can be reduced to a linear differential equation as follows: where is the thermal diffusivity of the workpiece. In the theory developed by Rosenthal [18], a constant heat source is located at the origin of a rectangular coordinate  system moving at a constant speed with respect to the stationary rectangular coordinate system as shown in Fig. 3. Using this assumption and the relations shown in Fig. 3, Eq. (3) can be simplified as follows: The solution of Eq. (4) must satisfy the following conditions: 1) Since the source is a point source, the heat flux around the source must tend to the value of the total heat generation caused by both frictional and plastic deformation delivered to the workpiece, 2) There is no heat dissipation from the workpiece to the surrounding atmosphere, 3) The temperature of the workpiece remains unchanged at very great distance from the source. In addition, Eq. (4) is more easily handled by replacing it with the following expression derived from Appendix A: Then, the final heat transfer governing equation becomes: Referring to Appendix B, the solution of Eq. (6) is provided for three-dimensional heat flows or thick plates in terms of the radial distance, R, from the heat source as follows: where Q is the total heat generation caused by both frictional and plastic deformation and also the plastic dissipation at a point source. For two-dimensional heat flows or thin plates, as described in Appendix C, the better results can be obtained by using the below equation: where H and K0 are the thickness of workpiece and the modified Bessel function of the second kind and zero order.

Heat generation
For modeling the heat generation in the FSW process, three contact states namely sticking, sliding, and partial sliding/sticking are considered. In the sticking or sliding condition, the material surface sticks, or slides to the moving tool surface, and the heat generation is produced by the plastic deformation or friction, respectively. In the mixed state, the heat is generated by both friction and plastic deformation in the shear layer. Using the definition of the tool rotational speed relative to surrounding materials (Fig. 4) and the mixed heat generation definition from the previous papers [14,19,20], the heat generation is formulated as follows:  where τ, P, ω, r, η, and µ are the shear yield stress of material, the normal pressure of the tool, the rotational speed of the tool, the radial distance of a point from the axis of the tool, the efficiency of process and the coefficient of friction, respectively; In Eq. (9), TP indicates the tool properties, and WP indicates the workpiece properties. In this paper, the shear yield stress (obtained by relation In addition, Ψ is the advancing retreating factor introduced by Darvazi et al. [24] and is about 20-30 in the present model, and δ is the sticking state variable and is defined as the velocity of the contact points at the workpiece surface relative to the tool point in contact (sticking, δ = 1; sliding, δ = 0; and sticking/sliding, 0<δ<1) [16]. Finally, by substituting the total heat generation in Eqs. (7) and (8), in the case of 2D and 3D heat flow, and using the superposition principle, the temperature of the FSW process under the steady-state condition can be calculated.

Temperature rise in workpiece caused by the tool shoulder and the pin tip
The surface of the tool shoulder is divided along the radial and circumferential direction, as indicated in Fig where dQ is heat generation for an element substituting instead of Q in Eqs. (7) and (8), and q can be obtained from Eq. (9). As shown in Fig. 6, for an arbitrary point C, the heat generated at point A plays the role of heat source; therefore, the temperature generated by each element is described as follows: where R and r' are the 3D and 2D distance from element center. The details of the calculation of these variables are given in Fig. 6. 6. The distance of one arbitrary point C from an element center.
Using the superposition principle and neglecting the nonlinearity effect of temperature in Eq.
(2) and considering the 3D based model, the temperature induced by the tool shoulder in the steady-state condition can be calculated as follows: where m and n are the number of divisions along the radial and circumferential directions, respectively. The temperature rise caused by the pin tip, pin tip T  , is calculated like the temperature rise caused by the shoulder.

Temperature rise in workpiece caused by the pin side
As shown in Fig. 7, the temperature generated by the pin side is considered as 2D heat flow, and the total heat generation is calculated as follows: 13 where r p is the radius of the pin of the tool. By substituting Eq. (14) in Eq. (12) for an element and integrating the effects of all elements, the temperature induced by the pin side in the steady-state condition can be calculated as follows: Definition of an element on the pin side.

Temperature rise in workpiece due to plastic dissipation
The heat generated due to the plastic deformation dissipation away from the interfaces in the thermo-mechanically affected zone (TMAZ) is about 4.4% (0.044) of the total heat generation [14,20,21]. Avila (2006) [22] used an inverse trigonometric function, in order to better approximate the domain of TMAZ as illustrated in Fig. 8. The inverse tangent range is arbitrarily restricted to (-π/2, π/2), and three constants are required to fit the function. The function used in this case is defined as follows: where all values of C1, C2, and C3 are constants depending on the tool dimensions in millimeters(mm). In the present study, the radius of the TMAZ is defined as 0 sin  TMAZ R r v  which is wider at the advancing side than the retreating side as shown in Figs. 8 and 9. The heat generation in the TMAZ is estimated using Eq. (9) and neglecting the term of frictional heating; therefore, the heat generation due to plastic dissipation can be expressed as: where S is a line element on the TMAZ curve shown in Fig. 9 and can be approximated as follows: Then, the temperature rise due to plastic dissipation in the TMAZ region in an arbitrary point is calculated as follows:

Temperature change in workpiece due to convection from surfaces
Since the previous analytical models had not considered the heat loss due to convection and radiation, the estimated temperature almost would be higher than the actual value. For solving this problem, two methods can be used. First, the temperature distribution can be multiplied by a factor to correlate the temperature results from the semi-analytical method and experimental results. The second method considers the heat loss due to convection and radiation as a heat sink and is defined as Eq. (20).
As shown in Fig. 10, the workpiece is divided along the x and y directions, and the distance between the heat source (the point B) and an arbitrary point (point A), R, is calculated; then the temperature drop of point A due to convection and radiation is calculated using Eq. (9) as follows:  where T0 is the ambient or initial temperature (22-25•C); h is the convection coefficient; ε is the surface radiation emissivity; and β is the Stefan-Boltzmann constant. Furthermore, the surface of the workpiece in contact with the backing plate is modeled as a convection condition with a higher convection coefficient, and the temperature drop due to it, bottom con T   , is calculated using Eqs. (20) and (21) without considering the radiation condition, and Z is replaced by H-Z. The implementation of Eq. (21) needs to use a routine for reaching a steadystate condition.

Finding total temperature distribution
In Eq.(2) the terms of nonlinearity are removed, and it is reduced to a linear differential equation. According to the superposition principle, which is correct for linear equations, the resultant effect is the sum of the effects produced by each component acting separately; therefore, the total temperature distribution can be calculated by Eq. (22). The temperature distributions for aluminum alloy 5083-H118 are shown in 3D views in Fig. 11 T  T  T  T  T  T  T  The maximum temperature occurs at the location with w equal to zero; in this point, the temperature distribution function is singular. In this paper, two methods are used for handling this singularity. In cases that the maximum temperature is unknown, only the part of the model is simulated in which R is not zero. In cases that the maximum temperature is known, R+ε is used instead of R, and the following relation is applied for finding the appropriate value of ε:

Materials and experimental procedure
In the present paper, the temperature-dependent thermal and mechanical properties of the material such as specific heat, thermal conductivity, density and yield shear stress are taken from the previous literature [13][14][15][16][17]23]. A numerical routine is constructed to determine the temperature-dependent thermal properties during the analysis. At the beginning of the analysis, the material properties at the ambient temperature are implemented to the model, and the temperature is calculated by Eq. (22); then using the previous temperature distribution, the material properties are updated and applied to the model estimating the new temperature distribution. This process is repeated sequentially until a steady-state temperature distribution is achieved. The schematic of temperature updating is shown in Fig.   12.  base material, supplied in 3 mm thick plates welded by the FSW process, the experimental process is installed. Friction stir welding is carried out on a CNC vertical milling machine. In the experiments, the thermal histories on surfaces are measured using K-type thermocouples located at different distances from the weld line, as shown in Fig. 13. The data used for the calculations such as tool geometry and welding conditions are presented in Table I.

Friction stir welding of 304L stainless steel
The temperature distribution of 304L stainless steel is calculated and compared with the experimental results of Zhu and Chao [13]. To maintain consistency, the material properties,

Friction stir welding of Mild steel 1018
The calculated temperature history on the advancing side at a distance of 12.7 mm from the weld line is compared with the experimental data [24] in Fig. 16. A good agreement between the calculated and measured results indicates that the developed model can be used for the prediction of temperature profiles and cooling rates.

Friction stir welding of aluminum alloy 5083-H118
Another verification of the proposed model is carried out by comparison with the experimental measurements. The welding parameters associated with the three different FSW conditions are presented in Table I. The thermocouples are inserted into holes with a depth of 1 mm from the top surface at a distance of 15, 21 and 32 from the welding centerline.
The comparison of the calculated temperature histories and the experimental data shows good agreement between them (Figs. 17, 18, and 19).

Friction stir welding of aluminum alloy 6061-T6
The calculated temperature histories at two locations are compared with the experimental data from the literature [25] in Fig. 20. Good agreements between the experimentally determined and the numerically calculated results at two locations indicate that the model can be used to examine the temperature profiles and cooling rates.

Aluminum alloy 7075-T7351
In order to verify the model, the calculated temperature variations are compared with the experimental data from Colegrove and Shercliff [23] that were obtained using thermocouples at various locations, as shown in Fig. 21. The calculated results match very closely with the experimental data.

Discussion
The temperature distribution has been calculated by a semi-analytical method, as shown in Fig. 22.
The heat-generating regions are divided into many elements that each element is considered as a heat source, and for an arbitrary point C, the heat generated at point An plays the role of a heat source (Fig.   22). Finally, the temperature distribution is calculated by integrating the effects of each heat source.
(a) (b) Fig.22. Schematic of any element of the shoulder region as a heat source.
As proposed by previous researchers [4] the heat source can be considered to be concentrated at the mid thickness of the workpiece in the center of the nugget and traveling with a constant linear velocity (v) (Fig. 23(a)), whereas in the presented model as shown in Fig. 23(b), the heat sources are distributed in all regions where the heat is generated. They also used the constant thermal properties at the ambient temperature, and did not consider the heat radiation or convection and backing plate in temperature calculation.  Fig. 24(a), the temperature distribution was calculated by a commercial CFD code, FLUENT, based on a 3D viscoplastic model [26], and in Figs. 24(b) and (c) the temperature was calculated by FE method based on solid mechanics heat transfer model [20,27]. As shown in Fig. 24(d), the present model can also predict the temperature distribution like the previous methods; however, it shows that the temperature distribution is asymmetric around the weld line, and the maximum temperatures are found in the back half of the shoulder/workpiece interface and toward advancing side. In comparison with the CFD and FE methods, the present method is time effective, and converges quickly, so it can be useful for the parametric study of the FSW process. Although the temperature-dependent thermal conductivity and diffusivity are implemented during the analysis, they are considered constant with temperature for each element (point heat source) in Eq. (2) and (3). This limitation causes some differences between the calculated results and the experimental data in the low-temperature zone, and these differences decrease in the high-temperature zone in which the metals exhibit an almost constant thermal conductivity and diffusivity behaviors (Figs. 13 through 20). Furthermore, this method can be used for additional studies about friction stir welding, which includes thermal cycles or thermal distribution and can be used for mechanical analysis such as thermal stress during the welding process.  This method predicts the asymmetric temperature at the advancing and retreating sides.

Conclusions
 Since this method considers the most parameters of the FSW process, it is suitable for studying the parameters affecting the temperature distribution.
 The present method is fast and accurate enough for predicting temperature distribution and is applicable for later mechanical analysis.

APPENDIX A: Finding the form of solution
By using ∅ = , Eq. (B.1) can be simplified as below: The solution of Eq. (B.2) is: By implementing the condition: I) The temperature of the workpiece remains unchanged at infinite distance from the source, i.e.: Then Eq. (5) becomes: II) The heat flux through the surface of the hemisphere drawn around the source must tend to the value of the total heat, Q, delivered to the workpiece, as the radius of sphere tends to zero: In addition, when the heat source is volumetric, such as that in the TMAZ region, the heat flux diffuses radially in a spherical shape (not hemisphere), and the Eq. (B.6) should be rewritten as follows:

APPENDIX C: Two-dimensional solution
Eq. (6) in cylindrical coordinate in the quasi-stationary state can be expressed as follows:

Declarations
Funding: Not applicable

Conflicts of interest/Competing interests:
The authors disclose that they do not have any competing interests.

Availability of data and material: Not applicable
Code availability: Not applicable

Authors' contributions: Not applicable
Ethics approval: This manuscript has not been simultaneously submitted to another journal for review. The work submitted is original and has never been published in any form or language before (partially or in full).                        [20] and b) the present method. Tables   Table I. The data used in the simulation of the FSW of aluminum alloys and steels.  Thermal conductivity and its average slop versus temperature [13][14][15][16][17].   The velocity of one point on the tool shoulder. De nition of an element on the surface of the tool shoulder.

Figure 6
The distance of one arbitrary point C from an element center.

Figure 7
De nition of an element on the pin side.

Figure 8
The suggested TMAZ pro le for 304L stainless steel at advancing and retreating side.

Figure 9
The shape of the TMAZ and de nition of its boundary.

Figure 10
Modeling of heat loss due to convection and radiation as a heat sink at the top surface Figure 11 The distribution of temperature for AA5083-H118 calculated via the present semi-analytical model induced by (a) the shoulder, (b) the pin, (c) the plastic dissipation in TMAZ and (d) total temperature for the welding condition of 800 rpm and 30 mm/min.

Figure 12
The algorithm of the numerical code for material updating written in Fortran language (in the present paper δT=1°C).

Figure 13
Experimental installation (a) workpiece and thermocouple installation and (b) the Tool geometry.

Figure 14
Comparison between the experimental data [13] and calculated temperature histories on the top of the retreating side of the workpiece at various distance from the weld line for the translational speed of 102 mm/min and the rotational speed of 300 rpm.

Figure 15
Comparison between the experimental data [13]and calculated temperature histories on the top of the retreating side of the workpiece at various distance from the weld line for the translational speed of 102 mm/min and the rotational speed of 500 rpm.

Figure 16
Comparison between the experimental [24] and calculated temperature histories on the advancing side at a point 12.7 mm from the weld line for the translational speed of 25 mm/min and the rotational speed of 450 rpm.

Figure 17
The temperature histories at various distance from the weld line for the rotational speed of 800 rpm and the translational speed of 100 mm/min.

Figure 18
The temperature histories at various distance from the weld line for the rotational speed of 1100 rpm and the translational speed of 100 mm/min.

Figure 19
The temperature histories at various distance from the weld line for the rotational speed of 1400 rpm and the translational speed of 100 mm/min.

Figure 20
Comparison between the experimental [25] and calculated temperature histories at points 2 mm below the top surface at a distance of ( a) 8 mm and ( b) 16 mm from the centerline for the rotational speed of 1400 rpm and the translational speed of 100 mm/min.

Figure 21
Comparison between the experimental [23] and calculated temperature histories at different locations for welding conditions of (a) 220 rpm and 220 mm/min, and ( b) 457 rpm and 457 mm/min.

Figure 22
Schematic of any element of the shoulder region as a heat source.

Figure 23
Modeling of heat source: (a) the Ferro and Bonollo model [4] and (b) the present model.

Figure 24
The Computed temperature pro les on the top surface of workpiece for the rotational speed of 300 rpm and the translational speed of 102 mm/min: (a) Nandan et al. [26], (b) Darvazi et al. [20], (c) He et al. [27] and (d) The present work.

Figure 25
The cross-sectional temperature distribution of the FSW in the transverse direction computed by: a) Darvazi et al. [20] and b) the present method.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. GraphicalAbstract.docx