Efficient analysis of welding thermal conduction using the Newton–Raphson method, implicit method, and their combination

Finite element analysis is commonly used to investigate the thermal-mechanical phenomena during welding. To improve the computing efficiency of finite element analysis for welding thermal conduction, a novel Newton–Raphson method (NRM) without the computation of inverse matrix and a hybrid method combing the NRM and conventional implicit method (IMP) were developed. Comparison of computing time between the hybrid method implemented in an in-house software JWRIAN and the IMP used in a commercial software ABAQUS indicated that the computing speed of the former was about 4.5 times faster than that of the latter. Additionally, compared to the conventional IMP, the NRM exhibited higher computing efficiency in the analysis of transient thermal conduction during the welding heating process. Meanwhile, a combined hybrid method of the NRM and IMP was verified to be more efficient in analyzing the welding thermal conduction throughout the heating and cooling processes. Moreover, the thermal cycles computed by the hybrid method were consistent with those from experimental measurement, indicating the high accuracy of the hybrid method. Furthermore, the hybrid method was used to predict the temperature field of the corner boxing fillet joint welded by a low transformation temperature weld metal for generation of compressive residual stress.


Introduction
Welding is one of the most practical techniques for assembling structures [1], such as ships, automobiles, and electronic goods [2]. For the effective design and analysis of the welding process, it is of great significance to accurately calculate the transient temperature field because the welding thermal cycles play a pivotal role in the prediction of the microstructure, residual stress, distortion, and their characteristics of welded joints [3][4][5][6].
To analyze the temperature field of a welded joint, most researchers have focused on the heat transfer during the welding process. Rosenthal [7] and Rykalin [8] were the first to use the principles of Fourier's theory of heat flow to understand the welding heat transfer. However, the theoretical solution of the temperature field is based on the following inaccurate assumptions: the material properties are independent of temperature; there is no heat exchange between welded components and the surroundings; the heat sources are concentrated at a point, line, or plane. To rectify the above shortcomings, researchers modified the heat source models and solutions. Eagar and Tsai [9] modified the theoretical solution and proposed a Gaussian-distributed two-dimensional-surface heat source to predict the temperature of the semi-infinite weldment. Jeong and Cho [10] successfully simulated the temperature field of fillet-welded joints by considering a bivariate Gaussian source. Additionally, Goldak et al. [11] developed a three-dimensional double-ellipsoidal moving heat source model. Based on the Goldak heat source model, Sabapathy et al. [12] considered the oscillation of the burner for a better modification of the weld pool geometry. There are some differences among welding processes in terms of physical characteristics and temperature distribution of the heat input. As a result, various heat source models were established utilizing different formulations [13][14][15][16][17], which showed a good agreement with experimental data in temperature fields of welded joints.
The temperature distribution during welding can be obtained according to the partial differential equation of transient heat conduction, which is usually written as follows: where ρ, c, T, t, λ and Q represent density (g/mm 3 ), specific heat (J/g°C), temperature (°C), time (s), thermal conductivity (J/mm s°C), and internal heat generation rate (W/mm 3 ), respectively. The x, y, and z coordinates represent the longitudinal, transverse, and thickness directions, respectively. Equation (1) can be solved analytically or numerically. The analytical method can only be applied to infinitely large geometries with simple boundary and initial conditions [18], while the numerical method is suitable for various boundary conditions and can also deal with nonlinear problems. Hence, the numerical method is the most commonly used in both academic and engineering applications. The main numerical method employed for welding simulations is the finite element method. With the development of simulation technology, several methods have been developed to calculate the thermal and mechanical behavior of welded joints, including the finite difference [19], finite volume [20], meshless [21], and boundary element [22] methods.
The transient heat conduction has been determined to be a time-dependent process [23], and considerable computational power is required to solve the ordinary differential equations in the time domain. To improve the computing efficiency of transient heat conduction simulations, a novel Newton-Raphson method (NRM) and its combination with the conventional implicit method (IMP) are proposed in this study to simulate temperature field. The accuracy and efficiency in computing temperature fields of a singlepass butt-welded pipe using the proposed methods are compared to those obtained with the conventional implicit method (IMP) in the commercial software ABAQUS. In addition, the computed temperature fields are compared to those obtained from experimental measurement. Furthermore, the proposed method is applied to predicting the thermal cycles of elongated weld in corner boxing fillet-welded joint.
2 Galerkin finite element method for three-dimensional temperature fields The finite element formulations of three-dimensional temperature fields can be obtained using the Galerkin's weighted residual method. The derivation of the finite element formulations has been demonstrated in several studies, so it is only briefly discussed in this paper.
The boundary conditions associated with Eq. (1) are listed as follows where k is tensor matrix of the thermal conductivity; ∇ is the gradient operator, ∇ ¼ ∂ ∂x i þ ∂ ∂y j þ ∂ ∂z k. n is normal unit vector, n ¼ n x i þ n y j þ n z k. S 1 , S 2 , and S 3 denote as the Dirichlet, Neuman, and Robin boundary conditions, respectively. The computational domain is V = S 1 ∪ S 2 ∪ S 3 , T , and q represent the given temperature and heat flux density, h is the convection heat transfer coefficient, n denotes the unit normal vector for the surface, and T a is the environment temperature.
Assuming that there are a total of m elements and n nodes owing to the discretization of the whole solution domain, the set of algebraic equations for Eq. (1) can be expressed as follows: ð5:2Þ ð5:1:1bÞ In the three-dimensional temperature field, the variable temperature T (x, y, z, t) is a function of space and time. Accordingly, the spatially discretized formulation with time differential elements ∂T/∂t (Eq. (5)) can be discretized suing the difference method as follows: where Δt is the time step, and θ ∈ [0, 1] is the weighting factor; with θ = 0, 0.5, 2/3, and 1, Eq. (6) corresponds to the forward Euler difference, Crank-Nicolson, Galerkin, and backward Euler methods, respectively [24]. In this study, we use the Galerkin difference method. By substituting Eq. (6) into Eq. (5), we obtain the following equations: Given the time level t + Δt, the temperature {T} t + Λt at all nodes can be computed as follows: With the above equation, corresponding to the conventional IMP, the three-dimensional temperature field can be easily obtained, but the computation can be time-consuming owing to the computation of the inverse matrix [A] −1 . Furthermore, if many degrees of freedom need to be analyzed, e.g., in the order of millions or more, a large amount of memory is also required, making a personal computer insufficient for many simulations.

