Numerical simulation of welding distortions in large structures with a simplified engineering approach

2.1 Thermal FE analysis

The thermal analysis of the welding process is based on a transient nonlinear heat transfer differential equation:

that incorporates the thermal material properties which are strongly temperature dependent. Here, k_x, k_y, and k_z are the thermal conductivities in the x, y and z directions, respectively; T denotes the temperature; Q is the heat input; ρ stands for the material density; C represents the specific heat capacity of the material; and t is the time. In order to solve Eq. (1) the following initial and boundary conditions are introduced

where N _x, N_y, and N_z are the direction cosine of the normal to the boundary; h_c and h _r are the convection and radiation heat transfer coefficients, respectively; q_s represents the heat flux over the body boundaries; T_r denotes the temperature of radiation; and T_∞ represents the ambient temperature. Radiation heat loses are expressed by the relation

where σ denotes the Stefan–Boltzmann constant, ϵ is the surface emissivity factor, and F is the configuration factor. The total heat input applied to the weld can be expressed with the following equation:

where η is the efficiency of the welding process, I is welding current, U is the arc voltage, and V_H is the volume of heat source, i.e. the volume of filler material.

According to the governing Eq. (1), the thermal FE equations for an element can be obtained by using the Galerkin method, as follows:

where c is the specific heat matrix, T represents the nodal temperature vector, k_T is the thermal conductivity matrix, while f is the nodal load vector in terms of internal heat source and boundary conditions. Matrices c, k and f are defined in terms of element shape functions matrix.

2.2 Mechanical FE analysis

Mechanical analyses in this work are conducted through standard elastoplastic procedure employing the implicit return mapping algorithm. Herein, it is assumed that both base and weld materials are homogeneous, isotropic, elastic-perfectly plastic solid that yields according to the von Mises criterion with the associated flow rule. The nonlinear material behavior is thus modeled via incremental plasticity while also assuming the geometrically nonlinear behavior of the welded plates. Furthermore, in the mechanical analysis of the present study, the influence of the steel phase transformation is not considered because its influence on the residual stress field and deformations is small in low-carbon steel [23].Additionally, the creep material behavior is neglected as the duration of the thermal cycles during the welding process is short. Accordingly, the total strain increment can then be additively decomposed into three components:

where the subscripts “e”, “p” and “th” denote reversible elastic, irreversible plastic and thermal strain increments, respectively. The governing equations used to carry out the mechanical FE analysis of the welding process are the well- known equilibrium equation (8), strain–displacement relationship (9) and constitutive equation for thermal elastoplastic material (10):

where L is the Cartesian coordinates partial derivatives matrix, σ is the stress tensor, b represents the body force vector, ϵ is the strain vector, u represents the displacement vector and C is the elastic constitutive material matrix. Using the virtual work principle and standard finite element method manipulations, the finite element equation can be derived in an incremental form as

where k₁ is the element stiffness matrix, k₂ is the element thermal stiffness matrix, v represents the vector of nodal degrees of freedom (DOF) and r is the nodal force vector. To solve the nonlinear equations (11), the iterative Newton-Raphson method is employed within the FE solver. Here, the updated values of the state variables at the end of the time step (ⁿ⁻¹t, ⁿ t) have to be calculated for the given values of the incremental strain tensor components Δε and the state variables at time ⁿ⁻¹t. More details about an algorithm for solving non-isothermal elastoplastic problems with the finite element method can be found in the authors’ earlier paper [24].

3.1 Small-scale T-joint model

