2D and 3D numerical simulation of fatigue crack growth path and life predictions of a linear elastic

The most successful engineering implementation of fracture mechanics is crack growth laws based on stress intensity factor ranges [43]. Irwin [44] defined the stress intensity factor, abbreviated as K, as a mechanical parameter that represents the state of stress at a crack point. For various common crack configurations, many analytical solutions for the SIFs were developed. On the other hand, these analytical solutions are restricted to simple crack geometries and loading conditions [2]. Small-sized yield requirements were implemented [45], which established a J-integral methodology for investigating the nonlinear behavior of the material [46]. The equivalent domain integral methodology replaces the finite-size domain with a divergence theorem by integration along the contour for FEM. The integral independent contour of the trajectory is as follows [46]: (5) J=∫C(Wn1−σijnj∂ui∂x)ds J = \int\limits_C {\left( {W{n_1} - {\sigma _{ij}}{n_j}{{\partial {u_i}} \over {\partial x}}} \right)ds} where W is the strain-energy density, σ_ij are the stresses, u_i is the local i axis displacement, s is the contour arc length, and n_j is the outward unit normal to the contour C around the crack tip (Figure 2A).

Fig. 2

For 2D structures, the contour integral i is replaced by area integral (Figure 2B). Then, Eq. (5) is modified as: (6) Jk=−∫A[W∂q∂x−σij∂ui∂x∂q∂x]dA−∫A[∂W∂x−∂∂x(σij∂ui∂x)]qdA−∫sti∂ui∂xqds \matrix{ {{J_k}} \hfill & = \hfill & { - \int\limits_A {\left[ {W{{\partial q} \over {\partial x}} - {\sigma _{ij}}{{\partial {u_i}} \over {\partial x}}{{\partial q} \over {\partial x}}} \right]dA} } \hfill \cr {} \hfill & {} \hfill & { - \int\limits_A {\left[ {{{\partial W} \over {\partial x}} - {\partial \over {\partial x}}\left( {{\sigma _{ij}}{{\partial {u_i}} \over {\partial x}}} \right)} \right]qdA} - \int\limits_s {{t_i}{{\partial {u_i}} \over {\partial x}}qds} } \hfill \cr }

The step by step procedure of the developed program software is illustrated in Figure 3. The details of these steps are explained in the literature [29, 45, 47,48,49,50]. The adaptive mesh refinement was explained in detail by Bashiri and Alshoaibi [45].

Fig. 3

In LEFM simulation, the J-integral by definition takes into account translational mechanical energy balance in the tip region of the crack along x axis. Eq. (6) simplifies the computation of mode I and II of SIFs but they do not compute K_I and K_II separately for multiple load conditions. In these cases, the following integral should be considered [50]: (7) Jk=−∫A[W∂q∂xk−σij∂ui∂xk∂q∂xj]dA−∫A[W∂q∂xk−σij∂∂xj(∂ui∂xk)]qdA−∫sti∂ui∂xkqds \matrix{ {{J_k}} \hfill & = \hfill & { - \int\limits_A {\left[ {W{{\partial q} \over {\partial {x_k}}} - {\sigma _{ij}}{{\partial {u_i}} \over {\partial {x_k}}}{{\partial q} \over {\partial {x_j}}}} \right]dA} } \hfill \cr {} \hfill & {} \hfill & { - \int\limits_A {\left[ {W{{\partial q} \over {\partial {x_k}}} - {\sigma _{ij}}{\partial \over {\partial {x_j}}}\left( {{{\partial {u_i}} \over {\partial {x_k}}}} \right)} \right]} qdA} \hfill \cr {} \hfill & {} \hfill & { - \int\limits_s {{t_i}{{\partial {u_i}} \over {\partial {x_k}}}qds} } \hfill \cr } where k represents the crack tip's index. The SIFs can be calculated using one of the following methods. The first method uses the J-integral and SIFs as: (8) J1=κ+18μ(KI2+KII2)J2=−κ+14μKIKII \matrix{ {{J_1}} \hfill & = \hfill & {{{\kappa + 1} \over {8\mu }}\left( {K_I^2 + K_{II}^2} \right)} \hfill \cr {{J_2}} \hfill & = \hfill & { - {{\kappa + 1} \over {4\mu }}{K_I}{K_{II}}} \hfill \cr }

