Numerical Modeling of Crack Growth under Mixed-Mode Loading

Abstract

The aim of this paper is to simulate the propagation of linear elastic crack in 3D structures using the latest innovation developed using Ansys software, which is the Separating Morphing and Adaptive Remeshing Technology (SMART), in order to enable automatic remeshing during a simulation of fracture behaviors. The ANSYS Mechanical APDL 19.2 (Ansys, Inc., Canonsburg, PA, USA), is used by employing a special mechanism in ANSYS, which is the smart crack growth method, to accurately predict the crack propagation paths and associated stress intensity factors. For accurate prediction of the mixed-mode stress intensity factors (SIFs), the interaction integral technique has been employed. This approach is used for the prediction of the mixed-mode SIFs in the three-point bending beam, which has six different configurations: three configurations with holes, and the other three without holes involving the linear elastic fracture mechanics (LEFM) assumption. The results indicated that the growth of the crack was attracted to the hole and changes its trajectory to reach the hole or floats by the hole and grows when the hole is missing. For verification, the data available in the open literature on experimental crack path trajectories and stress intensity factors were compared with computational study results, and very good agreement was found.

1. Introduction

Fracture mechanics is an efficient tool to optimize mechanical component performance, as failure and fracture are frequent industrial equipment problems. In reality, with regards to the ultimate load capability and progression of pre-existing cracks, computational methods are commonly used to assess the durability of cracked structures. The arbitrariness of the crack trajectory, however, requires appropriate computational methodologies, since the computational instability mechanisms influence the structural dynamic by decreasing strength and speeding up the mechanisms of crack propagation. Many of the simulations frameworks available in the literature are utilized using extended finite element method (XFEM) because of its capability to simulate complicated components, providing accuracy in the estimation of interfacial parameters among various materials. Numerical analyses such as finite element analysis have become an incredibly useful method for modeling the propagation of a crack and for determining associated parameters such as the energy release rate and stress intensity factors (SIFs). Over the last few decades, improvements to the finite element method (FEM) have continuously been made to improve the handling of complex material models and experimental conditions. During the meantime, new approaches and methods in several areas of study have been suggested and developed rapidly, including the FEM, Discrete Element Method (DEM) [1,2,3], Element Free Galerkin (EFG) method [4], XFEM [5,6,7], cohesive element method [8,9], and phase-field method [10]. Similarly to the XFEM, [11] proposes a novel implementation of the phase-field regularized cohesive zone model (PF-CZM), in which the additional degrees of freedom characterizing crack behavior are only associated with nodes within a small sub-domain. Similarities and differences between the gradient-enhanced damage models and phase-field models were discussed in details by [12]. In several applications, understanding interface mechanics is essential for determining and predicting overall material behavior. It is possible to capture the complex behavior of any interface by introducing a quantity called traction density proposed by [12]. Cracks created in structural parts by cyclic loads can move directly or change their direction according to the loading mode. A variety of commercial tools for pre-defined cracks have been developed in either constant or variable amplitude load conditions to evaluate the planar load mode of cracks, for example, AFGROW software for aerospace applications [13] and NASGRO software for fracture and fatigue crack growth analysis [14], which employ handbook solution for the SIFs and suitable fatigue crack growth (FCG) rate equations in the literature. Industrial structures, however, are mainly affected by mixed-mode load conditions that generate non-planar cracking path pathways. The above-mentioned commercial tools do not apply to the specimen experiencing mixed-mode loading. This means that the software modeling the realistic crack growth path and predicting fatigue life in a mixed-mode are in need for the industrial sector. In the case of loading in a mixed-mode, a better approximation of cyclic life requires a cyclic summation method based on an equivalent SIF, which is able to cover both the formation of arbitrary crack types and the effects of load interactions. A Finite Element methodology called XFEM has been introduced by Belytschko and Black [15] to address complications in FEM crack initiation and propagation simulations. In this process, during the crack growth stages, there is no need to re-mesh, reducing the required time for computation and associated error while updating the mesh and mapping the structure. In several works, XFEM was validated and can be accurately designed to model discontinuities such as cracks [16,17]. The implementation of XFEM using ANSYS can accurately model stationary cracks as well as growing cracks. The XFEM has been one of the most recent fields of study in computational mechanics and has attracted many researchers. Its major benefit is that the crack will move through an element without being remeshed. To do this, the information from a crack must be described and updated by a so-called Level Set process. Convenient structures are often subjected to various loading mechanisms, such as stress, shear, and torsion, which lead to mixed-mode interactions. The stress status in front of a crack was therefore highly dependent on the form of interactions in mixed-mode I/II, demonstrating the amount of the stresses at the crack tip. Mode I is the emphasis of scientific literature on the activity of the FCG, while cracks and defects are usually mixed-mode cracking, whereas crack growth in practical engineering structures such as the vessel pressure, pipelines, and fan blades appear to be mixed-mode too (I–II) [18]. The initiation and propagation of crack will usually refer in a complicated state to the governed stress intensity factors (SIFs) [19,20,21,22]. Fatigue cracks in various fields of engineering are important because they are one of the key causes of catastrophic failure. The sustainability of main components and launch systems for operation must be ensured. The reliability of critical components must be ensured in the case of launching systems for service. Consequently, it is primarily important in terms of reliability requirements to precisely predict and assess fatigue life in many industries. In certain applications, for example in the aerospace industry, fatigue studies require experimentation that leads to high costs, thus precise computational processes for crack propagation analysis are essential to predict the direction of crack growth and durability in both static and dynamic loads. Comprehensive works have been performed to establish an efficient model to simulate the propagation of fatigue cracks and fatigue life in order to reduce fatigue failures. Numerical analyses using XFEM implemented by ANSYS APDL.19.2 are an effective method to reduce the computational time and expenses involved with performing experimental studies. SIFs are being used to characterize displacement and stress around the front of the crack and analyze the LEFM crack growth under dynamic loading for isotropic material. According to the parametric study, material anisotropy can have a significant impact on the direction of crack growth [23,24,25,26]. Another study examined the influence of weak spots in composite structures that attract cracks and thereby prevent them from growing, such as the holes effect [27]. The SIF increases to a critical range as the crack grows, following which the structure might start deforming by starting a fracture. The general types of loads that depend on those used in fatigue life predictions attributed to the complex mode of applied loads and the geometry specifications are regularly mixed modes (I/II) [28,29,30,31]. The major motivations for this research were to contribute significantly to the use of ANSYS as an alternative approach for simulating fatigue crack propagation problems under mixed-mode loading, and monitoring the crack growth trajectory in the presence of holes in the geometry [32,33,34,35]. This paper aims to study the effect of the crack tip location related to the presented holes in the trajectory of the crack propagation, SIFs, and fatigue life of a three-point-bending specimen with three holes and six different configurations using the XFEM implemented by ANSYS APDL 19.2. There are several different articles that were studied for this geometry, but no comprehensive analysis implements ANSYS software using the SMART CRACK feature (e.g., [36]).