The geometry of the small-scale T-joint model, including the welding sequence and direction, is presented in Figure 1a. The size of the horizontal base plate (flange) is 500 × 500 mm and the height of the stiffener (web) is 140 mm. The thickness of the base plate and the stiffener are 9 mm and 12 mm, respectively. It is assumed that two steel plates are welded by fillet welds with 6 mm leg lengths in two single passes with no time in between the passes using the metal active gas (MAG) welding procedure. Since the initial gap between the plates was very small and the plates were tack-welded prior to the start of welding process, the influence of the small gap on the final residual stress and deflection field can be neglected [20]. Therefore, both the horizontal and vertical plates were modeled as a single unit. Moreover, due to the negligible heat input to the work piece, the tack welds’ heat influence has not been taken into the account in the numerical analyses conducted in this study. Furthermore, the small initial gap between the horizontal and the vertical plate is not taken into account. The welding parameters are provided in Table 1. The plates are free welded, without any restraints, but in the mechanical analysis the boundary conditions are imposed to the horizontal plate to prevent its rigid body motion as it is schematically shown in Figure 1. The material of the plates is a low-carbon shipbuilding steel SM400A whose chemical composition is given in Table 2. Its thermo-mechanical temperature dependent properties are given in Figures 2 and 3. Due to the lack of thermo-mechanical data for the weld filler material, the material properties were assumed to be the same as the base metal. The dimensions of the solid and shell parts for the shell/3D small-scale T-joint model are given in Figure 1b. Since the small-scale T-joint model is geometrically very similar to the model analyzed in [10] and is subjected to the same heat input, the size of the 3D zone in the shell/3D model is estimated in accordance with the recommendations set out in that authors’ previous study. Therefore, to obtain a continued stress transfer across the shell-to-solid interface, the size of the 3D zone is set to approximately 3 wall thicknesses of the thicker plate.

Figure 1

Geometry of small-scale T-joint model: (a) full 3D model and (b) shell/3D model

Figure 2

Thermal properties of SM400A steel [26]

Figure 3

Mechanical properties of SM400A steel [26]

The heat generation rate approach based on the heat flux prescribed on a weld element volume is employed for modeling the welding heat input in the thermal analyses. The uniformly distributed heat flux per weld volume Q = 5.22 × 10¹⁰ Jm⁻³s⁻¹ is applied [27, 28] according to Eq. (5) and Table 1, where the weld volume V_H is determined by the weld cross-sectional area and the welding speed. On the outer body surfaces, the following boundary conditions are assumed: convective heat transfer coefficient k = 10 Wm⁻²K⁻¹ and emissivity ϵ = 0.9. The welding process efficiency η = 80% is taken in accordance with the EN ISO 1011-1 norm.

Finite element meshes for the small-scale T-joint model are shown in Figure 4. To discretize the full 3D model (Figure 4a), 3D linear hexahedral elements with 8 nodes and full integration scheme are used: DC3D8 elements for the thermal analysis and corresponding C3D8R elements with reduced integration scheme for the mechanical analysis. As may be seen from Figure 4a, a fine mesh is used in the fillet weld region for both longitudinal and transverse directions to capture very high gradients of temperatures, strains and stresses caused by the heat, while the regions far away from the welding zone are discretized with a relatively coarse mesh to reduce the overall model size. The finite element mesh of the full 3D model consists of 16,068 elements. For the discretization of the shell/3D model (Figure 4b), DC3D8 elements and DS4 shell elements are applied for the thermal analysis and corresponding C3D8R and S4 shell elements are applied during the mechanical analysis to improve computational efficiency and reduce storage space. Here, the solid part of the shell/3D model is discretized with 3D continuum elements using the same mesh density as in the full 3D model (Figure 4a). The shell/3D FE mesh (Figure 4b) consists of 7,332 solid and 2,184 shell elements. The solid to shell transition is easily achieved by constraining the interface nodes to couple the motion/ temperature of a shell edge to the motion/ temperature of an adjacent solid face. The Abaqus “tie” constraint in the thermal analysis and the “shell to solid coupling” option in the mechanical analysis are used for this purpose. To check the mesh sensitivity for the convergence of solutions, the submodeling technique [29] is employed on a small volume of the full 3D small-scale T-joint model (Figure 1a). Its ability to implement it on a small part of the model containing a much finer finite element mesh, without the need of remeshing the whole model, is what makes it more efficient than the conventional method, and is the reason why it used in this work. As presented in the authors’ former work [30], the submodel boundary must lie far enough from the high temperature and stress gradients where nodal temperature and displacement solutions for the driven nodes of the submodel will be correctly interpolated from the results of the global model. Accordingly, the dimensions of the submodel volume were 78 × 32 × 20 mm³ (Figure 1a). The submodel consisted of about 6,120 solid finite elements while the same volume of the global model had only 564 ones.