The second method uses the relation between SIFs and J₁ and J₂ as: (9) KI=0.58μκ+1+(J1−J2+J1+J2)KII=0.58μκ+1+(J1−J2−J1+J2) \matrix{ {{K_I}} \hfill & = \hfill & {0.5\sqrt {{{8\mu } \over {\kappa + 1}}} + \left( {\sqrt {{J_1} - {J_2}} + \sqrt {{J_1} + {J_2}} } \right)} \hfill \cr {{K_{II}}} \hfill & = \hfill & {0.5\sqrt {{{8\mu } \over {\kappa + 1}}} + \left( {\sqrt {{J_1} - {J_2}} - \sqrt {{J_1} + {J_2}} } \right)} \hfill \cr }

Where κ is the elastic parameter defined as: κ=|(3−4ν)(3−ν)/(1+ν)plane stressplane strain \kappa = \left| {\matrix{ {\left( {3 - 4\nu } \right)} \cr {\left( {3 - \nu } \right)/\left( {1 + \nu } \right)} \cr } \matrix{ {{\rm{plane}}\,{\rm{stress}}} \hfill \cr {{\rm{plane}}\,{\rm{strain}}} \hfill \cr } } \right.

Adaptive mesh refinement is an optimization method used in the field of FE mesh. The a posteriori error estimation used in this scheme is derived from the previous mesh generation. The stress error norm, also known as relative stress norm error, is used to estimate the mesh refinement error. The ratio of element norm stress error to the average standard stress error of the entire region was calculated using h-type adaptive mesh optimization. The mesh size for each element in this scheme is provided as: (10) he=2Ae, {h_e} = \sqrt {2{A_e}}, where A_e is the triangle element area. For each element, the norm stress error is represented by: (11) ‖e‖e2=∫Ωe(σ−σ*)T(σ−σ*)2dΩ=∫Ωe({σxσyτxyσz}−{σx*σy*τxy*σz*})T({σxσyτxyσz}−{σx*σy*τxy*σz*})2dΩ \matrix{ {\left\| e \right\|_e^2} \hfill & = \hfill & {\int_{{\Omega ^e}} {{{\left( {\sigma - {\sigma ^*}} \right)}^T}\left( {\sigma - {\sigma ^*}} \right)2{\rm{d}}\Omega } } \hfill \cr {} \hfill & = \hfill & {\int_{{\Omega ^e}} {{{\left( {\left\{ {\matrix{ {{\sigma _x}} \cr {{\sigma _y}} \cr {{\tau _{xy}}} \cr {{\sigma _z}} \cr } } \right\} - \left\{ {\matrix{ {\sigma _x^*} \cr {\sigma _y^*} \cr {\tau _{xy}^*} \cr {\sigma _z^*} \cr } } \right\}} \right)}^T}\left( {\left\{ {\matrix{ {{\sigma _x}} \cr {{\sigma _y}} \cr {{\tau _{xy}}} \cr {{\sigma _z}} \cr } } \right\} - \left\{ {\matrix{ {\sigma _x^*} \cr {\sigma _y^*} \cr {\tau _{xy}^*} \cr {\sigma _z^*} \cr } } \right\}} \right)2{\rm{d}}\Omega } } \hfill \cr }

