1   Introduction

In a hostile environment, the prediction of fatigue life of structures such as aircraft fuselages, pressure vessels, pipes, and modern ships becomes more and more important (Dake et al. 2012; Zhang et al. 2016; Alessi et al. 2019; Li et al. 2020). It is estimated that approximately 90% of the structural failures in engineering practice during the service life is related to fatigue failure (Shipley and Becker 2002; Bergara et al. 2017). Existing surface or near-surface cracks largely increase the possibility of fatigue failure of structures. In the case of existing cracks, the crack propagation determines the rest life of a structure. The growth of a crack depends on factors such as the geometrical profile of initial cracks, loading conditions, welding effects, and stress status around a crack tip or crack front (Silva et al. 2017). With regard to all these factors, SIF is usually used for the quantification of a crack severity, and the estimation of the remaining life of a structure through Paris' law (Paris et al. 1961). Therefore, computation of the SIF in a proper and effective way is crucial for the prediction of a structural life.

Typically, there are two different ways for the estimation of SIF: analytical or empirical methods (Newman Jr and Raju 1981; Bowness and Lee 2000; BSI 2015) and numerical methods (Moës et al. 1999; Shi et al. 2010; Giner et al. 2008; Sukumar et al. 2000; Singh et al. 2012). The analytical methods only exist for simple cracked geometries subjected to simple loading conditions. For instance, Newman Jr and Raju (1981) proposed an empirical solution for a simple cracked plate with relative crack depth between 0 and 0.8. Lin and Smith (1999) and Wang (2002) investigated the SIF of a simple plate with the crack depth to thickness ratio larger than 0.8.

For arbitrary-shaped cracks in complex structures, the numerical method is a reliable resource for the modeling and calculation of SIF. A considerable amount of the numerical research based on traditional FEM has been found from researchers such as Gifford Jr and Hilton 1978; Shiratori and Miyoshi 1986; Bowness and Lee 1996; Wang 2002; Cai et al. 2017; Cai et al. 2018c; 2018a, 2018b. Shiratori and Miyoshi (1986) used the superposition of the "unit distributed load" on the cracked surface to calculate the SIF under arbitrarily distributed stresses. Bowness and Lee (1996) derived the weld magnification factors to account for the SIF for weld-toe cracks in T-butt joints. Gifford Jr and Hilton (1978) simulated the SIFs around the crack in a simple plate, directly integrating the SIF into the stiffness matrix as an extra unknown. In traditional FEM, the so-called spider web mesh around the crack tip/front has to be utilized for the sake of accuracy (Newman Jr and Raju 1979; 1981). By using such method, the displacement and stresses in the vicinity of the crack tip/front are first fitted to the elastic singular solution. Then, the SIFs can be calculated as a function of the distance between the crack tip/front and the selected vicinity point. Such points should be as close as possible to the crack tip/front. Finally, the approximate values of SIF at crack tip/front are extrapolated through fitted curve when the distance is set to be 0. However, this method by using traditional FEM is cumbersome, and the mesh has to conform to the complex geometrical topology of a crack (such as a crack on a pipe surface). Highly refined element meshes are therefore needed to guarantee the accuracy of simulation. In addition, re-meshing is required for crack growth simulation, which makes using traditional FEM complex and challenge.

In order to circumvent these limitations of the traditional FEM, the extended finite element method (XFEM) was introduced by Belytschko and Black (1999) based on the concept of unity partition by Melenk and Babuška (1996). It is an extension of the traditional FEM, with an essential updating of adding enrichment functions for the approximation of crack discontinuity and singularity. Therefore, the presence of a crack is ensured by the special enriched functions in conjunction with additional degrees of freedom. This technique allows the entire crack to be represented independently of the mesh. The accuracy of crack simulation can be guaranteed through a relative coarse mesh of a structure. No re-meshing is needed for the crack growth simulation. So far, the XFEM has been widely used to solve the problems of fracture mechanics. For instance, Shi et al. (2010) used the XFEM to predict the fatigue crack growth in a complex helicopter structure. Giner et al. (2008) applied the XFEM to fretting fatigue problems. The use of XFEM provided a relative accurate SIF calculation based on relative coarse meshes. Few modeling efforts are required to integrate a crack to complex structural models. Other relevant researches that have been found include the problem of crack growth with frictional contact (Dolbow et al. 2001), quasi-static crack growth (Sukumar and Prévost 2003), cracks in shells (Areias and Belytschko 2005b; da Silva et al. 2019), stationary and growing crack (Ventura et al. 2003), and 3D crack propagation (Areias and Belytschko 2005a).

