Fast 3D simulation of a single-pass steel girth weld

– The prediction of welding residual stresses requires accurate account to be taken of the couplings between heat transfer, metallurgy and stresses-strains in the heat aﬀected zone. To begin with, simulations of residual stresses for welding processes were performed at the beginning of the 1970s. Since then, calculations have been compared and validated with experimental measurements, the major problem remaining is the calculation time. Despite the technological evolution of computers, a 3D calculation can last several days. To avoid this diﬃculty, a 3D simpliﬁed approach is proposed in this article. It consists of decreasing the computation time for the current zone of the welds. To show the relevance of such an approach, the methodology developed is presented through the application of a single-pass steel girth weld. For this application, the computing time can be reduced by more than 67% compared to a standard “step by step” simulation.


Introduction
Welding has been in existence since the beginning of the 20th century and it has been widely used in many industrial fields. This assembling process induces alterations of structural and mechanical properties. Indeed the heat input causes heterogeneous cooling cycles in the material; there follow residual stresses in the heat affected zone (HAZ) and the seam weld. Measuring residual stresses inside the weld to study the propagation path of a crack, for example, is utterly impossible [1]. The numerical simulation of the welding process is therefore the only way to access the residual stress field by modeling heat transfer, metallurgy, and mechanics, as summarized in Figure 1.
For the moment, 3D numerical calculations of residual stresses are limited to several industrial fields such as the aeronautic or the nuclear industry. This is essentially due to the calculation times which can last several days. A large amount of studies have been carried out so as to decrease the computing times. For example, dynamic remeshing and later on, adaptive meshing were studied in welding simulations by Lindgren et al. [2], Runnemalm and Hyun [3], Duranton et al. [4]. Nastrom et al. [5] used a combination of solid and shell elements to model welded a Corresponding author: eric.feulvarch@enise.fr structures in order to reduce the number of finite elements.
Another approach, the first historically documented, consists of simplifying the problem with a 2D axisymmetric assumption. For example, Sarkani et al. [6] have shown that the thermo-mechanical 2D model of the single-pass weld of a T-joint can give fairly good results compared with those computed by a standard 3D model for the current section, except for the hoop stress which is 30% overestimated. All computations have been carried out without considering metallurgical phase changes. Another similar recent study of Xu et al. [7] shows that the axisymmetric 2D thermo-mechanical model of a butt-welded steel pipe leads to satisfactory results compared to the distribution of stresses obtained with a standard 3D model.  However, a significant difference in stress values can be observed for the hoop stress which is also overestimated. This is also observed by Duranton et al. [4] for a multipass weld where the metallurgical couplings have been taken into account to include tempering effects induced by each pass on the heat affected zone created by the previous passes. Tempering effects result in lowering the yield stress of the material. The consequences are very important to the residual stresses, and the temper bead welding technique is sometimes used to reduce the stress values intentionally. For a case study similar to the one of Duranton et al. [4], Feulvarch et al. [9] have shown that the 2D numerical results are in good agreement with stress measurements obtained by neutron diffraction and the deep hole drilling technique in the current region.

Article published by EDP Sciences
In this paper, we suggest a 3D simplified approach which can capture the stress distribution in the current zone but also on the start/stop regions. Given that the 2D assumption generally gives quite good results, the principle is to treat conventionally the start and stop zones with an equivalent 3D heat source and compute all the current region at the same time with an equivalent 2D source. The gain in computation time is thus directly related to the length of each weld and can be very significant. This methodology is presented through the application of a single-pass steel girth weld. The approach is validated in two steps for the application of a single-pass steel girth weld. First, the thermo-metallurgical results obtained are discussed according to those resulting from the standard "step by step" 3D analysis. The residual stresses are then compared showing the relevance of the simplified approach.

Geometry, mesh and manufacturing process
A T-joint links a tube with a 20 mm internal radius and a thickness of 12.5 mm with a 8 mm thickness plate (see Fig. 2). A hole with a 65 mm radius is drilled into the plate, the pipe then goes through this hole. The weld is made using the MAG (Metal Active Gas) process and a single-pass technique with a welding speed of 660 mm.min −1 . The welding stage starts and stops at the same point corresponding to the cross-section located at the angle value of 0 • . All components, including the filler metal, are made of ferritic steel of type S355. For the computations, all the data about the ferritic steel S355 are given in reference [8]. The 3D mesh used for the finite element simulations is plotted in Figure 3. The geometry of the seam weld has been design by means of a macrograph (see Fig. 4).
where ρ is the material volumetric mass, H is the specific enthalpy; T is the temperature; λ is the thermal conductivity; T (p) is the temperature of air and n is the unit outward normal vector to the boundary ∂Ω of the domain Ω representing the volume of the work pieces. h is a coefficient which represents the heat exchanged between the mechanical components and the air supposed to be at the temperature T (p) . The volume source term Q involved in the heat problem corresponds to the power provided by the welding process within the material per volume unit. The weak formulation of the heat equation is classically obtained by multiplying Equation (1) 1 by weighing function T * and integrating over the domain Ω. Integrating by parts and accounting for the boundary condition (1) 2 , one thus obtains the following weak formulation: Following the usual procedure, the weak formulation (2) is applied to the finite element approximation of the function T given by In this expression, N denotes the number of nodes, T i the value of the function T at node i and N i the shape function associated to this node. Following Galerkin's approach, the weighting function T * takes the same form.

