Numerical simulation of laser welded joints: modern fatigue analysis methodology

Welding always induces discontinuities and imperfections in the structure that allows for potential fatigue cracks. Welding effects thermal strains, which yield to residual stresses of the structure that have a noticeable effect on the fatigue behaviour of the structure. Welding inexorably leads to microstructure and geometry changes in the welding region. Material internal changes, residual stresses and microstructural changes can be simulated numerically, and the simulation results can be used in cyclic loading analysis in FEA.


Introduction
Welding always has a deterioration effect on the fatigue durability of structure [1][2][3][4][5][6]. The reason for the fatigue strength decrease is the weld joint-induced discontinuity. Welding inevitably leads to a discontinuity in geometry, microstructure and stress equilibrium imbalance, which have an effect on the material behaviour in static and in cyclic loading [5][6][7][8][9][10][11][12]. Residual stresses raise the stress levels in the welding region, but also alter material yield strength due to hardening. Therefore, simply adding residual stress on a weld joint and using base material properties may lead to a false outcome.
With a comprehensive simulation of welded joint residual stresses, microstructure, base material properties variation and deformations can be approximated [7,[13][14][15][16][17][18]. Deformations in structure and surface of a weld play an important role in the fatigue of welded joints. Angular misalignment, caused by deformation, includes additional stresses from a secondary bending stress [2,16]. A misalignment of a welded joint causes stresses arising from an additional bending component and the decrease of the nominal surface in the normal direction. The effects of misalignment are most evident in butt-joint welds and cruciform welds, where the increase in stress can be 30 % to 45 % [3].
Microstructure changes affect the fatigue behaviour [2,19,20]]. The terminal load needed to weld inexorably leads to microstructure changes in the heat-affected zone. Generally, laserwelding leads to an advantageous microstructure compared with traditional welding [2]. The material hardness rises in heat-affected zone in most metals [2,21,22]. The rise in hardness is due to the formation of a high-yield strength martensite. When cyclic loading is applied to the welded joint, the residual stress state results in material hardening and the rising of yield strength.
The welding simulation have been under study for decades as well as the fatigue behaviour of material. However, compounding of aforementioned is not done in many studies. The effect of residual stresses in fatigue analysis have been taken into account in study of Van Do et al. [23] and Shen et al. [6] with continuum damage mechanism approach. The residual stresses were taken into account in V-notch based analysis in [9,10,24]. In this study, the effect of residual stresses was studied by simulating welding processes and including the residual state in cyclic load analysis. The fatigue evaluated with convectional Smith-Watson-Topper equation with base material material properties.

Modelling welding
In welding, four states of material are present: solid, liquid, gas and plasma [15]. For a comprehensive simulation, multiple material variables and phenomena need to be taken into account. Welding simulation is temperature-dependent, as the material behaviour and state is dependent on temperature. Material properties also vary between states and temperature. For example, the steel modulus of elasticity constantly drops towards the melting point, as the steel becomes more elastic. In Figure 1, the temperature dependency of Yield strength, modulus of elasticity and thermal expansion coefficient are shown for S355 steel (properties are taken from [25]). In the liquid state, material properties vary significantly due to alloying elements. For example, the sulphur content for AISI 304 can change the flow direction of the molten pool surface, which affects the weld bead profile [15]. The simulation of welding can be divided into two different main approaches: the heat and fluid flow approach, HFF, and the thermo-mechanical and metallurgical approach, TMM [15]. The HFF approach is more comprehensive than TMM because it considers more phenomena; the HFF approach takes heat transfer, fluid flow, surface deformations and electromagnetic fields into consideration. In the HFF approach, a wide range of dynamic factors, such as surface tension, multiphase flows, electromagnetism and vaporisation, are taken into account.  Figure 2. Coupling of physical phenomenas according to Traidia, 2011 [15] The thermo-mechanical and metallurgical approach, TMM, is coupling heat transfer, metallurgical transformations and solid mechanics [15] and is more widely used in practice, as it easier to obtain [15,17,26].
Welding defects, such as porosity and undercuts, can be covered using the comprehensive heat and fluid flow approach [17]. Porosity is a significant problem in laser welding of aluminium and high-strength steels [17,26]. Porosity can be reduced with optimised welding parameters that can be obtained with numerical simulation. With the HFF approach, a wide range of physical phenomena, in addition to the weld pool dynamics, can be obtained, for example, residual stresses. As the HFF approach is widely multi-physical and takes into account multiple phenomena, the calculation is slower than the TMM approach. The TMM approach in its simplest is just coupling with the solid mechanism and heat transfer [15].

