High-Cycle Fatigue Crack Growth in T-Shaped Tubular Joints Based on Extended Finite Element Method

Abstract

High fatigue load, which exists widely in steel building structures, likely leads to brittle failure at the joints, supports, and so on. This can lead to the partial or total damage of the structure and even to cause the collapse of the whole structure. This article aims to provide a method to simulate high-cycle crack propagation in tubular joints, which is one of the most common types occurring in steel structures. Firstly, sixteen T-shaped tubular joint models under different load conditions and initial crack dimensions were built through the coordinate mapping method. Secondly, based on the extended finite element method (XFEM), an algorithm was developed by combining the secondary development in Abaqus and a quasistatic simulation method to simulate high-cycle crack growth in tubular joints under a constant amplitude. The results of the simulations were compared with experimental data. The study found that the surface stress calculated from the tubular joint models using the coordinate mapping method was close to the experimental data. Through the comparison of the crack propagation rate and the crack growth process between the simulation and experiment results, the simulation method was validated. When a crack penetrated the tube wall, the difference in the load cycles between the simulations and the experiment was 9.5%. The initial crack dimension had an impact on the crack propagation, with the decrease in the a/c and KII generally becoming the dominant factor with respect to the crack growth, while the fatigue life of the joints tended to increase.

1. Introduction

Due to their good mechanical and construction performance, tubular joints are widely used in large-span spatial structures, high-rise buildings, and many other engineering fields. Because of the sudden changes in geometry happening at welds, initial defects, and residual stresses points caused by welding, the welds become the weak points in tubular joints [1].

For steel tube structures in service under dynamic loads, such as wind loads, equipment vibrations, etc., stress concentrations that happen at the welds due to geometric discontinuity and initial defects can induce the initiation and propagation of fatigue cracks. When fatigue stress reaches its limit, the unstable propagation of cracks and brittle failure will occur, which may cause partial structural damage or integral destruction. Several examples can prove this point. For a high-rise building, it is likely that the curtain wall supporting structures will fall off under wind-induced fatigue loads, thereby causing casualties [2]. The same things take place in steel structures, where wind-induced loads may lead to fatigue damage at the beam-to-column connections, thus reducing the structure safety [3]. Equipment vibration, such as suspension cranes, can also generate high-cycle loading in large-span grid structures, thereby making the bolt joint a weak point [4,5]. When fatigue stress reaches its limit, the unstable propagation of cracks and brittle failure will occur, which may cause partial structural damage or integral destruction. In 1980, a brace of Alexander L. Kielland’s offshore platform suffered fatigue failure, thereby making the whole structure collapse in 15 min. A total of 123 people died in this accident [6]. Thus, in order to avoid fatigue failure accidents and ensure structural safety, a study of high-cycle fatigue crack propagation in structures, especially at the joints, needs to be conducted.

Currently, there are three main methods used in the study of high-cycle fatigue behavior. The first one is the S-N curve method, which is mainly used in fatigue designs [7]. In the second method, a predictable model is established to reveal how the fatigue damage is accumulated and to predict the behavior of the fatigue crack growth [8]. And the third one is the fracture mechanics method, which is mostly used to predict the residual fatigue life of tubular joints with initial cracks [9].

The high-cycle fatigue life of tubular joints can be divided into four stages from N₁ to N₄. N₁ is the recorded fatigue life of initial cracks using any detection techniques. N₂ represents the fatigue life of the first visual crack. N₃ denotes the fatigue life of the cracks that penetrate the wall thickness of the chord tubes. N₄ is the fatigue life at the point of component failure [10]. Typically, N₃ or N₄ have been used to describe the fatigue life of tubular joints. With the hot pot stress amplitude obtained using tests or numerical simulation, the S-N curve can be drawn. The S-N curve can show the influence of different dimensionless geometric parameters on the fatigue life of tubular joints [11]. However, it is hard to distinguish between the crack initiation and crack propagation [12] using the S-N curve; thus, the S-N curve cannot be used to study the process of fatigue crack propagation. A predictable model is based on probability statistics and the theory of fracture mechanics, thereby predicting the structural reliability at the theoretical level. But, it also cannot present the process of crack growth.

