Fatigue Crack Propagation Study of Bridge Steel Q345qD Based on XFEM Considering the Influence of the Stress Ratio

Abstract

In the past, fatigue cracks have appeared in the orthotropic steel decks of bridges shortly after they opened to traffic. Previous studies have shown that high tensile welding residual stress exists in welded joints of steel bridges, which significantly changes the average stress and stress ratio of the joints. However, traditional fatigue crack propagation (FCP) calculations based on the Paris equation do not consider the influence of the stress ratio. Steel Q345qD is a common material used in bridges. Therefore, it is meaningful to study the influence of the stress ratio on the FCP life of steel Q345qD. In this paper, an FCP equation based on the energy release rate considering the influence of the stress ratio is first derived and named the da/dN-ΔG-R equation. Next, three material parameters in the equation are determined based on the results from tests of steel Q345qD under different stress ratios. Then, a user subroutine based on the extended finite element method (XFEM) that considers the influence of the stress ratio is defined and the effects of mesh size are analyzed. Finally, the effects of the stress ratio on FCP are discussed and the adaptability of the da/dN-ΔG-R equation is verified by comparing the values obtained from the equation with experimental results. The results show that: with a 95% guarantee rate, three material parameters in the da/dN-ΔG-R equation are: log(C) = −10.71, m = 2.780, and γ = 0.957; in the numerical simulation, a mesh size of 1 mm is more appropriate than other mesh sizes as it shows better accuracy and efficiency; under the same energy release rate range, the crack growth rate decreases as the stress ratio increases; under the same loading amplitude and cycles, the fatigue life decreases as the stress ratio increases; and finally, the numerical results considering the influence of stress ratio based on the da/dN-ΔG-R equation are close to the test results, while the results without considering the stress ratio based on the Paris equation are inaccurate.

1. Introduction

Steel box girder bridges with orthotropic steel decks are widely used in long-span highway and railway bridges because of their light weight, convenient construction, and high load-carrying capacity [1,2]. However, in one case, serious fatigue cracks appeared in the steel deck shortly after the bridge opened to traffic [3,4]. The fatigue life of some steel box girder bridges is less than 10 years, which is far lower than the bridge design life [4,5]. Previous studies have shown that welding residual stresses significantly changed the mean stress level of welded joints of steel box girder bridges, that is, the stress ratio changed significantly under certain stress amplitudes. The high mean stress (or stress ratio) consequently reduced the fatigue life of the steel bridges [5]. Based on the test data regarding 14MnNbq steel used in bridges, Liu et al. found that the fatigue crack propagation (FCP) rate increased as the stress ratio under the same stress intensity factor amplitude increased [6]. However, traditional FCP calculations based on the Paris equation only reflect the relationship between FCP and the stress intensity factor amplitude, without considering the influence of the stress ratio [7].

Scholars have done extensive research regarding the influence of the stress ratio on FCP. Forman et al. [8] assumed that the crack growth rate would tend to infinity when the maximum stress intensity factor approaches the fracture toughness, and subsequently established a fatigue crack propagation equation that considered the influence of the stress ratio and fracture toughness. Walker [9] used experimental data to propose an FCP equation that considered the influence of the stress ratio. Donald and Lados [10] introduced the adjusted compliance ratio parameter into the Walker two-parameter formula and considered the influence of residual stress. Kwofie et al. [11,12] deduced the stress intensity factor under different stress ratios according to the transformation relationship between stress amplitudes and mean stresses, and created a driving force model that considered the influence of stress ratios. Li et al. [13] proposed a modified Walker equation to characterize the FCP of a cracked body under constant amplitude cyclic loading. Zhan et al. [14] presented an exponential equation that considered the stress ratio effect. Cano et al. [15] compared several different driving force models based on the energy release rate while also considering the effects of the stress ratio. The term $\Delta \sqrt{G}$ seems to be a valid crack driving force for adhesives, composites and polymers such as copolyester polyethylene terephthalate glycol. Xin et al. [16] incorporated the mean stress effects by employing the Walker equation to fit the fatigue crack propagation rate of steel S355 and S690 based on experimental data. Ribeiro et al. [17] suggested deterministic quadratic relations based on experimental results regarding the relationship between the crack closure quantitative parameter U and the applied stress ratio R and between the effective stress ratio R_eff and the applied stress ratio R in comparison with various crack closure models. An FCP model for two different stress ratios was generated based on the strain energy density approach and considering the residual stress effects [18]. Moarrefzadeh et al. [19] modified the Walker equation by considering the influence of welding residual stress on the effective stress intensity factor and effective stress ratio. Therefore, it is necessary to derive an equation considering the effect of the stress ratio to describe fatigue crack growth more accurately; experiments under different stress ratios are sometimes needed to obtain such an equation.