Simulation according to the HFF approach
The heat and fluid flow approach are used for welding parameter optimisation, in a defect, like porosity formation, and scientific studies [15,17]]. The fluid flow in a molten pool effects the shape and thus the quality of weld [15,17,27].

Fluid flow
The HFF approach takes heat transfer in the solid and fluid flow and the fluid into account. Molten metal flows in the weld pool are usually considered laminar due to the small size of the molten pool [15]. To evaluate the temperature and the hydrodynamics phenomena as a function of time, the classical incompressible Navier-Stokes equations for a conservation of mass (eq. 1), conservation of momentum (eq. 2) and energy conservation (eq. 3) are used [15,27,28].
where ∇ is a gradient operator, − → v is the velocity field in a weld pool, ρ is the density, p is the pressure, µ is the dynamic viscosity, − → F v is the sum of effective forces like gravity force, C p is the heat capacity, T is the temperature,Φ is the volumetric heat source in the weld pool, L f is the latent heat of fusion, and f l is the fraction liquid. The liquid fraction f l is dependent on the liquidus temperature, T l , of the material.
Unevenly distributed heat in a molten pool creates fluid flows between temperature zones [15,17]. The density of a molten metal depends on the temperature, and thus high terminal gradients induce natural a convection flow (Buoyancy effect) [15]. Buoyancy forces tend to create an outward flow from a heat source, increasing the width of a weld pool. Velocity fields induced by surface forces, such as, the Marangoni effect are stronger than velocity fields induced by natural convection. Also, forces such as gravity force are greater than natural convection force, but it does not contribute to flow creation as it is a constant. The buoyancy effect can be modelled with additional force [15]: where ρ is the density, k is the thermal expansion and T ref is the reference temperature. The Marangoni effect is one of the main components in molten pool flow forces [15]. The Marangoni effect on free surface can be described with equation 5; Marangoni surface tension determines the flow direction of the Marangoni flows. With a negative coefficient, the flow is outward, as the surface tension gradient is at its highest at the edge of the weld pool, and with a positive coefficient, the flow is inward. For pure metal, the coefficient is negative, but the presence of sulphur or oxygen can alternate it to positive. However, the coefficient is also dependent on temperature.
where µ is the dynamic viscosity, ∇ n is the normal gradient operator, v s is the tangential velocity, ∂γ ∂T is the surface tension gradient, ∇ s is the tangential gradient operator, and T is the surface temperature.

Laser beam absorption
The laser beam energy is absorbed with Fresnell absorption in the keyhole. The Fresnell absorption can be represented by a function of the incident angle, laser dependent coefficients and material-dependent coefficients [17]. In the model the beam is distributed into a discrete form. The absorption rate is calculated for each reflection point, and the absorbed energy is reduced from the total beam energy.

Simulation according to the TMM approach
The thermo-mechanical and metallurgical approach is usually used in solving weld bead (molten area) dimensions, metallurgic distributions and residual stresses and distortions [13,14,16,29,30]. The most influential variable in the TMM approach is the heat source. The heat source is given as a constant form heat source that is dragged through the welding line. To obtain good results in the form of heat, the course needs to be right.

Heat source
In the TMM approach, the most commonly used heat source in laser welding simulation is the Gaussian distribution or a combination of multiple Gaussian distributions [13,27,30].  The heat source is modelled as a moving heat source [14,29,30,32]. The Gaussian heat flux can be given in 2D in the following form: where P int is the power conducted to the part, r 0 is the initial radius, and r is the current distance from the weld centre.
The temperature evolution can be evaluated in the entire domain with a modified classical conservation equation of energy by excluding terms related to flow [33]: where ρ is the density, C p is the heat capacity, ∂T ∂t is the temporal gradient of temperature, k is the thermal conductivity, − → ∇T is the spatial gradient of temperature, and Φ is the source term. The course term Φ contains boundary conditions like heat flux, convection, radiation and contact conductivity [33].