Newton-Raphson method for computing the welding thermal conduction
To efficiently compute the transient temperature field due to the thermal conduction of the welded joint, a novel solving method based on the Newton-Raphson scheme is proposed.
The computed temperature T f g cal tþΛt is an approximate solution, and the residual error can be written as follows: Therefore, reducing the residual error results in a large difference in the simulation results. Here, a function related to the residual error is proposed as follows: where Q i is an element of the heat flow matrix {Q} n × 1 , A i, k is an element of the ith row vector of the conductivity matrix [A] n × n , T k is the nodal temperature, and NF is the number of nodes with unknown temperature. The purpose of our analysis is to make each residual error term B i ({T}) in Eq. (10) converge to zero, i.e., to solve the nonlinear algebraic equation B i ({T}) = 0. To this end, we differentiate Eq. (10) and obtain Subsequently, the modified NRM is applied to the nonlinear equation to seek the most reliable solution; each nodal temperature can be expressed as follows: where μ is the correction factor. Only the diagonal components a ii of the matrix [A] are employed to compute the temperature field. Therefore, there is no need to save the full components of matrix [A] or to solve the inverse matrix [A] −1 , making the proposed method more efficient in terms of computing memory and power consumption. Figure 1 shows flow charts of the IMP and NRM with the following parameters: where {Q 1 } is the heat flux vector at first iteration; {T} N is the temperature matrix at the Nth iteration; {ΔQ N } is the residual vector at the Nth iteration; ‖{Q 1 }‖ is the norm of the heat flux vector, whose value ; ‖{ΔQ N }‖ is the norm of the residual vector at the Nth iteration, whose value is tolaa is the computational tolerance; and NJ is the number of nodes with unknown temperature. The hybrid method is the combination of the IMP and NRM, as discussed later. Generally, for a higher computing efficiency, the NRM with small time steps is used for the welding heating process, while the IMP with large time steps is used for the cooling process.

