J-Integral Analysis of Crack Behavior in Composite Material Geometries

In the field of fracture mechanics, the accurate assessment of crack behavior in structural materials plays a critical role in ensuring safety and reliability, particularly in components subjected to complex mechanical loads. This is especially true for fiber-reinforced composite materials, which, due to their anisotropic and heterogeneous nature, exhibit fracture responses that significantly deviate from classical homogeneous materials like metals (Talreja et al., 2012; Judt et al., 2019). These advanced materials are increasingly used in aerospace, automotive, marine, and civil engineering applications, where lightweight structures must maintain high damage tolerance under service conditions.

Traditional fracture criteria such as the Stress Intensity Factor (SIF), derived from Linear Elastic Fracture Mechanics (LEFM), are often inadequate when applied to composite laminates. This limitation arises from the assumptions of isotropy and linear elasticity inherent in the SIF formulation, which fail to capture the complex interactions between fiber orientations, matrix behavior, and interface delamination Ousset (1999). To overcome these limitations, more generalized fracture parameters such as the J-integral, introduced by Rice (1968), have been employed to evaluate energy release rates under both linear and nonlinear conditions, even in anisotropic materials.

The J-integral is a path-independent contour integral that measures the energy available for crack propagation per unit of crack extension. It accounts for nonlinear material behavior and is thus applicable to elastic–plastic and anisotropic materials, including composite laminates with multiple failure mechanisms (Anderson, 2017; Agarwal et al., 1984). When implemented within the finite element method (FEM), particularly through commercial codes like Abaqus, the J-integral enables detailed evaluation of stress and displacement fields around crack tips without requiring mesh refinement at singularities (Murakami, 1987; Rybicki & Kanninen, 1977).

Numerical strategies based on contour integration have proven to be efficient and accurate, especially in problems involving complex crack geometries and heterogeneous media. For instance, Moës et al. (1999) introduced the Extended Finite Element Method (XFEM), which allows for cracks to be represented independently of the mesh, enhancing flexibility and convergence in fracture simulations. Later developments by Sukumar et al. (2000) further extended XFEM to laminated composites, enabling better modeling of interlaminar and intralaminar cracking.

Recent investigations have demonstrated the utility of the J-integral in analyzing cracked composite structures with geometric discontinuities, such as holes or notches. These features often act as stress concentrators that influence crack initiation and growth depending on their size, shape, and proximity to existing cracks (Abdullah et al., 2017; Fageehi & Alshoaibi, 2024). The interaction between a crack and a nearby hole alters the local stress field and thus the J-integral distribution, which can either amplify or shield the crack-tip driving force depending on the configuration (Zhao et al., 2017; Goushegir et al., 2015).

Experimental and numerical studies have shown that circular holes positioned close to edge cracks can lead to local redistribution of mechanical stresses, modifying the crack-driving energy landscape (Mouritz & Gibson, 2006). For example, employed a hybrid FEM-experimental approach to assess J-integral variation in perforated plates, confirming the sensitivity of energy release rates to geometric discontinuities. Similarly, Khatri et al (2018) analyzed the impact of hole spacing on fracture resistance using an XFEM-based simulation, revealing critical distances beyond which the shielding effect of the hole becomes negligible.

This study investigates the effect of geometric modifications specifically, circular holes placed near an edge crack on the distribution of the J-integral and the crack-driving force in fiber-reinforced polymer composite plates. These materials are widely used in advanced engineering applications due to their high strength-to-weight ratio. However, their anisotropic and heterogeneous nature leads to complex fracture behavior, especially in the presence of geometric discontinuities such as holes, which act as stress concentrators.

To evaluate this behavior, a dual approach was adopted, combining analytical modeling based on linear elastic fracture mechanics (LEFM) with three-dimensional numerical simulations using the Finite Element Method (FEM). The analytical calculations provided a reference framework using classical energy-based formulations suitable for homogeneous materials, while the numerical simulations enabled a more detailed and flexible investigation of stress and displacement fields in complex geometries and composite configurations.

The results show that the presence of nearby holes significantly affects the energy release rate at the crack tip, with the potential to either amplify or reduce the crackdriving force depending on the hole’s size, shape, and distance from the crack. Furthermore, the J-integral proved to be an effective parameter not only for assessing fracture resistance, but also for the early detection of crack initiation and propagation in composite materials. This makes it a valuable tool for the structural design and damage monitoring of advanced composite systems.

