Abstract

Translate

In the field of crack mechanics, predicting the direction of a crack is important because this will evaluate whether, when the crack propagates, it penetrates into important areas and whether the structure is dangerous or not. This paper will refer to three theories that predict the propagation direction of cracks: a theory of maximum tangential normal stress, a theory of maximum energy release, and a theory of minimum strain energy density. At the same time, the finite element method (FEM)–ANSYS program will be used to calculate stress intensity factors (SIFs), strain energy release rate (J-integral), stress field, displacement near a crack tip, and crack propagation phenomenon based on the above theories. The calculated results were compared with the results in other scientific papers and experimental results. This research used ANSYS program, an experimental method combined with FEM based on the above energy theories to simulate the J-integral, the SIFs, and the crack propagation. The errors of the SIFs of the FGM rectangular plate has a through-thickness center crack of 1.77%, J-integral of 4.49%, and crack propagation angle $θ_{c}$ of 0.15%. The FEM gave good errors compared to experimental and exact methods.

Full text

Turn on search term navigation

Translate

1. Introduction

FGM is a new-generation composite which is applied in aeronautical engineering (engineering fuselage), in medicine (making teeth, artificial bones), in national defense (bulletproof armor), and in the energy industry (insulation panels, turbines, reactors). The FGM material is a combination of two materials in which the volume ratio of each component changes smoothly and continuously [1,2,3,4,5].

In recent years, FEM is a numerical method that is widely used in mechanics to predict and model the mechanical behavior of structural materials. The FEM has emerged as an effective technique for analyzing crack problems. It is increasingly being used as a viable method in crack modeling that develops under the assumption of linear elastic fracture mechanics [6,7,8].

Furthermore, the FEM is one of the most popular numerical tools for finding approximate solutions to differential equations and partial differential equations. It has been successfully applied in many fields of science and technology to study, simulate, and predict the behavior of structures. Simulating cracks in textures requires the FEM to mesh in accordance with the pore geometry. This method is computationally expensive and complicated and also generates less accurate results to problems when updating the old mesh to a newer version, especially when there are many cracks [9].

For certain layers of FGM, the stress and displacement of cracks asymptotically come from fields of the same form as those in a homogeneous material (Eischen, 1983 [5]). Characteristic variations in the effects of materials manifest in near-SIFs and in higher-order conditions in asymptotic expansion. Erdogan (1995) [6] wrote an excellent review article on this topic. Although several analytic expressions for the stress strength factor for cracks in the FGM have been derived, the investigation is limited to semi-infinite or infinite domains and simple load cases. Analytical expressions for mixed-mode SIFs were obtained recently by Gu and Asaro (1997) [9]. In this case, the crack tip is oriented perpendicular to the material gradient. Although an infinite domain, Nonda and Erdogan (1994) [10] obtained a more general case of material gradients related to crack orientation. In order to correlate the fracture strength data and determine the critical crack length in more general samples, it is clear that a numerical method is necessary to determine the mixed-mode SIF in the FGM [10,11].

Computational fracture mechanics has progressed significantly in last decades, with diverse new methods for the stress evaluation, allowing us to evaluate complex fracture mechanics problems with relatively low computational costs. For metallic structures, as well as aluminum aircraft structures, typical approaches are based on linear-elastic fracture mechanics to evaluate the fatigue crack propagation in order to compute the fatigue life of these structures under different loading and boundary conditions. The greatest computational cost in these analyses is the characterization of stress and displacements fields, required for the determination of stress intensity factors (SIFs) [7].

The SIF is an extremely important parameter in crack mechanics, indicating a degree of stress concentration at a crack tip. In three-dimensional space, the SIFs, $K_{I}$ , $K_{I I}$ , and $K_{I I I}$ , characterize the three independent displacements of a crack, including mode I, mode II, and mode III. When predicting the propagation direction of two-dimensional cracks, the three energy methods, $σ_{θ θ m a x}$ , $S_{m i n}$ , and $G_{m a x}$ , all use two important parameters, $K_{I}$ and $K_{I I}$ , to calculate the bending angle of the crack [12,13,14].

This paper will present the theories based on these methods and some simple propagation crack models referenced from other references. We use ANSYS software 19.2 based on above three energy theories to calculate SIFs, strain energy release rate (J-integral), and crack propagation angle $θ_{c}$ of a rectangular, isotropic FGM plate having a through-thickness center crack and an internal crack. These mechanisms were also employed to check the experimental and the exact methods of calculating the J-integral, SIFs, and crack propagation phenomenon of a rectangular, isotropic FGM plate.

2. Method to Calculate the Propagation Direction of the Crack

2.1. The Theory of Maximum Tangential Normal Stress $σ_{θ θ m a x}$

The mixed expressions of an elastic stress field around a crack tip when expressed in polar coordinates are as follows [15]:

(1) $σ_{r r} = \frac{1}{\sqrt{2 π r}} c o s \frac{θ}{2} [K_{I} (1 + {s i n}^{2} \frac{θ}{2}) + \frac{3}{2} K_{I I} s i n θ - 2 K_{I I} t a n \frac{θ}{2}]$

(2) $σ_{θ θ} = \frac{1}{\sqrt{2 π r}} c o s \frac{θ}{2} [K_{I} {c o s}^{2} \frac{θ}{2} - \frac{3}{2} K_{I I} s i n θ]$

(3) $σ_{r θ} = \frac{1}{2 \sqrt{2 π r}} c o s \frac{θ}{2} [K_{I} s i n θ + K_{I I} (3 c o s θ - 1)]$

In which, $K_{I}$ and $K_{I I}$ are the two SIFs that characterize the two independent displacement forms of the crack, which are mode I and mode II.

The theory of maximum tangential normal stress $σ_{θ θ m a x}$ (first order) for materials confirms that a crack will develop in the direction perpendicular to the theory of maximum tangential normal stress. This theory was proposed by Sih and Erdogan in 1963 [15] (see Figure 1).

Derive Equation (2) in terms of the variable $θ$ and assign it to $0$ .

(4) $\frac{\partial σ_{θ θ}}{\partial θ} = 0$

After rearranging and setting $θ = ∆ θ_{c}$ , Equation (2) will have the following form:

(5) $\frac{K_{I I}}{K_{I}} = \frac{- \sin ∆ θ_{c}}{3 \cos ∆ θ_{c} - 1}$

Solving Equation (5) according to the variable $∆ θ_{c}$ , The bending angle of the crack will be calculated.

According to [16], based on the theory of maximum tangential normal stress, the bending angle of the crack can also be calculated from the following formula:

(6) $∆ θ_{c} = \arccos [\frac{3 K_{I I}^{2} + \sqrt{K_{I}^{4} + 8 K_{I}^{2} K_{I I}^{2}}}{K_{I}^{2} + 9 K_{I I}^{2}}]$

In addition, the bending angle of the crack can also be calculated according to the formula from reference [17] as follows:

(7) $∆ θ_{c} = 2 \arctan [\frac{- 2 K_{I I} / K_{I}}{1 + \sqrt{1 + 8 {(K_{I I} / K_{I})}^{2}}}]$

According to Formula (7), if $K_{I I} = 0$ , then $∆ θ_{c} = 0$ (purely opening form (mode I only)). If $K_{I I} > 0$ , then the crack angle $∆ θ_{c} < 0$ . If $K_{I I} < 0$ , then the crack angle $∆ θ_{c} > 0$ .

2.2. The Theory of Minimum Strain Energy Density $S_{m i n}$

This theory was proposed by Sih in 1974 [18]. Sih developed the formula for calculating a strain energy density $S$ according to stress intensity factors $K_{I}$ and $K_{I I}$ as follows:

(8) $S = a_{11} K_{I}^{2} + 2 a_{12} K_{I} K_{I I} + a_{22} K_{I I}^{2}$

with

(9) $a_{11} = \frac{1}{16 μ} [(1 + c o s θ) (κ - c o s θ)]$

(10) $a_{12} = \frac{1}{16 μ} s i n θ [2 c o s θ - (κ - 1)]$

(11) $a_{22} = \frac{1}{16 μ} [(κ + 1) (1 - c o s θ) + (1 + c o s θ) (3 c o s θ - 1)]$

where

E

is elastic modulus, and

ν

is Poisson’s ratio; shear module

μ = \frac{E}{2 (1 + ν)}

;

κ = \frac{3 - ν}{1 + ν}

for plane stress;

κ = 3 - 4 ν

for the plane strain.

A crack will grow in a direction $θ = ∆ θ_{c}$ , where the strain energy density is minimal.

(12) $\frac{d S}{d θ} = 0 and \frac{d^{2} S}{d θ^{2}} > 0$

The crack starts to propagate when the strain energy density reaches its maximum value $S = S_{c r}$ .

According to reference [19], the maximum value $S_{c r}$ is calculated according to the following formula:

(13) $S_{c r} = (1 - 2 ν) (1 + ν) K_{I C}^{2} / 2 E$

where

K_{I C}

is the critical SIF.

2.3. Theory of Maximum Energy Release $G_{m a x}$

This theory is based on the calculation of Hussain in 1974 [20]. These are the SIFs, $K_{I} (θ)$ and $K_{I I} (θ)$ , of an initial major crack with a bending part with a very small angle $θ$ at the tip, calculated based on the SIFs, $K_{I}$ and $K_{I I}$ , of normal cracks (see Figure 2).

(14) $K_{I} (θ) = g (θ) (K_{I} \cos θ + \frac{3}{2} K_{I I} \sin θ)$

(15) $K_{I I} (θ) = g (θ) (K_{I I} \cos θ - \frac{3}{2} K_{I} \sin θ)$

(16) $g (θ) = (\frac{4}{3 + \cos^{2} θ}) {(\frac{1 - θ / π}{1 + θ / π})}^{\frac{θ}{2 π}}$