Metallurgical transformations and thermo-metallurgical couplings
Three types of interaction between thermal and metallurgical phenomena must be considered [9]: metallurgical transformations depend on the thermal history, metallurgical transformations are accompanied by latent heat effects which modify temperature distributions, thermal properties are phase dependent. Different approaches can be used to describe the transformation kinetics in steels [10][11][12]. These approaches rest upon the modeling of transformations under either isothermal conditions (IT diagrams) or anisothermal conditions. In the latter case, parameters come from CCT diagrams. The model we use in this study belongs to the second category. Generally, those models can be expressed, in the case of one transformation, with a differential equation. Leblond and Devaux [11] have suggested the use of a simple first order differential equation of the form: where p is the proportion of the created phase. p (T ) denotes the phase proportion obtained after an infinite time at temperature T ; τ (T ) a time delay depending on temperature T and F (Ṫ ) a function representing the dependency of the transformation rate on the temperature rate. These parameters are adjusted for each transformation in order to represent the CCT diagram. One can note that if the F (Ṫ ) is chosen proportional toṪ , then the resulting phase proportion only depends on temperature as in the case for martensitic transformation which follows the Koistinen-Marburger law [13]: where M s and b respectively characterize the temperature at the beginning of transformation and the evolution of transformation with temperature. The effect of stresses on transformation kinetics has been widely examined [14], but, in practice, the lack of material data is a strong limitation to the application of models so that this effect is neglected.
In the present case study, the following metallurgical transformations are considered for the ferritic steel S355: Once p is known, the enthalpy H of the mixture may be obtained from the simple "linear mixture rule": where H 1 (T ) denotes the enthalpy of the parent phase 1, the proportion of which is (1−p) and H 2 (T ) that of the product phase 2, the proportion of which is p. This expression disregards possible thermodynamic interactions between the phases, it is widely considered to provide an acceptable approximation of the enthalpy of alloys during phase changes. The material volumetric mass ρ and the thermal conductivity λ of the mixture may be obtained in the same way. The model has been generalized in order to consider several transformations between several phases [11]. The material thermal properties are then obtained through mixture laws between the thermal properties of all the phases.