Several methodologies can be used for modeling FCP, including the mesh-free method [20], boundary element method [21,22], discrete element method [23], traditional finite element method (FEM) [24,25], and extended finite element method (XFEM) [26,27,28,29]. The most used numerical simulation methods are traditional FEMs, such as ANSYS [24,25], the specialized crack propagation simulation software ZENCRACK [30], and the XFEM in ABAQUS [27,28,29]. In the traditional finite element method ANSYS, a wedge singular element should be established at the crack tip because of stress singularity, as shown in Figure 1a [25]. The stress intensity factor in ANSYS can be calculated, but the energy release rate cannot be extracted directly. ZENCRACK is a specialized finite element software that is used to simulate crack propagation. A special crack block is set at the crack tip, which can move with the crack propagation, as shown in Figure 1b [30]. This software is good at crack simulation but has poor universality and often needs to be used in combination with other conventional software. The XFEM was first proposed by Belytschko in 1999 [26]. The enrichment function is introduced into the conventional finite element displacement mode to reflect the crack discontinuity, and the level set function is used to trace the crack surface. Therefore, as shown in Figure 1c [28], the crack can penetrate the element directly without the use of refined mesh at the crack tip, which is easier in terms of meshing than ANSYS and ZENCRACK. It can also clearly show the continuous crack propagation process. At present, the XFEM in ABAQUS generally uses the virtual crack closure technique (VCCT) to calculate the energy release rate, and its accuracy is higher than that of the stress intensity factor calculated by the node displacement extrapolation method. It also has a good ability to solve fatigue fracture problems involving linear, nonlinear, and even anisotropic materials [29].

In this paper, an FCP equation based on the energy release rate considering the influence of the stress ratio is first derived. Then, three material parameters in the equation are determined based on tests of bridge steel Q345qD under different stress ratios. Finally, a user subroutine based on the XFEM that considers the influence of the stress ratio is defined. The effects of the stress ratio on FCP are discussed. The numerical simulation results considering and not considering the influence of the stress ratio are compared with the experimental results.

2. Deriving an FCP Equation Considering the Stress Ratio R Based on the Energy Release Rate G

An FCP calculation based on the Paris equation without considering the stress ratio R is shown in Equation (1):

$$\frac{da}{dN}=C_0\cdot \Delta K^{m_0}$$

(1)

where C₀, m₀ are material parameters, ΔK is the stress intensity factor range, a is the crack length, N is the number of cycles, and da/dN is the FCP rate.

Based on the Paris equation, Walker [9] proposed a fatigue crack growth equation that considered the stress ratio R based on ΔK, as shown in Equation (2):

$$\frac{da}{dN}=C{\left(\frac{\Delta K}{{(1-R)}^{1-\gamma }}\right)}^m$$

(2)

where C, γ, m are material parameters and R is the stress ratio.

For opening mode cracks, the relationship between the energy release rate G and stress intensity factor K is shown in Equation (3) [31]:

$$G=\frac{1}{E^{\prime }}{K}_{}^2$$

(3)

where E’ = E for plane stress, E’ = E/(1 – ν²) for plane strain, E is material elasticity modulus, and ν is Poisson’s ratio.

The stress ratio R and intensity factor range can be calculated as Equation (4):

$$R=K_{\min }/K_{\max },\Delta K=K_{\max }-K_{\min }$$

(4)

Thus, the energy release rate range is expressed in Equation (5) based on the relationship between G and K in Equation (3).

$$\Delta G=G_{\max }-G_{\min }=\frac{K_{\max }^2-K_{\min }^2}{E^{\prime }}$$

(5)

Substituting Equation (4) into Equation (5), Equation (6) is obtained.

$$\Delta G=\frac{1}{E^{\prime }}\cdot \frac{1+R}{1-R}\Delta K^2$$

(6)

ΔK can be solved from Equation (6) and substituted into Equation (2). Thus, as shown in Equation (7), an FCP equation considering the stress ratio R based on the energy release rate G, which is named the da/dN-ΔG-R equation, is obtained.

$$\frac{da}{dN}=C{\left(E^{\prime }\cdot \Delta G\cdot \frac{{(1-R)}^{2\gamma -1}}{(1+R)}\right)}^{\frac{m}{2}}$$

(7)

It can be seen from the derivation process that the material parameters C, γ, m in Equation (7) are the same as those in Equation (2). Therefore, a standard compact tensile test using the specifications outlined in [32] can be used to obtain the material parameters in the da/dN-ΔG-R equation.