According to Irwin’s general expression [21], the energy release rate $G$ for the initial crack with a bent part with an angle $θ$ will be given as follows:

(17) $G (θ) = \frac{1}{E^{'}} (K_{I}^{2} (θ) + K_{I I}^{2} (θ))$

with

E^{'} = E

for plane stress, and

E^{'} = \frac{E}{(1 - ν^{2})}

for plane strain.

Combined with Equations (14)–(16), Equation (17) becomes the following:

(18) $G (θ) = \frac{1}{4 E^{'}} g^{'} (θ) [(1 + 3 \cos^{2} θ) K_{I}^{2} - 8 \sin θ \cos θ K_{I} K_{I I} + (9 - 5 \cos^{2} θ) K_{I I}^{2}]$

The propagation angle of the crack is found by minimizing $G (θ)$ .

(19) $\frac{\partial G (θ)}{\partial θ} = 0$

And the following stability conditions must be satisfied:

(20) $\frac{\partial^{2} G (θ)}{\partial θ^{2}} < 0$

The general form of Equation (18) can be abbreviated as follows:

(21) $G (θ) = \frac{1}{4 E^{'}} [A_{11} K_{I}^{2} (θ) + A_{22} K_{I I}^{2} (θ) + 2 A_{12} K_{I} (θ) K_{I I} (θ)]$

(22) $[\begin{matrix} A_{11} \\ A_{12} \\ A_{22} \end{matrix}] = g^{2} (θ) [\begin{matrix} 4 - 3 \sin^{2} θ \\ - 2 \sin 2 θ \\ 4 + 5 \sin^{2} θ \end{matrix}]$

2.4. Comparison Between Three Methods of Calculating Crack Propagation Direction

The following is a formula comparing the results between the theory of maximum tangential normal stress with the theory of maximum energy release and the theory of minimum strain energy density as referenced from [22,23,24,25,26]. For the convenience of comparison, we put the following:

(23) $M^{e} = \frac{2}{π} \tan^{- 1} (\frac{K_{I}}{K_{I I}})$

In addition, this research would like to present another comparison between the three methods of calculating the propagation direction of the crack. This result is referenced from the literature [15,18,20,27,28,29,30,31,32].

Rice [33] and Eshelby [34] introduced a contour integral surrounding a crack front, as illustrated in Figure 3:

(24) $J = \int_{Γ} (W d y - \vec{T} \frac{\partial \vec{μ}}{\partial x} d s)$

Here, J denotes the energy release rate (J/m²); W represents the elastic strain energy density (J/m³); $\vec{μ}$ is the displacement vector evaluated at the differential element $d s$ ; n is the outward unit normal to the contour $Γ$ ; $d s$ is a differential segment along the contour; $\vec{T} (\partial \vec{u} / \partial x) d s$ corresponds to the input work; $s$ denotes the arc length; $Γ$ is an arbitrary counterclockwise contour; $\vec{T}$ refers to the traction (tensile) vector acting on the body enclosed by $Γ$ .

Rice [33] and Riveros [35] certified a path independent concepts and found that stress energy release rate $G$ is equal to $J$ -integral:

(25) $J = G = \frac{K_{I}^{2} + K_{I I}^{2}}{E^{'}}$

where

E^{'} = E

for plane stress,

E^{'} = \frac{E}{(1 - ν^{2})}

for plane strain.

3. A Results and Discussions: A FGM Plate with the Through-Thickness Center Crack

Figure 4 illustrates the FGM plate along with its boundary conditions and dimensions. The material was characterized by a Poisson’s ratio of $ν$ = 0.25, with geometric parameters a = 0.5 cm, H = 10 cm, and W = 2 cm. An applied stress of σ = 10⁶ N/m² was considered. The Young’s modulus E varied exponentially across the plate from the left to the right edge. The analysis was conducted under the assumption of a plane strain condition.

(26) $E (x) = E_{0} e^{λ x}$

where

E_{0} = 10^{3}

Pa, and

λ = 2

Experimental tests were performed on specimens with dimensions H = 10 cm, W = 2 cm, a = 0.5 cm, and a thickness of t = 0.5 cm, fabricated using 3D printing technology. The specimens were made of nylon ABS, with a Poisson’s ratio of ν = 0.25, and prepared in accordance with the ASTM E399 standard; ASTM International, West Conshohocken, PA, USA, 11 July 2023. Tests were conducted using an LF2373 Friction Tester under varying loading conditions. Two testing speeds were applied: v₁ =2 mm/min for Case 1 and v₂ =3 mm/min for Case 2. In both cases, a constant load of F = 100 N was applied, based on the relationship F = σ/A, where A = W × t is the cross-sectional area. The experimental setup is shown in Figure 5.

The vertical displacements $u_{y}$ obtained from the finite element method (FEM) are presented in Figure 6, corresponding to the model with dimensions H = 10 cm, a = 0.5 cm, and W = 2 cm (see model in Figure 4). The maximum displacement predicted by FEM is Max( $u_{y}$ ) = 0.15882 mm. In comparison, experimental measurements yielded Max( $u_{y}$ ) = 0.16335 mm for Case 1, with a percentage error of 5.55%, and Max( $u_{y}$ ) = 0.17032 mm for Case 2, corresponding to a percentage error of 6.42%. These results are summarized in Figure 6 and Figure 7 and Table 1. Figure 7 illustrates the stress–displacement curves obtained from the FEM and the two experimental cases. The maximum stress was observed to be Max( $σ_{y y}) = 1$ MPa (or 10⁶ N/m²), with all three curves showing close agreement. These findings confirm that the FEM approach provides reliable predictions when validated against experimental data.

The strain results $ε_{y}$ , obtained using the finite element method (FEM), are shown in Figure 8 for a model with dimensions H = 10 cm, a = 0.5 cm, and W = 2 cm, consistent with the geometry presented in Figure 4. The maximum strain calculated by the FEM is Max( $ε_{y}$ ) = 0.001531 = 0.001531 mm/mm. Experimental results yielded Max( $ε_{y}$ ) = 0.0014833 mm/mm for Case 1, corresponding to a percentage error of 3.12%, and Max( $ε_{y}$ ) = 0.001579 mm/mm for Case 2, with an error of 3.13%. These findings are summarized in Figure 8 and Figure 9 and Table 1. Figure 9 presents the stress–strain curves obtained from FEM and both experimental cases. The maximum stress was observed to be Max( $σ_{y y}) = 1$ MPa (or 10⁶ N/m²), and the three curves show a high degree of consistency. These results further confirm the accuracy and reliability of the FEM in predicting mechanical behavior compared to experimental methods.

Figure 10 presents the delamination propagation results along the y-direction, obtained from FEM simulations and two experimental cases. The maximum delamination predicted by the FEM is Max( $ε_{y}$ ) = 0.1543 mm. In comparison, experimental results show Max( $ε_{y}$ ) = 0.1648 mm for Case 1, with a percentage error of 6.37%, and Max( $ε_{y}$ ) = 0.1585 mm for Case 2, with an error of 2.65%. The load–delamination propagation curves for FEM and both experimental cases are also illustrated in Figure 10. The maximum applied load was F = 100 N, calculated based on the relation F = σ/A, where σ = 10⁶ N/m² and the cross-sectional area A = 2 × 0.5 × 10⁻⁴ m². The three curves exhibit close agreement, confirming the consistency between the FEM predictions and experimental observations. These results further validate FEM as a dependable tool for analyzing delamination behavior.

Results of the FEM and two case experiments show that the extension of crack length $∆ a$ of a FGM plate are illustrated in Figure 11. The maximum extension of the crack length by the FEM at the time $t =$ 13 (s) was Max( $∆ a$ ) $=$ 0.6439 (mm), by the experiment for Case 1 was Max( $∆ a$ ) $=$ 0.6423 mm, % Error $=$ 0.25%; and for Case 2 was Max( $∆ a$ ) = 0.6286 mm, % Error = 2.43%. Figure 11 illustrates the crack length extension versus time curves obtained from the FEM analysis and two experimental cases. The three curves exhibit nearly identical trends, indicating a strong correlation between numerical predictions and experimental observations. These findings further confirm the reliability and accuracy of the FEM approach in modeling crack propagation behavior.

As shown in the newly included Figure 6a, the finite element mesh was refined in the vicinity of the crack to capture the high stress gradients accurately. A mapped mesh was employed in the crack region with a higher element density, while a coarser mesh was used in regions far from the crack to optimize computational efficiency. The element type used was PLANE183 for 2D modeling, which is suitable for simulating linear elastic behavior in fracture mechanics.

The specimen was modeled with symmetry in both geometry and loading conditions. Displacement boundary conditions were applied as follows: the bottom edge of the specimen was fixed in the Y-direction (to prevent rigid body motion), and a uniform tensile displacement was applied at the top edge in the Y-direction. The lateral edges were left free to simulate uniaxial tensile loading. These settings replicate a mode I crack-opening scenario under plane stress conditions. These updates are now described in Figure 6b,c. The general model used in the numerical analysis, including crack geometry, boundary conditions, and material gradient direction, is presented in Figure 6d.

The stress intensity factor (SIF) analysis results, based on the model parameters W = 2 cm, σ = 10⁶ N/m², and a = 0.5 cm, are taken from reference [30] as follows:

(27) $K_{I} = α σ \sqrt{π a} = 4.2152 MPa \sqrt{m m}$

The crack geometry correction factor is adopted from Ref. [30] as follows:

(28) $α \approx {[1 + 0.5 {(a / w)}^{2} + 20.46 {(a / w)}^{4} + 81.72 {(a / w)}^{6}]}^{1 / 2} = 1.0635$

Based on Equation (25), the following results were obtained:

(29) $E^{'} = \frac{E}{(1 - ν^{2})} = 1009.0817 P a$

(30) $J = \frac{K_{I}^{2}}{E^{'}} = 0.0176 \times 10^{12} Pa \cdot mm = 0.0176 mJ / {mm}^{2}$

The mode I stress intensity factors (SIFs), $K_{I}$ , obtained from the FEM analysis are presented in Figure 12. The FEM yielded a value of $K_{I}$ = 4.2896 MPa· $\sqrt{m m}$ , while the exact method produced $K_{I}$ = 4.2152 MPa· $\sqrt{m m}$ , corresponding to a percentage error of 1.77%, as reported in Table 2. Figure 13 provides a comparison of SIF versus crack length curves obtained from FEM, the exact method, and other reference results. All five curves demonstrate close agreement. These findings confirm that FEM provides accurate and reliable predictions for stress intensity factors.

Figure 14 presents the J-integral results obtained using the finite element method (FEM) and the exact method (as defined in Equation (30)). The FEM yielded a value of J = 0.01681 mJ/mm², while the exact method produced J = 0.0176 mJ/mm², resulting in a percentage error of 4.49%, as reported in Table 3. Figure 15 shows a comparison of the J-integral versus crack length curves obtained from FEM, the exact method, and other reference data. The five curves are in close agreement, further confirming that the FEM provides reliable and accurate results for fracture analysis.

The results $M^{e}$ (see Equation (23)) and $θ_{c}$ (see Equation (7)) are shown in Figure 16—by the FEM, exact, and reference energy methods. $M^{e}$ by the FEM: $M^{e} =$ 0.31, by the exact method: $M^{e} =$ 0.32; % Error $M^{e}$ by the FEM compared to the exact method: Error $=$ 3.12%; Values $θ_{c}$ by FEM: $θ_{c} =$ −64.5°; by the exact method $θ =$ −64.6°; % Error $θ_{c}$ by the FEM compared to the exact method: Error $=$ 0.15% (Table 4). Because values $K_{I I} > 0$ , values $θ_{c} < 0$ —by theory (Equation (7)). Figure 16 compares the curves $θ_{c}$ versus $M^{e}$ of the FEM, exact method, and reference energy methods; the five curves lie the same. The values $M^{e}$ calculated at each step according to the three theories are approximately equal to 1 (between 0–1). Therefore, the bending angle of the crack calculated according to the three theories at each step has nearly equal values. This is consistent with the graph comparing the results between the three theories that predict the propagation direction of the crack referenced from the literature. Therefore, the crack path modeled according to the three theories has almost same form. The results illustrate that the FEM is the reliable method.

The following criteria were compared:

a. Maximum Tangential Stress (MTS) Criterion:

Advantage: Simple and widely used, especially for brittle materials; provides good agreement with experimental data under mode I-dominated loading. Disadvantage: Less accurate for mixed-mode loading or ductile fracture. Applicability: Effective in predicting initial crack propagation in brittle materials, such as ceramics or high-strength steels.

b. Maximum Energy Release Rate (G-max) Criterion:

Advantage: Energy-based and theoretically sound; considers both mode I and mode II contributions. Disadvantage: Requires accurate energy release rate calculation; less intuitive than stress-based methods. Applicability: Suitable for materials with significant energy dissipation, such as composites or polymers.

c. Minimum Strain Energy Density (SED) Criterion:

Advantage: Accounts for the energy density around the crack tip, providing a more localized perspective. Disadvantage: Highly sensitive to mesh refinement and may yield ambiguous directions in certain cases. Applicability: Useful in complex geometries or when stress fields are not symmetric.

The influence of these criteria on the predicted crack path and the observed deviations in the numerical results is also discussed. This comparative analysis helps identify the most suitable criterion for different material types and loading conditions in practical fracture mechanics applications.

A comparative analysis between our proposed method and other relevant studies in the field of crack mechanics and finite element simulation is provided. The proposed approach was compared with studies by Sih, 1974 [18], Erdogan and Sih, 1963 [19], and Hussain et al., 1974 [20], which utilize similar FEM-based techniques for crack propagation analyses. The comparison highlights the following key points of novelty and practical contribution in this work:

1. Novelty: Unlike prior studies that primarily focus on stress intensity factor evaluation or static crack modeling, the proposed approach integrates dynamic crack propagation modeling, multi-theory comparison, and refined meshing near crack tips within the ANSYS environment, providing enhanced accuracy and broader applicability.

2. Practical Application: The methodology demonstrated in this paper can be directly applied to structural integrity assessment in components such as pressure vessels, pipelines, and aircraft panels, where understanding crack direction and growth behavior is critical for failure prevention.

In the revised manuscript, we have included a new subsection (Section 3, lines 329–350), presenting a parametric study that evaluates the influence of different gradient indices (n = 0.5, 1, 2, and 5) on the stress intensity factors (SIFs) and crack propagation direction (Figure 17). The results demonstrate that a higher gradient index generally leads to a reduction in SIF values due to the increased contribution of the tougher material near the crack tip. Crack paths tend to curve toward the region with lower stiffness, consistent with previous observations in FGM fracture mechanics. Additionally, we briefly discuss how alternative material combinations (e.g., ceramic–metal vs. metal–metal FGMs) could further impact crack growth behavior, and we outline this as a direction for future work to keep the scope focused.

4. Conclusions

In this research, we use the FEM to simulate the crack tip behavior in two dimensions, and stress–displacement fields near a crack tip could be more exact compared with the experimental method. The calculation of the SIFs and J-integral will be faster than the analytical theory when encountering a complex crack model. Concurrently, it was possible to relatively calculate the growth and shape of the two-dimensional crack, as it propagates under different conditions of the model, material, and loading method.

The SIFs, J-integral, and crack propagation angle mentioned in Section 3 were calculated based on three energy methods. The results obtained via the FEM were better than those obtained using the experimental method, when compared to the exact method. Findings can be introduced as a novel approach to calculate the SIF, J-integral, and crack propagation phenomenon of the FGM plate.

Further, this research presented the FEM based on the three energy methods above to calculate the SIFs, J-integral, and crack propagation phenomenon of the FGM plate based on the reference models. Such mechanisms have widespread applications in industrial practice and real-world engineering.

Author Contributions

Conceptualization, S.-C.H. and H.-D.N.; methodology, H.-D.N.; software, H.-D.N.; validation, S.-C.H. and H.-D.N.; formal analysis, H.-D.N.; investigation, S.-C.H.; resources, H.-D.N.; data curation, H.-D.N.; writing—original draft preparation, data curation; writing—review and editing, data curation; supervision, S.-C.H. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Maximum tangential normal stress in polar coordinate system.

Figure 2 Initial primary crack with a bending part with angle $θ$ : (a) original crack; (b) kinked crack.

Figure 3 $J$ -integral contours around the crack surfaces.

Figure 4 A through-thickness center crack model.

Figure 5 The experimental setup.

Figure 6 (a) The finite element mesh was refined in the vicinity of the crack. (b) The bottom edge of the specimen was fixed. (c) A uniform tensile displacement was applied at the top edge. (d) The general model used in the numerical analysis. (e) The results of displacement $u_{y}$ .

Figure 7 Comparison of displacement $u_{y}$ of FEM and experimental results.

Figure 8 The results of strain $ε_{y} .$

Figure 9 Comparison of strain $ε_{y}$ of FEM and experimental results.

Figure 10 Comparison of delamination propagation along the y-direction between FEM simulation and experimental results.

Figure 11 Crack length extension $∆ a$ : FEM vs. experimental results.

Figure 12 The result of SIF.

Figure 13 Comparison of SIFs versus crack length of FEM, exact, and other reference results (Sih, 1974 [14]; Erdogan and Sih, 1963 [16]; Hussain, et al., 1974 [20]).

Figure 14 J-integral results obtained from FEM.

Figure 15 Comparison of J-integral versus crack length of FEM, exact, and other reference results (Sih, 1974 [14]; Erdogan and Sih, 1963 [16]; Hussain, et al., 1974 [20]).

Figure 16 Comparison of modeling of $θ_{c}$ versus $M^{e}$ of FEM, exact, and reference energy methods (Sih, 1974 [14]; Erdogan and Sih, 1963 [16]; Hussain, et al., 1974 [20]).

Figure 17 Effect of gradient index on SIF $K_{I}$ in FGM.

Table 1

Comparison of a maximum strain and displacement between FEM and experimental method.

	FEM	Exper.		Error (%)
		Exper. 1	Exper. 2	Exper. 1	Exper. 2
Displacement $u_{y}$ (mm)	0.15882	0.16335	0.17032	5.55	6.42
Strain $ε_{y}$ (mm/mm)	0.001531	0.0014833	0.001579	3.12	3.13

Table 2

Comparison of SIFs $K_{I}$ between the exact method and the FEM.

Exact (MPa· $\sqrt{m m}$ )	4.2152
FEM (MPa· $\sqrt{m m}$ )	4.2896
Error (%)	1.77

Table 3

J-integral comparison: exact method vs. FEM.

FEM mJ/mm²	0.01760
Exact mJ/mm²	0.01681
Error (%)	4.49

Table 4

Comparison of the $M^{e}$ and $θ_{c}$ between the exact method and the FEM.

	$M^{e}$	$θ_{c}$ (Degree)
Exact	0.32	−64.6
FEM	0.31	−64.5
Error (%)	3.12	0.15

References

1. Mutnbak, M.; Abbadi, A.; Mousa, S.; Abd-Elhady, A.A.; Sallam, H.E.M.; Reda, R.M. Effects of specimen geometry and size on mode I and mixed mode fracture behavior of high strength fiber reinforced concrete. Sci. Rep.; 2025; 15, 15286. [DOI: https://dx.doi.org/10.1038/s41598-025-99013-9] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/40312553]

2. Nguyen, H.-D.; Huang, S.-C. Use of XTFEM based on the consecutive interpolation procedure of quadrilateral element to calculate J-integral and SIFs of an FGM plate. Theor. Appl. Fract. Mech.; 2023; 127, 103985. [DOI: https://dx.doi.org/10.1016/j.tafmec.2023.103985]

3. Nguyen, H.-D. Using Finite Element Method to Calculate Strain Energy Release Rate, Stress Intensity Factor and Crack Propagation of an FGM Plate Based on Energy Methods. Int. J. Mech. Energy Eng. Appl. Sci.; 2025; 17, 10345.

4. Nguyen, H.-D.; Huang, S.-C. Using the Extended Finite Element Method to Integrate the Level-Set Method to Simulate the Stress Concentration Factor at the Circular Holes Near the Material Boundary of a Functionally-Graded Material Plate. J. Mater. Res. Technol.; 2022; 21, pp. 4658-4673. [DOI: https://dx.doi.org/10.1016/j.jmrt.2022.11.062]

5. Ma, L.; Wu, L.-Z.; Guo, L.-C.; Zhou, Z.-G. Dynamic behavior of a finite crack in the functionally graded materials. Mech. Mater.; 2005; 37, pp. 1153-1165. [DOI: https://dx.doi.org/10.1016/j.mechmat.2005.05.004]

6. Naebe, M.; Shirvanimoghaddam, K. Functionally graded materials: A review of fabrication and properties. Appl. Mater. Today; 2016; 5, pp. 223-245. [DOI: https://dx.doi.org/10.1016/j.apmt.2016.10.001]

7. Nguyen, H.-D. Analysis for Stress Concentration Factor at the Circular Holes Near Materials Border by Extended Finite Element Method. Proceedings of the 4th International Online Conference on Materials; Online, 3–5 November 2025.

8. Nguyen, H.D.; Huang, S.-C. Applying Quarter-Point Singular Element in Finite Element Method to Calculate J-Integral and Stress Intensity Factors of a Rectangular Isotropic FGM Plate. Materials; 2025; 27, 12085.

9. Sladek, J.; Sladek, V.; Zhang, C. An advanced numerical method for computing elastodynamic fracture parameters in functionally graded materials. Comput. Mater. Sci.; 2005; 32, pp. 532-543. [DOI: https://dx.doi.org/10.1016/j.commatsci.2004.09.011]

10. Nguyen, H.-D.; Huang, S.-C. Extend Finite Element Method for A Plate with Two Holes and Multiple Cracks. Proceedings of the Taiwan International Conference on Ocean Governance; Kaohsiung, China, 23 September 2021.

11. Nguyen, H.-D.; Huang, S.-C. Extended Finite Element Method for Modeling Arbitrary Discontinuities with Center Edge Crack in Finite Dimensional Plate under Tension and Shear. Proceedings of the International Conference on Advanced Technology Innovation; Okinawa, Japan, 1 September 2020.

12. Sérgio, E.R.; Gonzáles, G.L.G.; Gonzáles, G.L.G.; Vasco-Olmo, J.M.; Antunes, F.V.; Prates, P.; Díaz, F.A.; Neto, D.M. A comparison between FEM predictions and DIC results of crack tip displacement field in AA2024-T3 CT specimens. Eng. Fract. Mech.; 2025; 318, 110964. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2025.110964]

13. Monfared, M.; Ayatollahi, M. Dynamic stress intensity factors of multiple cracks in an orthotropic strip with FGM coating. Eng. Fract. Mech.; 2013; 109, pp. 45-57. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2013.07.002]

14. Sih, G.C.; Cha, B.C.K. A fracture criterion for three-dimensional crack problems. Eng. Fract. Mech.; 1974; 6, pp. 699-723. [DOI: https://dx.doi.org/10.1016/0013-7944(74)90068-X]

15. Nguyen, H.D. Using the eXtended Finite Element Method (XFEM) to Simulate Own Frequency under External Influences of a Closed System based on Dynamic Compensation Method. J. Int. Multidiscip. Res.; 2025; 27, 985.

16. Erdogan, F.; Sih, G.C. On the Crack Extension in Plates Under Plane Loading and Transverse Shea. J. Basic Eng. Dec.; 1963; 85, pp. 519-525. [DOI: https://dx.doi.org/10.1115/1.3656897]

17. Gao, X.W.; Zhang, C.; Sladek, J.; Sladek, V. Fracture analysis of functionally graded materials by a BEM. Compos. Sci. Technol.; 2007; 68, pp. 1209-1215. [DOI: https://dx.doi.org/10.1016/j.compscitech.2007.08.029]

18. Nguyen, H.-D.; Huang, S.-C. Designing and Calculating the Nonlinear Elastic Characteristic of Longitudinal–Transverse Transducers of an Ultrasonic Medical Instrument Based on the Method of Successive Loadings. Materials; 2022; 15, 4002. [DOI: https://dx.doi.org/10.3390/ma15114002]

19. Nguyen, H.-D. Extended Finite Element Approach for Simulating Arbitrary Openings in Functionally Graded Plates. Int. J. Mech. Energy Eng. Appl. Sci.; 2025; 17, 885.

20. Hussain, M.A.; Pu, S.L.; Underwood, J. Strain energy release rate for a crack under combined mode-I and -II. Natl. Symp. Fract. Mech.; 1974; 560, pp. 2-28.

21. Nguyen, H.D. Application of the eXtended Finite Element Method for Simulating Random Holes in Functionally Graded Plates. Mater. Proc.; 2025; 27, 975.

22. Kim, J.-H.; Paulino, G.H. Finite element evaluation of mixed mode stress intensity factors in functionally graded materials. Int. J. Numer. Meth. Eng.; 2002; 53, pp. 1903-1935. [DOI: https://dx.doi.org/10.1002/nme.364]

23. Wu, C.; He, P.; Li, Z. Extension of J integral to dynamic fracture of functional graded material and numerical analysis. Comput. Struct.; 2002; 80, pp. 411-416.

24. Song, S.H.; Paulino, G.H. Dynamic stress intensity factors for homogeneous and smoothly heterogeneous materials using the interaction integral method. Int. J.Solids Struct.; 2006; 43, pp. 4830-4866. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2005.06.102]

25. Dai, M.-J.; Xie, M.-Y. A novel inverse extended finite element method for structural health monitoring of cracked structures. Ocean. Eng.; 2025; 325, 120786. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2025.120786]

26. Nguyen, H.-D.; Huang, S.-C. Rubber’s crack of a platen in papermaking machines at a crack tip of hyperelastic materials with finite element method. Proceedings of the IEEE 7th International Conference on Applied System Innovation; Chiayi, China, 24–25 September 2021.

27. Bayesteh, H.; Mohammadi, S. FEM fracture analysis of orthotropic functionally graded materials. Compos. Part B; 2013; 44, pp. 8-25. [DOI: https://dx.doi.org/10.1016/j.compositesb.2012.07.055]

28. Nguyen, H.-D.; Huang, S.-C. Calculating Strain Energy Release Rate, Stress Intensity Factor and Crack Propagation of an FGM Plate by Finite Element Method Based on Energy Methods. Materials; 2025; 67, 8885.

29. Singh, I.V.; Mishra, B.K.; Bhattacharya, S. FEM simulation of cracks, holes and inclusions in functionally graded materials. Int. J. Mech. Mater. Des.; 2011; 7, pp. 199-218. [DOI: https://dx.doi.org/10.1007/s10999-011-9159-1]

30. Nestor, P. Fracture Mechanics; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2004.

31. Nguyen, H.-D.; Huang, S.-C. The Uniaxial Stress Strain Relationship of Hyperelastic Material Models of Rubber Cracks in the Platens of Papermaking Machines Based on Nonlinear Strain and Stress Measurements with the Finite Element Method. Materials; 2021; 14, 7534. [DOI: https://dx.doi.org/10.3390/ma14247534]

32. Nguyen, H.-D. Using Matlab to Study of the Response of the System to Oscillatory Influences be an Oscillatory Effect on the Resting Antenna. Int. J. Multidiscip. Res. Growth Eval.; 2025; 27, 885.

33. Rice, J.R. A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks. J. Appl. Mech.; 1968; 35, pp. 379-386. [DOI: https://dx.doi.org/10.1115/1.3601206]

34. Eshelby, J.D. Calculation of Energy Release Rate. Prospects of Fracture Mechanics; Sih, G.C.; Von Elst, H.C.; Broek, D. Noordhoff: Groningen, The Netherlands, 1974; pp. 69-84.

35. Riveros, G.A.; Gopalaratnam, V.S. Post-cracking behavior of reinforced concrete deep beams: A numerical fracture investigation of concrete strength and beam size. ASCE Struct. J.; 2005; 47, 3885.

© 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Calculating Strain Energy Release Rate, Stress Intensity Factor and Crack Propagation of an FGM Plate by Finite Element Method Based on Energy Methods

Content area

Abstract

Full text

1. Introduction

2. Method to Calculate the Propagation Direction of the Crack

2.1. The Theory of Maximum Tangential Normal Stress σθθmax

2.2. The Theory of Minimum Strain Energy Density Smin

2.3. Theory of Maximum Energy Release Gmax

2.4. Comparison Between Three Methods of Calculating Crack Propagation Direction

3. A Results and Discussions: A FGM Plate with the Through-Thickness Center Crack

4. Conclusions

2.1. The Theory of Maximum Tangential Normal Stress $σ_{θ θ m a x}$

2.2. The Theory of Minimum Strain Energy Density $S_{m i n}$

2.3. Theory of Maximum Energy Release $G_{m a x}$