Extension of the average strain energy density method for crack initiation prediction in functionally graded materials and identification of a toughening mechanism

Abstract

This study extended the average strain energy density theory, traditionally applied for notch analysis in isotropic materials, to predict mixed mode I/II fracture in functionally graded materials. The effectiveness of existing experimental criteria designed for orthotropic materials is examined for their capacity to predict crack initiation in graded materials. A toughening mechanism based on fiber bypassing during crack propagation is introduced and investigated through experimental methods. Scanning electron microscopy images are employed to elucidate the role of this mechanism in enhancing the fracture resistance of chopped/short fiber-reinforced composites. The findings reveal that the average strain energy density method demonstrates strong correlation and agreement with experimental–numerical data, affirming its robustness for predicting crack initiation in mixed mode fracture scenarios.

Graphical abstract

Keywords

crack / functionally graded materials / average strain energy density / mixed mode I/II fracture criterion / toughening mechanism / crack propagation

1 Introduction

Functionally graded materials (FGMs) represent a remarkable class of advanced composite materials characterized by a continuous, spatial variation in composition and microstructure. This unique feature enables FGMs to exhibit tailored mechanical and thermal properties, making them particularly advantageous in high-performance applications across various industries, including aerospace, biomedical engineering, and electronics. However, as with many advanced materials, FGMs are not immune to structural imperfections. The presence of cracks can pose significant challenges, compromising both their integrity and functionality under operational stresses [1,2].

The study of crack behavior in FGMs is essential to understanding their overall performance and reliability. Cracks in these materials can propagate in complex manners due to the non-homogeneous distribution of material properties, leading to intricate fracture behavior that deviates from traditional, homogeneous materials. Consequently, the development of accurate fracture criteria for FGMs is crucial for predicting crack behavior and designing more resilient materials.

Wang and Zhou [3] employed the specific finite element method to simulate crack propagation paths in FGMs, tuning model parameters at the crack tip. They used the maximum tangential stress (MTS) criterion to estimate crack initiation directions and assumed that the stress field around cracks in FGMs behaves similarly to that in isotropic materials.

Hossein Izadi et al. [4] introduced an innovative technique for producing FGMs utilizing randomly distributed chopped fibers, primarily designed for experimental investigations. They characterized the mechanical behavior of these FGMs, including properties such as elastic modulus and fracture toughness. In subsequent studies [5,6], the authors proposed a new mixed mode I/II fracture criteria for the fabricated FGM, incorporating the fracture process zone (FPZ) and using strain energy release rate theory. The authors utilized the concept of effective fracture toughness and damage factor to account for the damage zone effect. Their results indicated that while the addition of chopped fibers leads to a consistent increase in elastic modulus, it does not necessarily enhance the mode I fracture toughness.

Pan et al. [7] carried out a probabilistic analysis of mixed mode fracture behavior in FGMs, focusing on an internal slant crack. Their study explored how stochastic variations in microstructural features influence bulk mechanical responses and produced analytical expressions for stress intensity factors (SIFs). The results underscored the crucial role of microstructural randomness in shaping the statistical distributions of shear modulus and SIFs, particularly in the presence of large slant cracks.

Pandey and Patel [8] analyzed the quasi-static crack growth in a metal-ceramic FGM, observing stable crack growth from the stiff side to the compliant side and unstable growth in the opposite direction. The evaluation of fracture toughness and crack tip opening displacement increased with distance from the stiff side, aligning with Raveendran’s micromechanics predictions based on alumina weight fractions. Bao et al. [9] examined crack propagation and damage mechanisms in functionally graded ultra-high-performance cementitious composites with varying layers and fiber types, finding that a three-layer design enhanced flexural strength and deflection capacity. Analysis via digital image correlation (DIC) revealed a nonlinear crack opening displacement profile and differences in acoustic emission energy growth rates among the layers, indicating improved performance with mixed fibers.

Ren et al. [10,11] developed a variational damage model to simulate fracture in solids, incorporating a threshold to prevent damage under low energy and enabling automatic crack evolution. The approach yields sharper crack interfaces than conventional phase-field methods and is compatible with standard finite element (FE) frameworks. Their model effectively captures complex crack behavior in both two and three dimensional (2D and 3D) static and dynamic fracture scenarios. In another study [12], the authors proposed a dual-horizon peridynamic model enhanced with variational damage to simulate dynamic brittle fracture without explicit bond-breaking, improving numerical stability. By linking damage to the positive strain energy density (SED) through spectral decomposition, the model ensures realistic crack behavior and avoids interpenetration upon crack closure. Its effectiveness is demonstrated in both 2D and 3D fracture simulations.

While research on crack initiation and propagation in FGMs remains limited, extensive research has been conducted on developing analytical methods for crack growth prediction in orthotropic composites, such as wood and composites made from continuous glass and carbon fibers.