Figure 4

Finite element for small-scale T-joint model: (a) full 3D model and (b) shell/3D model

Thus, the submodel mesh had more than ten times greater density than the equivalent part of the global model. The same element types are used for both global and submodel. In the thermal analysis, the addition of weld filler is simulated by the element birth and death technique [31, 32] available in the finite element software package Abaqus. The subsequent mechanical analysis is then done without applying the element birth and death technique [33] to further accelerate the calculation speed.

3.2 Large panel model

To investigate the deformations in the large welded panel using the shell/3D modeling technique, an experimentally tested model taken from the literature [34] is considered. Figure 5 shows the geometry of the large panel model with restraint conditions and the weld leg dimension of the fillet weld. As can be seen from Figures 1 and 5, the large panel model size is six times greater than that of the small-scale T-joint model. On a horizontal plate of dimensions 1000 mm × 1500 mm, two vertical plates of dimensions 140 × 1500 mm are MAG welded without any gap between the plates in four single weld passes. Since the welding sequences are not given in [34], they are assumed as shown in Figure 5. The welding conditions, material properties and boundary conditions are the same as in the small-scale T-joint model. Additionally, the same finite element types and mesh densities are used to discretize the shell/3D model geometry of the large panel in thermal and mechanical analysis, according to the shell/3D small-scale T-joint model. The finite element mesh consisting of 48,336 solid and 12,312 shell elements of the large panel model is shown in Figure 6. All the numerical simulations in this study are performed using the Abaqus/Standard software.

Figure 5

Large panel model geometry

Figure 6

Finite element mesh of large panel model

4.1 Small-scale T-joint model

Figure 7 shows a comparison of results for the temperature profiles predicted by the full 3D and shell/3D models at the middle surface of the horizontal plate along line A-B shown in Figure 1, 150 s after the beginning of the welding process. It is obvious that the temperature profiles are almost the same along the observed path. The peak temperature reaches about 350^∘C for the full 3D model and the deviation from the shell/3D model is negligible. The temperature history of the arbitrarily selected point N1 located at the bottom surface of the horizontal plate (Figure 1a) near the weld bead for the first 300 s after the beginning of the welding process is given in Figure 8. Again, the temperatures of both models during the welding process and cooling stage to ambient temperature are almost the same with negligible deviations. Finally, it can be concluded that the use of a shell/3D technique gives almost identical spatial-temporal temperature distribution during the welding and cooling process as with the application of a full 3D model.

Figure 7

Temperature profiles for full 3D and shell/3D models 150 s after beginning of welding process at middle surface of horizontal plate along line A-B in Figure 1

Figure 8

Temperature-time history of node N1 (Figure 1a) for full 3D and shell/3D models for first 300 s after beginning of welding process

Figure 9 shows the horizontal plate deflections in y-direction at middle surface along line A-B (Figure 1) after the cooling of the welded T-joint to ambient temperature (20^∘C). The maximum deflection in y-direction (vertical deflection) of the full 3D small-scale T-joint model is about 5.4 mm for both full 3D and shell/3D models. The full field vertical deflection distributions after the cooling process to ambient temperature obtained by the full 3D model and shell/3D modeling techniques are provided in Figures 10a and 10b.

Figure 9

Horizontal plate deflection in y-direction at middle surface after welding process along line A–B in Figure 1

Figure 10

Full field vertical deflection distributions after cooling process: (a) full 3D model and (b) shell/3D model

These results are in agreement with the previous literature reports and with the experimental results obtained by the authors in [10]. As seen in the figures, almost the same results are obtained with both the shell/3D and full 3D models. Thus, it is shown that the recommendations set out by the authors of this contribution in [10] for the estimation of the minimum 3D zone size in the shell/3D model can be applied not only to the T-joint models with the same thicknesses of flange and web plates but also to the T-joint models consisting of plates with different thicknesses. This is very useful in engineering practice because the fillet welded stiffened panel structures are generally composed of longitudinal stiffeners with attached plates of different thicknesses.