This study investigates the fracture behavior of a finite, rectangular composite plate containing an edge crack, subjected to uniform uniaxial tensile loading applied along the top and bottom edges. This loading configuration is designed to simulate pure Mode I crack opening conditions, which are commonly encountered in structural components under tensile stresses such as in aircraft fuselage skins, wing panels, or lightweight structural assemblies. Figure 1 provides a schematic representation of the cracked composite plate, including its geometry, loading, and boundary conditions.

Figure 1.

Schematic of the composite plate with an edge crack.

Figure 2.

Crack-tip geometry showing the coordinate system and the J-integral contour.

Figure 3.

Local coordinate system defined ahead of the crack tip, illustrating the stress components (σ_x, σ_y, τ_xy, τ_yx), radial distance r, and angle θ.

The plate is made of a carbon fiber-reinforced polymer (CFRP) composite, where high-strength carbon fibers are embedded within an epoxy resin matrix. This class of advanced composite materials is widely recognized for its superior strength-to-weight ratio, engineered anisotropic behavior, and excellent resistance to fatigue, corrosion, and environmental degradation. These attributes make CFRP composites highly suitable for demanding applications in the aerospace industry as well as in high-performance structural engineering contexts.

The laminate stacking sequence and fiber orientations were carefully selected to reflect practical configurations commonly used in real-world engineering designs. This ensures that the numerical and theoretical assessments closely mirror actual service conditions. By adopting a realistic stacking arrangement, the study enables a comprehensive evaluation of the fracture response of the composite plate under tensile loading. Such investigations are essential for understanding crack initiation and propagation mechanisms in composite structures, ultimately contributing to the development of more robust and reliable design methodologies.

The stress field equations near the crack tip are shown below (1-4) below for a mode I crack. For a crack as shown in Fig. 1, the stress fields are: 1φ(Z),avecZ=x+iy\varphi (Z),\;avec\;Z = x + iy 2σx+σy=2Re[φ(Z)]{\sigma _x} + {\sigma _y} = 2{\mathop{\rm Re}\nolimits} [\varphi (Z)]

In Equation (1), the function φ (Z) denotes the Westergaard stress function, where Z is a complex variable composed of real and imaginary parts corresponding to the spatial coordinates. This function allows for the compact expression of the stress field in the vicinity of a crack in linear elastic materials.

As shown in Equation (2), the sum of normal stresses is equal to twice the real part of φ (Z), which is an analytic function within the elastic domain. This analytical approach, initially introduced in the context of linear elasticity, has been widely used in fracture mechanics due to its ability to represent singular stress fields near crack tips (Westergaard, 1939). 3σx+σy+2irxy=2Zφ′(Z){\sigma _x} + {\sigma _y} + 2i{r_{xy}} = 2Z\varphi '(Z)

In Eq. (3), φ′(Z) represents the derivative of the Westergaard stress function. This expression establishes the relationship between stress components in the vicinity of a crack and the complex variable Z, offering a mathematically efficient formulation for evaluating stress fields based on complex potentials (Westergaard, 1939). 4φ′(Z)=KI2π·1Z1/2\varphi '(Z) = {{{K_I}} \over {\sqrt {2\pi } }}\cdot{1 \over {{Z^{1/2}}}}

In Eq. (4), K_I represents the stress intensity factor. This expression relates the derivative of the Westergaard stress function to K_I and the complex variable Z, thereby enabling the calculation of stress fields in the vicinity of the crack tip (Westergaard, 1939). 5Y=1.12−0.23(aw)+10.55(aw)2−21.72(aw)3+30.39(aw)4Y = 1.12 - 0.23\left( {{a \over w}} \right) + 10.55{\left( {{a \over w}} \right)^2} - 21.72{\left( {{a \over w}} \right)^3} + 30.39{\left( {{a \over w}} \right)^4} 6KI=Y·σ·πa{K_I} = Y\cdot\sigma \cdot\sqrt {\pi a}