2. ANSYS Mechanical APDL 19.2

ANSYS Mechanical APDL 19.2 XFEM was applied to model the crack propagation direction for the three-point bending beam with three holes for various configurations. There are three types of cracks that can be modeled within ANSYS: arbitrary, seminal, and pre-meshed cracks. The pre-meshed crack process includes the front of the crack which is used by the “Smart Crack Growth” simulation tool, in which the failure criterion is the SIFs. Instead of using the enrichment area (splitting) of the XFEM, SMART automatically updates the mesh from crack geometry changes due to crack growth on each solution stage, thereby eliminating long pre-processing sessions. The previously generated node sets were assigned to the front of the crack and the top and bottom sides of the crack. The new feature introduced in ANSYS is the tetrahedron mesh-based crack called the smart crack growth, which introduces the “pre-meshed crack” requirement after completion in which the user can choose the crack growth option type. The mesh around the crack tip should be refined using the sphere of influence methodology around the geometrical edge moving across the thickness. SMART can be used with a higher order tetrahedron (entire domain of the model) mesh only. The named geometric regions to be defined are the edge, top surface, and the bottom surface of the crack. Three procedures were primarily used to illustrate the fatigue assessment of materials, which include the fracture mechanics method developed by [37], the method of strain-life introduced by [38], and the stress-life method introduced by [39]. For mixed-mode fatigue crack growth analysis, the maximum circumferential stress criterion was implemented as a crack propagation criterion by the ANSYS software. The following are the equations for the direction of crack propagation in ANSYS [40,41]:

$$\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 ANSYS Mechanical APDL 19.2, it can be found that when using XFEM, the crack growth simulation is restricted to region II of the typical fatigue crack growth behavior graph, which can be defined in the modified formula of Paris law equation as:

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

(2)

From Equation (2), the number of life cycles for fatigue can be estimated with a crack growth increment 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 SIF formula can be formulated as follows [40]:

$$\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$ = range of SIF in mode I loading and $\Delta K_{II}$ = range of SIF in mode II loading.

There are many approaches formulated for evaluating the SIFs based on numerical analysis, including the use of XFEM and other numerical methods. The approaches utilizing the XFEM include the method of displacement extrapolation (DEM), interaction integral method (IIM), stress extrapolation method (SEM), and node displacement method (NDM) [42]. The interaction integral technique is usually the most accurate method and it is able to estimate K_I and K_II separately. ANSYS provides two approaches for the determination of SIFs, which are the DEM and IIM. The second approach was utilized because it is numerically simpler to execute and thus offers greater preciseness and less mesh requirements. This approach depends on the domain integral approach [43] where an auxiliary field is used to separate K_I from K_II, as this ability is missing in the domain integral itself. According to the mixed-mode stress intensity factors K_I, K_II, and K_III, the energy release rate was formulated as [43,44]:

$$J(s)=\frac{K_I^2+K_{II}^2}{E^{*}}+\frac{1+\nu }{E}{K}_{III}^2\hspace{1em}\hspace{1em}{E}^{*}=[\begin{array}{cc}\frac{E}{(1-{\nu }^2)}& \mathrm{Plane} \mathrm{strain}\\ E& \mathrm{Plane} \mathrm{stress}\end{array}.$$

(5)

For the superimposed state, Equation (5) becomes:

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

(6)

where superscript (S) represents the superimposed case, J(s) is the actual state domain integral, J^aux(s) is the auxiliary state domain integral, and I(s) is an integral part of the actual and auxiliary terms that connect.

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

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

(7)

By setting $K_{II}^{aux}=1$ and $K_I^{aux}=K_{III}^{aux}=0$, and the selection $K_{III}^{aux}=1,\hspace{1em}{K}_I^{aux}=K_{II}^{aux}=0$ gives the relationships between K_II, K_III:

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

(8)

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

(9)

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

3. Numerical Results and Discussions

3.1. Three-Point Bending Beam with Different Configurations

The ANSYS Mechanical APDL 19.2’s ability to predict crack paths was verified for various initial crack lengths and crack locations in six different configurations of the three-point bending beam shown in Figure 1. Ingraffea and Grigoriu [45] investigated this example for the first time. The results showed high sensitivity to small changes in the initial position in the crack path. This example has been commonly used to test the predictive ability of the system as a benchmark for numerical models. The considered geometry has 2L = 508 mm, width = 203.2 mm, and thickness = 12.7 mm with three different configurations according to the initial crack length and its position as shown in

Table 1. Configurations of the three-point bending beam specimen.

Specimen No.	Length of Crack, a (mm)	Crack Position, b (mm)	Hole Layout
1	25.4	152.4	With hole
2	25.4	152.4	Without hole
3	38.1	127	With hole
4	38.1	127	Without hole
5	63.5	152.4	With hole
6	63.5	152.4	Without hole

.

For all specimens, the dimensions, positions of the hole, and sizes of the hole are the same. The main differences are the initial length of the crack and its locations as displayed in Table 1. For the first configurations, the length of the crack is a_o = 25.4 mm and is positioned at a distance of d = 152.4 mm from the beam’s mid-length. The second configurations varies with an initial crack length of a_o = 63.5 mm from the first one. The third configurations vary from the first length to the initial crack length of a_o = 38.1 and also include the initial location of the crack with a distance of d = 127 mm from the beam’s mid-length. The specimen was subjected to a cyclic point load acting at a value of 4448 kN at the top mid-span position. The material properties were taken as elasticity modulus, E = 205 GPa, yield stress ${\sigma }_y=516$ MPa, threshold stress intensity factor $\Delta k_{th}=80\hspace{0.17em}\mathrm{MPa}\sqrt{\mathrm{mm}},$ $\Delta K_{IC}=730 \mathrm{MPa}\sqrt{\mathrm{mm}}$, coefficient of Paris’ law, C = 1.2 × 10⁻¹¹, Paris law exponent m = 3, and Poisson’s ratio ν = 0.3. The initial mesh with a moderate mesh density for this geometry is shown in Figure 2 with a number of nodes = 119,977 and a number of elements = 78,959. Each specimen has been simulated as illustrated in Table 1 with and without hole configurations. As in the experimental specimen, specimens are constrained to all degrees of freedom from the position of the left pin, and the right pin only has a first degree of freedom of translation Ux. This case study aims to investigate the Ansys 19.2 APDL’s capabilities for crack propagation prediction and stress intensity factors evaluation with the new smart crack growth features in the mixed-mode loading environment. For each crack growth increment, the angles of the crack extension and the locations of the crack tip were predicted. The initial mesh for the first specimen before starting the simulation is shown in Figure 2.

3.1.1. Specimen 1

The simulated crack propagation in this specimen moves between the lower and middle hole and approaches the middle hole from its right side. This indicates that the shear stress intensity factor (K_II) around the holes is dramatically increased, which forces the step-sizes of crack to be decreased. During propagation, the results of the crack trajectory agree excellently with the experimental data given by [45,46] and with the numerical results when using the coupled extended meshfree–smoothed meshfree method presented by Ma et al. (2020) [47], XFEM results using ABAQUS software obtained by Dirik and Yalçinkaya [36], and numerical results obtained by Huynh et al. (2019) [5] using A polygonal XFEM with numerical integration for linear elastic fracture mechanics as presented in Figure 3a–d, and e, respectively. As shown in this figure, the predicted crack growth path in the present study is closer to the experimental trajectory compared to other numerical results.

The crack propagation trajectory predicted in the present study was also matched to those results for the same geometry and boundary condition using the Dual Boundary Element Method (DBEM) introduced by Price and Trevelyan [48] and to the numerical result reported by [49] using the Boundary element method (BEM), as shown in Figure 4 (left). The predicted crack growth trajectory in the present study as shown in Figure 4 (right) was more closely matching the experimental crack growth path predicted by Bittencourt et al. (1996) [46] than to those predicted in Figure 4 (left).

The findings of this simulation for dimensionless stress factors were compared with those from XFEM using the smooth nodal stress technique by [50], as shown in Figure 5 with identical results. It is observed that as the crack approaches the hole, SIFs tend to change at a greater amplitude. Figure 6 shows the equivalent Von-Mises stress distribution for specimen 1 with enlargement of the area around the crack tip.

3.1.2. Specimen 2

Figure 7 displays the final stage of crack growth trajectory modeled by ANSYS for specimen 2, which has the same configurations as specimen 1 except for the presence of holes. As opposed to the geometry with holes, the influence of the holes was observed. Therefore as demonstrated in the experiments conducted by [45,46], holes act as crack stoppers and attract cracks to grow towards themselves. Figure 8 presented the comparison for the crack growth trajectory with the experimental path predicted by Ingraffea and Grigoriu (1990), [45] with an identical crack trajectory.

3.1.3. Specimen 3

The crack in this specimen Grow above the lower hole and comes to an end on the left side of the middle hole. The direction of the crack growth for specimen 3 is seen in Figure 9 and compared to the experimental crack trajectory conducted by Bittencourt et al. (1996) [46], as shown in Figure 10 (left) and with the numerical results using the coupled expanded meshfree-smoothed meshfree method provided by Ma et al. (2020), [47] as seen in Figure 10 (right) with close agreement.

The results of the present study for the dimensionless stress intensity factors were compared with those results obtained by Peng et al. (2017), [50] using XFEM with a smooth nodal stress method, as shown in Figure 11 with comparable results. Furthermore, Figure 12 shows the equivalent Von-Mises stress distribution for specimen 3.

3.1.4. Specimen 4

This specimen has the same outside dimensions as specimen 3 except for the presence of the three holes inside the geometry. Figure 13 shows the predicted crack growth trajectory using ANSYS. This crack trajectory was compared to the experimental path predicted by Ingraffea and Grigoriu (1990) [45] as seen in Figure 14. The predicted crack growth trajectory was identical to the experimental trajectory.

3.1.5. Specimen 5

The geometry in this specimen is the same as specimen 2, but with a different crack length of 63.5 mm. Figure 15 demonstrates the final stage of the crack propagation trajectory from ANSYS compared with numerical results obtained by Peng et al. (2017), [50] using smooth nodal stresses in the XFEM and a double interpolation method with XFEM (XDFEM), as well as with XFEM results using ABAQUS software obtained by Dirik and Yalçinkaya [36] as shown in Figure 15a–c, respectively. In Figure 16a, comparison of the simulated crack tip coordinates and experimental path [46] was presented for this specimen configurations with a closer crack path. Figure 17 also displayed the last stage of the crack propagation for this specimen from ANSYS simulation along with Von-misses stress contours.

For this specimen, the predicted SIFs K_I and K_II were compared with the numerical results using the dual boundary method obtained by Andrade and Leonel (2020), [51] as shown in Figure 18 and Figure 19 with excellent agreement. At first, mixed-mode crack growth was developed by the crack due to the higher value of the second mode of SIF (K_II). After that, mode-I propagation took precedence after the initial increments, before the crack reached the holes. In these instants, the second mode of SIF (K_II) differed substantially, resulting in changes in the direction of the crack, as shown by Figure 19.

3.1.6. Specimen 6

The exterior dimensions of this specimen are the same as those of specimen 5 without the insider holes. Figure 20 displays the predicted crack growth path using ANSYS. This crack trajectory was compared to the experimental path predicted by [45], as shown in Figure 21. The predicted crack growth trajectory was similar to the trajectory of the experiment.

In these simulation sequences, holes act as a crack stopper and attract a crack trajectory to growth, as reported in the experimental study of the specimens by [45].

4. Conclusions

A numerical analysis of incremental crack growth was conducted using ANSYS Mechanical APDL for the three-point bending beam. For various specimens of the three-point bending beam with different configurations based on the original crack length and its position from the mid-spam, the results of the ANSYS simulation were compared with experimental and numerical data, with excellent agreement found for all cases. Six sample configurations with and without hole configurations were simulated in total. The ability to determine the crack direction has been shown through this sequence of simulations. The hole acts as a crack stopper in these simulation sequences and attracts a crack path to growth. Such findings reinforce that the algorithm can be used in damage tolerance designs to identify crack-stopping holes.