The main advantage of the shell/3D model can be clearly seen in Table 3, which shows that the shell/3D model shortens the simulation time by 52% compared with the full 3D model. Additional savings in computing time for the shell/3D model can be achieved by parametrically varying the sizes of the 3D zone on the flange and web as well as the weld leg dimensions for each specific geometry of the T-joint model. However, the time required to perform the parametric analyses and define the optimal 3D zone size could be greater than the savings in CPU time in that case. The numerical simulations of the small scale T-joint are performed on a PC with an Intel® Core^™i3 processor with 8 GB of RAM memory.

Furthermore, the results of the mesh sensitivity study are presented here. Figure 11 shows a temperature profile comparison of the global model and submodel along line A-B (Figure 1) 150 s after the beginning of the welding process. As already mentioned, the submodel has the mesh of extreme density, i.e. more than ten times denser than the global mesh, so it can be concluded that the global mesh of the T-joint small-scale joint is properly designed. The same conclusion can be drawn from Figure 12,where it can be seen that the vertical deflection of the horizontal plate of the global model converges to the submodel values. A slight deviation of the results is observed only in the weld region where the large stress gradients occur.

Figure 11

Temperature profiles of global model and submodel 150 s after beginning of welding process at middle surface of horizontal plate along line A-B in Figure 1

Figure 12

Horizontal plate deflection of global model and submodel in y-direction at middle surface after welding process along line A–B in Figure 1

Finally, it can be concluded that the use of a shell/3D model has negligible influence on the temperature and deformation distribution in comparison with the full 3D model. Moreover, the mesh sensitivity analysis verified that the finite element mesh of the small-scale model is sufficiently dense and thus the mesh density can be successfully expanded on a larger-size model, i.e. a large panel model.

4.2 Large panel model

The next step is to test the shell/3D modeling technique to predict the out-of-plane deformations and residual stresses in a large panel model with two longitudinal stiffeners. The results of the mechanical analysis are presented in comparison with the experimental results from the literature [34]. Figure 13 shows the deflection profile of the base plate along line A-B shown in Figure 5. The maximum deflection computed by the shell/3D model reaches 12.6 mm which is very close to the experimental measurement and by using the ISM method. Furthermore, it is apparent that a nearly full numerically obtained curve matches well with the experimental measurements. A similar conclusion can be drawn for the vertical deflections along line C-D (Figure 5) shown in Figure 14, where the ISM method misses the experimental curve trend completely as it is not able to take the welding process history into the account, while the proposed method shows a very good matching with the experimental data on the right half of the figure. On the other hand, the partial deviation of the numerical results from the experimental measurements could be attributed to the welding parameters used in this work, mainly the order of welding passes, time between the passes and mechanical boundary conditions, which were not provided in the literature [34] where the model was taken from. Further investigation into the welding conditions is needed to address these discrepancies, which is out of the scope of this paper. Additionally, from Figures 13 and 14, it can be seen that the panel deflection profiles obtained by the ISM are strongly dependent on the input data. As they are obviously applied symmetrically, the calculated deflections are symmetric along paths A-B and C-D. This is sometimes very far from reality so that more realistic solutions can be obtained using the presented shell/3D TEP model. Generally, the peak deflections for both shell/3D and ISM methods agree well, but the trend of the shell/3D model is more realistic compared with the experimental measurements. A full field of vertical deflections can be seen in Figure 15. The numerical simulations of the welding process for a large panel model are performed on a Sun three-node cluster with 48 cores in total, 32 GB of RAM, and 500 GB of storage each.

Figure 13

Horizontal plate deflection in y-direction at middle surface after welding process along line A–B in Figure 5

Figure 14

Horizontal plate deflection in y-direction at middle surface after welding process along line C–D in Figure 5

Figure 15

Full field deflection distribution of large panel model