Comparison of accuracy and computing time with commercial software
To highlight the advantages of the newly proposed hybrid method, accuracy and computing time between in-house software JWRIAN and commercial software ABAQUS were compared firstly. Here, a pipe-welded joint was modeled by a total of 432,000 solid elements and 577,152 nodes, as shown in Fig. 2. The time step used in both heating process during welding and cooling process after welding was 0.1 s. Both computations were performed in the same computer m-Book K690XN-M2SH2   Figure 3 shows the temperature-dependent The density is assumed to be 7850 kg/m 3 independently of the temperature. In addition, the room temperature and heat transfer coefficient are set to 20°C and 10 W/m 2 K, respectively. Figure 4 shows the computing times of the proposed hybrid method in the software JWRIAN and IMP method in ABAQUS. Evidently, there exists a large difference in computing time between these two methods. The total computing time of the ABAQUS is about 4.5 times longer than that of the JWRIAN. Specifically, the main difference between these two software is the computing time needed for heating process. The hybrid method using the NRM method in heating process accelerates the computation to a great extent, fully highlighting the advantage of the newly developed method. Additionally, Fig. 5 shows the thermal cycles of specified nodes computed with JWRIAN and ABAQUS. It can be seen that the welding thermal cycles of the three specified nodes computed with JWRIAN and ABAQUS are the same. Figure 6 a and b show the temperature distribution comparison using JWRIAN and ABAQUS, respectively. The overall distributions are the same. The maximum temperatures are 1835 and 1838°C, respectively, indicating little difference between the two software and confirming the high accuracy of JWRIAN.

Comparison of computing time and accuracy of different methods and models
Today, the IMP (Eq. (5)) is the most common method in computing welding thermal conduction. However, the hybrid method as revealed in the above section is superior to the IMP method in terms of computing efficiency. Here, the relationship between unknown freedoms of models (numbers of nodes) and computing time was investigated for quantitative verification of the computing efficiency of the NRM, IMP, and their hybrid methods. The butt-welded joint of the pipe shown in Fig. 2 was used again for computation using the in-house software JWRIAN. In addition, six types of simulation models with different numbers of nodes (type Hexa8 element) are presented in Table 1. The time step used in heating process was 0.1 s and the time step used in the cooling process was 0.5, 0.1, and 0.5 s, respectively for the IMP, NRM, and hybrid methods. The computations were performed in the same computer Dell Precision 5820 (Xeon W-2145/32GB memory) with 4 cores. The welding parameters and material properties for the butt-welded joint of the pipe were set as the same as that defined in the above section. Figure 7 a shows the simulation results. From the figure, it can be seen that the computing time of both methods Generally, the IMP is suitable for a time-consuming simulation of a physical phenomenon because a large time increment is acceptable [25]. In the NRM, the initial iteration may deviate from the solution with a large time increment. Therefore, small-time increments are often adopted to ensure the robustness of the NRM [26]. Interestingly, as shown in Fig. 7b, the simulation results of model-6 confirm the advantages and disadvantages of the IMP and NRM. Hence, the hybrid method was proposed to combine the advantages of the short-time heating process of the NRM and the long-time cooling process with the IMP to analyze the heat conduction in welded joints. The results of the hybrid method are also shown in Fig. 7, which confirms that the hybrid method can accelerate the computation by a considerable amount.
To verify the accuracy of the proposed methods, the thermal cycles of five specified nodes of model-6 were computed with the different computational methods, and the results were compared. As shown in Fig. 8, the thermal cycles of the three specified nodes computed with the three computational methods during heating process are very similar. It is worth noting that the thermal cycles computed with the hybrid method closely follow the trend of those computed with the IMP and NRM. Figure 9 compares the maximum temperature of the model-6 obtained with the three methods. The maximum temperature obtained using the IMP, NRM, and hybrid method is 1822, 1835, and 1835°C, respectively, indicating little difference among the three methods and confirming their high accuracy.