Material model
The material model has a key role in the simulation of welding. The material model is composed of the thermal behaviour and the elasticity of the material. The total material strain composition is shown in equation 8 [18]. The elastic part of the strain can be expected to follow Hooke's law; the plastic strain part can be defined with hardening laws such as the power law or equivalent (equation 9): where ε a is the total strain, ε el is the elastic strain, ε pl is the plastic strain, ε th is the thermal strain, and ε ph is the phase transformation strain due to the volume change between various phases.
where σ is the total stress, σ y is the initial yield stress, H is the hardening law coefficient, and N is the hardening law exponent, α th is the temperature-dependent coefficient of thermal expansion, and I is the identity tensor.

Phase transformation
In steel welding, phase transformation is always present due to the high thermal input. Phase transformation involves notable changes in the kinetics of a metal, which are not always associated solely with the external thermal impact, but also with the rate of thermal changes [34,35]. The phase fractions in a solid state after cooling can be obtained through mathematical models such as Jonson-Mehl-Avrami-Kolmikorov, Leblond-Devaux and Koistinen-Marburger. If more than two phases are calculated, the transformations of phases are included with CCT diagrams.

Material model affected cyclic loading
In cyclic fatigue analysis, the nonlinear behaviour needs to be taken into account. The transformation from the elastic to plastic part is defined with a yield function. The yield surface definition by von Mises is given in equation 11 [36]. In the theory of plasticity strain increments are rate-dependent. The rate and direction of plastic flow are determined with the flow rule (equation 12) [36]. The direction of plastic flow is in the normal direction about the yield surface.
ε pl =λ ∂f ∂σ (12) where S is the deviatoric part of the stress tensor, : is the double contraction, σ y is the yield stress, tr is a trace operator andλ is the plastic multiplier.
Exceeding the material yield strength develops plasticization. Plasticization affects material properties, like yield strength, by increasing or decreasing it. The hardening is accounted by two models: isotropic and kinematic hardening. Isotropic hardening determines the change of the yield surface size and kinematic hardening the change of the yield surface location. Isotropic hardening is accounted for using exponential laws after Voce or equivalent (equation 13) [36]. Kinematic hardening is modelled by taking the plastic behaviour in respect of "kinematic shift", α [36] into account. The function defines the pressure-independent yield surface in respect of backstresses given in equation 14.
where σ| 0 is the yield surface at zero point plastic strains, Q ∞ is the maximum hardening, b is the hardening slope, ε pl is the equivalent plastic strain, σ is the stress tensor, α is the backstresses and σ 0 is the yield surface.

Material fatigue model
Fatigue behaviour can be considered for parts and structures under cyclic loading. Fatigue propagation can be calculated through the relation of the number of cycles to the initiation of a macrocrack [37]. This relation is called the damage parameter. The Smith-Watson-Topper, SWT, is the most used damage parameter [5,37,38]. In Smith-Watson-Topper, the stress, σ n,max , is calculated in the maximum strain plane direction [39]. (15) where σ n,max is the maximum stress, σ f is the fatigue strength coefficient, N f is the total number of cycles, b is the fatigue strength exponent,ε f is the fatigue ductility coefficient and c is the fatigue ductility exponent.