Jernkvist [13] formulated a mixed mode I/II fracture model for wood, considering crack propagation parallel and perpendicular to the grain direction, within the framework of linear elastic fracture mechanics. It established that cracks generally propagate along the fibers, allowing a unified approach regardless of the crack’s initial orientation or mode mixity. Daneshjoo et al. [14] introduced a mixed mode I/II failure model tailored for orthotropic materials, which incorporates the energy absorbed within the FPZ and builds upon the principles of SED approach. The new criterion demonstrated better compatibility with failure phenomena in laminated composites and wood species, showing strong agreement with experimental data and outperforming existing criteria. Kumar et al. [15] examined the influence of plate curvature on the mixed mode SIF of multiple cracks in riveted joints under fatigue loading, using both numerical and experimental methods. It was found that higher SIF values occurred in curved panels compared to flat panels, particularly for lower crack depth ratios, due to the “flattening” effect under loading. Experiments were conducted on Al 2024-T3 specimens, confirmed a higher crack growth rate in curved panels attributed to additional bending stresses at the rivet holes, demonstrating a strong correlation with the numerical findings.

Pitti et al. [16] introduced a mixed mode fracture specimen designed to ensure stable crack growth during creep loading while accounting for viscoelastic behavior under mixed mode conditions. Utilizing the M-integral approach within FE software, the authors optimized the specimen’s geometry to effectively evaluate fracture parameters, viscoelastic properties, and creep crack growth across various mixed mode ratios. Golewski [17] investigated the fracture toughness of concrete incorporating 0%, 20%, and 30% class F fly ash under mode I loading, determining critical SIFs and crack tip opening displacements using DIC. Results indicated that concrete with 20% fly ash exhibited high fracture toughness, while that with 30% showed low toughness, highlighting the effectiveness of the DIC method in analyzing crack initiation behavior.

Hua et al. [18] investigated the fracture behavior of cracked materials through modified mixed mode I-II fracture criteria that include T-stress considerations. A comprehensive comparison of experimental results from five different cracked configurations demonstrated notable variations in predictive accuracy due to discrepancies in T-stress magnitudes and signs around the crack tips. The study identified specific criteria that effectively predicted fracture behavior for various specimen types, highlighting that generalized MTS and related criteria provided reliable predictions depending on the loading conditions and T-stress characteristics.

Dong et al. [19] used DIC and numerical simulation to study mixed mode fracture in concrete beams, analyzing FPZ evolution. Their findings showed that the crack opening-to-sliding displacement ratio remained relatively constant before peak load, and FPZ length was influenced by specimen size and mode I/II SIF ratio. Specifically, shorter ligament lengths or lower mode I/II ratios prevented full FPZ development.

Studying crack propagation and trajectory is crucial for understanding the behavior of materials under different loading conditions [20]. By investigating the mechanisms that occur during crack path, researchers can develop strategies to prevent catastrophic failures and improve material durability [21]. Numerous studies have been conducted in various fields to address this challenge.

Braun and Ariza [22] enhanced an existing linear elastic lattice framework for anisotropic materials by incorporating damage progression through a softening constitutive relationship, enabling its application to composite materials. They verified the approach by analyzing dynamic fracture behavior in unidirectional carbon fiber composites, focusing on impact-loaded rectangular notched beams. Their findings, reported through the evolution of crack patterns, propagation velocity, and crack length, demonstrated good consistency with earlier experimental and computational results. Ricoeur et al. [23] analyzed the anisotropy in elastic properties and crack growth resistance of short fiber reinforced composites, noting its significant effect on crack deflection and path predictions. It determined that cracks in the transverse direction exhibit the highest toughness, influenced by fiber orientations shaped during injection molding and described using Gaussian random field models. The research also explored factors affecting crack path predictions from FE simulations and stochastic behaviors related to bifurcation phenomena under mode I loading.

Lin and Li [24] and Li et al. 25] developed a new bridging model for cement-based composites reinforced with randomly oriented discontinuous flexible fibers, incorporating slip-hardening interfacial shear stress. They theoretically derived the bridging stress against crack opening diagram and discussed its implications for properties such as tensile strength and ductility, using a 2% polyethylene fiber reinforced cement composite as an example. This model effectively captured the slip-hardening behavior and accurately predicted increased toughness and fracture energy, addressing limitations in existing constant interface shear stress models.

Tiu et al. [26] explored the anisotropic fracture behavior of layered structures in short-fiber particulate-reinforced composites compared to conventional flowable particulate-reinforced composites. Through three-point bending (3PB) tests on monolithic and bilayer specimens, the researchers developed R-curves, revealing that crack deflection and fiber bridging substantially enhanced fracture resistance, especially in specimens with aligned fibers.

Fakoor et al. [27] developed a theoretical model based on the reinforcement isotropic solid approach to evaluate fracture parameters in fiber-reinforced composites at 0° and 90° fiber-crack orientations. The model was validated experimentally, and a new expression was introduced to relate fracture toughness to the damage volume. Their results showed higher crack resistance when the crack runs perpendicular to the fiber direction.

Jiao et al. [28] investigated crack initiation and energy evolution in fractured sandstones under uniaxial compression, finding that internal fractures significantly reduce rock strength. They introduced a new parameter, the initiation factor Z, and established its relationship to crack propagation and fracture toughness. Their findings enhance the understanding of mechanical responses in fractured rock and provide insights for underground engineering applications.