Multipass butt welding
The microstructural evolution, mechanical properties, and residual stress in the welds and heat-affected zones of multipass-welded joint are strongly affected by welding thermal cycle and material properties of base metal. Therefore, it is meaningful and necessary to accurately compute the welding thermal cycles in welded joints. Figure 10 displays a multipass butt weld and its three-dimensional finite element model [27] for welding thermal conduction analysis using the hybrid method. In addition, measurements of the temperature field in the welded joint were experimentally conducted to validate the computed results. The welding parameters for the butt-welded joint were set as follows: a current of 220 A, voltage of 26 V, welding speed of 3.3 mm/ s, and thermal efficiency of 0.85. The room temperature in performing this experiment is 12°C. Figure 11 shows the temperature-dependent thermal properties of the SM490 steel and heat transfer coefficient. The density is assumed to be 7850 kg/m 3 at all temperatures, and the melting point of the welded components is set to 1450°C.
The thermal cycles during the first pass welding at two positions were measured by a thermocouple, as shown in Fig. 12a. Figure 12b shows the simulated and measured thermal cycles at the two positions during the first pass welding. It can be seen that the computed and measured thermal cycles agree well, demonstrating the high accuracy of the results of the hybrid method.
As shown in Fig. 13, the shape of the computed molten zone and the experimental weld cross sections are consistent. This indicates that the hybrid method is efficient and accurate for welding thermal conduction analysis.

Thermal cycles in elongated weld with low transformation temperature welding material
It is well-known that tensile welding residual stress is detrimental to the fatigue strength of welded joint [28,29]. To minimize the unfavorable tensile stress, the low transformation temperature (LTT) alloys, which take advantage of martensitic transformation, have been developed during the past two decades [30][31][32]. Furthermore, it has been investigated that an elongated weld method overlaying LTT weld metal on the corner boxing fillet-welded joints can greatly enhance fatigue lives to at least 4 times [33]. In practice, both the martensite start (Ms) temperature of LTT alloys and interpass temperature during welding process strongly affect the residual stress development in multipass-welded joints, which has been verified by a number of studies [34][35][36]. Generally speaking, much beneficial compressive residual stress will be left in the weld zone if its temperature during welding can be controlled above the Ms temperature. Since welding processes can greatly affect the temperature distribution and welding thermal cycles in both the molten zone and heataffected zone, it is necessary to predict the thermal cycles of the LTT bead before practical welding. Figure 14 shows the half model of corner boxing fillet-welded joint, where the mesh size around weld zone is 2 mm. The welding parameters for the LTT welds were set as follows: a current of 200 A, voltage of 32 V, swing speed of 12 mm/s, and thermal efficiency of 0.85. The total welding time is 33 s. The density is assumed to be 7500 kg/m3 at all temperatures, while the other thermal properties like thermal conductivity and specific heat are shows in Fig. 15. The root temperature is set to 25°C, while the heat transfer coefficient is the same as that shown in Fig. 11b. In addition, the melting point and Ms temperature of LTT weld metal are 1450°C and 200°C, respectively. Figure 16 shows the temperature distributions of welded joints with different welding directions before cooling process. It can be seen in Fig. 16a that the temperature of the whole LTT weld exceeds the Ms temperature (200°C) in the downward welding. Conversely, as displayed in Fig. 16b, the temperature at the front LTT weld toe is lower than the Ms temperature in the upward welding. Furthermore, the thermal cycles of LTT passes shown in Fig. 17 also demonstrate that heat dissipation is faster in the welded joint with upward welding compared to that using downward welding. This suggests that downward welding is more suitable for the fabrication of LTT weld using the same welding parameters. On the other hand, large heat input should be used in the upward welding.

Conclusions
To efficiently analyze welding thermal conduction, a novel NRM and a hybrid method combining the NRM and the conventional IMP were developed. Additionally, the computing time and accuracy obtained with the hybrid method were compared with those computed with the IMP using the commercial software ABAQUS. Furthermore, the newly developed hybrid method was applied to the analysis of welding thermal cycles in practical welded joints. Based on the analyses, the following conclusions can be drawn: 1. The computing speed of the hybrid method was about 4.5 times faster than that computed with the IMP using the commercial software ABAQUS. 2. Compared with the IMP, the NRM can significantly accelerate the computing speed while ensuring high accuracy; this is because the NRM does not require the calculation of the inverse matrix. 3. By combining the advantages of the IMP and NRM, the hybrid method can further shorten the computing time while ensuring high accuracy. 4. The thermal cycles computed by the hybrid method were in good agreement with the experimentally measured ones, indicating the high accuracy of the hybrid method. 5. The thermal cycles of elongated LTT weld computed by the hybrid method indicated that the temperature in all weld zone due to downward welding was higher than  the Ms temperature, which was a necessary condition for generation of compressive residual stress.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.