However, the average norm stress error over the whole domain is represented as: (12) ‖e^‖2=1m∑e=1m∫ΩeσTσ 2dΩ=1m∑e=1m∫Ωe{σxσyτxyσz}T{σxσyτxyσz}2dΩ \matrix{ {{{\left\| {\hat {\rm e}} \right\|}^2}} \hfill & = \hfill & {{1 \over m}\sum\limits_{e = 1}^m {\int_{{\Omega ^e}} {{\sigma ^T}\sigma \,2{\rm{d}}\Omega } } } \hfill \cr {} \hfill & = \hfill & {{1 \over m}\sum\limits_{e = 1}^m {\int_{{\Omega ^e}} {{{\left\{ {\matrix{ {{\sigma _x}} \cr {{\sigma _y}} \cr {{\tau _{xy}}} \cr {{\sigma _z}} \cr } } \right\}}^{\rm{T}}}\left\{ {\matrix{ {{\sigma _x}} \cr {{\sigma _y}} \cr {{\tau _{xy}}} \cr {{\sigma _z}} \cr } } \right\}2{\rm{d}}\Omega } } } \hfill \cr } where m is the overall number of elements in the entire domain and σ^* is the smoothed stress vector. The integration with the triangular isoparametric dimension will be converted in the FEM by the summation of quadratics according to the Radau rules as follows: (13) ‖e‖e2=∫−11∫−11({σ(ξ, η)xσ(ξ, η)yτ(ξ, η)xyσ(ξ, η)z}−{σ(ξ, η)x*σ(ξ, η)y*τ(ξ, η)xy*σ(ξ, η)z*})T({σ(ξ, η)xσ(ξ, η)yτ(ξ, η)xyσ(ξ, η)z}−{σ(ξ, η)x*σ(ξ, η)y*τ(ξ, η)xy*σ(ξ, η)z*})te det Jedξdη=∑p=13({σ(ξp, ηp)xσ(ξp, ηp)yτ(ξp, ηp)xyσ(ξp, ηp)z}−{σ(ξp, ηp)x*σ(ξp, ηp)y*τ(ξp, ηp)xy*σ(ξp, ηp)z*})T({σ(ξp, ηp)xσ(ξp, ηp)yτ(ξp, ηp)xyσ(ξp, ηp)z}−{σ(ξp, ηp)x*σ(ξp, ηp)y*τ(ξp, ηp)xy*σ(ξp, ηp)z*})te det JeWp \matrix{ {\left\| {\rm{e}} \right\|_e^2} \hfill & = \hfill & {\int\limits_{ - 1}^1 {\int\limits_{ - 1}^1 {{{\left( {\left\{ {\matrix{ {\sigma {{\left( {\xi ,\,\eta } \right)}_x}} \cr {\sigma {{\left( {\xi ,\,\eta } \right)}_y}} \cr {\tau {{\left( {\xi ,\,\eta } \right)}_{xy}}} \cr {\sigma {{\left( {\xi ,\,\eta } \right)}_z}} \cr } } \right\} - \left\{ {\matrix{ {\sigma \left( {\xi ,\,\eta } \right)_x^*} \cr {\sigma \left( {\xi ,\,\eta } \right)_y^*} \cr {\tau \left( {\xi ,\,\eta } \right)_{xy}^*} \cr {\sigma \left( {\xi ,\,\eta } \right)_z^*} \cr } } \right\}} \right)}^{\rm{T}}}} } } \hfill \cr {} \hfill & {} \hfill & {\left( {\left\{ {\matrix{ {\sigma {{\left( {\xi ,\,\eta } \right)}_x}} \cr {\sigma {{\left( {\xi ,\,\eta } \right)}_y}} \cr {\tau {{\left( {\xi ,\,\eta } \right)}_{xy}}} \cr {\sigma {{\left( {\xi ,\,\eta } \right)}_z}} \cr } } \right\} - \left\{ {\matrix{ {\sigma \left( {\xi ,\,\eta } \right)_x^*} \cr {\sigma \left( {\xi ,\,\eta } \right)_y^*} \cr {\tau \left( {\xi ,\,\eta } \right)_{xy}^*} \cr {\sigma \left( {\xi ,\,\eta } \right)_z^*} \cr } } \right\}} \right){t^e}\,\det \,{J^e}d\xi d\eta } \hfill \cr {} \hfill & = \hfill & {{{\sum\limits_{p = 1}^3 {\left( {\left\{ {\matrix{ {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_x}} \cr {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_y}} \cr {\tau {{\left( {{\xi _p},\,{\eta _p}} \right)}_{xy}}} \cr {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_z}} \cr } } \right\} - \left\{ {\matrix{ {\sigma \left( {{\xi _p},\,{\eta _p}} \right)_x^*} \cr {\sigma \left( {{\xi _p},\,{\eta _p}} \right)_y^*} \cr {\tau \left( {{\xi _p},\,{\eta _p}} \right)_{xy}^*} \cr {\sigma \left( {{\xi _p},\,{\eta _p}} \right)_z^*} \cr } } \right\}} \right)} }^{\rm{T}}}} \hfill \cr {} \hfill & {} \hfill & {\left( {\left\{ {\matrix{ {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_x}} \cr {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_y}} \cr {\tau {{\left( {{\xi _p},\,{\eta _p}} \right)}_{xy}}} \cr {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_z}} \cr } } \right\} - \left\{ {\matrix{ {\sigma \left( {{\xi _p},\,{\eta _p}} \right)_x^*} \cr {\sigma \left( {{\xi _p},\,{\eta _p}} \right)_y^*} \cr {\tau \left( {{\xi _p},\,{\eta _p}} \right)_{xy}^*} \cr {\sigma \left( {{\xi _p},\,{\eta _p}} \right)_z^*} \cr } } \right\}} \right){t^e}\,\det \,{J^e}{W_p}} \hfill \cr } and similarly (14) ‖e^‖2=1m∑e=1m∑p=13({σ(ξp, ηp)xσ(ξp, ηp)yτ(ξp, ηp)xyσ(ξp, ηp)z})T({σ(ξp, ηp)xσ(ξp, ηp)yτ(ξp, ηp)xyσ(ξp, ηp)z})tedetJeWp \matrix{ {{{\left\| {{\rm{\hat e}}} \right\|}^2}} \hfill & = \hfill & {{1 \over m}\sum\limits_{e = 1}^m {\sum\limits_{p = 1}^3 {{{\left( {\left\{ {\matrix{ {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_x}} \cr {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_y}} \cr {\tau {{\left( {{\xi _p},\,{\eta _p}} \right)}_{xy}}} \cr {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_z}} \cr } } \right\}} \right)}^{\rm{T}}}} } } \hfill \cr {} \hfill & {} \hfill & {\left( {\left\{ {\matrix{ {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_x}} \cr {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_y}} \cr {\tau {{\left( {{\xi _p},\,{\eta _p}} \right)}_{xy}}} \cr {\sigma {{\left( {{\xi _p},\,{\eta _p}} \right)}_z}} \cr } } \right\}} \right){t^e}\det {J^e}{W_p}} \hfill \cr } where t^e is the thickness of the element for a plane stress case and t^e = 1 is the plane strain case. Consequently, the relative stress norm error ζ_e for each item is considerably <5%, which is acceptable for many applications in engineering. Thus, (15) ζe=‖e‖e‖e^‖≤ζ {\zeta _e} = {{{{\left\| e \right\|}_e}} \over {\left\| {\hat e} \right\|}} \le \zeta and the new element relative stress error level is identified as a permissible error by: (16) εe=‖e‖eζ‖e^‖≤1 {\varepsilon _e} = {{{{\left\| e \right\|}_e}} \over {\zeta \left\| {\hat e} \right\|}} \le 1