Despite the significant role of FGMs in various industries and the clear necessity to investigate their fracture behavior, in this study the application of the average strain energy density (ASED) theory, commonly used for notch analysis in isotropic materials, has been extended to predict mixed mode I/II fracture in FGMs. The validity of the proposed criterion has been examined through comparison with experimental data reported by Hossein Izadi et al. [5]. Additionally, the effectiveness of experimental criteria developed for orthotropic materials has been examined for predicting crack initiation in FGMs. Moreover, this study introduces a mechanism involving bypassing fibers, which serves as a toughening mechanism during crack propagation. This mechanism has been investigated through experimental methods, employing images acquired from scanning electron microscopy (SEM) to elucidate its role in enhancing the material’s fracture resistance.

2 Theoretical studies

2.1 Main assumptions

The position of the crack tip in FGMs is one of the most important parameters in predicting the fracture behavior of them during crack initiation. In accordance with the application of FGMs and the fabricated graded composite described by Hossein Izadi et al. [5], the crack tip is assumed to be located on the most brittle side of the material along the gradation profile (Fig. 1). The fabricated FGM consisted of an epoxy matrix reinforced with chopped glass fibers, and the gradation was achieved by varying the fiber volume fraction across the specimen, specifically at 0%, 20%, 40%, and 50%. Since the epoxy matrix exhibits lower fracture resistance compared to the reinforcing fibers, crack propagation following initiation predominantly occurs within the matrix (Fig. 2). Therefore, fracture theories originally developed for isotropic materials, such as the ASED method, can be appropriately applied to predict crack initiation, as the crack tends to grow in the isotropic phase of the material.

2.2 Average strain energy density criterion for functionally graded materials

The ASED method is traditionally employed for predicting mixed mode I/II fracture in notched specimens. However, given its robustness, simplicity, and proven predictive accuracy, extending this well-established approach to cracked materials, is valuable.

According to the ASED theory, fracture occurs when the average value of SED over a given control volume \overline{W} reaches a critical value W_{\mathrm{C}}. This critical value is a property of the material and does not depend on the geometry of the notched specimen.

The general form of stress distribution ahead of an U-shaped notch is as follows (Fig. 3(a)):

(1){\sigma }_{ij}=\frac{r^{{\lambda }_1-1}}{\sqrt{2\pi }}{K}_{\mathrm{I}}^{\mathrm{N}}{f}_{ij}\left(\theta \right)+\frac{r^{{\lambda }_2-1}}{\sqrt{2\pi }}{K}_{\mathrm{I}\mathrm{I}}^{\mathrm{N}}{g}_{ij}\left(\theta \right),

where the f_{ij}\left(\theta \right) and g_{ij}\left(\theta \right) are angular functions, and are defined in the following way [29,30]:

(2)\begin{array}{rl}{f}_{\theta \theta }=& \frac{1}{\left(1+{\lambda }_1\right)+{\chi }_1\left(1-{\lambda }_1\right)}\\ & \cdot \left[\left(1+{\lambda }_1\right)\cos \left(1-{\lambda }_1\right)\theta +{\chi }_1\left(1-{\lambda }_1\right)\cos \left(1+{\lambda }_1\right)\theta \right],\\ f_{rr}=& \frac{1}{\left(1+{\lambda }_1\right)+{\chi }_1\left(1-{\lambda }_1\right)}\\ & \cdot \left[\left(3-{\lambda }_1\right)\cos \left(1-{\lambda }_1\right)\theta +{\chi }_1\left(1-{\lambda }_1\right)\left(-\cos \left(1+{\lambda }_1\right)\theta \right)\right],\\ f_{r\theta }=& \frac{1}{\left(1+{\lambda }_1\right)+{\chi }_1\left(1-{\lambda }_1\right)}\\ & \cdot \left[\left(1-{\lambda }_1\right)\sin \left(1-{\lambda }_1\right)\theta +{\chi }_1\left(1-{\lambda }_1\right)\sin \left(1+{\lambda }_1\right)\theta \right],\\ g_{\theta \theta }=& \frac{1}{\left(1-{\lambda }_2\right)+{\chi }_2\left(1+{\lambda }_2\right)}\\ & \cdot \left[\left(-1+{\lambda }_2\right)\sin \left(1-{\lambda }_2\right)\theta +{\chi }_2\left(1+{\lambda }_2\right)\left(-\sin \left(1+{\lambda }_2\right)\theta \right)\right],\\ g_{rr}=& \frac{1}{\left(1-{\lambda }_2\right)+{\chi }_2\left(1+{\lambda }_2\right)}\\ & \cdot \left[\left(-3-{\lambda }_2\right)\sin \left(1-{\lambda }_2\right)\theta +{\chi }_2\left(1+{\lambda }_2\right)\sin \left(1+{\lambda }_2\right)\theta \right],\\ g_{r\theta }=& \frac{1}{\left(1-{\lambda }_2\right)+{\chi }_2\left(1+{\lambda }_2\right)}\\ & \cdot \left[\left(1-{\lambda }_2\right)\cos \left(1-{\lambda }_2\right)\theta +{\chi }_2\left(1+{\lambda }_2\right)\cos \left(1+{\lambda }_2\right)\theta \right].\end{array}