3. Determining Material Parameters in the da/dN-ΔG-R Equation by FCP Tests

3.1. FCP Tests

The chemical composition of steel Q345qD is shown in

Table 1. Chemical composition of steel Q345qD [33].

Chemical Composition	C	Si	Mn	P	S	Nb	V
Mass proportion (%)	≤0.18	≤0.55	0.90~1.70	≤0.025	≤0.020	≤0.06	≤0.08
Chemical composition	Ti	Cr	Ni	Cu	Mo	N	Als
Mass proportion (%)	≤0.03	≤0.80	≤0.50	≤0.55	≤0.20	≤0.012	≥0.015

[33]. To obtain the material parameters C, γ, m of steel Q345qD, a standard FCP test was designed according to the metallic materials fatigue test specifications [32]. The size of the standard compact tensile specimen (CT specimen) is shown in Figure 2. The thickness of the specimen is 10 mm, the machining notch length a_n is 9 mm, and the initial crack length a₀ is 10 mm to eliminate the effect of machining notch on crack growth [32]. The material is steel Q345qD.

Table 2. Test loading scheme.

Specimen Case	Stress Ratio R	Minimum Load Pmin (kN)	Maximum Load Pmax (kN)
CT-1	0.1	1.67	16.67
CT-2	0.3	6.43	21.43
CT-3	0.5	15	30
CT-4	0.7	35	50

shows the test loading scheme. FCP tests in which the loading amplitude ΔP was kept at 15kN were conducted under four different stress ratios: 0.1, 0.3, 0.5, and 0.7.

As shown in Figure 3, a QBG-100 fatigue testing machine with a maximum load of 100 kN and maximum loading frequency of 50 Hz was used as the loading equipment. The axial tensile cyclic loading mode was used and the experimental maximum load and loading frequency were 50 kN and 20 Hz, respectively. The stress-controlled constant amplitude loading method was used until the crack reached the critical value and the specimen fractured. The initial crack length was 10 mm. During the test, the crack length was measured by visual inspection. The data regarding the crack length a and the number of cycles N can be recorded every 1 mm at the beginning of crack growth and recorded every 2~3 mm when the specimen is close to being fractured. The crack length a was measured on the front and back sides of the specimen, and the average value was taken as the actual crack length. For each case, three specimens were taken to avoid experimental error. The number of cycles N is the arithmetic mean value of the three specimens under the same crack length a. Figure 4 shows the FCP process of the CT specimen.

3.2. Analysis of Test Results

The a-N curves under different stress ratios measured in the tests are shown in Figure 5a, indicating that the specimen crack propagates faster as the stress ratio increases. For example, at N = 6 × 10⁴, the crack lengths are about 12 mm, 13 mm, 16 mm and 26 mm when R equals to 0.1, 0.3, 0.5 and 0.7, respectively. In order to obtain the FCP rate da/dN-ΔK curve, a seven-point incremental polynomial method recommended by GB/T6398-2017 is used to process the data [32]. That is, for each data point (a_i, N_i), the previous three points and next three points are taken, and these seven data points are used to fit a quadratic polynomial. The slope of the fitted curve at this data point is regarded as the FCP rate da/dN. ΔK under different crack lengths can be calculated using Equation (8). The FCP rate curves in Figure 5b show that the FCP rate increases as the stress ratio increases. Since the slopes and intercepts of these fitted curves vary with the stress ratio, this means that the Paris equation including only two material parameters and without considering the influence of the stress ratio is inaccurate. Equation (8) is shown below:

$$\Delta K=\frac{\Delta P}{B\sqrt{W}}\cdot \frac{2+\alpha }{{(1-\alpha )}^{3/2}}\left(0.886+4.64\alpha -13.32{\alpha }^2+14.72{\alpha }^3-5.6{\alpha }^4\right)$$

(8)

where ΔP is the load amplitude, B is the specimen thickness, W is the horizontal distance from the loading point to the edge of the specimen (as shown in Figure 2), a′ is the calculated crack length which is the horizontal distance from the loading point to the crack tip, and α = a′/W.

Equations (1) and (2) are expressed in logarithmic form, as shown in Equations (9) and (10), respectively. The material parameters C₀, m₀ in the Paris Equation (9) and C, γ, m in the Walker Equation (10) of steel Q345qD can be obtained by experimental data fitting, as shown in Figure 6. The material parameters in the da/dN-ΔG-R equation are the same as those in the Walker Equation (10) and are obtained according to the derivation process described in the last section. The material parameters of steel Q345 are shown in