In spite of the fact that engineering standards try to provide accurate SIF solutions, inadequate estimations are generally produced in engineering practice. For instance, in BS 7910 (BSI 2015), an abnormal phenomenon of SIF in a cracked T-butt joint accounting for welding effect has been observed during our former research work. Compared with BS 7910, the calculation results of SIF at the surface points of welded specimens cannot be well predicted, with a large discrepancy appearing. Since these empirical formulas of SIF in BS 7910 are originally based on results of traditional FEM, it is reasonable to doubt the confident utilization of traditional FEM for the computation of the SIF at such surface points of a cracked welded specimen. In order to explore such abnormal increase at the surface points of cracked welded specimens, the objective of this paper is to numerically investigate the SIF of a cracked T-butt joint using different methods including BS 7910, XFEM, and FEM, revealing the possible causes of such discrepancies through parametric study. Hence, the research work in this paper will facilitate the improvement of the analytical prediction of SIF in industry design standards BS 7910.

In order to obtain confidence of further simulation, the comparison work based on a single cracked plate without welding will be firstly conducted. Afterwards, investigation of the SIF in a cracked T-butt joint with welding which is widely used in engineering domain is deployed. Comparisons in term of SIF among analytical predictions in BS 7910, traditional FEM, and XFEM are carried out.

2   Methodologies and Models of SIF Calculation

In this section, the methodologies including analytical solutions and numerical methods for calculation of SIF are described. Specifically, the general geometrical descriptions of both a simple plate and a T-butt joint have been presented in Section 2.1. The analytical solutions for the estimation of SIF in these types of structures are presented in Section 2.2, while the numerical models through both traditional FEM and XFEM are presented in Section 2.3. The partition strategy of structures and the mesh of the cracked region are described. A convergence analysis for the selection of contours in XFEM is performed. It should be noted that only J-integral is used to calculate the SIF of cracks during simulation in this paper. Alternatively, the numerical method in terms of VCCT (virtual crack closure technique) can be used to evaluate the SIF of cracks, especially for the structures under multiple loads (Silva et al. 2017). However, we do not take it into account in current research.

2.1   Geometrical Description

Two typical types of structures are used in this investigation. One is a simple plate, and the other one is a T-butt joint, as illustrated in Figure 1. The general geometrical information of the plate is 2h in plate length, 2b in plate width, and T in thickness. The crack is located at the central of the plate with the half-crack length of c, the crack depth of a, and the crack angle of ϕ. For the T-butt joint, h_att is the height of stiffener. In this paper, the variation of the crack dimension is listed in Table 1, with a fixed thickness T of 16 mm. The specific geometrical variation for further parametric study is based on the research by Bowness and Lee (2000). The crack dimension varies, with crack depth (a) changing from 0.005T to 0.9T, and crack half-length (c) changing from a to 10a. In addition, the half-width of the structures (b) is set to changeable, equal to 5T when the crack length c is less than the plate thickness T, but otherwise is set to the length of 5c. A half-elliptical partial-through surface crack is put at the center of both structures. Only a transverse stiffener in the T-butt joint is adopted due to the lack of data.

Figure  1  The sketch of two typical types of structures used in this research

Download: Full-Size Img

2.2   Analytical Formulas

The British Standard 7910 (BSI 2015) suggests an empirical method to assess SIF and the crack growth of cracked structures. For the calculation of SIF on a simple finite plate with a center-crack subjected to either tension or bending loads, Eq. (1) proposed by Newman Jr and Raju (1981) is used. The equation is expressed as the function of the crack angle (ϕ), the crack depth (a), the crack length (2c), the plate thickness (T), and the plate width (2b). The crack shape is semi-elliptical shape, as shown in Figure 1. Two types of loads may be applied to the surface-cracked plate: remote uniform tension (σ_app) and/or remote bending (S_b).