Parameters K_{\mathrm{I}}^{\mathrm{N}} and K_{\mathrm{I}\mathrm{I}}^{\mathrm{N}} are notch stress intensity factors, {\lambda }_1 and {\lambda }_2 represent the Williams’ eigenvalues and also, {\chi }_1 and {\chi }_2 are dependent parameters on the opening angle of notch. The SED has a linear form and is defined as follows, where E\left(x\right) and v\left(x\right) represent the elastic modulus and Poisson’s ratio at the crack tip and varies depending on the position of the crack in the FGM [31,32].

(3)W\left(r,\theta \right)=W_1\left(r,\theta \right)+W_{12}\left(r,\theta \right)+W_2\left(r,\theta \right),

(4)\begin{array}{rl}W\left(r,\theta \right)=& \frac{1}{2E\left(x\right)}[{\sigma }_{11}^2+{\sigma }_{22}^2+{\sigma }_{33}^2-2v\left(x\right)\\ & \cdot \left({\sigma }_{11}{\sigma }_{22}+{\sigma }_{11}{\sigma }_{33}+{\sigma }_{22}{\sigma }_{33}\right)+2\left(1+v\left(x\right)\right){\sigma }_{12}^2].\end{array}

In 2D and polar condition, Eq. (4) is simplified to the following form:

(5)W\left(r,\theta \right)=\frac{1}{2E\left(x\right)}\left[{\sigma }_{\theta \theta }^2+{\sigma }_{rr}^2-2v\left(x\right){\sigma }_{\theta \theta }{\sigma }_{rr}+2\left(1+v\left(x\right)\right){\sigma }_{r\theta }^2\right].

By substituting the stress equations in Eq. (5), the SED relationship can be organized and rewritten in the following way:

(6)W\left(r,\theta \right)=\frac{1}{2E\left(x\right)}\left[A_{11}{\left(K_{\mathrm{I}}^{\mathrm{N}}\right)}^2+A_{12}{K}_{\mathrm{I}}^{\mathrm{N}}{K}_{\mathrm{I}\mathrm{I}}^{\mathrm{N}}+A_{22}{\left(K_{\mathrm{I}\mathrm{I}}^{\mathrm{N}}\right)}^2\right],

where the A_{ij}(\theta ) coefficients are defined as follows:

(7)A_{11}=\frac{r^{2\left({\lambda }_1-1\right)}}{2\pi }\left[f_{\theta \theta }^2+f_{rr}^2-2v\left(x\right)f_{rr}{f}_{\theta \theta }+2\left(1+v\left(x\right)\right)f_{r\theta }^2\right],

(8)\begin{array}{rl}{A}_{12}=& \frac{r^{{\lambda }_1+{\lambda }_2-2}}{\pi }[f_{\theta \theta }{g}_{\theta \theta }+f_{rr}{g}_{rr}-2v\left(x\right)\left(f_{\theta \theta }{g}_{rr}+f_{rr}{g}_{\theta \theta }\right)\\ & +2\left(1+v\left(x\right)\right)f_{r\theta }{g}_{r\theta }],\end{array}

(9)A_{22}=\frac{r^{2\left({\lambda }_2-1\right)}}{2\pi }\left[g_{\theta \theta }^2+g_{rr}^2-2v\left(x\right)g_{\theta \theta }{g}_{rr}+2\left(1+v\left(x\right)\right)g_{r\theta }^2\right].

The elastic deformation energy within a region with a radius of R surrounding the notch tip is:

(10)\begin{array}{rl}G(R)=& {\displaystyle {\int }_{A}W\mathrm{d}A={\displaystyle {\int }_0^R{\displaystyle {\int }_{-\gamma }^{+\gamma }[W_1\left(r,\theta \right)+W_{12}\left(r,\theta \right)}}}\\ & +W_2\left(r,\theta \right)]r\mathrm{d}r\mathrm{d}\theta .\end{array}

Due to the symmetry of the integration field around the notch bisector, the contribution of W_{12}\left(r,\theta \right) would vanish and the elastic deformation can be rewritten as follows:

(11)\begin{array}{rl}G\left(R\right)=& G_1\left(R\right)+G_2\left(R\right)\\ =& \frac{1}{E\left(x\right)}\left[\frac{I_1\left(\gamma \right)}{4{\lambda }_1}{\left(K_{\mathrm{I}}^{\mathrm{N}}\right)}^2{R}^{2{\lambda }_1}+\frac{I_2\left(\gamma \right)}{4{\lambda }_2}{\left(K_{\mathrm{I}\mathrm{I}}^{\mathrm{N}}\right)}^2{R}^{2{\lambda }_2}\right],\end{array}

where I_1\left(\gamma \right) and I_2\left(\gamma \right) are:

(12)I_1\left(\gamma \right)={\int }_{-\gamma }^{+\gamma }\frac{1}{2\pi }\left[f_{\theta \theta }^2+f_{rr}^2-2v{f}_{rr}{f}_{\theta \theta }+2\left(1+v\right)f_{r\theta }^2\right]\mathrm{d}\theta ,