This means each element with ɛ_e is required to be refined and the new size of the mesh needs to be expected. Here the asymptotic convergence rate criteria are implemented, which presume: (17) ‖e‖e∝hep {\left\| e \right\|_e} \propto h_e^p where p is the approximation of polynomial order. In this study, p = 2 since the quadratic polynomial is used for the calculation of FEs. Thus, the new element estimated size is: (18) hN=1εehe {h_N} = {1 \over {\sqrt {{\varepsilon _e}} }}{h_e}

Depending on the user numbers of refinement, the existing mesh will be considered as the new background mesh and the advancing front method will be replicated. The element type was triangular with six nodes.

As illustrated in Figure 4, a four holes’ rectangular plate of 10 mm diameter and a 6 mm long edge crack is analyzed. The dimensions of the plate are 100 mm × 100 mm × 10 mm. The plate is supported on the bottom and is subjected to a tension force of 10 MPa. The plate is made of aluminium 7075-T6, whose properties are shown in Table 1.

Fig. 4

The geometry is fixed in the x and y directions at the bottom face in the developed program, while in Ansys the geometry is fixed in the x, y, and z directions.

The meshed geometry created in Ansys workbench has 14,737 nodes and 7,177 elements, while the meshed geometry generated in the developed program has 3,933 nodes and 1,911 elements, as represented in Figure 5.

Fig. 5

The predicted crack growth paths in both software is shown in Figure 6, and these were also in good agreement with that path predicted by Knowles and Sternberg [51] using the fast multipole method (FMM) boundary element method.

Fig. 6