Table 3. Material parameters of steel Q345qD.

FCP Equation	Guarantee Rate	Material Parameters
		log(C0) or log©	m0 or m	γ
Paris	Mean	−10.80	2.687
	With 95% guarantee	−10.54	2.847
Walker, da/dN-ΔG-R	Mean	−10.85	2.696	0.947
	With 95% guarantee	−10.71	2.780	0.957

.

$$\mathrm{Paris}:\log (da/dN)=\log C_0^{}+m_0\log (\Delta K)$$

(9)

$$\mathrm{Walker}:\log (da/dN)=\log C+m\log (\Delta K)+m(\gamma -1)\log (1-R)$$

(10)

4. Numerical Simulation Considering the Stress Ratio R Based on the XFEM

4.1. FEM of the CT Specimen

Figure 7 shows the FEM of the CT specimen and the model size that is consistent with the test. The mesh is fine near the crack and coarser away from the crack. Two reference points are defined such that they are coupled to the inner surface of the circular hole. This is done so that the boundary conditions and loads can be applied directly to the reference points. Reference point (RP)-1 constrains translation in the x, y, and z directions and rotation in the x and z directions. Translation and rotation in the x, y and z directions of RP-2 are constrained. The cyclic load P is applied on RP-1 according to the loading scheme in Table 2. The elastic modulus and Poisson′s ratio of steel Q345qD are 2.06 × 10⁵ MPa and 0.3, respectively [34]. As shown in Figure 7, the cyclic load P is only in the XY plane. Therefore, the CT specimen under the cyclic load P experiences plane stress because the stress is only in the XY plane and the stress perpendicular to the XY plane can be ignored. Under these conditions, E’ in Equations (3) and (7) equals to E, that is 2.06 × 10⁵ MPa.

4.2. FCP User-Defined Subroutine Considering R Based on the XFEM

When the energy release rate of a crack tip satisfies G_th < ΔG < G_pl, a fatigue crack begins to grow. G_pl is the upper threshold value of the energy release rate, which is close to the steel fracture toughness G_C. G_th is the lower threshold value of the energy release rate. The stress intensity factor K_th of steel Q345qD is 92 MPa·mm^1/2 [35]. The energy release rate threshold G_th can be obtained through Equation (3). When the maximum energy release rate of a crack tip is greater than the lower threshold value G_th, a fatigue crack propagates along the direction of the maximum principal stress in ABAQUS. Assuming that the ith element at the crack tip is the crack growth direction, the virtual crack closure technique (VCCT) is used to calculate the energy release rate ΔG_i. ΔN_i is calculated using Equation (7). The element with the lowest number of cycles is the final fracture element, and the number of cycles increment ΔN and crack propagation increment Δɑ are obtained. The cycle continues until the crack length reaches the final crack length ɑ_f. To consider the influence of the stress ratio R, a user-defined subroutine is used in the simulation. A flow chart of FCP considering R effects based on the XFEM is shown in Figure 8.

4.3. Effect of Mesh Size

Figure 9a shows four different mesh sizes: 0.5 mm, 1 mm, 1.5 mm, and 2 mm. Calculation results of these four different meshing models under ΔP = 15 kN and R = 0.1 are shown in Figure 9b. Under the same cycles, the crack length decreases as the mesh size increases. For example, under N = 2 × 10⁵, the crack lengths are 22.5 mm, 22.3 mm, 19 mm, and 18 mm for mesh sizes of 0.5 mm, 1 mm, 1.5 mm and 2 mm, respectively. The results are towards insecurity with the increase of mesh size. However, the calculation results obtained when using a mesh size of 1mm are close to those obtained when using a mesh size of 0.5 mm, while the calculation efficiency is obviously improved. Therefore, a 1 mm mesh size is more appropriate to use in the simulation.

The maximum principal stress distributions of the FCP process with a mesh size of 1mm under ΔP = 15 kN and R = 0.1 are shown in Figure 10. The stress is large near the crack tip and decreases rapidly away from it. The stress at the crack tip increases with the degree of crack propagation. As shown in Figure 10b–d, the stress at the crack tip is about 500~2200 MPa at the steady crack propagation stage. When the model is approaching the fracture point, the stress reaches 1.2 × 10⁵ MPa, as shown in Figure 10f. The increment number of cycles decreases while the crack propagation length remains constant, indicating that the growth rate of the fatigue crack increases.

4.4. Effect of the Stress Ratio R