The fracture mechanics method was established based on the Paris law [13]. According to this method, the process of crack propagation from its initial stage to its critical stage can be shown, and the curves of the fatigue crack propagation rate and residual fatigue life can be calculated [14]. The drawback of this method mainly lies in how to develop a precise model to ensure the accuracy of the stress intensity factor (SIF) at the crack front.

In the early years, the shell element was used in tube modeling, and the linear spring was used to simulate cracks [15], or the tube model was established using the shell element, while the weld was set up using solid elements. Bownees and Lee [16] found that shell elements and linear springs failed to consider the size effect of tubular joints and cracks; as a result, it could only be used when the ratio of the crack depth to the tube thickness ranged from 0.35 to 0.8. Cao et al. [17] reported that simulation results acquired from their model of tubular joints with shell elements and solid elements were less accurate than their model using solid elements only. Thus, in order to acquire more accurate results, using solid elements for modeling has become the broad consensus. However, weld solid modeling and initial cracks are difficult to derive due to the limit of finite element software (FE software (Abaqus/CAE V6.1)). Therefore, Lie et al. [18] proposed a coordinate mapping method. In this method, an FE mesh was initially generated at the plate component; then, the transition from the plate component to the tube component was conducted using coordinate mapping. Based on this transition, the tube and initial cracks could be modeled and meshed more accurately. The literature [18,19,20] has used this coordinate mapping method to establish a numerical model of T-shaped and Y-shaped tubular joints with initial cracks.

With the development of computer technology, the numerical study of high-cycle fatigue propagation has transformed from the SIF in a single stage to the whole process of crack propagation. Fatigue crack propagation criteria have been programmed in FE software, but they are only applicable to simulate low-cycle fatigue crack propagation due to the correlation between the crack growth rate and the energy release rate at the crack tip [21], which is hard to calculate for high-cycle fatigue crack propagation simulations. Some scholars use the secondary development function of FE software to develop packages that realize high-cycle simulation. Yagi et al. [22] developed a fatigue crack propagation algorithm, and Adhikari [23] used Abaqus (CAE V6.1) and Fracn3D (V7.5.5) to model the process of high-cycle fatigue crack growth in tubular joints. These simulations were based on the quasistatic simulation method, in which the crack propagation process was divided into several steps. The SIF at the crack front was obtained first. The space location of the crack was then calculated according to fatigue crack propagation criteria [24]. The reason for using the quasistatic simulation method is that the secondary development can only develop the software’s function, but it cannot change its operation logic. Taking Abaqus as an example, the dynamic simulation method combined with fatigue growth needs the damage parameters of the material, which can only be used in ultra-low cycle or low-cycle fatigue problems [25]. When using the finite element method (FEM), remeshing must be conducted at each step. With an increase in the crack dimension, the mesh of the crack becomes more complicated to derive, and it is harder to update the crack model due to huge computing costs. Yang et al. [9] replaced the FEM with the boundary element method to avoid remeshing the whole component, but it was inevitably necessary to remesh the elements at the crack fronts, the intersections of the cracks, and the outer surfaces of the tubular joints.

In 1999, Belytschko [26] provided an extended finite element method (XFEM) to simulate crack propagation. The XFEM retains the advantages of the FEM and overcomes its drawbacks, without requiring remeshing for the crack propagation and high-density mesh at the crack front. The XFEM can also be used in secondary development techniques. Nikfam et al. [27] simulated the high-cycle crack propagation of a T-shaped butted joint via the XFEM and secondary development in Abaqus. Dirik et al. [28] developed an algorithm based on the XFEM to simulate the overloaded retardation effect of a plate component.

However, at present, few studies take the XFEM into consideration when simulating the high-cycle fatigue crack propagation of tubular joints. Therefore, this paper aims to combine the XFEM and high-cycle fatigue crack propagation to study the fatigue crack growth in T-shaped tubular joints via Abaqus and secondary development. This paper is organized as follows: Section 2 briefly introduces the basic theory of the XFEM. Section 3 develops two algorithms—one to model T-shaped tubular joints and other to realize high-cycle fatigue crack propagation. Section 4, based on the tests given in the literature [22], validates the accuracy of the simulation method. In Section 5, based on the algorithms developed, a parametric study is conducted to obtain the influence of the initial crack dimension on the high-cycle fatigue crack propagation.