The maximum and minimum principal stresses were compared side by side in both programs, as represented in Figure 7. In Figure 7, the crack starts to grow in a straight line as there is no influence of any holes. As it approaches the hole, the path propagates towards the hole. The crack deviated its direction because the first upper hole was not close enough to attract the crack to sink in the hole. It continued to propagate in a straight line up until it reached near the second upper hole, close enough to the crack trajectory so that the crack sinks into this hole. The predicted estimation of the mixed-mode stress intensity parameters K_I and K_II are depicted in Figure 8. These estimation values accurately represented the crack growth trajectory. The values of the first mode of SIFs decreased in the interval of the crack length between 15 mm and 20 mm consequent to the influence of the first upper hole. In contrast, the second mode of stress intensity factor increased slightly for the same interval, as Figure 8 shows. As the crack trajectory missed the first upper hole, the value of K_I starts to increase again until the crack reached the second upper hole.

Fig. 7

Fig. 8

Figure 9 shows the comparison for the stress contour distribution in y direction between the developed program results, Ansys results, and the results obtained by Knowles and Sternberg [51] using FMM boundary element method.

Fig. 9

A more complex mixed-mode crack growth problem was investigated to illustrate the performance of the developed program to represent the propagation of cracks from beginning notches in a three holes single edge cracked plate and precisely prognosticate the values of SIFs. The geometry of three holes single edge cracked plate, and its final adaptive mesh, is presented in Figure 10. This case was investigated for the first time by Ingraffea and Grigoriu [53]. The results revealed that small changes in the crack path's initial position significantly impacted the results. As a benchmark for numerical models, this example is commonly used to assess the system's predictive abilities. As indicated in Table 2, the geometry analyzed has 2L = 508 mm, width = 203.2 mm, and thickness = 12.7 mm, with two alternative configurations depending on the initial crack length and position. The initial meshes for this geometry, generated by the developed program and Ansys before starting the simulation, are shown in Figure 11.

Fig. 10

Fig. 11

The length of the crack in the first configurations is a_o = 25.4 mm, and it is positioned at a distance of b = 152.4 mm from the mid-length of the beam. The second configuration differs from the first by having an initial crack length of a_o = 38.1 mm, and it is positioned at a distance of b = 127 mm. The specimen is subjected to a cyclic point load at the upper mid-span position with a force of 4,448 kN. The material properties are shown in Table 3.

Specimens are constrained to all degrees of freedom from the left pin position, similar to the experimental specimen. In contrast, the correct pin only has a first degree of freedom of translation UX. For the first configuration of specimen number one, the simulated crack propagation using the developed program and Ansys software is shown in Figure 12 with an identical crack growth trajectory. The crack growth trajectory goes between the lower and middle holes in this specimen and approaches the middle hole from the right side. It implies that the shear stress intensity factor (KII) surrounding the holes has significantly increased, forcing the crack step-sizes to be decreased.

Fig. 12

The crack trajectory using the developed program and Ansys ends up in excellent agreement with experimental data presented in the literature [47, 48], numerical results using the coupled extended meshfree-smoothed meshfree method presented by Knowles and Sternberg [51], and XFEM results obtained by Zhang and Tabiei [21]. Numerical results obtained by Kanth et al. [14] using a polygonal XFEM with numerical integration as shown in Figure 13 A, B, C, D, E, and F, respectively.

Fig. 13

The results of this simulation for dimensionless stress factors by using the developed program and Ansys are shown in Figure 14, which are almost identical to each other. As can be seen in this figure, the SIFs increase in magnitude as the crack approaches the hole, demonstrating the hole's influence on the crack trajectory. For the second configuration of specimen number two, the simulated crack propagation using the developed program and Ansys software is shown in Figure 15 with an almost identical crack growth trajectory. The crack in this specimen propagates above the lower hole to the middle hole's left side.

Fig. 14

Fig. 15

Figure 16 depicts the direction of crack growth for specimen 2, which was compared to the experimental crack trajectory performed in the literature [47, 48] and the numerical findings utilizing the coupled expanded meshfree-smoothed mesh-free technique provided by Knowles and Sternberg [51], as illustrated in Figure 16 A, B, C, and D respectively.

Fig. 16

The predicted values of the two modes of SIFs in the developed program and Ansys are depicted in Figure 17 with excellent agreement. As observed in the figure, the lower hole has no influence on the crack growth direction, since the initial crack was far from the lower hole, while the middle hole has a stronger influence, attracting the crack to sink into the hole.

Fig. 17