Mechanics
The mechanical analysis is based on the momentum equation where inertial effects are neglected.
where σ is the Cauchy stress tensor. It is performed with the infinitesimal strain theory. Therefore, σ is assumed to depend linearly on the elastic strain tensor e by means of the fourth rank elastic tensor C as follows As the effect of plastic dissipation on heat transfer and the influence of stresses on metallurgical transformations are neglected, the mechanical analysis can be uncoupled from the thermo-metallurgical simulation. In this way, heat transfer and metallurgy are involved in the mechanical analysis through four effects [9]: the thermal strains; the volume changes due to the transformations (contraction during heating and expansion during cooling as shown in Fig. 5); the influence of temperature and of the phases on the behavior law (generally, the yield strength of phases is higher if they have been created with high cooling rates); the transformation induced plasticity.
The strain rate tensor˙ is thus decomposed into a thermal and metallurgical part˙ th , an elastic part˙ e and a plastic part˙ p that includes both classical plasticity and transformation induced plasticity: The volume changes due to transformations are a major cause for the appearance of residual stresses and strains. To include those changes, a simple solution consists of changing the standard thermal strain into the following thermal and metallurgical strain: where th 1 (T ) denotes the thermal strain of the parent phase 1 and th 2 (T ) that of the product phase 2. th reflects both the thermal strains of the phases and the volume change between the latter as shown in Figure 5 [15]. The stress-strain relation must be temperaturedependent, representative of the phase mix, and reproduce the transformation induced plasticity phenomenon. The material behavior law during transformation is usually supposed to be elastoplastic. Plastic strain rate˙ p is expressed as the sum of three terms, proportional to stress variationσ, temperature variationθ and phase proportion variationsṗ respectively: The first two terms represent the conventional plastic strain rate, while the third one represents the transformation-induced plastic strain rate. Transformation plasticity can be associated with two physical phenomena: -Greenwood-Johnson mechanism [16]: the volume difference between the phases can generate microscopic internal stresses which are sufficient to induce plastic strain in the weaker phase (austenite). -Magee mechanism (martensitic transformation only) [17]: in the presence of an external stress martensite needles are formed in a preferential direction.
The plasticity induced by transformation has been studied by many authors [16][17][18][19]. Different models of transformation plasticity have been proposed, but most deal with the Greenwood-Johnson mechanism only. Among all the models related to the behavior of steels during phase transformation, the model proposed by Leblond and coworkers [20], based on a micro-mechanical analysis, is widely used. The model requires the temperature dependent stress-strain relations of all the phases which is quite difficult to obtain [21]. The hardening can be isotropic or kinematic, and viscoplastic effects can also be considered [22]. In the present study, the behavior of the different constitutive phases of S355 steel is assumed to be elastoplastic with isotropic hardening. All details on the elastoplastic behavior including the mixed phase relations used in this work are detailed by Devaux et al. [23].
To apply the finite element method, the weak formulation for the mechanical problem is obtained by multiplying Equation (7) by weighing function U * and integrating over the domain Ω. The boundary conditions can be prescribed displacements or stress vectors imposed. Integrating by parts and considering that the surfaces of the mechanical components are stress free, one thus gets: Following the usual procedure, the weak formulation (12) is applied to the finite element approximation of the displacement U which is similar to the one used for the temperature approximation:

Thermal simulation
As mentioned in the previous section, the internal heating due to plastic dissipation is always neglected considering the small strain rates generated by a welding operation [24]. Therefore, the simulation can be achieved through a staggered method which consists of splitting the analysis into two steps. First, the temperature and phase variations are determined as a function of time. Then, the mechanical computation uses the previous results to get displacements, stresses and distortions.

Equivalent heat source modeling
To compute the thermal cycles, several input data are necessary: the geometries of the different parts, the material properties and the equivalent heat source. Some methods are proposed in literature in order to model the heat input. For example, Goldak et al. [25] propose a two half-ellipsoid heat source shape which needs to be readjusted by means of multiple coefficients. This method is well adapted to processes with filling materials. The readjustment can take place using thermocouple sensors to measure the temperature at different locations but this method cannot always capture the fast temperature changes near the welding pool because of stability problems [7]. A second approach consists of preparing the macrograph of a cross-section of the weld; and fitting heat source parameters to make sure that simulated geometric  appearance of and the molten zone matches the geometric dimensions observed on the weld macrograph as shown in Figure 4.
In this study, we have chosen the second method for its simplicity by considering a simple heat source shape plotted in Figure 6. The heat input is represented by two cylinders of power per volume unit Q(x, t) (see Eq. (1)): a small cylinder of high power (Fig. 6A) in order to simulate the penetration observed on the macrograph and a second bigger one of lower power (Fig. 6B) for modeling the shape of the molten zone. Figure 7 shows the comparison between the molten zone identified on the macrograph, and the distribution of maximum computed temperatures on the cross-section located at an angle of 180 • in the current region (Fig. 2). The tridimensional simulation is carried out on the whole geometry of the weld using the 3D mesh plotted in Figure 3 and a standard implicit step  Fig. 9. Comparison of the temperature evolutions at node E in 3D at 180 • and in 2D (see Fig. 3).
by step approach. In this paper, all the simulations have been performed using the computer code Sysweld r [8].

