Computational Simulation of 3D Fatigue Crack Growth under Mixed-Mode Loading

Abstract

The purpose of this research was to present a simulation modelling of a crack propagation trajectory in linear elastic material subjected to mixed-mode loadings and investigate the effects of the existence of a hole and geometrical thickness on fatigue crack growth and fatigue life under constant amplitude loading. For various geometry thickness, mixed-mode (I/II) fatigue crack growth studies were carried out to utilize a single edge cracked plate with three holes and compact tension shear specimens with various loading angles. Smart Crack Growth Technology, a new feature in ANSYS, was used in ANSYS Mechanical APDL 19.2 to predict the cracks’ propagation trajectory and their consequent fatigue life associated with evaluating the stress intensity factors. The maximum circumferential stress criterion is implemented as a direction criterion under linear elastic fracture mechanics (LEFM). According to the hole position, the results demonstrate that the fatigue crack grows towards the hole due to the unbalanced stresses on the hole induced crack tip. The results of this simulation are verified in terms of crack growth paths, stress intensity factors, and fatigue life under mixed-mode load conditions, with several crack growth studies published in the literature showing consistent results.

1. Introduction

Most structures are subjected to cyclic loading, tension, and shear loads, known as mixed-mode fatigue loading. Fatigue failure is the most common type of failure in such configurations, caused by cracks and other defects in the components. The primary objective of fracture mechanics is to evaluate whether or not a structure will fail based on the presence of a crack. Crack growth is crucial in engineering structures, because it significantly impacts the quality and stability of engineering structures. Thus, the safety or reliability of engineering structures is vital to predicting the crack propagation path. As a result, in many industries, the accurate estimation of the crack path and fatigue life is essential in terms of reliability. Fracture mechanics is an essential tool in current materials science for enhancing the mechanical performance of mechanical components. A significant parameter for estimating a cracked structure’s lifetime is the stress intensity factor (SIFs). The stress intensity factor is defined physically as the intensity of load transmitted through the crack tip area due to introducing a crack into the component [1].

The SIFs evaluate the severity of the stress induced by remote loading around the crack tip. The associated instantaneous value of SIFs would follow the changes in crack geometry and stresses during crack growth. The stress intensity factor has a complicated function of the load, boundary conditions, crack propagation, geometry, and material characteristics. The fatigue crack propagation is evaluated using the equivalent stress intensity factor in Paris’ law. Experimental investigations are essential for fatigue assessment in several applications, such as the aerospace industry and the aviation industries. However, accurate calculation methods are needed to analyze crack propagation to prevent crack propagation and fatigue life in both static and dynamic loading [2]. The failure was caused by (a) faults such as interfaces and cracks and (b) the nature of fluctuating loads. When exposed to variable loads, cracks tend to initiate and grow until the structure no longer bears the load, leading to catastrophic failure. These cracks are characterized as fatigue cracks, and the predicted life is one of the most critical factors in determining the structure’s reliability. It can be calculated by adding the required number of loading cycles to nucleate the fatigue crack and cause the failure. The crack propagation rate is generally determined using the relationship between SIFs and the crack geometry. Crack initiation and propagation are always associated with SIFs in a complicated state in general [3,4,5,6]. In recent studies, the extended finite element technique (XFEM) introduced by [7] has frequently been adopted to study the fracture due to crack growth. Based on a standard finite element framework, it uses a unique displacement feature to allow for discontinuities, eliminating the need to re-mesh continuously during the crack tip growth process.