Fatigue crack propagation da/dN-ΔG curves at R ratios of 0.1, 0.3, 0.5 and 0.7 are shown in Figure 11. The simulation results were calculated based on the XFEM using the da/dN-ΔG-R equation, as shown in Figure 8, and were close to the experimental results. Figure 11 shows that under the same energy release rate range ΔG, the crack growth rate da/dN decreases as the stress ratio R increases. However, under the same number of cycles N, the crack length a increases along with R, as shown in Figure 5a. This means that in the experiment, the fatigue crack growth rate increases as R increases, which seems to contradict the results in Figure 11. The reason for this is that ΔG is not a constant value under the same load amplitude ΔP and crack length a, and instead increases along with R, as shown in Equation (6). ΔK is a constant value under the same ΔP and crack length a according to Equation (8). For example, when a = 20 mm, ΔG is 8334 N/m, 12,664 N/m, 20,458 N/m, and 38,642 N/m and da/dN is 2.45 × 10⁻⁷, 2.54 × 10⁻⁷, 2.67 × 10⁻⁷, and 2.87 × 10⁻⁷ for R = 0.1, 0.3, 0.5, and 0.7, respectively, as shown in the Figure 11. Therefore, under the same crack length a and load amplitude ΔP, the fatigue crack growth rate increases as R increases. In the experiment, under the same number of cycles N, the crack length a increases as R increases.

Figure 12 compares the experimental and numerical simulation results under different stress ratios. In Figure 12, the initial crack length a₀ is 10 mm, which is given in Section 3.1. Numerical simulation results considering the stress ratio R were calculated based on the XFEM using the da/dN-ΔG-R equation, as shown in Figure 8. The simulation results without considering R were calculated using the Paris equation. Figure 11 shows that the simulation results considering R from this paper are closer to the experimental data. However, the simulation crack length a considering R is slightly smaller than the test results obtained under the same number of cycles N. This may be because in the tests, the fixture and specimen are not completely fixed constraints, which is slightly different from the simulation. The simulation results without considering R calculated by the Paris equation, which are close to the mean values of different stress ratios, are inaccurate; when R is equal to 0.1 or 0.3, the crack length without considering R is larger than the test value, while it is smaller than the test value when R is equal to 0.5 or 0.7. Therefore, these findings suggest that the method based on the da/dN-ΔG-R equation proposed in this paper should be used to calculate the FCP life of steel Q345qD.

5. Conclusions

In this paper, an energy release rate-based fatigue crack propagation equation was derived and used to describe the FCP rate of bridge steel Q345qD using the experimental data obtained from various stress ratio tests. Then, a user subroutine based on the XFEM that considers the influence of the stress ratio was defined and the effects of mesh size were analyzed. Finally, the effects of the stress ratio were discussed and the adaptability of the equation was verified by comparing the values obtained from the equation with the experimental results. The main conclusions of this research are the following:

• A fatigue crack propagation equation based on the energy release rate considering the influence of the stress ratio was derived and named the da/dN-ΔG-R equation.
• Three material parameters in the da/dN-ΔG-R equation were determined with a 95% guarantee rate based on the standard compact tensile tests of steel Q345qD under different stress ratios. These parameters are log(C) = −10.71, m = 2.780, and γ = 0.957, respectively.
• A fatigue crack propagation subroutine using the da/dN-ΔG-R equation based on the XFEM was defined. In the numerical simulation, a CT specimen mesh size of 1 mm showed better accuracy and efficiency than other mesh sizes.
• The effects of the stress ratio on fatigue crack propagation were discussed. The results showed that under the same energy release rate range, the fatigue crack growth rate decreases as the stress ratio increases.
• Under the same loading amplitude and number of cycles, the fatigue life decreases as the stress ratio increases. The numerical results considering the influence of stress ratio using the da/dN-ΔG-R equation based on the XFEM were close to the test results, while the results without considering the stress ratio based on the Paris equation were inaccurate.

6. Future Work

• ∙ When considering the deck-rib welding joints of steel bridges, the residual stress is much greater than the external vehicle load stress [5]. Therefore, the real effective stress ratio will change significantly and not be equal to the external load stress ratio. The method proposed in this paper that considers the influence of the stress ratio will help guide future studies regarding the effects of welding residual stress on the fatigue life of steel bridges.
• ∙ In practical engineering, the load forms, such as tension, bending, shearing, and some combination forms, are complicated than those used in this study. Thus, the influence of load forms on FCP should be studied further.
• ∙ In practical engineering, the initial crack forms are random and require special flaw detection equipment to be detected. Two-dimensional surface cracks are also common. Therefore, further research regarding methods for describing the propagation of two-dimensional surface cracks is also needed.