where the application domains are $ 0 \leqslant a/c \leqslant 1.0 $, $ 0 \leqslant a/T \lt 0.8 $, 2c/b <  0.5, and $ 0 \leqslant \phi \leqslant \pi $. Q is the shape factor, expressed as Eq. (2). Function F is the tension factor, expressed as Eq. (3). S_b is the remote outer-fiber bending stress, which can be calculated by S_b = 6M/bT².

where

The function f_ϕ is an angular function which is used for the correction of the crack tip plasticity, assuming a very small plastic crack tip.

The correction factor f_w for a finite-width plate is:

The function H has the form of Eq. (10), which is developed by curve fitting and engineering judgment. The expressions of parameters H₁, H₂, and p are not listed here for clarity reason, referring to Newman Jr and Raju (1981).

For the analytical solution of T-butt joint with a transverse stiffener, the used empirical formula is proposed by Bowness and Lee (2000) based on a parametric study on FEM simulation. A weld magnification factor (M_k) is used to take into account the effect of welding and the geometrical factors during the calculation of SIF, as expressed in Eq. (11). Where Y_attachment is the geometrical factor which accounts for the presence of welding and attachment; Y_noattachment is the geometrical factor of a cracked simple plate without attachment. The specific formula is expressed as Eq. (12) for the deepest point in a crack under membrane loading. Equation (13) is for the crack surface point under membrane loading. The application domain of these formulas is $ 0.005 \leqslant a/T \leqslant 0.9 $. For the specific expressions of f₁, f₂, and f₃, they are from Table 7 and Table 9 in the paper of Bowness and Lee (2000). For clarity reason, they are not presented here.

For the comparison of FEA results in Section 3, normalized values are used. For the simple plate, the calculated SIF is normalized by K_ref ($ K_{\text {ref}}=\sigma _{\text {app}} \sqrt {\pi a/Q} $). Where Q is expressed as Eq. (2), and σ_app is the remote uniform tension stress. For the T-butt joint, only the weld magnification factor is used for comparison, which is expressed as Eq. (11). The geometry factor Y (both Y_attachment and Y_noattachment) is calculated by $ Y=K/(\sigma _{\text {app}} \sqrt {\pi a}) $ in the structures with and/or without attachments. It should be noted that J-integral is, generally, first provided during a FEA simulation. Under this situation, the SIF is expressed as $ K=\sqrt {JE/(1-\nu ^{2})} $ for plane strain condition (the deepest point of a crack in this paper), while $ K=\sqrt {JE} $ for plane stress condition (the surface point of a crack in this paper). Where J is the calculated J-integral, E is the Young modulus, and ν is the Poisson ratio.

2.3   Numerical Models in FEM

The numerical models for the calculation of SIF through traditional FEM have been developed through ABAQUS (Simulia 2016), as seen in Figures 2 and 3. Only the structure of a simple plate is deployed. The mesh strategy is the same as the work by Newman Jr and Raju(1979, 1981). Figure 2 illustrates the geometrical partition and the general distribution of mesh seed in the plate. A bias seed strategy is then used to guarantee a refined mesh of the cracked region. The arrow is the direction of seed increase at the respective structural edge. Figure 3 shows the mesh with a "spider web" in the cracked region in detail. The utilization of "spider web" is to fashion uniform integral paths for J-integral. At least four layers of element are required around the crack front for the accuracy of simulation according to the research of Cai et al. (2017). The SIF will be derived through J-integral in simulation by software. A special collapsed element is used in the crack front to account for the crack singularity.

Figure  2  Geometrical partition and the distribution of mesh seed in a cracked plate (half model displayed)

Download: Full-Size Img

Figure  3  Mesh of a plate with a "spider web" mesh in the cracked region for SIF calculation by traditional FEA

Download: Full-Size Img