XFEM was utilized to perform crack growth analysis without updating the mesh to evaluate SIFs [8,9,10]. Nowadays, numerous codes could be used to model cracks propagating through structures, e.g., FEM [11] and BEM [11,12], both in mixed-mode [13,14,15,16] or in-plane [17]. In addition, the methods used for calculating K values, kink angles, and crack-growth rates can be numerous [6,7,8]. Significant work was done to establish appropriate models for evaluating the fatigue crack growth (FCG) and fatigue life to avoid fatigue failures. Numerous experimental methods were proposed; nevertheless, the methods are frequently time-consuming and expensive to perform. Using a simulation strategy that involves numerical analysis and the ANSYS APDL.19.2 finite element approach to decrease laboratory effort, time, and expenses is a great way to save time and money. Meanwhile, there are several of modern finite element simulation tools for fatigue crack propagation problems, including ANSYS, ABAQUS, FRANC3D, etc. Many fatigue crack problems identified to date by the literature use different computational approaches in simulation for simple and complex structures in 2D and 3D [2,18,19,20,21,22,23,24]. An energy-based criterion was developed by [25] using additively manufactured materials to analyze the effects of material anisotropy on the behavior of fatigue creak growth. The criteria are established based on the concept of strain energy density in the anisotropic domain by including T-stress into the solution and considering the effects of the specimen geometry. Therefore, this work uses the ANSYS APDL 19.2 to predict the mixed-mode SIFs accurately and the associated fatigue life for the single edge cracked plate and compact tension shear specimens. In particular, three methods have been widely used to illustrate the fatigue analysis of materials: the method of fracture mechanics developed by Paris and Erdogan [26], the method of strain-life independently proposed by Coffin [27], and the method of stress-life proposed by Wöhler [28]. The first approach was employed in this study for predicting fatigue life, whereby the crack tip could be described separately by the stress intensity factors. The primary motivation for this study was to make a significant contribution to the use of ANSYS as an efficient tool for modelling FCG under the condition of mixed-mode loading and monitoring the influence of the holes on the crack growth trajectory.

2. Mixed Mode Fatigue Life Evaluation Procedure Using ANSYS

ANSYS can model three kinds of cracks: arbitrary, semi-elliptical, and pre-meshed. The smart crack growth analysis tool uses the crack front in the pre-meshed crack approach, and the stress intensity factor is the failure criteria. The rendered nodes will distribute to the front, bottom, and top of the crack, and the front of the crack is included in the pre-meshed crack process, which is used by the “Smart Crack Growth” simulation tool. This is used instead of employing the XFEM’s enrichment area (splitting), which eliminates the re-mesh process around the crack tip and assumes that the discontinuities cut the element entirely. In this case, the displacement formulation does not account for the presence of singularity. As the crack propagates, the newly introduced crack segments are always assumed to have cohesive zone behavior. SMART automatically updates the mesh from crack geometry changes due to crack propagation on each solution stage, reducing the need for long pre-processing sessions. The sphere of influence was used to refined the crack tip mesh and around the geometric edge that passes through the thickness. The geometric regions ‘contain the crack tip, the crack’s top surface, and the bottom surface of the crack. However, each of these regions is connected to a set of nodes that will be utilized for analysis. In order to analyze a mixed-mode fracture, it is essential to have a proper understanding of the crack growth direction. The direction of crack extension is defined by an angle θ measured from the initial crack plane [29,30,31]. The direction of a mixed-mode crack is definitely determined by the ratio of modes stress intensity at the crack tip. The ANSYS software considers mixed-mode loading where the maximum circumferential stress criterion is implemented. The following are the formulas for the direction of crack propagation in ANSYS [32,33]:

$$\theta ={\cos }^{-1}\hspace{0.17em}\left(\frac{3K_{II}^2+K_I\sqrt{K_I^2+8K_{II}^2}}{K_I^2+9K_{II}^2}\right)$$

(1)

where:

K_I = Max K_I during cyclic loading and K_II = Max K_II during cyclic loading.

In the present simulation using ANSYS, the simulation of crack propagation is restricted to region II of the typical crack propagation under fatigue loading that can be represented as:

$$\frac{da}{dN}=C{(\Delta K_{eq})}^m$$

(2)

From Equation (2), for a crack increment, the number of life cycles of fatigue may be predicted as:

$${\displaystyle \underset{0}{\overset{\Delta a}{\int }}\frac{da}{C{(\Delta K_{eq})}^m}}={\displaystyle \underset{0}{\overset{\Delta N}{\int }}dN=\Delta N}$$

(3)

The equivalent range for the stress intensity factor formula can be found as follows [19]:

$$\Delta K_{eq}=\frac{1}{2}\cos \frac{\theta }{2}\left[(\Delta K_I(1+\cos \theta ))-3\Delta K_{II}\sin \theta \right]$$

(4)

where:

$\Delta K_I$ = the range of SIF in mode I loading and $\Delta K_{II}$ = the range of SIF in mode II loading.

Based on numerical analysis, there are many methods formulated for estimating the SIFs. The integral interaction method is usually the most accurate method that can estimate K_I and K_II separately. ANSYS proposes two methods, the Displacement Extrapolation Method (DEM) and the Integral Interaction Method (IIM), to determine SIFs. The second method was adopted because it is numerically easier to implement and has better precision and less mesh requirement. This approach depends on the domain integral approach [34], where an auxiliary field uses to separate K_I from K_II, as this ability is missing in the domain integral itself. The energy release rate is expressed in terms of mixed-mode stress intensity factors K_I, K_II, and K_III, which were proposed by [34,35] as follows:

Equation (5) becomes for the superimposed state:

$$J(s)=\frac{K_I^2+K_{II}^2}{E^{\ast }}+\frac{1+\nu }{E}{K}_{III}^2\hspace{1em}\hspace{1em}{E}^{\ast }=[\begin{array}{c}\frac{E}{(1-{\nu }^2)}\\ E\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\end{array}\begin{array}{l}\mathrm{Plane} \mathrm{strain}\hfill \\ \mathrm{Plane} \mathrm{stress}\hfill \end{array}$$

(5)

where

$$\begin{array}{ll}\hfill J^S(s)& =\frac{1}{E^{\ast }}\left[{(K_I+K_I^{aux})}^2+{(K_{II}+K_{II}^{aux})}^2\right]+\frac{1+\nu }{E}{(K_{III}+K_{III}^{aux})}^2\hfill \\ & =J(s)+J^{aux}(s)+I(s)\hfill \\ \\ \hfill \hspace{1em}I(s)& =\frac{1}{E^{\ast }}(2K_I{K}_I^{aux}+2K_{II}{K}_{II}^{aux})+\frac{1+\nu }{E}(2K_{III}{K}_{III}^{aux})\hfill \end{array}$$

(6)

where superscript (S) denotes the superimposed state, J(s) is the actual state domain integral, J^aux(s) is the auxiliary state domain integral, and I(s) is integral with interacting auxiliary and actual terms.

By setting $K_I^{aux}=1$ and $K_{II}^{aux}=K_{III}^{aux}=0$, Equation (6) yields:

$$K_I=\frac{E^{\ast }}{2}I(s)$$

(7)

By setting $K_{II}^{aux}=1$ and $K_I^{aux}=K_{III}^{aux}=0$, that gives the relationships between K_II, K_III as follow:

$$K_{II}=\frac{E^{\ast }}{2}I(s)$$

(8)

$$K_{III}=\mu \hspace{0.17em}I(s)$$

(9)

where E and $\mu$ are the modulus of elasticity and modulus of rigidity, respectively.

3. Numerical Results and Discussion

3.1. Single Edge Cracked Plate with a Three-Hole

A rectangular plate with a dimension of 120 mm × 65 mm × 16 mm contains two holes with a diameter of 13 mm near both ends of the plate and a 20-mm hole near the center of the plate, as depicted in Figure 1a. An initial crack of 10 mm is located at the center edge of the plate. The plate is made from Aluminum 7075-T6 with a Young’s Modulus of 71.7 GPa, yield strength of 469 MPa, the ultimate strength of 538 MPa, Paris’ law coefficient, C = 5.27 × 10⁻¹⁰, Paris law exponent, m = 2.947, fracture toughness of 938.25 MPa mm 0.5, and a Poisson’s ratio of 0.33. A 100 N load was applied at the top and bottom holes. The applied force is P = 20 kN with a stress ratio R = 0.1, and the behavior of the material is assumed to linear elastic. The initial mesh generated by Ansys is shown in Figure 1b, with an element size of 1 mm, generating 495,354 nodes and 348,892 elements.

As seen in Figure 2, a straightforward comparison is made of the crack path simulated with ANSYS software on the present study, one from experiment [36], the other three from XFEM [36], XFEM with ABAQUS software [37], and by using XFEM with a controllable crack propagation strategy [38] as displayed in Figure 2a–d. It can be shown that the numerically obtained crack path agrees well with the predicted experimental path.

The effect of thickness on fatigue crack growth parameters was investigated for this geometry using three different thicknesses, 4, 8, and 16. It is noticeable that there no influence of the geometry thickness on the crack extension trajectory. In contrast, there is a more significant influence on the stress and strain distribution as the thickness increased or decreased. Therefore, as thickness increases, the stress–strain field in the crack tip tends to change from the plane–stress state to the plane–strain state, meaning that the crack tip is in the tension state in all three directions and small scope the plastic zone will be limited. As the geometry thickness increases in a plane strain state, critical fracture toughness decreases, and brittle fracture is more likely to occur, which is much more dangerous than ductile fracture.

Consequently, as the geometry thickness is increased, the equivalent Von-Mises stress and the maximum principal stress were decreased with approximately the same ratio as shown in Figure 3. Figure 4 illustrates the equivalent elastic strain distribution for the three different thicknesses and gives additional insight into the effect of thickness on mechanical properties. The equivalent elastic strains decreased as the thickness increased, approximately in the same proportion as the equivalent Von-Mises and maximum principal stresses, as seen in this figure. The stress concentration is greater at the crack tip on the right top of the hole, as shown in Figure 3. The other stress value, which is somewhat lower than the previous value, is concentrated on the mid-right of the hole, providing the continuity of the crack to grow until the right end of the specimen, similar to the experimental study.

Figure 5 and Figure 6 display the estimated values of the two SIFs modes, K_I and K_II, respectively. As it can be seen in these two figures, the expected values for both K_I and K_II, while still under plane stress for geometry thicknesses of 4 mm and 8 mm, were very close to each other and much higher than those for 16-mm thickness, which was considered a plane strain condition. The values of K_II were increased with negative values, which only affect the crack growth direction to deviate toward the hole.

Since the fatigue life cycles are associated with the equivalent stress intensity factor, the number of cycles changes as $\Delta K_{eq}$ changes. The fatigue life cycles for different specimens is as follows: 16 mm start at 9109 and go up to 11,031, 8 mm start at 18,228 and go up to 20,284, 4 mm start at 20,766 and go up to 22,734, as shown in Figure 7.

3.2. Compact Tension Shear Specimen (CTS)

Figure 8a shows CTS geometry, while Figure 8b shows its proposed loading angles. The distribution of forces according to the load angles considered in ANSYS is shown in Figure 9. The applied load F is compared with the corresponding loads on other holes of the six holes, based on the given formulations [39]:

$$F_1=F_6=F(0.5\cos \alpha +\frac{c}{b}\sin \alpha )$$

(10)

$$F_2=F_5=F\sin \alpha$$

(11)

$$F_3=F_4=F(0.5\cos \alpha -(c/b)\sin \alpha )$$

(12)

where c and b are length parameters (c = b = 54 mm) and α is the loading angles as represented by Figure 2.

Table 1. Forces F₁ to F₆ as a function of loading angle α.

α	F2 = F5	F1 = F6	F3 = F4
30	0.5 F	0.933 F	–0.067 F
45	0.707 F	1.061 F	–0.354 F
60	0.866 F	1.116 F	−0.616 F

shows the results for all loading forces at different loading angles.

The considered material was Aluminum alloy having the following material properties as shown in

Table 2. Properties of the considered material.

Properties	Metric Units Value
Elasticity Modulus, E	74 GPa
Poisson’s ratio, υ	0.33
Yield strength, σy	517 MPa
Ultimate strength, σu	579 MPa
Threshold SIF, Kth	3.15 MPa $\sqrt{m}$
Fracture toughness, KIC	32.95 MPa $\sqrt{m}$
Paris’ law coefficient, C	4.3378 × 10−7
Paris law exponent, m	2.6183

:

The simulation is conducted under fatigue loading ratio R = 0.1, applied load F = 14 kN, and the ratio of the initial crack length to the width of the specimen (a/W = 0.5). This simulation considered three different thicknesses, which are 3 mm, 6 mm, and 12 mm, for three angles of loading 30°, 45°, and 60°. The Ansys simulated model shown in Figure 7 is entirely compatible with the experimental work conducted by [40,41,42,43].

Figure 10 shows the typical finite element mesh for each specimen with different thickness associated with the number of nodes and elements simulated by ANSYS. The element size in all simulation was 1 mm, and the total number of nodes and element for the three different thicknesses are shown in

Table 3. Total number of nodes and element for each specimen with different thicknesses.

	t = 12 mm	t = 6 mm	t = 3 mm
No. of nodes	897,045	623,867	443,766
No. of elements	615,479	419,579	288,809

. The geometry is meshed using (Q8) eight-noded isoparametric quadrilateral elements. At each step of the crack propagation, the finite element meshes are completely re-meshed, and the mixed-mode SIFs, K_I and K_II, are estimated.

The predicted crack growth path for the three different angles and three different thicknesses are shown in Figure 9 through the Von-Mises stress distribution. The effects of loading angles on the crack growth trajectory were demonstrated clearly. Three different thicknesses, 3, 6, and 12, were used to investigate the effect of thickness on fatigue crack propagation for this geometry. The crack growth trajectory was unaffected by the geometrical thickness, whether the thickness increased or decreased. However, there was a significant impact on the stress and strain distribution, as seen in Figure 10 for the corresponding equivalent Von-Mises stress, which decreased with approximately the same ratio as the specimen thicknesses were increased.

Comparisons were made for the predicted crack growth trajectories with different loading angles, as illustrated in Figure 11 with the experimental results obtained by [42] for a 30° loading angle, in Figure 12 with the experimental data reported by [44] for a 45° loading angle, and finally in Figure 13 with the experimental results obtained by [42] for a 60° loading angle. As illustrated in these three figures, the predicted crack extension trajectories were similar to the experimental pathways.

The number of cycles estimated in this investigation for different loading angles of 30°, 45°, and 60° is compared to experimental data reported by [42,44] under mixed-mode loading (I/II), resulting in significant agreements, as shown in Figure 14, Figure 15 and Figure 16.

4. Conclusions

In this study, mixed-mode (I/II) fatigue crack growth investigations were performed on a single edge cracked plate with three holes and different thicknesses as well as for a CTS specimen for different thicknesses and loading angles. It was observed that the geometry thickness did not affect the crack growth trajectory. In contrast, when the thickness of the specimen increased or decreased, it had a more significant impact on the stress and strain distribution. As a result, the stress–strain field at the crack tip tended to move from the plane stress state to the plane strain state as thickness increases, asserting that the crack tip is under tension in all three dimensions and the plastic zone will be constrained in a small scope. Accordingly, as the geometry thickness is increased, the equivalent Von-Mises stress and the maximum principal stress were decreased with approximately the same ratio of increasing the thickness. The fatigue crack grew toward the hole, as predicted, due to unbalanced stresses at the crack tip induced by the hole. In these simulation sequences, holes act as a crack stopper and attract a crack trajectory to grow. These results confirm the algorithm’s ability to identify crack-stopping holes in damage tolerance designs. The predicted fatigue crack propagation paths, stress intensity factors, and fatigue life cycles were validated by other researchers’ experimental and numerical results.