(13)I_2\left(\gamma \right)={\int }_{-\gamma }^{+\gamma }\frac{1}{2\pi }\left[g_{\theta \theta }^2+g_{rr}^2-2v{g}_{rr}{g}_{\theta \theta }+2\left(1+v\right)g_{r\theta }^2\right]\mathrm{d}\theta .

In the aforementioned relationships, the variable \gamma is in radians and the area on which the integration is carried out is:

(14)A\left(R\right)={{\displaystyle \int }}_0^R{{\displaystyle \int }}_{-\gamma }^{+\gamma }r\mathrm{d}r\mathrm{d}\theta =R^2\gamma .

Upon averaging over the area A\left(R\right), the resulting value of elastic deformation energy is found to be:

(15)\overline{W}=\frac{G\left(R\right)}{A\left(R\right)}=\frac{1}{E\left(x\right)}\left[e_1{\left(K_{\mathrm{I}}^{\mathrm{N}}\right)}^2{R}^{2\left({\lambda }_1-1\right)}+e_2{\left(K_{\mathrm{I}\mathrm{I}}^{\mathrm{N}}\right)}^2{R}^{2\left({\lambda }_2-1\right)}\right].

e_1 and e_2 represent the geometric constants and are defined as follows:

(16)e_1(2\alpha )=\frac{I_1(\gamma )}{4{\lambda }_1\gamma },

(17)e_2(2\alpha )=\frac{I_2(\gamma )}{4{\lambda }_2\gamma }.

In a cracked specimen, the parameters K_{\mathrm{I}}^{\mathrm{N}} and K_{\mathrm{I}\mathrm{I}}^{\mathrm{N}}, which represent the mode I/II SIFs of the notch, are converted to K_{\mathrm{I}}\left(x\right) and K_{\mathrm{I}\mathrm{I}}\left(x\right), which denote the mode I/II SIFs of the crack at point x. So, the ASED in a control volume around the crack tip can be written as follows:

(18)\overline{W}=\frac{1}{E\left(x\right)}\left[\frac{e_1{K}_{\mathrm{I}}^2\left(x\right)}{R_{\mathrm{C}}^{2\left(1-{\lambda }_1\right)}}+\frac{e_2{K}_{\mathrm{I}\mathrm{I}}^2\left(x\right)}{R_{\mathrm{C}}^{2\left(1-{\lambda }_2\right)}}\right].

To extend the applicability of the ASED method to cracked specimens, the specific case is considered in which the notch opening angle approaches zero. In this limiting condition, a sharp crack can be viewed as a special case of a notch with a vanishing opening angle, allowing the theoretical framework of ASED, originally formulated for notched geometries, to be adapted for crack tip analysis. In the specific case of cracked specimens, i.e., 2\alpha =0 (Fig. 3(b)), the values of geometric constants and Williams’ eigenvalues are equal to e_1=0.134, e_2=0.341 and {\lambda }_1={\lambda }_2=0.5 [33]. Therefore, Eq. (18) will be changed in the following way:

(19)\overline{W}=\frac{1}{E\left(x\right)}\left[\frac{\left(0.134\right)K_{\mathrm{I}}^2\left(x\right)}{R_{\mathrm{C}}\left(x\right)}+\frac{\left(0.341\right)K_{\mathrm{I}\mathrm{I}}^2\left(x\right)}{R_{\mathrm{C}}\left(x\right)}\right].

A circular control volume is considered around the crack tip, having a radius R_{\mathrm{C}}\left(x\right) and can be defined using the following expression [34]:

(20)R_{\mathrm{C}}\left(x\right)=\frac{\left(1+v\left(x\right)\right)\left(5-8v\left(x\right)\right)}{4\pi }{\left(\frac{K_{\mathrm{I}\mathrm{C}}\left(x\right)}{{\sigma }_{\mathrm{t}}\left(x\right)}\right)}^2.

As mentioned earlier, according to the ASED theory, crack initiation occurs when the average value of SED over a given control volume reaches a critical value. This critical value can be determined in terms of tensile strength and elastic modulus of the material at the crack tip:

(21)W=W_{\mathrm{C}}.

(22)W_{\mathrm{C}}=\frac{{\sigma }_t^2\left(x\right)}{2E\left(x\right)}.

By substituting Eq. (22) into Eq. (21) and considering Eqs. (19) and (20), the ASED criterion for FGMs, will be obtained as follows:

(23)\begin{array}{rl}& \frac{8\pi }{\left(1+v\left(x\right)\right)\left(5-8v\left(x\right)\right)}\left(\left(0.134\right)K_{\mathrm{I}}^2\left(x\right)+\left(0.341\right)K_{\mathrm{I}\mathrm{I}}^2\left(x\right)\right)\\ =& K_{\mathrm{I}\mathrm{C}}^2\left(x\right).\end{array}