2. Basic Theory of XFEM

The XFEM enriches the displacement field of the FEM to solve the weak discontinuities represented by inclusions and the strong discontinuities represented by cracks. Equation (1) is the general formula of the XFEM.

$$u^h\left(x\right)=u^{\mathrm{FE}}\left(x\right)+u^{\mathrm{enr}}\left(x\right)={\sum }_{i=1}^n{N}_i\left(x\right)u_i+{\sum }_{j=1}^p{N}_j\left(x\right)\phi \left(x\right)q_j$$

(1)

where u^FE(x) is the displacement field for the FEM, u^enr(x) is the enrichment displacement field for the XFEM, N_i(x) and N_j(x) are shape functions for nodes i and j, respectively, u_i is the displacement field for the FEM nodes, φ(x) is the enrichment function (which is changed for different problems), q_j is the displacement field for the XFEM nodes, and n and p denote the number of nodes for the FEM and XFEM, respectively.

For crack propagation problems, based on the state of the different elements, the XFEM divides all of the elements into three sections: a conventional continuous section, a crack penetration section, and a crack tip section. Then, Equation (1) can be changed into Equation (2).

$$u_{\mathrm{xfem}}\left(x\right)=u^{\mathrm{FE}}\left(x\right)+u^{\mathrm{H}}\left(x\right)+u^{\mathrm{tip}}\left(x\right)$$

(2)

where, u^H(x) is the approximate displacement field for the elements being penetrated by cracks, u^tip(x) is the approximate displacement field for the elements which have a crack tip.

In order to precisely describe the displacement fields, two enrichment functions are introduced, which are F(x) and H(x) (See Figure 1). The former is used to solve the stress singularity at the crack tip, and the latter is used to solve the jump in the displacement field across the crack surfaces in the fully cut elements. Equation (3) is the displacement field for F(x) and H(x).

$$u_{\mathrm{xfem}}\left(x\right)={\displaystyle \sum _{i\in I}{N}_i\left(x\right)}{u}_i+{\displaystyle \sum _{j\in I_{\mathrm{step}}}{N}_j\left(x\right)}H\left(x\right)a_j+{\displaystyle \sum _{k\in I_{\mathrm{tip}}}{N}_k\left(x\right)}F\left(x\right)b_k$$

(3)

where, I, I_step, and I_tip represent the sets of conventional continuous elements, elements enriched by H(x), and elements enriched by F(x), respectively. a_j and b_k denote the displacement fields of the elements enriched by H(x) and F(x), respectively. N_j(x) and N_k(x) are the shape functions of node j and k, respectively.

Equation (4) is the formula of H(x) [29].

