A semi-analytical method for the calculation of double-ellipsoidal heat source parameters in welding simulation

Structural integrity assessment of welded components requires the determination of residual stress and deformations induced by the welding process. The implementation of both a numerical static mechanical analysis and a thermal transient analysis is needed in order to accurately determine stress and strain fields. The most faithful representation of the heat flow provided by the welding torch in thermal simulations is realized with Goldak’s double ellipsoid, especially in case of MIG/MAG and combined TIG-MIG/MAG welding. The correct determination of the temperature trends during the welding process simulation is essential for the subsequent mechanical FEA. In order to obtain the above-mentioned temperature trends, some parameters must be determined ex ante, e.g., the geometrical characteristics of Goldak’s double ellipsoid. Generally, a large number of trial simulations must be performed in order to define these parameters. However, manual fitting is often complicated. To overcome such difficulties, in this work, it is reported a combined analytical and numerical method to determine the parameters using the genetic algorithm NSGA-II. First-attempt parameters are obtained through the solution of the Fachinotti’s analytical formulation. This first stage requires a restricted set of experimental data, such as weld pool extension and one temperature trend in the proximity of the seam. Afterwards, first-attempt parameters and experimental data are used to determine the final parameters by means of a set of numerical thermal simulations guided through the genetic algorithm optimization. To confirm the good-functioning of the method, a numerical analysis for the simulation of the experimental data and another FE simulation performed with the final results of the method are presented.


Introduction
Following the welding process, residual stresses and distorsions generally occur in the joint. Residual stress can be added to those acting on the joint during the normal operation, resulting in unforeseen stress conditions. In addition, in structures fatigue cracks often occur in welded joints. The correct determination of residual stresses is therefore essential for static and fatigue assessments. The implementation of advanced structural simulations of welded joints requires adequate thermal simulations. These must be able to faithfully replicate the temperature trends that occur in the welding process. Rate temperature changes have a significant influence on the metallurgical phase transformations that have a pivotal role in thermal dilatation and in material IOP Publishing doi: 10.1088/1757-899X/1214/1/012023 2 behaviour, especially for ferritic steels [1,2,3]. Goldak's double ellipsoid [4] is the state-of-art model for thermal simulations of a moving source of heat flux. The TIG and MIG/MAG welding can be simulated with outstanding results using the double ellipsoid modelling. This model has been used for the simulation of many types of welded joints, such as butt-welded plates, girth joints, fillet and cruciform welds.
The calibration of Goldak's model is achieved by varying the parameters of the double ellipsoid in the transient FE simulations. Correct calibration is obtained by replicating the the weld pool extension of the joint and the temperature trends provided by the thermocouples in the experimental tests. This procedure is rarely straightforward and is generally carried out manually by trial and error. In the case of multi-pass welds, it is difficult to determine the parameters. Welding parameters such as voltage, current and welding speed can vary with each pass, requiring as many calibrations as passes [5]. Manual calibration is challenging even in the case of special weld pool shapes, that can be achieved by combining several heat flux models. This work has been carried out in order to facilitate the implementation of numerical thermal simulations by eliminating the manual calibration of the heat flux. The method outlined below employs the genetic algorithm to guide the process of Goldak's parameter determination through 2D numerical models. When parametrically configured these models allow many automated simulations to be performed with a low computational burden. These simulations can be used with negligible errors when conduction in the longitudinal direction can be neglected. Such modelling may lead to increasing errors and deviations from the experimental results, where longitudinal conduction becomes important, due to the material properties or low welding speed.
An analytical approach is employed before the parameter definition by using the genetic algorithm and 2D parametric simulations. This preliminary approach enables to limit the population of 2D analyses to be performed reducing the range of research of the heat flux parameters.