In Eq. (5) provides an empirical formula for the geometry correction factor y, which depends on the relative crack length. Eq. (6), the Mode I stress intensity factor K_I is expressed as a function of the applied stress. These formulations are commonly used to predict fracture behavior in cracked plates under linear elastic conditions (Tada et al., 2000; Anderson, 2017). 7σxx(r,θ)=KI2πrcos(θ2)[ 1−sin(θ2)sin(3θ2) ]{\sigma _{xx}}(r,\theta ) = {{{K_I}} \over {\sqrt {2\pi r} }}\cos \left( {{\theta \over 2}} \right)\left[ {1 - \sin \left( {{\theta \over 2}} \right)\sin \left( {{{3\theta } \over 2}} \right)} \right] 8σyy(r,θ)=KI2πrcos(θ2)[ 1+sin(θ2)sin(3θ2) ]{\sigma _{yy}}(r,\theta ) = {{{K_I}} \over {\sqrt {2\pi r} }}\cos \left( {{\theta \over 2}} \right)\left[ {1 + \sin \left( {{\theta \over 2}} \right)\sin \left( {{{3\theta } \over 2}} \right)} \right] 9τxy(r,θ)=KI2πrcos(θ2)sin(θ2)cos(3θ2){\tau _{xy}}(r,\theta ) = {{{K_I}} \over {\sqrt {2\pi r} }}\cos \left( {{\theta \over 2}} \right)\sin \left( {{\theta \over 2}} \right)\cos \left( {{{3\theta } \over 2}} \right)

In Eq. (7), the normal stress σ_xx is expressed in terms of the radial distance r and the angular coordinate θ, with K_I being the stress intensity factor. Similarly, Eq. (8) gives σ_yy which describes the normal stress in the crack region. Eq. (9) represents the shear stress τ_xy (r, θ) all three equations providing essential expressions for stress components in a cracked body under Mode I loading (Westergaard, 1939).

3.1.

Application of the J-Integral to Composite Materials

10J=∫Γ(wδ1j−σij∂μi∂x1)njdsJ = \int\limits_\Gamma {\left( {w{\delta _{1j}} - {\sigma _{ij}}{{\partial {\mu _i}} \over {\partial {x_1}}}} \right)} \;{n_j}ds 11(w=∫0εijσijdεij)\left( {w = \int_0^{{\varepsilon _{ij}}} {{\sigma _{ij}}} d{\varepsilon _{ij}}} \right)