$$H\left(x\right)=\left\{\begin{array}{l}1\hspace{1em}\hspace{1em}\mathrm{if} \left(x-x^{*}\right)·n\ge 0\\ -1\hspace{1em} \mathrm{other} \mathrm{conditions}\end{array}\right.$$

(4)

where x denotes the integral point, x* is the point that is closest to x at the crack surface, and n is the normal vector of x*.

Equation (5) is the formula of F(x), which applies to isotropic material [30].

$$F\left(x\right)=\left[\sqrt{\rho }\sin \frac{\psi }{2}\sqrt{\rho }\cos \frac{\psi }{2}\sqrt{\rho }\sin \psi \sin \frac{\psi }{2}\sqrt{\rho }\sin \psi \cos \frac{\psi }{2}\right]$$

(5)

where ρ and ψ are the local polar coordinates at the crack tip.

The most obvious difference between the FEM and the XFEM is that the XFEM does not need to remesh, which is realized using a level set function. The enrichment functions of the XFEM are functions of the level set function. Because of this, when the crack grows, updating the level set function can ensure that the enrichment function and enrichment region are simultaneously updated [31]. Two level set functions are used to describe the location of the crack front. The first one is $\varphi \left(x,t\right)$, which describes the crack surface, where x represents the mesh nodes in the model, and t denotes a point in time. When a point is above the crack surface, the value of $\varphi \left(x,t\right)$ is greater than zero and vice versa (See Figure 1a). The second one is Φ(x,t), which describes an orthogonal surface, where x and t have the same meaning as in $\varphi \left(x,t\right)$ (See Figure 1b). Equations (6) and (7) are the formulas of these two level set functions.

$$\varphi \left(x,t\right)=\left(x-x_{\mathrm{tip}}\right)·t$$

(6)

$$\Phi \left(x,t\right)=\pm \min ‖x-x_i‖$$

(7)

where x_tip is the coordinate of the crack tip, and x_i is the point on the crack path.

3. Algorithm of High-Cycle Fatigue Crack Propagation

3.1. Procedure

The high-cycle fatigue crack propagation algorithm was developed in Python and is operated using the command line interface in Abaqus. The quasistatic numerical simulation method was used to establish the framework, and the results are shown in Abaqus/CAE. The algorithm has 7 steps:

Step 1: Establish the component with the proper mesh size and use the coordinate mapping method to create the T-shaped tubular joint with its initial cracks. Then, insert it into a new .cae file.

Step 2: Define the material property, boundary, and loading conditions, as well as the geometry of discontinuity, which refers to the initial crack and the region the crack may be across. After defining these parameters, submit the job.

Step 3: Obtain the SIF and coordinates of the crack front from the .odb file. Then, calculate the crack length increment based on the fatigue crack propagation criteria.

Step 4: Judge whether the crack length reaches the critical length. Then, judge whether the SIF at the crack front reaches the fracture toughness. As long as one shows as true, the loop is terminated.

Step 5: Update the crack model. Insert the updated crack model into the file created in Step 1 to replace the previous crack.

Step 6: Redefine the geometry of discontinuity and then submit the job.

Repeat Step 2 to Step 6 until the loop ends. Figure 2 is the flowchart of the algorithm.

3.2. Establishment of T-Shaped Tubular Joint with Initial Crack

Figure 3 shows the principle of the coordinate mapping method. In Figure 3, x–y–z define the local coordinate system of the tube, where the z axis overlaps with the axis of the tube and is perpendicular to the x–y plane. X–Y–Z define the local coordinate system of the plate where the X axis overlaps with the x axis and Y axis, and the Z axis is parallel to the y axis and z axis. R_cir denotes the outer diameter of the tube. t is the wall thickness. In this method, a plate model with the proper mesh is initially built in the X–Y–Z coordinate system, then point A is mapped on the plate to the point a on the tube in the x–y–z coordinate system via Equation (8). For the other points, we used the same method until the plate model was transformed into the tube model.

$$\begin{array}{l}x=s\cos \left(Y/s\right)\\ y=s\sin \left(Y/s\right)\\ z=Z\end{array}$$

(8)

where s is a parameter that is related to the position of the nodes in the X–Y–Z coordinate system. For the nodes on the upper surface of the plate, s is equal to R_cir. For the nodes on the lower surface of the plate, s is equal to R_cir-t. For those nodes between the upper and lower surfaces, the value s needs linear interpolation.

For a T-shaped tubular joint, a specific process for mapping the coordinate method is as follows: Firstly, the tubular joint was divided into five parts, including the chord, brace and weld, intersections of the chord and brace, crack, and the fixed end, which were named as Part 1 to Part 5. Secondly, we built a plate model of Part 1 to 4, and we then transformed them into a tube or curved plate model using the mapping coordinate method. Thirdly, we assembled Part 1 to 5. The above modeling process is shown in Figure 4.

3.3. Calculation of SIF and Crack Growth Increment

The calculation of the SIF in Abaqus was realized by using the interaction integral method, through which the mode-I, -II, and -III SIFs could be indirectly acquired using the J integral (See Equation (9)) [32].

$${\left[K_{\mathrm{I}},K_{\mathrm{II}},K_{\mathrm{III}}\right]}^{\mathrm{T}}=4\mathsf{\pi }B{\left[J_{\mathrm{int}}^{\mathrm{I}},J_{\mathrm{int}}^{\mathrm{II}},J_{\mathrm{int}}^{\mathrm{III}}\right]}^{\mathrm{T}}$$

(9)

where, K_I, K_II, and K_III represent the mode-I, -II, and -III SIFs, respectively, which denote the interaction J integral corresponding to the mode-I, -II, and -III cracks, respectively. B is a logarithmic energy system matrix.

When acquiring the three types of the SIF, the equivalent stress intensity factor can be calculated using Equation (10).

$$\Delta K_{\mathrm{eq}}^2={\left(\Delta K_I+k\left|\Delta K_{\mathrm{III}}\right|\right)}^2+2\Delta K_{\mathrm{II}}^2$$

(10)

where ΔK_I, ΔK_II, and ΔK_III denote the increments of the mode-I, -II, and -III SIFs, respectively. k is an empirical parameter, which is 1.0.

In order to calculate the SIF of the crack front, Abaqus discretizes it into several points by using mesh cutting. Then, the distribution of the SIF at the crack front is characterized by the SIFs of these points (See Figure 5). The same method is also applied when calculating the increments of the crack front. By calculating each point’s increment, the increments of the crack front can be achieved.

Figure 6 shows a local coordinate system established with a certain discrete point (named M) at the crack front as the original point. In Figure 6, x–y–z is a global coordinate, where the x–y plane overlaps with the lower surface of the crack. X–Y–Z define the local coordinate system set up at the point M where the X axis and Y axis overlap with the normal and tangent direction of the point, respectively. Two lines across point M are set up. One line is parallel to the x axis and is named x*, and the other line is parallel to the y axis and is named y*. The angle between X axis and line y* is α. Δa is the crack growth increment. θ is the direction of crack growth.

Δa is calculated using the Paris law. The other points’ increments can be calculated through Equation (11).

The maximum crack growth increment Δa_max is predefined to determine the increment for the point with a maximum ΔK_eq at the crack front in a single step. The other points’ increments can be calculated through Equation (11).

$$\Delta a_i=\left(\frac{\Delta K_{\mathrm{eq},i}^m}{\Delta K_{\mathrm{eq},\max }^m}\right)\Delta a_{\max }$$

(11)

where ΔK_eq,i is the equivalent SIF range at a discrete point i, and ΔK_eq,max is the maximum SIF range at the crack front.

A maximum tangential stress criterion [33] was adopted to calculate θ; see Equation (12):

$$\theta =2{\tan }^{-1}\left[\frac{-2\Delta K_{\mathrm{II}}}{\Delta K_{\mathrm{Ieq}}+\sqrt{{\left(\Delta K_{\mathrm{Ieq}}\right)}^2+8{\left(\Delta K_{\mathrm{II}}\right)}^2}}\right]$$

(12)

where ΔK_Ieq is an equivalent SIF range with respect to the mode-III effect, which is determined through Equation (13).

$$\Delta K_{\mathrm{Ieq}}=\Delta K_{\mathrm{I}}+B\left|\Delta K_{\mathrm{III}}\right|$$

(13)

where B is an empirical parameter, which is 1.0.

Equation (14) gives the formula to calculate the crack growth increment under the global coordinate system.

$$\begin{array}{l}\Delta x=\Delta a\cos \theta \cos \alpha \\ \Delta y=\Delta a\cos \theta \sin \alpha \\ \Delta z=\Delta a\sin \alpha \end{array}$$

(14)

3.4. Update of Crack Model

During crack propagation, the number of discrete points at the crack front remains constant or increases. However, the increase in the number of points is uncertain, which makes crack model updating more difficult. Abaqus offers a function to determine whether the mesh is associated with the geometry model. When decoupling them, the mesh turns into an orphan mesh, and any modification made to the orphan mesh cannot be influenced by the geometry model. A method for updating the crack model can be proposed if the coordinates of the discrete point of the crack front for step n and step n + 1 are known in advance. Firstly, we build and mesh a 2D plane model based on the x and y coordinates of the two crack fronts; see Figure 7a. Secondly, we transform the mesh into an orphan mesh to modify the nodes coordinates. Thirdly, we modify the nodes’ z coordinates based on the z coordinates of the two crack fronts; see Figure 7b, Finally, we import the crack model of step n to assemble the modified and imported models; see Figure 7c. The updated model is named as the crack model of step n + 1. We import it into the file where the tubular joint is, transform the orphan mesh into geometry dimensions, and submit the job.

However, when updating, a special condition may occur. Generally, only a single-layer mesh is generated in the crack depth direction, as shown in Figure 7. But, with the increase in the crack dimension, a multilayer mesh may be generated; see Figure 8a. Besides the nodes on the crack fronts, other nodes’ coordinates are unknown, which are called redundant nodes. After modifying, the redundant nodes are not changed in their space position, thereby causing a huge deformation of the crack mesh; see Figure 8b. In order to solve this problem, a solution was proposed. Firstly, we identified the nodes subjected to the crack fronts of step n and n + 1, which were recorded as set Ω₁. Secondly, we recorded the rest nodes as set Ω₂. Thirdly, we calculated the distance between the nodes in Ω₁ and the nodes in Ω₂. Finally, we made the nodes in Ω₂ move in the direction of the minimum distance until they overlapped with the nodes in Ω₁.

4. Numerical Simulation

4.1. Comparative Test

This research used the T-shaped axial tensile test and the high-cycle fatigue test described in the literature [22] as an example to verify the accuracy of the algorithm developed in Section 3. Figure 9 is the diagram of the tubular joint model. The geometry dimensions of the joint were as follows: for the chord, its outer diameter D, length L, and thickness T were 406.4 mm, 1386 mm, and 12.6 mm, respectively. For the brace, its outer diameter d, length l, and thickness t were 139.7 mm, 430 mm, and 5 mm, respectively. G1~G4 were the geometric parameters of the fixed end, which were 450 mm, 32 mm, 150 mm, and 110 mm, respectively.

The Young’s modulus E of the steel used in the tubular joint was 206 Gpa. The Poisson’s ratio ν, yield strength, and ultimate strength were 0.3, 410 MPa, and 500 MPa, respectively. For the weld, the E and ν were 210 GPa and 0.3, respectively. The material constants used in the Paris law, C and m, were 2.66704 × 10⁻¹¹ m/cycle and 2.75, respectively. The load was applied to the end of the brace, and a fixed constraint was added at the end of the chord. For the axial tensile test, the load was 70 kN; for the high-cycle fatigue test, the maximum load was 140 kN, and the stress ratio was 0.05. In order to obtain the stress distribution around the weld, three directions were set on the chord and brace surfaces. The locations of the strain gauges are illustrated in Figure 10.

4.2. Numerical Simulation

Seven models were built and designated as A1 and B1. For A1, the model was used to examine the accuracy of the T-shaped tubular joint modeled using the mapping coordinate method. The accuracy of the high-cycle fatigue crack propagation algorithm was verified by comparing the results obtained from B1 to those of the tests. The models’ geometric dimensions and material parameters were the same as the ones mentioned in Section 4.1. The initial length and depth of the cracks were set as 2c₀ = 2 mm and a₀ = 1 mm, respectively.

In order to increase computing efficiency, different meshing schemes were used in different regions. For the chord, fixed end, and the intersection of the chord and brace, the mesh size was 10 mm. For the brace and weld, they were divided into nine regions. Figure 11 shows the meshing schemes of the brace and weld in the form of the plate model. Regions 1 and 5 were connected with other parts of the tubular joint. Region 2 was the mesh of the weld. Region 3 was where the crack could grow, whose dimensions were 0.25 mm × 0.7 mm. The type of all of the meshes was C3D8R.

4.3. Simulation Results

4.3.1. Results of Axial Tensile Numerical Simulation

Figure 12 shows the comparison of the surface stresses obtained from the test and simulation. For Line A, the two curves showed the same trend. Due to the influence of the weld polishing [22], the differences in stress on Line A were relatively significant from 0 mm to 5 mm. When the distance was beyond 10 mm, the simulated results tended to be steady, and the difference from the test result was about 20 MPa. For Line B and Line C, the simulation results showed great agreement with the test ones. Therefore, based on the result comparison, the proposed algorithm used for the T-shaped tubular joint modeling was validated.

4.3.2. Results of Fatigue Crack Propagation Simulation

Equation (12) was used to calculate crack growth increment, where Δa_max = 0.8 mm. The whole growth process had 80 steps; Figure 13a shows the crack model of the ultimate step, which grew in three-dimensional space and finally formed a double-curved crack with a length of 100.6 mm. Figure 13b shows the angle of the crack at the direction of depth. During propagation, the crack growth direction was not perpendicular to the surface of the tubular joint all of the time. For convergence, an assumption was adopted that the crack grew at a constant angle consistently, which was calculated using Equation (13). The angle obtained using the simulation was 9.88°, and the test result was 14° [22].

Because a predefined crack was adopted in the simulation, the number of load cycles obtained using the simulation only included the stage when the crack grew from the initial dimension to the critical dimension. For the stage from the crack initiation to the initial dimension, the simulation did not consider it. Thus, the simulated number of load cycles was lower than the test result, and a modified method should be taken. Under the initial stage, the depth and SIFs at the deepest point of the crack were substituted into Equation (10) and Equation (15) to calculate the equivalent stress intensity factor K_eq,0 and the number of load cycles N₀, which were used to represent the life of the crack initiation.

$$N_0=\frac{a_0}{C{\left(\Delta K_{\mathrm{eq},0}\right)}^m}$$

(15)

After calculating, N₀ was equal to 275,585. Thus, the number of load cycles added N₀ for each step. Figure 14 shows the crack propagation path of the simulation and the test. The results of the test and simulation show that, before growing to 50 mm, the crack grew around the weld toe. Then, the crack grew in the direction of deviation from the weld toe. In the test, when the length of the crack was 12.2 mm, the number of load cycles was 4.1 × 10⁵. In the simulation, the corresponding result was about 408,651, which was close to the test result. However, the critical length obtained from the test and simulation had a difference. This was because the material constants C and m adopted in the simulation were abstracted from the literature [22], which were relatively conservative. As a result, in each step, the fatigue crack growth rate calculated in the simulation was higher than that in the test. And with the increase in the crack dimension, the influence of the material constants was more obvious. In addition, the simulation did not take residual stress of the weld and weld initial defects into consideration. The ideal condition assumption may have decreased the simulated critical length of the crack.

In the test, benchmarking was used to trace the process of crack growth at the rupture surface; see Figure 15a. BM1~BM5 are the five outlines of cracks. The corresponding result is shown in Figure 15b, where the simulated crack outlines showed great agreement with the test ones. Figure 16 shows the curves of the fatigue crack propagation rates. When the number of load cycles was less than 500 thousand, the simulation results closely matched with the test results. When the number of load cycles was beyond 500 thousand, the simulated curves were slightly higher than the test ones; the general trend was the same. When the crack penetrated the wall thickness, the simulated fatigue life was 573,635, and the counterpart in the test was 634 thousand. The difference was 9.5%. Figure 17 shows the changes in crack aspect ratio a/c during the crack propagation. Comparison of the simulated curve with the seven groups of data recorded in the test showed the same trend, and the simulated curve agreed well with the experimental results.

5. Influence of Initial Crack Dimensions on Fatigue Crack Growth

The initial crack dimension was predefined in Section 4 according to the referenced literature [22]. However, these predefined crack dimensions were not based on experimental data. To explore the influence of the initial crack dimension on the fatigue crack growth, nine T-shaped tubular joints with different initial crack dimensions were built; see

Table 1. Models with different initial crack dimensions.

No.	L/mm	D/mm	T/mm	l/mm	d/mm	t/mm	2c/mm	a/mm	a/c
P1	1000	200	10	400	100	5	1.2	0.75	1.25
P2	1000	200	10	400	100	5	1.5	0.75	1
P3	1000	200	10	400	100	5	1.8	0.75	0.833
P4	1000	200	10	400	100	5	2.1	0.75	0.714
P5	1000	200	10	400	100	5	2.4	0.75	0.625
P6	1000	200	10	400	100	5	3	0.75	0.5
P7	1000	200	10	400	100	5	3.75	0.75	0.4
P8	1000	200	10	400	100	5	5	0.75	0.3
P9	1000	200	10	400	100	5	7.5	0.75	0.2

. The material parameters and meshing method were the same as for the models in Section 4. An axial tensile force of 30 MPa was applied at the end of the brace, and the displacement and rotation of the x axis and z axis were constrained at the end of the chord.

As shown in Figure 18, the initial crack dimension affected the SIFs at the crack front under the initial stage. When the a/c was no more than 0.625, the distribution curves of K_eq were similar to those of K_I, thus meaning that K_I played a dominant role in the crack growth. However, with the decrease in the a/c, the curves of K_I and K_II gradually fit together. When the a/c was lower than 0.625, the curve shape of K_eq was close to that of K_II, thus indicating that K_II replaced K_I as the main factor to control the crack propagation.

The θ of the P1~P9 models under the initial stage was calculated through Equation (9); see Figure 19. The result shows that when the a/c ranged from 0.625 to 1.25, the variation in θ was small, while when the a/c was lower than 0.625, the θ increased sharply. This indicates that the a/c is related to the θ. According to related tests [34], when there was no change in the fatigue loads, boundary conditions, geometric dimensions of the joints, and materials, only slight differences were observed in the crack angles and crack dimensions. Meanwhile, the location where the cracks were initiated was also similar. In other words, the location where the crack initiates, the crack dimension, and the crack angle have a relationship with the fatigue loads and boundary conditions. This relationship exists within a certain range. Thus, based on this result, an assumption was proposed as follows: When the a/c is in a certain range named Ω, if the change in θ of is small, then each a/c value in Ω can be used to determine the initial crack size. Then, by taking a relatively small value of a, the corresponding value of c is obtained. Therefore, models P1~P5 could be used for further study.

Figure 20 shows the crack propagation rates with different initial crack dimensions. The result shows that the crack propagation rates obtained from the five models P1~P5 were similar. The fatigue lives simulated were 160,416, 166,542, 171,412, 190,818, and 188,903, respectively. Except for P4, in which the fatigue life was a little higher than that of P5, the fatigue life tended to gradually increase with the decrease in the a/c, and the largest difference was 30,402.

6. Discussion

Fatigue damage is a cumulative process. When initial fatigue damage is detected, experimental study is a useful method to acquire the fatigue crack propagation and fatigue life. However, the whole process may be time-consuming, and the cost may be relatively high. Thus, numerical simulation is an alternative.

Based on the XFEM and Paris law and combined with Abaqus secondary development, this study provides a novel method to realize crack propagation under high-cycle fatigue load, and it acquires fatigue life via this method. We further studied the influence of the initial crack size on the crack growth and fatigue life, which can be acquired by detection.

It is known that fatigue simulation is available now. The reason why we conducted this study can be summarized as follows. For previous study, the theory of simulation is mainly based on the FEM, which needs to be remeshed when a crack grows [14,34]. It means that in each incremental step, the remeshing scheme is significant in order to acquire accurate results, which undoubtedly decreases computational power and increases technology difficulty. But the XFEM avoids remeshing due to its enrichment function, thereby simplifying the process of computation. Nikfam M.R. et al. [27] used this theory to realize crack propagation in a weld T-joint. And we adopted the same method to study the welded tubular T-joint, which is largely used in steel building structures.

Based on the current study, numerical simulation is deemed to be relatively preferable without considering experimental error, which causes differences in the fatigue life and crack growth path. Further study will be conducted to close these differences. And this study focuses on amplitude high-cycle fatigue load; however, in engineering practice, random amplitude fatigue is more common, such as wind loads and equipment vibration. And under many circumstances, the initial fatigue crack is not a single one. How to simulate several crack propagations and analyze the mechanism of crack combinations needs to be further studied.

7. Conclusions

The XFEM and secondary development techniques were used to realize high-cycle fatigue crack propagation simulations in T-shaped tubular joints. And, the main conclusions are summarized as follows:

(1) This study provided a method to realize high-cycle fatigue crack propagation simulations for T-tubular joint based on the XFEM and Abaqus secondary development. The crack propagation path, crack growth rate, and fatigue life were acquired, and these results show that the simulation results were in great agreement with the test results, thus indicating that this novel method is effective.

(2) Further study was conducted in terms of the initial crack size, which can be used for simulation when a lack of initial crack size is acquired by detection. When the a/c was in a certain range, the difference in the crack angle θ was not significant, thus indicating that the a/c had a relatively small influence on the crack growth. This can be used as a sensitive analysis for the crack size.

(3) Simulations of five tubular joints with different initial crack dimensions were conducted. The fatigue lives of these joints ranged from 160,416 to 190,818 cycles. With the decrease in the a/c, the fatigue life tended to increase. From an engineering perspective, when a lack of initial crack size is detected, the smaller a/c can acquire a safer result.