Fatigue behaviour analysis including residual stresses
The effect of residual stresses of the welded joint in fatigue behaviour can be studied by including them in cyclic loading analysis. Residual stress decreasing the effect of fatigue life has also been a point of discussion in several studies [1,6,9,[40][41][42]. Residual stresses effect material behaviour in cyclic loading and, thus, fatigue behaviour. Residual stress state can be simulated and included in the cyclic analysis.
The analysis was done to butt joint. The test rod geometry was according to the test rod geometry used by Liinalampi et al. (2016) [43] and Lillemae et al. (2017) [32] in fatigue tests. The geometry of the test rod was ideal: it did not include any initial offset, misalignment or other dimensional or manufacturing defects. The geometry defects affect the fatigue behavior of the welded joint. With ideal geometry, residual stress effect will accent.
The dimensions and mesh of the test rod used in the simulation are shown in Figure  5. The element size was 0.3 mm in the welding region and increased steadily to the end of the test specimen. One end of the simulation test rod was constrained in all directions and another end in Y (22) and Z (33) directions (see Figure 6). The thermal convection 20 W/(m 2 K) was set on all free surfaces. The welding heat source was a combination of a double ellipsoid and a double cone shape with widths were 1.8 mm and 1.2 mm. The power of the weld used was 2000 W with an efficiency rate of 65 %. The welding speed was 0.022 m/s. The temperature depended material model of S355 was used in welding simulation with an annealing temperature of 1400 o C. Material temperature dependency is shown in Figure 1. The calculated residual stress field in the test rod is shown in Figure 6. In the literature, it is stated that the residual stresses are in material yield strength. The residual stresses based on welding simulation were partly exceeding the original pulling test yield strength. However, because the yield surface has changed due to the hardening, the residual stresses are within yield strength. Transverse stresses, σ 22 , are the highest in value, but the longitudinal stresses, σ 11 , have the strongest effect on fatigue strength [14]. The welding was simulated with Geon-X's VIRFAC program. VIRFAC uses the linear element in the welding simulation, which caused some discontinuity in the stress field. The discontinuity would have been smaller with parabolic elements. In a cyclic analysis, S355 base material properties were used in the fatigue analysis. The material phase change effect on the material properties was obsolete. The cyclic analysis was done with Abaqus -program. The simulated residual stresses and hardening resulting from stresses over material yield strength were calculated as an initial condition to the Abaqus. The base material properties used (taken from [37]) are shown in Table 1.  In Figures 7 and 8, the difference in the gradient fields between the test rod with introduced model, including residual stresses and blanc test rod, is shown. In Figure 7, the gradient field of strain and stress when the test rod with residual stresses was affected to tensile stress. The gradient field does not have a clear trend in direction. In a normal case (Figure 8), the gradient field has a clear trend in the direction of the load.  Because of the stress gradient direction SWT parameter as it does not work. Therefore, the stresses were calculated according to the stress definition of Sines -a multiaxial fatigue strength criteria method [44]. In the method, the stress criteria are the square root of the second invariant of the Cauchy stress tensor. Stress amplitudes determine the stress tensor. The term of mean hydrostatic stress is added as shown in equation 16 [44]. In the calculations, the stresses were defined as √ where a s and b s are the coefficients associated with fatigue criteria.
The calculated fatigue strength with the stress ratio R = 0 (Sim. SWT (R=0)) and R=-1 (Sim. SWT (R = -1)) are shown in Figure 9. The least square method was used in curve fitting. For comparison in figure laser welding fatigue test results from the literature are also plotted. Test results from the literature use S355 as the base material, and the stress ratio is close to R=0.  [43]. A butt joint presented in Figure 5 was used in the above-mentioned. The test data from the KeKeRa project (Ke1) [45] was shear data and, thus, multiplied by √ 3. The Ke1 tests were done with lap-joint. Cyclic loading altered the material stress-strain behaviour. The hysteresis loop of stress and strain moved as a result of hardening. The stresses were levelled at cyclic loading.
The poor mesh caused stress concentrations on some particular elements. In these elements, the number of cycles for failure was close to zero. This effect was corrected by excluding 5 % of the N f results in the calculation of the number of cycles for the S-N curve.
In Figure 10, the slope of S-N curves is compared; the calculated slope corresponded well with the slopes from literature. The lines are fitted to each data point sets with the least-squares curve fitting. The comparison is done with tests with stress ratios close to zero. It can be concluded from the results that residual stresses have a significant effect on fatigue behaviour. The slope of the curves coincides with the curves fitted from the literature test data.

Conclusions
With numerical simulation, the as-welded state of a structure can be computed, and the results can be used in structure behavior analysis. A wide range of phenomena such as residual stresses, phase changes, liquid flows, etc. can obtain through simulation. But by including an excessive amount of physics, the analysis is more challenging and timeconsuming. Thus some model-reduction is needed in practice. The scope of a welding simulation can be optimized to correspond to the results needed, e.g. residual stresses. It is stated in the literature [14,16,32] that the distortions have a deteriorating effect on fatigue behaviour of welded joint. The distortions resulting from the welding are constituted in residual strains. In this study, heat transfer and solid mechanics were taken into account to simulate residual stresses. The effect of residual stresses in cyclic loading was studied, and the results showed that they have a clear effect on welded joint fatigue behavior. The residual stress state affects stress equilibrium and raises the effect of residual distortion. The residual stresses also effect on material behavior through material hardening, which has an effect on the material cyclic stress state. The results of a fatigue analysis made to test rod with simulated residual stress state, corresponded well with the test results from the literature.