This energy-based criterion offers a method to forecast the moment when cracks initiate in FGMs. The ASED criterion’s validity for FGMs was assessed by comparing its predicted fracture limit curve to experimental data, including critical mixed mode SIFs and mode II fracture toughness, obtained from Hossein Izadi et al.’s study [5]. In this study, the critical mixed mode I/II SIFs (for a 3 mm crack in the no-fibers section at angles of 0°, 15°, 30°, 45°, and 60°) and the pure mode II fracture toughness of the fabricated fiber-reinforced FGM were obtained using 3PB and four-point bending (4PB) tests, respectively. In this study, after experimentally determining the fracture loads through 3PB and 4PB tests, these loads were applied as boundary conditions in a FE model of the cracked FGM. The FE model accurately represented the geometry, material gradation, and loading configuration of the tested specimens. Using this model, the critical mixed mode SIFs were calculated at the crack tip corresponding to the onset of fracture. This combined experimental–numerical approach enabled a precise determination of fracture parameters, which were subsequently used to evaluate the accuracy of the proposed ASED-based fracture criterion. The good agreement between the ASED predictions and experimental results suggests that the ASED method effectively predicts fracture behavior in FGMs (Fig. 4).

3 Experimental criteria

Experimental criteria are highly valuable despite being costly. One of the key advantages of these criteria is their ability to provide accurate predictions for specific types of materials, regardless of the complex failure mechanisms involved, thus reducing computational requirements. Many of these criteria, developed decades ago, rely on curve fitting methods that often fix the initial and end points of the curves, which can be considered a limitation. However, their proven ability to predict fractures accurately underscores their value in materials research. Also, all of these criteria are developed based on the nonlinear elastic fracture mechanics as follows:

(24)f_{\mathrm{C}}\left(K_{\mathrm{I}},K_{\mathrm{I}\mathrm{I}},\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{i}\mathrm{e}\mathrm{s}\right)=0.

Due to the complex failure mechanisms of orthotropic materials, there is a greater interest in developing empirical criteria specific to this material type. Table 1 lists some of the most significant and well-known experimental fracture criteria for orthotropic materials.

To assess the applicability of these criteria to the FGMs, it is sufficient to substitute the extracted values of K_IC and K_IIC into these criteria, plot them, and then compare the mixed mode behavior of the material obtained from experimental tests. Figure 5 illustrates the experimental fracture limit curves in comparison with the experimental results of critical SIFs used in the previous section.

Referring to Fig. 5, it is important to note that the Hunt and Crouger criterion coincides with the Wu criterion, both of which demonstrate relatively good predictions of the material’s mixed mode fracture behavior. The difference between the coefficients used in their respective formulations is minimal (approximately 0.005), leading to nearly identical fracture limit curves. The Spencer and Barnby criterion, with its adjustable coefficients, effectively aligns with experimental critical SIFs for low mode II values. The Mall and Murphy criterion accurately predicts material mixed mode fracture behavior, but its accuracy decreases slightly for higher mode II values. On the contrary, the Leicester criterion is conservative and despite providing correct predictions, does not accurately forecast the material’s behavior, potentially leading to increased design costs. The Leicester criterion suggest that the material is weak against crack initiation, while literally it exhibits greater resistance against mixed mode fracture. In comparison to the other experimental criteria, the Spencer and Barnby criterion shows a greater material’s resistance to crack initiation, implying that this criterion is less conservative than the others.

The differences observed between the experimental pure mode II results and the plotted diagrams, as well as other experimental results, can be attributed to the change in test method from 3PB (for mixed mode I/II tests) to 4PB (for pure mode II test). The material’s fracture behavior against mixed mode fracture, which is dependent on the type of loading, changes with the test method alteration. As a result, the experimental criterion diagrams plotted based on the 3PB results may not accurately predict the pure mode II SIF obtained from the 4PB test.

It can be predicted that for the results obtained from the 3PB tests, the critical mode II SIF remains relatively constant without significant changes once the crack angle reaches 60°. It can be inferred that the K_II value obtained at a 60° angle is relatively equivalent to K_IIC in this material, indicating that under 3PB loading the material’s fracture behavior with a 60° crack is comparable to its behavior under pure shear loading conditions.

4 Microstructure study

To investigate the fracture mechanism of the fabricated FGM in Hossein Izadi et al.’s study [5], SEM images were taken from the mixed mode bending specimens tested under 3PB loading. SEM is a powerful imaging technique widely utilized in various scientific disciplines for high-resolution surface analysis of materials, offering exceptional magnification capabilities to visualize structures at the nanoscale level with remarkable detail. To enhance the conductivity of the non-conductive fabricated FGM, specimens with different crack angles were coated with gold. Gold is commonly used as a coating material in SEM imaging to address the non-conductivity of samples, preventing charge accumulation that may cause image distortions or sample damage. The application of a thin layer of gold coating enables better imaging and analysis in SEM, enhancing resolution and contrast in the captured images. The SEM images were captured using the QUANTA 450 to study crack propagation behavior from various crack angles. The images obtained indicate that the bonding between the components of the graded material is highly effective, suggesting that the manufacturing process of the FGM was carried out successfully. Furthermore, no defects were observed in the structure of the specimens.