Description of the semi-analytical method
Goldak's double-ellipsoidal heat source is defined by nine parameters [4], as shown in Eq. (1) for the front ellipsoid and in Eq. (2) for the rear ellipsoid. Generally, welding process parameters are known a priori. Welding Procedure Specification (WPS) collects data including: Voltage V o, current I and welding speed v. The efficiency η depends on the welding process used. Geometrically, the double ellipsoid is defined by four parameters: a, b, c f and c r . As shown in Figure 1, a parameter is the width semi-axis, b is the depth of the double-ellipsoid while c f and c r parameters are the front and the rear semi-axis of the ellipsoid. The last two factors represent the fraction of the total heat deposited in the front and rear ellipsoid. These factors are related to c f and c r geometrical parameters through the Eqs. 3 and 4. In these relations, it is also present the α parameter that is equal to two to guarantee continuity of the q functions at z = 0 ( Figure 1) [6].
(2)  Ultimately, the geometrical parameters are unknown in the simulation of a weld using the double ellipsoid schematisation of the heat flux (a, b, c f and c r ). Often these parameters can be related to each other, as reported in [7,8], with c r = 4c f and c f = b. In this way, the only unknowns to determine are a and b.
The method presented in this work makes it possible to determine the two unknown geometric parameters by knowing WPS parameters, weld pool extension and peak temperature T peak reached in the welded joint in a point as close as possible to the seam. The extension of the molten pool is characterised by measuring the depth of penetration (WPD) and the width (WPW) of the weld, as shown in Figure 2.
The method presented is based on the use of the NSGA-II genetic algorithm [9] applied to a 2D parametric transient thermal model. In order to determine the geometric parameters, these factors are varied by the NSGA-II algorithm. For each value pair, a and b, a 2D analysis is performed and the results in terms of T peak and weld pool extension (WPW, WPD) are analysed. Subsequent analyses are guided by the genetic algorithm according to the results and their deviation from the experimental data. At the end of the analyses, it is possible to determine the optimal parameters a and b by imposing a criterion for choosing between those solutions satisfying the above constraints, i.e. on Pareto front.
Through the use of Fachinotti's analytical method [6], have been obtained the starting values of the geometric variables to be set in the genetic algorithm. This analytical model solves the thermal field induced by a moving heat source in a semi-infinite body. However, this method neglects convective and radiative heat exchanges. Furthermore the thermal properties of the material are considered constant as the temperature changes. In any case, the preliminary determination of the geometrical parameters of the double ellipsoid makes it possible to reduce the optimisation burden in the subsequent phase by limiting the range of research of the a and b semi-axis. A MATLAB [10] code has been developed to determine the best set of geometric starting values solving the Fachinotti's method. In input, only the WPS data, the weld pool extension and the T peak in a point near the seam are needed, even at this early stage.

Reference case study
In order to test the novel semi-analytical method presented in this work, a numerical reference case study has been realised. The reference case (REF) is the thermal and structural simulation of a bead-on-plate weld. Geometry and position of the seam and the dimensions of the plate are reported in Figure 3. The plate is 10 mm thick.
Transient thermal analysis of the bead on plate welding has been realized with 3D linear brick elements with 8 nodes. In order to streamline the simulations, only half the plate was IOP Publishing doi:10.1088/1757-899X/1214/1/012023 4 simulated using symmetry. In the thermally altered zone the mesh is very dense (Figure 4) to better handle the high thermal gradients and to adequately fit Goldak's double ellipsoid equation. The elements in the center line of the seam have a side size of 0.5 mm, and the whole model is composed of 140,400 elements and 155,952 nodes. Convection and radiation are imposed on all external faces of the plate, except the symmetry face, that is considered adiabatic. The convective exchange coefficient is 15W /m 2• C and the emissivity coefficient is 0.8, while the ambient temperature, which coincides with the initial temperature of the plate, has been set to 20 • C. Table 1 shows the temperature-dependent thermal and structural properties of 304 stainless steel [11]. The solidus temperature of the 304 steel is 1,340 • C while the liquidus temperature is 1,390 • C. The welding parameters are resumed in Table. In the Table 2 are reported the double ellipsoid parameters (a, b, c f and c r ), that have been chosen in order to obtain thermal results with physical significance, using typical values that can be found in literature for this type of welding and material. In the same Table are shown: the speed of welding v, the voltage V o, the current I and the efficiency η.   [11].  Figure 4. Reference 3D FE model mesh. Simulation time step has been imposed very short, 0.05 seconds, to best simulate the movement of the welding torch. The reference model during the transient thermal simulation is shown in Figure 5. The grey zone in the figure represents the molten weld pool. The resulting extension of the pool and the peak temperature reached in the center line of the welding are reported in Table 3. These values represent the experimental parameters that has been used to test the method in object.