3D simplified approach
As mentioned in the introduction, many authors have shown that results obtained with a standard 3D FEM approach can be approximated in the current region by means of 2D axisymmetric modeling. This configuration is less time consuming and less complicated. It corresponds to the assumption that all the filling material is deposited at the same time. Unfortunately, the 3D equivalent heat source can not be directly implemented in the 2D simulation, because the associated heat transfers are also 3D.
Let Q(x, t) be the power per volume unit representing the 3D equivalent heat source defined in the previous subsection; Sarkani et al. [6] have shown that a 2D equivalent heat source can be obtained by multipying the 3D source by a coefficient γ. Let us consider a cross-section, the authors explained that the coefficient γ for this section is due to the power transferred by the neighbour cross-section where the welding torch is located before its arrival. In this work, this coefficient has been fitted using a 2D axisymmetric model based on the 2D mesh plotted in Figure 3. The value of γ has been obtained by comparisons of the maximum temperatures (see Fig. 8). It is taken equal to 0.96. As shown in Figure 9, the temperature evolutions are very close. Unfortunately, 2D modeling does not allow to simulate heat transfer hence residual stresses in the start and stop zones.
The principle of the technique developed in this paper is to use the 3D equivalent heat source for the start and stop areas and the 2D equivalent heat source for the current region as shown in Figure 10. Let Q s (x, t) be the new 3D heat source; the simplified 3D approach consists of replacing Q(x, t) by Q s (x, t) in Equation (1). Considering a weld composed of n sections with a current region containing m sections, the position of each section k is denoted x k . As shown in Figure 10, the new heat source is identical to the initial equivalent heat source in the start zone. For the current region, the heat input is the same in all sections. It is the same as the one received by the section numbered i + 1 for a standard simulation, multiplied by the coefficient γ. γ is necessary because there is very little heat transfer from one section to another in the current region. When the heating starts in the current region, at the same time, the section j receives a heat input identical to the one of the section numbered i + 1 for a standard simulation. In this way, the time for the heating of the current region does not explicitly appear in the simulation. Then, the source continues its path along the welded joint in the stop area. As a result, the gain of computing time can be approximated by the ratio m/n for the proposed case study. Figures 11-13 show the change in temperature for different position of the heat source. The first position corresponds to an angle of 22.5 • of the axis of the heat source. This angle is common to the three simulations since it is upstream of the welding current region. The second position corresponds to an angle of 180 • for the step by step simulation. For the simplified simulations, the second angle is the one of the beginning of the current region. The third position is out of the current region just before coming back to the start point. One can note that for the third case (see Fig. 13), the start area has not enough time to cool down and the temperature is still very important. This can be observed in Figure 14.
A comparison of the distribution of the maximum temperatures is plotted in Figure 15 for the standard simulation (case a) and the new simplified method with 2 lengths of the current zone: 90 • → 270 • (case b) and 45 • → 315 • (case c). One can note that the results obtained with the simplified approach are in good agreement with the ones of the standard simulation. As shown in Figure 16, the temperature evolutions are again very close to the standard simulation. In Figure 16, the time has been readjusted to correspond to an identical instant of heat input through the cross-section located at 180 • .

Mechanical analysis
Once the thermo-metallurgical simulation has been carried out, the mechanical computation uses the temperature and phase distribution in order to get displacements and stresses. For practical reasons, the mesh used for the thermo-metallurgical computation is also the one used for the 3D mechanical simulation. This approach avoids the problem of the internal variables projection from one mesh to another. The finite elements corresponding to the filling material are mechanically activated under the solidus temperature to model the addition of material according to the welding thermal kinetics [9]. Figure 17 shows the distribution of the von Mises residual stresses for the 3 thermo-metallurgical simulations developed in the previous section.
As shown in Figures 17 and 18, a slight difference in the current region appears between the reference (a) and the simplified simulation (b). This can be explained by the fact that the weld is relatively short and that the calculation (a) does not exhibit a real steady state corresponding to the current region. Indeed, the hypothesis of a current region appears clearly  for case (b) on half the welding (90 • → 270 • ). The difference with the standard simulation will be mitigated in the case of longer welds. The difference appears more clearly with case c. This is explained by the fact that the starting and stopping areas are halved compared to case b. This last case shows that the lengths of the start and stop zones must be chosen carefully in order to obtain satisfactory residual stresses in the whole weld. As shown in Figure 14, the start area in case (c) has not enough time to cool down as quickly as for cases (a) and (b). Therefore, the thermo-metallurgical state of the start zone differs significantly compared to the reference simulation (a) when the heat source reaches the end of the welding sequence.  allow it to cool down sufficiently. Nevertheless, the hoop stress is well computed in the current region for cases (b) and (c) as shown in Figure 19, compared to the reference case (a).
The computing times on a standard 3.4 GHz PC with 4Go memory are given in Table 1 by using an identical time step for the time integration.

Conclusion
Welding simulations are all the more difficult as they need refined 3D models in order to capture multiphysical couplings in the HAZ but the calculation times can last several days. In this paper, a 3D simplified approach has been proposed. It consists of decreasing the computation time for the current zone of the weld. To show the   relevance of such an approach, the methodology developed is presented through the application of a single-pass steel girth weld. For this application, the computing time can be reduced by more than 67% compared to a standard "step by step" simulation. In a general way, the gain can be estimated by the ratio of the current region length to the weld total length. It can be very significant for very long welds. This approach can be easily extended to multi-pass welding or to other continuous deposit processes such as additive manufacturing.