As illustrated in Fig. 6, when the crack extension reaches a fiber colony, it exhibits a distinct behavior: it rotates and bypasses the fibers before continuing in its main direction. This deviation can be attributed to the significantly higher strength of the fibers compared to the matrix material. The crack initially tends to propagate along the weaker matrix bed. However, upon encountering a cluster of fibers, it faces substantial resistance, preventing further movement in that direction. This resistance arises from the strong interfacial bonding between the fibers and the matrix, as well as the load transfer mechanisms that occur at the fiber-matrix interface.

Consequently, the crack shifts to a path of lower resistance, deviating from its original trajectory and continuing its propagation. Once past the fiber cluster, the crack reverts to its original direction and resumes its progression (Fig. 7), demonstrating the influence of fiber strength and interface characteristics on crack path deviation. This phenomenon highlights the intricate interaction between the crack and the fiber-matrix interface, showcasing the role of fiber strength in guiding crack propagation within the material. The crack deflection mechanism observed in Figs. 6 and 7, is a crucial aspect of the toughening behavior in fiber-reinforced composites, contributing to the increased fracture toughness and resistance to crack growth. This action dissipates the energy applied to the crack, which, in turn, slowdowns crack propagation and increases the material’s overall strength.

5 Conclusions

The ASED theory, traditionally used for notch fracture analysis, was extended in this study to investigate crack initiation in FGMs. The ability of other experimental criteria to predict crack initiation in FGMs was also examined. Key findings of this study are as follows.

1) The ASED theory demonstrates good agreement with experimental results and shows potential for predicting mixed mode I/II fracture in FGMs.

2) The Hunt and Crouger criterion and the Wu criterion demonstrate appropriate mixed mode fracture prediction for the investigated FGM.

3) The Spencer and Barnby criterion is the least conservative, while the Leicester criterion is the most conservative among the investigated experimental criteria.

4) The accuracy of the Spencer and Barnby criterion’s prediction decreases as the mode II SIF increases.

5) The crack tip parameters, particularly critical SIF values, are dependent on the specimen geometry and loading conditions.

6) Cracks tend to propagate through the matrix rather than the fibers due to the higher strength of the fibers compared to the matrix.

7) When the crack reaches a fiber cluster, it rotates and bypasses the fibers by shifting to a path of lower resistance, deviating from its original trajectory and then reverts to its initial direction.

References

[1] Tran T V , Nguyen-Xuan H , Zhuang X . Investigation of crack segmentation and fast evaluation of crack propagation, based on deep learning. Frontiers of Structural and Civil Engineering, 2024, 18(4): 516–535

[2] Suo Y , Chen Z , Rahman S . Mixed-mode fracture behaviour of semicircular bend shale with bedding layer. Arabian Journal for Science and Engineering, 2021, 46(7): 6967–6978

[3] Wang L F , Zhou X P . Fracture analysis of functionally graded materials by the field-enriched finite element method. Engineering Fracture Mechanics, 2021, 253: 107875

[4] Hossein Izadi S M , Fakoor M , Mirzavand B . Characterization of elastic modulus and fracture toughness of randomly oriented chopped glass fibers functionally graded materials. Steel and Composite Structures, 2024, 53: 91–101

[5] Hossein Izadi S M , Fakoor M , Mirzavand B . A novel mixed mode fracture criterion for functionally graded materials considering fracture process zone. Theoretical and Applied Fracture Mechanics, 2024, 134: 104710

[6] Hossein Izadi S M , Fakoor M , Mirzavand B . Developing a new mixed mode I/II fracture criterion for cracked functionally graded materials considering the fracture process zone effect. Acta Mechanica, 2024, 236: 21–35

[7] Pan H , Song T , Ge H . A probabilistic study on the mixed-mode fracture in functionally graded materials. Engineering Failure Analysis, 2021, 120: 105038

[8] Pandey V K , Patel B P . Study of quasi-static crack growth in aluminum/alumina functionally graded material using crack gauge. Theoretical and Applied Fracture Mechanics, 2020, 110: 102777

[9] Bao S , Zhang Y , Liu H , Niu Y , Zhang W , Zeng K . Study on crack propagation of functionally graded ultra-high performance cementitious composite. Engineering Fracture Mechanics, 2023, 282: 109143

[10] Ren H , Zhuang X , Zhu H , Rabczuk T . Variational damage model: A novel consistent approach to fracture. Computers & Structures, 2024, 305: 107518

[11] Ren H , Rabczuk T , Zhuang X . Variational damage model: A new paradigm for fractures. Frontiers of Structural and Civil Engineering, 2025, 19(1): 1–21

[12] Ren H , Zhuang X , Bie Y , Rabczuk T , Zhu H . Dual-horizon peridynamics-based variational damage modeling for complex dynamic fractures. Theoretical and Applied Fracture Mechanics, 2025, 138: 104974

[13] Jernkvist L O . Fracture of wood under mixed mode loading. Engineering Fracture Mechanics, 2001, 68(5): 565–576

[14] Daneshjoo Z , Shokrieh M M , Fakoor M , Alderliesten R C . A new mixed mode I/II failure criterion for laminated composites considering fracture process zone. Theoretical and Applied Fracture Mechanics, 2018, 98: 48–58