In Eq. (10), the expression represents a general form of the work integral used to calculate crack-tip parameters in the context of fracture mechanics. The strain energy release rate w, defined in Eq. (11), quantifies the energy associated with material deformation under applied loading. These formulations are fundamental for evaluating energy balance and fracture criteria in stressed materials (Rice, 1972). 12J=KI2E′J = {{K_I^2} \over {E'}} 13E′={ EE1−v2E' = \left\{ {\matrix{ E \cr {{E \over {1 - {v^2}}}} \cr } } \right. 14J=KI2E′=KI2(1−v2)EJ = {{K_I^2} \over {E'}} = {{K_I^2\left( {1 - {v^2}} \right)} \over E}

Equation (12) presents the expression of the J-integral, where E′ denotes the effective Young’s modulus, which accounts for the state of stress in the material. Its value depends on whether the problem is considered under plane stress or plane strain conditions, as defined in Equation (13). By substituting the relevant parameters, this leads to the formulation given in Equation (14), which is widely used in fracture analysis under plane strain conditions. These expressions form a fundamental basis of linear elastic fracture mechanics (LEFM), initially introduced by Rice (1968) and later elaborated upon in comprehensive references such as Anderson (2017).

Table 1 provides a detailed overview of the computed values of the normal stress component σ_yy and the corresponding J-integral as a function of the radial distance r from the crack tip. The calculations were carried out under a constant applied stress of σ = 20 MPa and at an angular orientation of θ = 0°, which corresponds to the direction perpendicular to the crack faces. The values were determined using a combination of analytical formulations given by Equations (5), (6), (8), (12), and (14). These expressions are derived within the framework of linear elastic fracture mechanics (LEFM) and are specifically adapted to characterize the near-tip stress and energy release rate behavior in cracked composite structures.

Table 1.

rradius (mm)	J (N/mm2)	σyy theoretical (MPa)
0.2	15.9	894.92
0.6	15.9	516.68
1	15.9	400.22
1.4	15.9	338.25
2	15.9	283.00
2.4	15.9	258.34
2.8	15.9	239.18
3	15.9	231.07
4	15.9	200.11
5	15.9	178.98

Figure 4 shows the theoretical distribution of the normal stress σ_yy, which decreases sharply and non-linearly as the radial distance r increases from the crack tip. The stress reaches a maximum of 894.92 MPa at r = 0.2 mm, then drops rapidly up to r = 0.6 mm. Beyond this point, the decrease becomes more gradual, reaching approximately 178.98 MPa at r = 5 mm. This behavior, characteristic of Linear Elastic Fracture Mechanics (LEFM), highlights the strong stress concentration near the crack tip.

Figure 4.

Stress vs. distance r from notch tip at θ = 0° (theoretical method).

Figure 5.

Assembly of composite materials.

Figure 6.

The Von Mises stress contour plot of the composite plate without a hole, presented at different radial distances (r) from the crack tip.

Finite element software was used to perform the numerical simulation, applying the contour integral technique to compute both the normal stress σ_yy and the J-integral at the crack tip, by varying the radial distance r and comparing the results with theoretical calculations. Initially, a plate containing a hole was modeled, and the influence of different hole radii and the distance between the hole center and the crack tip was investigated in relation to the stress behavior in the fracture zone.

Property	Symbol	Value	Unit
Longitudinal Young’s modulus	E1	140000	MPa
Transverse Young’s modulus	E2 = E3	9000	MPa
In-plane shear modulus	G12 = G13	5000	MPa
Out-of-plane shear modulus	G23	5000	MPa
Poisson’s ratio	V12 = V13	0.30	–
Poisson’s ratio	V23	0.40	–
Density	ρ	1.5 × 10−6	kg/mm3

Table 2.

r (mm)	σyy theoretical (MPa)	σyy numerical (MPa)
0.2	894.92	884.9
0.6	516.68	454.6
1	400.22	345.9
1.4	338.25	290.6
2	283.00	242.2
2.4	258.34	220.7
2.8	239.18	190.1
3	231.07	197.1
4	200.11	170.2
5	178.98	152

Figure 7 shows the variation of normal stress σ_yy with radial distance r from the crack tip, comparing theoretical and numerical results. Both curves decrease rapidly and nonlinearly as r increases, with peak stress near r = 0, indicating a stress concentration. Theoretical values are consistently higher than numerical ones, especially close to the crack tip. This difference is expected due to numerical approximations that may smooth stress peaks. Despite this, the similar trends confirm that the numerical model represents the physical stress behavior effectively.

Figure 7.

Theoretical vs. numerical σ_yy stress as a function of r (mm) from the notch tip at θ = 0°.

This study analyzes the effect of a hole located near the crack tip on both the normal stress distribution σ_yy and the J-integral value. The investigation focuses on a composite plate containing a hole with varying radii, positioned along the X-axis at different distances S from the crack tip. The analysis is performed using the finite element method (FEM), with all simulation parameters including mesh element sizes kept constant to ensure a direct comparison with the results of the plate without a hole.

5.1.

Numerical results

The data in Table 3 present the values of the normal stress σ_yy for different distances S (15 mm, 20 mm, 25 mm, 30 mm, 35 mm), while Figure 8 illustrates the influence of this distance on the stress distribution. It is clearly shown that as the distance S increases – in other words, as the hole is positioned farther from the crack tip – the stress in the immediate vicinity of the crack decreases. This behavior highlights a local “stress relief” effect resulting from the geometric discontinuity (the hole) being positioned farther away. Furthermore, it is observed that increasing the radial distance (r) from the crack tip leads to a rapid and nonlinear decrease in the normal stress σ_yy, reflecting the gradual dissipation of the stress field around this critical zone.

Figure 8.

Dimensions of perforated plate.

Figure 9.

Perforated plate mesh.

Figure 10.

The Von Mises stress contour plot of the composite plate hole of varying radii (r) and distances (S) from the notch tip.

Table. 3.

r (mm)	σyy (MPa)	σyy (MPa)	σyy (MPa)	σyy (MPa)	σyy (MPa)	σyy (MPa)
0.2	885.1	885	846.3	831.7	811.5	884.9
0.6	455.4	455	454.6	442.2	427	454.6
1	387.2	371.1	364.5	333.7	333.4	345.9
1.4	281.8	279.9	251.6	242.6	229.8	290.6
2	237.7	228.3	228.1	226.5	218.8	241.9
2.4	202.5	201.4	200.9	199.5	199	204.2
2.8	189.1	186.9	186.1	185.5	184,7	190.1
3	190.6	187.4	179.5	179.1	177.7	196.9
4	169.4	162.9	157.3	155.7	154.9	170.2
5	150.7	144.4	142.1	140.3	139.3	152

Table 4 presents the normal stress values σ_yy for a perforated composite plate, showing the average stress values corresponding to different distances S (15 mm, 20 mm, 25 mm, 30 mm, 35 mm) from the crack tip. These values are compared with the normal stress measured at the crack tip as the radial distance (r) increases. Figure 11 illustrates a gradual reduction in the average stress values for increasing S, which remain lower than the σ_yy values observed in the perforated plate. This indicates that as the distance (S) increases, meaning the hole is located farther from the crack tip, the average normal stress σ_yy decreases. Moreover, as the radial distance r increases, moving away from the stress concentration zone, the σ_yy value progressively declines, reflecting the behavior typical of stress attenuation away from concentration zones.

Figure 11.

The plots of the stress σ_yy versus for different values of distance S (in the plate with one hole).

Figure 12.

Plots of the σ_yy stress for different average values of the distance S (in the plate with and without a hole).

Table 4.

r (mm)	σyy (MPa) – Without Hole	σyy (MPa) – With Hole
0.2	884.9	851.92
0.6	454.6	444.84
1.0	345.9	357.98
1.4	290.6	257.14
2.0	241.9	227.88
2.4	204.2	199.6
2.8	190.1	186.46
3.0	196.9	182.92
4.0	170.2	160.44
5.0	152	143.36

Table 5 presents the J-integral values for a perforated plate with the hole positioned at various distances from the crack tip, compared to an unperforated plate. When the hole is located close to the crack tip, the J-integral values are significantly higher than in the case without a hole (for instance, J = 20.0 N/mm² with a hole compared to J = 17.86 N/mm² without a hole). This behavior indicates a strong stress concentration resulting from the interaction between the hole and the crack tip, which increases the energy required for crack propagation. As the distance between the hole and the crack tip increases, this effect gradually diminishes. The J values for both configurations (with and without a hole) begin to converge, and the difference between them becomes marginal at larger distances. This suggests that the influence of the hole is highly localized and weakens as it moves further from the crack front. These results align well with the known behavior of composite materials under load, where geometric discontinuities play a crucial role in stress distribution. The hole acts as a stress concentrator when it is close to the notch, but its effect fades with distance, and the material’s response tends to return to its undisturbed state. This trend is also illustrated in Figure 13, where the plot is limited to a maximum distance of S = 35 mm to support our conclusions.

Figure 13.

J-integral plots for S = 35 mm (in the plate with and without hole).

Table 5.

r (mm)	J (N/mm2) S = 15 mm	J (N/mm2) S = 20 mm	J (N/mm2) S = 25 mm	J (N/mm2) S = 30 mm	J (N/mm2) S = 35	J (N/mm2) Without holes
0.2	17.85	17.84	18.78	19.41	20.01	17.86
0.6	14.29	14.27	14.25	14.39	14.46	14.22
1	15.05	15.31	14.85	14.51	14.75	13.80
1.4	13.78	13.88	13.80	13.94	13.88	13.67
2	14.01	13.92	13.90	13.71	13.64	13.57
2.4	14.19	14.10	13.92	13.72	13.66	13.59
2.8	14.45	14.09	13.90	13.88	13.72	13.59
3	14.61	14.12	13.93	13.88	13.62	13.52
4	15.41	15.35	14.27	13.99	13.86	13.49
5	16.46	15.06	14.62	14.24	13.98	13.48

This study has highlighted the critical role that geometric discontinuities such as holes can play in the fracture behavior of composites plats. The analysis, based on linear Elastic Fracture Mechanics (LEFM) and supported by numerical simulations using the Finite Element Method (FEM), allowed for a detailed understanding of how the position and geometry of holes influence the evolution of the stress field around the crack tip.

Using the J-integral as an energy-based fracture criterion proved to be highly effective in quantifying the material’s sensitivity to crack propagation. It enabled the evaluation of the energy required for crack growth while accounting for the local effects of geometrical singularities. The results showed that the presence of a hole near the crack front leads to a significant increase in the J-integral value, indicating a strong local stress concentration and a higher fracture energy. However, this effect gradually diminishes as the hole moves farther from the crack tip, demonstrating that its influence is highly localized.

Interestingly, it was also observed that holes can be strategically placed to ast as stress-relief elements, which help to slow down or even arrest crack propagation when properly placed. This behavior opens up promising perspectives for passive crack control in composite structures, particularly in applications where durability and structural integrity are critical, such as aerospace and marine industries.

In conclusion, this study confirms that the J-integral is a robust and reliable tool for assessing fracture toughness in composite materials, especially in the presence of geometric discontinuities. It provides both researchers and engineers with a theoretical and practical framework for anticipating fracture mechanisms and designing structures that are more tolerant to damage and less susceptible to catastrophic failure.