Python script is used to develop models for parametric study, accounting for the change of crack parameters. The variation of crack parameters is listed in Table 1. The isotropic steel is used with the Young modulus of E = 2.1 × 10⁵ MPa and the Poisson ratio of ν = 0.3. The three-dimensional element C3D8R is utilized. Only a quarter of a plate is used due to the symmetry of structure. The symmetrical boundary conditions in both X plane and Y plane are used, as shown in Figure 3. A specific constraint on Z displacement is used for the vertex of the cracked region. A distributive tensile force (1 MPa) is exerted on the end surface of the plate in order to compare with data from Wang (2002) and Newman Jr and Raju (1979).

2.4   Numerical Models in XFEM

Structures in terms of a simple plate and a T-butt joint are used for the modeling and calculation of SIF through XFEM. The numerical models are shown in Figure 4. As explained in Section 1, there is no special mesh strategy required for the XFEM model of a simple plate. General partition of geometry is the same as the one in Section 2.3.

Figure  4  Mesh of structures for SIF calculation by XFEM

Download: Full-Size Img

Due to the complexity of a T-butt structure, a special partition of geometry is proposed for the cracked region. A partition region in terms of a rectangular box (W × H × L in Figure 4) is utilized in the weld-toe area in order to produce a regular mesh around the crack. Equation (14) shows the dimension of the varied box according to crack parameters, where X is the number of contours (a is the crack depth, c is the half-crack length, I is a parameter to make sure the crack is always inside the specified box, and T is the thickness of plate. seed_crack_a and seed_crack_c indicate the characteristic mesh sizes at H edge and L edge, respectively. The applied tension is σ_app = 1 MPa on the surface of structures along Y direction. The surface X = 0 is set as a symmetry boundary. The other boundary conditions are the same with the models described in Section 2.3.

In order to obtain an accurate result through XFEM simulation, a convergence analysis in terms of both the enriched radius and the number of contours (X) in the enriched region is conducted in this section. The number of the contours in the enriched region will affect the calculation accuracy of crack parameters such as J-integral and SIF. The enriched region (or the enriched radius) is set by the utilization of the characteristic mesh size (l_c) of structures, which is equal to seed_Crack_a. In this research, it is set to l_c, 2l_c, 3l_c, 6l_c, 10l_c, 20l_c, and 40l_c, respectively. The specific dimension of the plate is designated as a/T = 0.3, a/c = 0.4, T = 16 mm, and b = h = 80 mm. A simple tension loading condition (1 MPa) is exerted. The empirical results are calculated through Eq. (1) proposed by Newman Jr and Raju (1981), which is explained in the analytical part in Section 2.2. As a result, the calculated SIF values are $ K_{Ia} =1.18 M\text {Pam}^{1/2} $ at the deepest point of a crack and $ K_{Ic} = 0.84 M\text {Pam}^{1/2} $ at the surface point of a crack, respectively.

According to the sensitivity study, we select 8 contours (X= 8), and select the length l_c for the enriched region for the following simulation. In order to further investigate the calculation accuracy of SIF, a comparison of SIF value from each contour layer is carried out. Figure 5 shows the simulation results of SIF on a T-butt joint. Here, X is set to 8, while enriched region is set to l_c. Eight different contours have been compared with each other. For the uncracked ligament in every specific crack front point (different crack angles as shown in Figure 1, contours 2 to 8), the K_I value has a very small variation. It reflects the path independence characteristic of SIF around a crack. However, we found that the K_I value along the crack front (contour 1) has increased a lot. This is the so-called oscillation behavior. Hence, the simulation result from one contour may not be reliable. The corresponding results should be intentionally removed. In addition, there is a sudden drop of K_I at ϕ = 90^∘, which is caused by an extrapolation error in software. Therefore, in order to obtain accurate simulation results in the following research, the crack surface data is calculated at ϕ = 87^∘, instead of ϕ = 90^∘.

Figure  5  A comparison of K_I with respect to the point along crack front ($ \phi \subseteq [0,90^{\circ }] $) from different layers of the contours (T-joint, X= 8)

Download: Full-Size Img

3   Results and Discussions

3.1   Comparison of a Simple Plate

In this section, the comparisons in terms of normalized SIF are made among the simulation results from traditional FEM, the simulation results from XFEM, and the results from empirical predictions (Section 2.2) for a simple cracked plate without welding effect.

Figures 6 and 7 show the comparison results in terms of normalized SIF between simulation results by FEM and by empirical predictions. Since the application domain of Eq. (1) is 0 ≤ a/T <  0.8, the known simulation results from Wang (2002) are used for comparison at a/T >  0.8 (a/T = 0.9 and 0.95). A relative small difference (1.78%) has been found for the SIF at the deepest point of a crack, while the average difference is 4.93% for the surface point of a crack. A possible cause of such discrepancy is the different numbers of elements. In the original numerical model of Newman Jr and Raju (1981), a maximum of 4000 elements were used, whereas the maximum element in the model of this paper is 240 000. Hence, a good agreement has been obtained, which provides the confidence to utilize FEM for the prediction of SIF for a cracked plate without welding effect.

Figure  6  Comparison diagrams of the normalized SIF with respect to a/T on a simple plate between traditional FEM and empirical formulas: SIF at the deepest point of a crack

Download: Full-Size Img

Figure  7  Comparison diagrams of the normalized SIF with respect to a/T on a simple plate between traditional FEM and empirical formulas: SIF at the surface point of a crack

Download: Full-Size Img

The SIF on the same cracked plate is also calculated based on XFEM. Figures 8 and 9 show the comparison results in terms of the normalized SIF between simulation results by XFEM and the results by empirical predictions. Likewise, for a/T >  0.8 (a/T = 0.9 and 0.95), the known simulation results from Wang (2002) are used for comparison due to the limited application domain of the formulas from Newman Jr and Raju (1981). The average difference for the deepest crack point is 5.12%, while the difference for the surface crack point is 5.83%. Compared with the analytical solutions, the SIF discrepancies introduced by both FEM and XFEM are comparable. Accounting for the factor that there is no need for special mesh refinement in the crack region, it can be concluded that XFEM is a good alternative for SIF prediction of cracked structures without welding effect. Such mutual verification indicates that the empirical formulas used for SIF of a cracked simple plate in BS 7910 are still satisfied with current engineering practice. Sufficient confidences are obtained to conduct the following investigation on cracked T-butt joints through the same simulation strategy.

Figure  8  Comparison diagrams of the normalized SIF with respect to a/T on a simple plate between XFEM and empirical formulas: SIF at the deepest point of a crack

Download: Full-Size Img

Figure  9  Comparison diagrams of the normalized SIF with respect to a/T on a simple plate between XFEM and empirical formulas: SIF at the surface point of a crack

Download: Full-Size Img

3.2   Comparison of a T-Butt Joint

In this section, the crack SIF of a T-butt joint is calculated by both XFEM and empirical formulas. The result in terms of weld-toe magnification factors (M_k) for semi-elliptical cracks (Figure 1) in T-butt joints is compared. Empirical formulas are from Bowness and Lee (2000), which have been further included in the British Standard BS 7910. Only half model is considered during simulation. The mesh strategy, applied force, and boundary conditions are described in Section 2. During simulation, the a/T ratio varies from 0.04 to 0.9. The corresponding a/c ratio varies by 0.1, 0.2, 0.4, 0.7, and 1.0.

Figure 10 shows a typical "butterfly" pattern (Schijve 2001) of the von Mises stress distribution in the crack tips (at both the deepest crack point and the surface point). Figures 11 and 12 show the comparison results in terms of the weld magnification factors (M_k) at both the deepest crack point and the surface crack point on a T-butt joint by XFEM. As shown in Figure 11, it is obvious that there is a good agreement between XFEM simulation and the empirical predictions from Bowness and Lee (2000) at the deepest crack point, i.e., an average difference of 2.98%. However, it is found that a relative large difference, say 27.31%, has appeared at the surface crack point (as seen in Figure 12). The simulation values from XFEM are much lower than the values obtained from Bowness and Lee (1996, 2000). Since the predictions of Bowness and Lee (2000) are majorly based on the data from traditional FEM, it is reasonable to consider that such a large discrepancy is similar to the discrepancy between FEM and XFEM for the computing of SIF accounting for the welding effect.

Figure  10  The typical "butterfly" pattern of the von Mises stress distribution in the crack tip regions

Download: Full-Size Img

Figure  11  Comparison diagrams of the normalized SIF with respect to a/T on a T-butt joint between XFEM and empirical formulas: SIF at the deepest point of a crack

Download: Full-Size Img

Figure  12  Comparison diagrams of the normalized SIF with respect to a/T on a T-butt joint between XFEM and empirical formulas: SIF at the surface point of a crack

Download: Full-Size Img

In order to clarify the possible causes of the large discrepancy of M_k at the surface crack point, an extra discussion is carried out in this section. As mentioned above, the empirical formulas were actually derived from a series of numerical simulation by traditional FEM. Due to the limitation of computation capacity at that time, the number of elements of their numerical models was restricted to 2500. Hence, its calculation accuracy may be doubted.

Another possible reason is the irregular distribution of the von Mises stress with an exceptional increase around the surface crack tip by traditional FEM. As described in Section 2.2 for the calculation of M_k, the accuracy of stress distribution determines the accuracy of SIF which is the determined factor of M_k. As demonstrated from the research of Zettlemoyer and Fisher (1977), a constant stress gradient exists in the weld-toe region of a welded stiffener. Hence, we also re-simulated the SIF of a T-butt joint by traditional FEM with the "spider web" strategy described in Section 2.3. It is found that a high-varied stress gradient exists in the weld-toe region with a much higher stress concentration (2.363 MPa) at the surface crack point through traditional FEM, as seen from the top image in Figure 13. Instead, a smooth and constant von Mises stress is produced at the same place by XFEM. The simulated value (1.268 MPa) is much lower than the one by FEM, as seen from the bottom image in Figure 13. Therefore, we consider that the "spider web" strategy in FEM is the major reason of such discrepancy in Figure 12, which does not properly produce such constant stress gradient along the weld-toe. The XFEM can be used for more accurate prediction of M_k of a T-butt joint. It is therefore suggested to update the empirical formulas proposed by Bowness and Lee (2000).

Figure  13  von Mises stress (range between 0 and 2.36 MPa with the color varying from blue to red) distribution around the weld-toe in T-butt joint: simulation results by traditional FEM (top image) and results by XFEM (bottom image)

Download: Full-Size Img

4   Conclusions

This paper aims to explore the discrepancy of SIF at the surface points of the welded T-butt joints with semi-elliptical cracks, and its possible causes. A numerical investigation of SIF is therefore performed, deploying both XFEM and FEM. Results from both numerical simulation and empirical predictions are used for comparison. Two types of structures are used, including a simple cracked plate without welding and a T-butt joint with welding. A series of parametric studies with respect to the variation of surface crack parameters (a/T and a/c) are carried out in order to compare SIF and/or M_k. The conclusions are presented as follows:

1) Compared with the empirical solutions in BS 7910, the computing accuracies of SIF through both FEM and XFEM are comparable in this paper. The XFEM is a good alternative for SIF prediction of a simple cracked plate and a cracked T-butt joint.

2) The mutual comparison results indicate that the empirical formulas used for the calculation of SIF in a cracked simple plate without welding in BS 7910 are still satisfied with current engineering utilize.

3) For the prediction of SIF in a cracked T-butt joint with a transverse stiffener, there is a good agreement between XFEM and empirical predictions at the deepest crack point (a difference of 2.98% in terms of weld magnification factor (M_k)). However, a large discrepancy (as large as 27.31%) appears at the crack surface point.

4) In FEM, the spider web strategy has affected the constant stress gradient in the weld-toe, and further produced large von Mises concentration at crack surface points. This phenomenon causes the large SIF/M_k discrepancy of a T-butt joint with a transverse stiffener. Accounting for welding effect, XFEM has effectively overcome this abnormal phenomenon since a smooth and constant von Mises stress gradient is successfully simulated.

5) The research in this paper indicates a possible inaccuracy estimation of SIF at the crack surface point in a T-butt joint used by BS 7910. As a result, a potential updating of the relevant empirical formulas in BS 7910 is recommended for a better prediction in engineering practice.