Analytical determination of the preliminary double ellipsoid parameters
As stated previously, the first step of the semi-analytical method is the determination by a MATLAB code of the starting value of the geometrical parameters a and b. This code solve numerically the Fachinotti's integral (Eq. 5), that does not have a closed-form solution. The index i in the Eqs. 6 and 7 can be referred to front f or rear r. The solution of the integral is implemented on a numerical mesh that has a limited extension around the centre line of the seam to limit the calculation time.
Implemented code solves the integral over the domain, by discretely varying the parameters a and b. For each set of a and b parameters, the peak temperature T peak and the T liquidus isotherm are calculated. These isotherms allow the determination of the WPD and WPW. Finally, these values are compared with the experimental ones. The pair of parameters a and b that allows the smallest deviation of peak temperature T peak and weld pool extension from the experimental data are automatically selected and become the starting parameters for the next step. The maximum and minimum values of the parameters a and b used in the calculation code can be related to the experimental dimensions of the molten pool. In this application: a min = 0.1 WPW, b min = 0.1 WPD, a max = 1.2 WPW, b max = 1.2 WPD.
The preliminary parameters obtained through the analytical method above described are given in the Table 4.

Final geometric parameters determination using the NSGA-II algorithm with 2D parametric analysis
The last phase of geometric parameters determination employs a 2D parametric model to perform transient thermal analyses guided by the NSGA-II genetic algorithm. The 2D model  thermally simulates the welding process by testing several sets of the a and b parameters. The FE model simulates only half plate by using the symmetry with quadrilateral 8-noded quadratic elements. The model is composed of 800 elements and 2,521 nodes, convection and radiation have been imposed in all external lines excluding the symmetry line, which was considered adiabatic. The mesh of the 2D model is shown in Figure 6. The area near the bead has square elements with 0.5 mm sides. Preliminary parameters enable to set the mean values of the variables (a mean = pa and b mean = pb) and the relative ranges of variation in the NSGA-II algorithm (a max,min = ±60% pa and b max,min = ±60% pb).  . By means of a selection criterion, the best solution will be identified within the group to which the final pair of geometric parameters a and b is associated. In this application, the solution that minimised the sum of the deviations of the T peak , WPW and WPD results was chosen, i.e. a linear choice method was set up with an equal weight function for the three results. The final Goldak's parameters, deduced through the method presented, are reported in Table  5.   The weld pool extension results and the comparison with the REF ones are reported in Table  6. WPW and WPD resulting from the thermal simulation with the final Goldak's parameters are in good agreement with the reference ones.
In Figure 8, it is showed the origin of Cartesian axes used to present the thermal results. The temperature trends are reported in Figure 9, evaluated at the center line of the seam, at middle length of plate, for both REF and FEM models. In this figure the experimental and the FEM T peak values can be observed. The latter has a percent deviation of −5.50% from the REF ones. In Figures 10 and 11 Figure 15. Stress on path D.

Conclusions
In conclusion, a semi-analytical method was developed for the determination of the geometrical parameters of the Goldak's ellipsoid knowing only a peak temperature near the welded seam and through the extension of the molten pool. The resolution of the Fachinotti's analytical model reduces the search range of the double ellipsoid parameters in the following step, which is implemented through the NSGA-II genetic algorithm that guides the 2D numerical simulations.
It is through such simulations that the final Goldak's ellipsoid parameters are determined. The thermal and structural results are encouraging even if obtained on preliminary numerical models of a bead on plate welding. It was decided to test the method presented on a numerical IOP Publishing doi:10.1088/1757-899X/1214/1/012023 11 reference model and not on an experimental test in order to exclude any unknown variables that might influence the preliminary tests.
The method can be easily automated through the use of fully parametric lists that can be adapted to the simulation of butt welds or fillet welds, including those realised with equivalent or hybrid models [12,13]. Moreover, this semi-analytical approach is well suited to determine the geometric parameters of multi-pass welds.