[15] Kumar S S , Clement H A , Karthik R . Mixed mode fracture analysis of multiple cracks in flat and curved stiffened panels of aircraft fuselage structures. Archive of Applied Mechanics, 2017, 87(11): 1815–1828

[16] Pitti M R , Dubois F , Pop O . A proposed mixed-mode fracture specimen for wood under creep loadings. International Journal of Fracture, 2011, 167(2): 195–209

[17] Golewski G L . Measurement of fracture mechanics parameters of concrete containing fly ash thanks to use of digital image correlation (DIC) method. Measurement, 2019, 135: 96–105

[18] Hua W , Zhu Z , Zhang W , Li J , Huang J , Dong S . Application of modified fracture criteria incorporating T-stress for various cracked specimens under mixed mode I-II loading. Theoretical and Applied Fracture Mechanics, 2024, 129: 104184

[19] Dong W , Wu Z , Zhou X , Dong L , Kastiukas G . FPZ evolution of mixed mode fracture in concrete: Experimental and numerical. Engineering Failure Analysis, 2017, 75: 54–70

[20] Sørensen B F , Gamstedt E K , Østergaard R C , Goutianos S . Micromechanical model of cross-over fibre bridging—Prediction of mixed mode bridging laws. Mechanics of Materials, 2008, 40(4–5): 220–234

[21] Gutkin R , Laffan M L , Pinho S T , Robinson P , Curtis P T . Modelling the R-curve effect and its specimen-dependence. International Journal of Solids and Structures, 2011, 48(11–12): 1767–1777

[22] Braun M , Ariza M P . A progressive damage based lattice model for dynamic fracture of composite materials. Composites Science and Technology, 2020, 200: 108335

[23] Ricoeur A , Lindner F , Zarjov K . Stochastic aspects of crack deflection and crack path prediction in short fiber reinforced polymer matrix composites. European Journal of Mechanics. A, Solids, 2022, 95: 104598

[24] Lin Z , Li V C . Crack bridging in fiber reinforced cementitious composites with slip-hardening interfaces. Journal of the Mechanics and Physics of Solids, 1997, 45(5): 763–787

[25] Li V C , Leung C K Y . Steady-state and multiple cracking of short random fiber composites. Journal of Engineering Mechanics, 1992, 118(11): 2246–2264

[26] Tiu J , Belli R , Lohbauer U . Rising R-curves in particulate/fiber-reinforced resin composite layered systems. Journal of the Mechanical Behavior of Biomedical Materials, 2020, 103: 103537

[27] FakoorMGhodsiSIzadiS M HDaneshjooZ. Exploring fracture toughness in the damage zone of fiber-reinforced composites: A study on crack orientation effects. Arabian Journal for Science and Engineering, 2025, 1–10

[28] Jiao Y , Zuo Y , Wen Z , Chen Q , Zheng L , Lin J , Chen B , Rong P , Jin K , Du S . Crack-tip propagation laws and energy evolution of fractured sandstone. Engineering Failure Analysis, 2024, 166: 108832

[29] Berto F , Lazzarin P . A review of the volume-based strain energy density approach applied to V-notches and welded structures. Theoretical and Applied Fracture Mechanics, 2009, 52(3): 183–194

[30] Livieri P , Lazzarin P . Fatigue strength of steel and aluminium welded joints based on generalised stress intensity factors and local strain energy values. International Journal of Fracture, 2005, 133(3): 247–276

[31] Lazzarin P , Zambardi R . A finite-volume-energy based approach to predict the static and fatigue behavior of components with sharp V-shaped notches. International Journal of Fracture, 2001, 112(3): 275–298

[32] Razavi N , Aliha M R M , Berto F . Application of an average strain energy density criterion to obtain the mixed mode fracture load of granite rock tested with the cracked asymmetric four-point bend specimens. Theoretical and Applied Fracture Mechanics, 2018, 97: 419–425

[33] Lazzarin P , Berto F . Some expressions for the strain energy in a finite volume surrounding the root of blunt V-notches. International Journal of Fracture, 2005, 135(1): 161–185

[34] Yosibash Z , Bussiba A , Gilad I . Failure criteria for brittle elastic materials. International Journal of Fracture, 2004, 125(3): 307–333

[35] Wu E M . Application of fracture mechanics to anisotropic plates. Journal of Applied Mechanics, 1967, 34(4): 967–974

[36] Leicester R H . Application of linear fracture mechanics to notched timber elements. Progress in Structural Engineering and Materials, 2006, 8(1): 29–37

[37] Spencer B , Barnby J T . The effects of notch and fibre angles on crack propagation in fibre-reinforced polymers. Journal of Materials Science, 1976, 11(1): 83–88

[38] Hunt D G , Croager W P . Mode II fracture toughness of wood measured by a mixed-mode test method. Journal of Materials Science Letters, 1982, 1(2): 77–79

[39] Mall S , Murphy J F , Shottafer J E . Criterion for mixed mode fracture in wood. Journal of Engineering Mechanics, 1983, 109(3): 680–690

RIGHTS & PERMISSIONS

Higher Education Press