Quasi-static damage evolution of carbon fiber-reinforced polymer laminates under three-point bending

Abstract

Carbon fiber reinforced polymer (CFRP) laminates inevitably undergo bending loads during service, but bending failure modes are complex and not fully understood. This study investigates the quasi-static progressive damage of CFRP laminates under three-point bending primarily through finite element analysis, with experiments serving to validate the overall trends and final failure modes. The influence of interlaminar shear strength (ILSS) and stacking thickness on the flexural properties is studied. Experimentally, as the stacking thickness increases from 0.6 mm to 1.2 mm, the flexural strength decreases from 2099 MPa to 1790 MPa. When ILSS increases from 20 MPa to 80 MPa, the peak load increases by approximately 1.4–1.9%, and the simulation indicates that the predominant failure mode tends to shift to fiber fracture. To analyze the damage evolution behavior of the interlaminar interface, an interlaminar interface catastrophic failure model based on Taylor expansion is established. The damage of the interlaminar interface exhibits a power-law characteristic in the simulation with a power exponent of 0.5. At ILSS20 MPa, the damage coefficient λ increases from 1.88 to 4.70 as the stacking thickness increases. As the displacement approaches the failure point, the damage rate rises rapidly, suggesting a power-law singularity (−0.5).

Keywords

progressive damage CFRP laminates stacking thickness interlaminar shear strength

Introduction

Carbon fiber reinforced polymer (CFRP) laminates are widely used in aerospace, automotive, wind energy, and other industrial sectors due to their high specific strength, specific stiffness, and excellent fatigue resistance.^1,2,3 These mechanical advantages arise from the anisotropic nature of the material, which allows engineers to tailor fiber orientations and stacking sequences to meet specific loading requirements.^4,5,6 However, this same anisotropy also introduces complex failure mechanisms—such as fiber fracture, matrix cracking, and delamination—that are not fully understood, particularly under bending loads.⁷ Among various service conditions, bending is a prevalent loading mode for laminates, and two key factors—stacking thickness and interlaminar shear strength (ILSS)—critically influence their damage evolution and failure behavior.⁸

The failure mode of composite laminates is influenced by numerous factors, encompassing interface characteristics, material thickness, fiber ply angle,⁹ stacking sequence,^10,11 interlaminar hybridization,^12,13,14,15 fiber treatments,¹⁶ weave patterns,¹⁷ surface scratches,¹⁸ temperature,¹⁹ and other service conditions. Carrillo et al.⁹ found that composite ply orientation had only a minor effect on the flexural properties of fiber-metal laminates. In contrast, Husain et al.¹⁰ and Gai et al.¹¹ demonstrated that interlaminar angle and stacking sequence differences significantly influence stress distribution, damage evolution, and flexural strength in pure CFRP composites, indicating a strong effect of material combination. Numerous scholars engage in the hybridization of each layer of materials and modify the stacking sequence to examine their impacts on flexural properties and failure modes. Studies on carbon/glass,¹² jute/madar/glass,¹³ carbon/basalt,¹⁴ and aluminum/basalt/pineapple¹⁵ hybrids collectively show that fiber hybridization improves failure modes and toughness, but flexural strength and modulus are still dominated by the longitudinal fiber type. Moreover, some studies suggest that the flexural performance of composite materials is controlled by the strength of the surface material, and by altering the stacking sequence, it is possible to delay the propagation of cracks in the thickness direction.^13,14 In addition to the surface material, surface scratches on CFRP laminates suppress delamination, promote trans-laminar fracture, and reduce load-bearing capacity in a geometry-dependent manner.¹⁸ Beyond these factors, the mechanical behavior of CFRP laminates is also influenced by interface characteristics and layer thickness considerations, which will be discussed below.

The interface of composite laminates comprises both interlaminar and in-plane interfaces, both of which have a significant impact on the damage modes of composite laminates. Numerous scholars improve the mechanical properties of composite laminates by enhancing the properties of interlaminar interface, employing methods such as reinforcement with carbon black,²⁰ toughening modification of epoxy resin²¹, matrix modification with carbon nanotubes,²² fiber surface modifications,^16,23 nanoparticle surface coatings,²⁴ and the use of electrospun nanofiber interleaves.²⁵ Nevertheless, under various loading conditions, the influence of the interlaminar interface on mechanical properties is different. Li et al.²⁶ discovered that the interlaminar interface properties had a greater impact on the flexural strength than on the tensile strength of glass fiber reinforced aluminum laminate. Additionally, composite laminates also encompass interlaminar interfaces. Numerous scholars have noted that weak interlaminar interface properties lead to delamination fracture in unidirectional laminates, whereas strong interlaminar interface properties result in a stepwise fracture behavior in unidirectional laminates.^23,27,28 The main cause of the damage in the composite laminates is a result of both in-plane and interlaminar damage mechanisms competing with each other. Under bending loads, the laminate is more susceptible to interlaminar damage.²⁹ Although many scholars have discussed the influence of interfaces on the mechanical properties of laminates, there are few studies on damage failure at the interface. The main reason for this is that the damage evolution process is challenging to capture experimentally.

Furthermore, material thickness is also a crucial factor influencing the mechanical properties of composite laminates.³⁰ Laux et al.³¹ investigated the impact of thickness on quasi-isotropic carbon/epoxy composites under multiaxial loading. The findings revealed that laminate thickness had a greater effect on the failure behavior under tensile shear load, whereas its influence is less significant under compressive shear load. Liu et al.³² examined how the thickness of carbon fiber reinforced thermoplastic unidirectional laminates affects their tensile and flexural properties. Upon reducing the thickness of the single layer from 0.105 mm to 0.045 mm, the tensile strength exhibited a 79.8% increase, while the flexural strength increased by 24.8%. The investigation of bending performance in the aforementioned study solely considered the impact of single-layer thickness, overlooking the effects of total layer count and overall thickness. Therefore, it is essential to investigate the influence of layer count and thickness on bending performance.

Despite the valuable contributions of the above-mentioned studies, certain aspects remain relatively underexplored. First, most existing studies on flexural behavior have focused on symmetric, quasi-isotropic, or surface-damaged laminates, while fewer investigations have systematically addressed pristine unidirectional laminates under three-point bending. Second, although the effects of hybridization and stacking sequence on mechanical properties have been widely studied, the effects of stacking thickness and interlaminar shear strength (ILSS) on damage evolution have received relatively less attention. Third, while experimental observations of final failure modes are common, quantitative characterization of the damage evolution process, particularly the progressive damage of the interlaminar interface, is still limited due to experimental difficulties. Finite element analysis (FEA) is considered to be helpful for the analysis of material properties.^33,34,35 It is crucial to develop numerical simulation methods to predict the mechanical response of laminates, particularly for conducting in-depth analyses of the damage and failure processes of the interlaminar interface.

Despite the valuable contributions of previous studies, three specific research gaps remain unaddressed. First, most existing flexural studies focus on symmetric or quasi-isotropic laminates; the progressive damage evolution of pristine unidirectional (UD) CFRP laminates under three-point bending has not been systematically quantified. Second, although the effects of hybridization and stacking sequence on final strength are well documented, how stacking thickness and interlaminar shear strength (ILSS) influence the damage evolution process (rather than just peak load) remains unclear. Third, and most importantly, no quantitative model exists to describe the catastrophic failure behavior of the interlaminar interface from damage initiation to complete fracture, limiting the mechanistic understanding of bending-induced delamination.

To fill these gaps, this work combines finite element simulations with three-point bending to achieve three objectives:

(1) To systematically quantify the effects of stacking thickness (4, 6, and 8 layers) and ILSS (10–80 MPa) on both the flexural properties and the progressive damage evolution of pristine UD CFRP laminates.

(2) To establish an interlaminar interface catastrophic failure model based on Taylor expansion, which reveals that interface damage follows a power-law characteristic with an exponent of 0.5.

(3) To quantitatively link the damage rate to stacking thickness, demonstrating that at ILSS20 MPa, thicker laminates exhibit faster damage progression and more brittle failure.

These findings provide new insights into the damage evolution mechanisms of CFRP laminates and offer practical guidance for structural design.

Numerical model

Three-point bending finite element model

Under three-point bending, the damage and failure of CFRP laminates occur rapidly, making it difficult to capture the progressive damage evolution process experimentally in real time. Therefore, finite element analysis provides an effective alternative to investigate the damage evolution mechanisms with high temporal and spatial resolution. In this study, the three-point bending models are established using Abaqus. Five different sizes of CFRP laminates are constructed according to Table 1. All fiber bundles are arranged axially along the X-direction, as illustrated in Figure 1. In the figure, “Ply” represents the fiber layer and “Coh” represents the interface between adjacent fiber layers. “Ply-1” is the lower layer, while “Ply-4” (C₄), “Ply-6” (C₆) and “Ply-8” (C₈) are the upper layers of each laminate. “Coh-2” (C₄), “Coh-3” (C₆) and “Coh-4” (C₈) are the interfaces of the neutral layers of each laminate. The loading nose and supports have a diameter, D, of 10 mm. A reference span-to-thickness ratio of 60:1 is selected in accordance with ASTM D7264,³⁶ which permits ratios of 16:1, 20:1, 32:1, 40:1, and 60:1 for flexural testing. To systematically investigate the individual and coupled effects of thickness and span, five specimen configurations are designed as listed in Table 1. These models are grouped into three categories: (1) constant span-to-thickness ratio group (S:T = 60:1) including C₄L₃₆, C₆L₅₄, and C₈L₇₂, which follow ASTM D7264³⁶ for simulation-experiment validation but involve simultaneous changes in thickness and span; (2) constant thickness group comparing C₄L₃₆ versus C₄L₅₄ and C₈L₅₄ versus C₈L₇₂ to isolate the pure span effect; and (3) constant span group comparing C₄L₅₄, C₆L₅₄, and C₈L₅₄ (span fixed at 54 mm) to isolate the pure thickness effect. This design decouples thickness effects from structural scale effects, directly addressing the concern that changing span length might bias the observed damage evolution.

Table 1.

Specimen model parameters.

Name	Length × width × Thickness,L × W × T (mm)	Span, S (mm)	Span-to-Thickness, S:T
C₄L₃₆	50 × 13 × 0.6	36	60:1
C₄L₅₄	65 × 13 × 0.6	54	90:1
C₆L₅₄	65 × 13 × 0.9	54	60:1
C₈L₅₄	65 × 13 × 1.2	54	45:1
C₈L₇₂	90 × 13 × 1.2	72	60:1

Figure 1.

Geometric models of laminates with different stacking thickness under three-point bending.

Damage constitutive of composite laminates

The damage modes of composite laminates are complex, including fiber damage, matrix damage and interface cracking.³⁷ Therefore, a constitutive model based on in-plane failure and interlaminar failure is introduced in this study. Damage in composite laminates arises from the interplay between interlaminar and in-plane mechanisms. Interlaminar damage, characterized by delamination at ply interfaces, is typically triggered by shear stresses and weak interfacial bonding. In contrast, in-plane damage, such as fiber breakage and matrix cracking within individual plies, is governed by the material’s tensile and compressive strength.

In-plane damage constitutive model

The 2D Hashin damage criteria are employed as the fracture and damage criterion for CFRP laminates, which is suitable for unidirectional fiber laminates. The required material parameters are outlined in Table 2. The volume content of epoxy resin in laminates is 40%. The Hashin damage initiation criteria are provided in Table 3.

Table 2.

Parameters of CFRP laminates for simulation⁴.

Properties	Value
Elastic modulus (GPa)	$E_{1} = 125$ , $E_{2} = 7.8$
Poisson’s ratio	$v_{12} = 0.32$
Shear modulus (GPa)	$G_{12} = G_{13} = 4.0$ , $G_{23} = 3.6$
Density (kg/m³)	$ρ = 2032$
Strength (MPa)	$X^{T} = 1800$ , $X^{C} = 1200$ , $Y^{T} = 50$ , $Y^{C} = 150$ , $S^{L} = S^{T} = 93$
Fracture energy (kJ/m²)	$G_{11}^{T} = 40$ , $G_{11}^{C} = 40$ , $G_{22}^{T} = 0.25$ , $G_{22}^{C} = 0.75$

Table 3.

Hashin damage initiation criteria³⁸.

Failure modes	Criteria
Fiber tension $(σ_{11} > 0)$	$F_{f}^{t} = {(\frac{σ_{11}}{X^{T}})}^{2} + α {(\frac{τ_{12}}{S^{L}})}^{2} \geq 1$
Fiber compression $(σ_{11} < 0)$	$F_{f}^{c} = {(\frac{σ_{11}}{X^{C}})}^{2} \geq 1$
Matrix tension $(σ_{22} > 0)$	$F_{m}^{t} = {(\frac{σ_{22}}{Y^{T}})}^{2} + {(\frac{τ_{12}}{S^{L}})}^{2} \geq 1$
Matrix compression $(σ_{22} < 0)$	$F_{m}^{c} = {(\frac{σ_{22}}{2 S^{T}})}^{2} + [{(\frac{Y^{C}}{2 S^{T}})}^{2} - 1] \frac{σ_{22}}{Y^{C}} + {(\frac{τ_{12}}{S^{L}})}^{2} \geq 1$

In Table 2 and Table 3, σ₁₁, σ₂₂ and τ₁₂ are the stress tensions, X^T, X^C represent the tensile and compressive strength in the x-direction, and Y^T and Y^C represent the tensile and compressive strength in the y-direction. Additionally, S^L and S^T represent the longitudinal and transverse shear strength. The coefficient α determines the contribution of the shear stress to the fiber tensile initiation criterion.

Before damage occurs, the material behaves linearly elastic. When damage reaches the initiation criterion, the material stiffness degrades, and the damage variable at this time can be expressed as d_i.

d_{i} = \frac{u_{i e q}^{f} (u_{i e q} - u_{i e q}^{0})}{u_{i e q} (u_{i e q}^{f} - u_{i e q}^{0})} (i = f, m)

(1)

where

u_{i e q}

is the equivalent displacement,

u_{i e q}^{0}

is the initial equivalent displacement and

u_{i e q}^{f}

is the displacement at which the material is completely damaged in each failure mode. The evolution of damage variable is illustrated in Figure 2. The fracture energy corresponds to the area of the triangle in Figure 2.

Figure 2.

Relationship between equivalent stress and equivalent displacement.

Interlaminar interface damage constitutive model

The cohesive zone model is widely used in the simulation of delamination in composite laminates.^39,40 The parameters utilized in this study are presented in Table 4. Interface damage follows the traction-separation law and the interface damage initiation criterion is as follows:

{\frac{〈 t_{n} 〉}{t_{n}^{0}}}^{2} + {\frac{t_{s}}{t_{s}^{0}}}^{2} + {\frac{t_{t}}{t_{t}^{0}}}^{2} = 1

(2)

where

t_{n}^{0}

t_{s}^{0}

, and

t_{t}^{0}

represent Mode I, Mode II and Mode III failure strength. The symbol

〈 〉

represents the Macaulay bracket operator, which is defined for every

〈 b 〉 = (b + | b |) / 2

. For linear softening, an evolution of the damage variable, D, is used that reduces to

D = \frac{δ_{m}^{f} (δ_{m}^{\max} - δ_{m}^{0})}{δ_{m}^{\max} (δ_{m}^{f} - δ_{m}^{0})}

(3)

where

δ_{m}^{f} = 2 G^{C} / T_{e f f}^{0}

with

T_{e f f}^{0}

as the effective traction at damage initiation.

δ_{m}^{\max}

refers to the maximum value of the effective displacement attained during the loading history. Damage evolution follows the Benzeggagh–Kenane (BK) law,⁴¹ which defines the fracture energy as a function of the mode mix.

G_{n}^{C} + (G_{s}^{C} - G_{n}^{C}) {\frac{G_{S}}{G_{T}}}^{η} = G^{C}

(4)

where

G_{S} = G_{s} + G_{t}

G_{T} = G_{n} + G_{S}

, and

η

is a BK exponent. The BK law is a semi-empirical criterion, and its BK exponent can be determined through fitting.⁴¹ In this study, the BK exponent

η = 1.45

is adopted from the literature.⁴ The linear softening law is selected for two main reasons. First, linear softening is one of the simplest and most widely used formulation for delamination simulations, balancing accuracy and computational efficiency.⁴² Second, under less favorable numerical conditions, discontinuities in the bilinear constitutive law caused instabilities, whereas the linear softening law does not suffer from this issue.⁴³ Therefore, this approach is appropriate for the comparative parametric study in this work. The penalty stiffness of the zero-thickness cohesive layer should be sufficiently high to prevent artificial compliance introduced by the cohesive layers, but not so high as to cause numerical difficulties or excessive computational cost.⁴⁴ The penalty stiffness value of 10⁶ N/mm³ adopted in this study follows the common practice in the literature,⁴ which falls within the typical range of 10⁴ to 10⁸ N/mm³⁴⁵ and is therefore considered reasonable.

Table 4.

Mechanical properties of cohesive layers⁴.

Properties	Symbol	Value	Unit
Density	ρ	1200	kg/m³
Normal penalty stiffness	$K_{n n}^{0}$	10⁶	N/mm³
First shear penalty stiffness	$K_{s s}^{0}$	10⁶	N/mm³
Second shear penalty stiffness	$K_{t t}^{0}$	10⁶	N/mm³
Mode I tensile strength	$t_{n}^{0}$	10, 20, 40, 50, 80	MPa
Mode II shear strength	$t_{s}^{0}$	10, 20, 40, 50, 80	MPa
Mode III shear strength	$t_{t}^{0}$	10, 20, 40, 50, 80	MPa
Mode I fracture toughness	$G_{I C}$	0.25	N/mm
Mode II fracture toughness	$G_{I I C}$	0.75	N/mm
Mode III fracture toughness	$G_{I I I C}$	0.75	N/mm
BK exponent	$η$	1.45	-
Initial thickness	$T_{c}$	1	mm

Additionally, considering the influence of ILSS on flexural properties, a parametric analysis of ILSS (

t_{n}^{0}

t_{s}^{0}

t_{t}^{0}

) is conducted. Assuming that the interlaminar strengths in different modes (I/II/III) are equal,

t_{n}^{0} = t_{s}^{0} = t_{t}^{0}

, the shear strength values are varied among 10, 20, 40, 50 and 80 MPa to investigate the effect of different shear strengths through simulation. This ILSS range is selected based on literature values for CFRP laminates. As shown in Table 5, reported ILSS values for CFRP composites range from approximately 20 MPa to 170 MPa. Therefore, the selected range of 10–80 MPa covers typical ILSS values for practical CFRP laminates^46,47,48,49 and extends to lower values to systematically investigate the transition from delamination-dominated failure (low ILSS) to fiber fracture-dominated failure (high ILSS).

Table 5.

Literature values of interlaminar shear strength (ILSS) for CFRP laminates.

Source	ILSS values (MPa)
Ren et al.⁴⁶	20.3, 28.5, 36.2, 40.3
Yi et al.⁴⁷	38.8, 45.2, 46.7
Perogamvros et al.⁴⁸	49, 55, 64, 66
Barbero⁴⁹	69.6, 76.8, 137.2, 170

Mesh and boundary condition

The boundary conditions are as follows. As shown in Figure 3, face1 and its opposite face: U2 = 0; Face2 and its opposite face: U2 = 0, UR1 = 0, UR3 = 0; Supports: U1 = U2 = U3 = UR1 = UR2 = UR3 = 0. The displacement load is applied in the U3 direction. U1, U2 and U3 are the axial displacements in 1, 2 and 3 directions, respectively. UR1, UR2, and UR3 are rotations on the 1, 2, and 3 axes, respectively. The contact properties between the loading nose and the laminate are set as face-to-face contact. The fiber layer utilizes an eight-node continuous shell element (SC8R), and the interlaminar interface employs an 8-node three-dimensional cohesive element (COH3D8) with a thickness of zero. As shown in Figure 4, simulations of the C₆ laminate with different mesh sizes indicate that the results stabilize when the mesh size is smaller than 2 mm. The peak load difference between the 0.5 mm and 0.8 mm meshes is less than 1‰. Therefore, balancing computational time and accuracy, a mesh size of 0.5 mm × 0.5 mm is adopted.

Figure 3.

Finite element mesh and boundary conditions.

Figure 4.

The variation of the peak load with mesh size for the C₆L₅₄ laminate with an ILSS of 50 MPa.

In order to ensure calculation accuracy and prevent non-convergence, an explicit dynamics solver is employed for quasi-static analysis. Mass scaling is not used in the simulations. A smooth step loading function is applied with a total loading time of 1 s. Furthermore, the quasi-static condition is verified by examining the kinetic energy relative to the internal energy. As presented in Figure 5, the kinetic energy constituted only 0.28%, 0.09%, and 0.46% of the internal energy at the peak load for the C₄L₃₆, C₆L₅₄, and C₈L₇₂ laminates, respectively. Such low ratios justify that the kinetic energy is negligible and the analysis accurately represents quasi-static behavior.

Figure 5.

Changes in internal energy and kinetic energy over time in simulation.

Experiment

Materials and sample preparation

The CFRP unidirectional laminates with three different stacking thicknesses used in this study are manufactured by Dezhou Blue Label Composite Materials Technology Co., Ltd. According to the manufacturer’s description, the laminates are produced using a hot-press consolidation process with wet lay-up techniques. The thickness of a single fiber layer is 0.15 mm. The manufacturing equipment is shown in Figure 6. The material parameters provided by the manufacturer are presented in Table 6.

Figure 6.

Manufacturing equipment: (a) Autoclave; (b) Cutting machine.

Table 6.

Material parameters of laminates.

Properties	Value
Fiber elastic modulus (GPa)	230
Fiber elongation (%)	2.1
Fiber tensile strength (MPa)	4900
Fiber density (g/cm³)	1.8
Weave pattern	Unidirectional pattern
Curing temperature (°C)	120–140
Glass transition temperature (Tg/°C)	120–130
Resin volume content (%)	33

The preparation process involves the following steps: Mixing the resin and hardener evenly, and setting it aside for later use. Placing a layer of T700 carbon fiber cloth in the mold and using a roller to evenly spread the mixed resin on the surface. Repeating the steps of stacking carbon fiber cloth and applying resin until the desired number of layers is achieved. Place the stacked materials in a vacuum box to remove residual air bubbles from the resin. Place the degassed material into an autoclave for solidification at a temperature of 150°C and a pressure of 1.5 MPa. After 10 h, when curing is complete, taking out the CFRP laminate and cutting it.

The fiber volume content of the CFRP laminate prepared in this experiment is 67%. The fiber volume content used in the simulation is 60%. The effect of fiber volume fraction is discussed in Section 4.1. The mechanical properties of the CFRP laminate in the experiment are higher than those given by the simulation data, with a difference value ≤12%. The objective of this study is to investigate the patterns of quasi-static progressive damage evolution, rather than to achieve accurate strength predictions. Therefore, despite this difference, and considering the experimental conditions and simulation convenience, no corrections are made to the simulation data. It’s important to note that this discrepancy does not affect the comparison of various parameters, as the comparison is conducted within the experimental group and the simulation group.

Flexural test

According to ASTM D7264 testing standards,³⁶ a three-point bending test is conducted. The diameter of the loading nose and supports is 10 mm, and the span-to-thickness ratio (S:T) of all the specimens is 60:1. The quasi-static loading rate is 1 mm/min. The thickness of the prepared laminates is as follows: C₄ (0.6 mm), C₆ (0.9 mm), and C₈ (1.2 mm), as shown in Figure 7. Each group tested 4 samples and the flexural strength is calculated as follows:

σ_{f} = \frac{3 P_{t} S}{2 W T^{2}}

(5)

where P_t represents the peak load.

Figure 7.

(a) Side view of the specimens; (b) Top view of the specimens; (c) Test equipment.

Results and discussion

Load-displacement curves and damage evolution

Firstly, a comparison is made between the simulated and experimental load-displacement curves under three different stacking thicknesses using the constant span-to-thickness ratio models (C₄L₃₆, C₆L₅₄, C₈L₇₂) that follow ASTM D7264. The additional constant span models (C₄L₅₄, C₆L₅₄, C₈L₅₄) are then used to isolate the pure thickness effect, as discussed in Section 4.4. Subsequently, numerical simulations are conducted to analyze the damage evolution process of C₄L₃₆ laminate. The notations in Figure 8 are defined as follows: fiber tensile initiation criterion (FT), fiber compression initiation criterion (FC), and interlaminar interface damage (SDEG). Their values range from 0 to 1. Figure 8(a) illustrates the simulated load-displacement curves of the laminates. It is evident that with increasing stacking thickness, the peak load and fracture displacement of the laminate also increase. Figure 8(b)-(d) present a comparison of the simulation and experimental load-displacement curves of C₄L₃₆, C₆L₅₄, and C₈L₇₂ laminates. The simulation and experiment curves exhibit a good fit. It is noteworthy that the peak load and fracture displacement of the simulation curves are smaller than those of the experiment. This discrepancy may be attributed to the inconsistency in fiber volume content between the simulation and the experimental setup, as mentioned in section 3.1. The fiber volume content is 67% in the experiment and 60% in the simulation. According to the rule of mixtures formula, the relationships between the longitudinal modulus E₁, the longitudinal strength σ₁, and the fiber volume fraction V_f can be obtained as follows⁴⁹:

E_{1} = V_{f} E_{f} + (1 - V_{f}) E_{m}

(6)

σ_{1} = V_{f} σ_{f} + (1 - V_{f}) σ_{m}

(7)

Where E_f is the elastic modulus of the fiber, E_m is the elastic modulus of the matrix, σ_f is the strength of the fiber, σ_m is the stress in the matrix at the strain corresponding to fiber failure. For fiber-reinforced polymer laminates, the modulus and strength of the fibers are significantly higher than those of the matrix, resulting in an approximately linear relationship between the longitudinal modulus and strength of the composite and the fiber volume fraction.⁴⁹ An 11.6% increase in the longitudinal modulus and strength of CFRP is expected when the fiber volume fraction rises from 60% to 67%. As shown in Table 7, increasing the fiber volume fraction from 60% to 67% reduces the maximum simulation error from −12.1% to −2.7%, confirming that the discrepancy is primarily attributable to the difference in fiber volume fraction. It should be emphasized that the primary objective of this study is to investigate the trends and mechanisms of quasi-static progressive damage evolution, rather than to achieve accurate absolute strength predictions. Therefore, despite the improved agreement at 67%, the material parameters corresponding to a fiber volume fraction of 60% are consistently used throughout all simulations. This choice does not affect the validity of the comparative parametric analyses, as all comparisons are conducted under consistent material properties.

Figure 8.

(a) Simulated load-displacement curves of the laminates; (b) Comparison of simulation and experimental load-displacement curves for C₄L₃₆ laminate; (c) Comparison of simulation and experimental load-displacement curves for C₆L₅₄ laminate; (d) Comparison of simulation and experimental load-displacement curves for C₈L₇₂ laminate; (e) Damage evolution cloud diagram of C₄L₃₆ laminate with ILSS of 50 MPa.

Table 7.

Effect of fiber volume fraction on simulated peak load with ILSS of 50 MPa.

Properties	ExperimentPeak load (N)	Simulation $V_{f} = 60 %$		Simulation $V_{f} = 67 %$
Properties	ExperimentPeak load (N)	Peak load (N)	Error	Peak load (N)	Error
C₄L₃₆	182 ± 2.9	160	−12.1%	178	−2.2%
C₄L₅₄	/	106	/	121	/
C₆L₅₄	254 ± 3.4	229	−9.8%	247	−2.7%
C₈L₅₄	/	406	/	442	/
C₈L₇₂	310 ± 8.1	286	−7.7%	309	−0.1%

Figure 8(e) illustrates the damage initiation and evolution process of each layer (Ply) and interface (Coh) of the C₄L₃₆ laminate under simulation conditions with an interlaminar shear strength (ILSS) of 50 MPa. The C₆L₅₄ and C₈L₇₂ laminate exhibit similar phenomena; hence, the damage initiation and evolution process of only C₄L₃₆ is presented. δ₀, δ₁, and δ₂ represent different displacement loadings, corresponding to Figure 8(a). When the displacement is loaded to δ₁, damage initiation is relatively large below the loading nose of the upper and lower layers (Ply-4, Ply-1), and damage occurs at the neutral plane Coh-2. After reaching displacement δ₂, fiber layers Ply-1 and Ply-2 below the neutral plane are fractured, resulting in interface Coh-1 cracking, and the damage evolution of interface Coh-1 is rapid. This may be attributed to the higher ILSS, which could lead to a less distinctive evolution process of interface damage. Therefore, the subsequent analysis in the following sections assesses the impact of ILSS on laminate damage evolution.

Influence of ILSS on damage and failure

Differences in interface strength can influence the fracture mode. Therefore, numerical simulations are employed to conduct a parametric analysis of the ILSS. In the cohesive zone model, $t_{n}^{0}$ , $t_{s}^{0}$ , and $t_{t}^{0}$ represent the ILSS. For ease of calculation, a total of five ILSS values (10, 20, 40, 50, and 80 MPa) are considered when the stacking number is 4 layers (C₄L₃₆). The fracture mode under an ILSS of 50 MPa has been introduced in the previous context, and the discussion of the other four ILSS values will follow.

Figure 9 illustrates the load-displacement curves of the C₄L₃₆ laminate under five different ILSS values. The curves demonstrate that as the ILSS increases from 10 MPa to 40 MPa, both fracture displacement and peak load increase. However, the fracture displacement and peak load remain essentially unchanged as the ILSS increases from 40 MPa to 80 MPa. Key positions are marked: δ₀ represents a specific moment in the rising section of the curve before the peak load, δ₁ represents the peak load of each curve, δ₂ represents the fracture, δ_1.4 and δ_1.7 represent two displacement loading moments between δ₁ and δ₂ in ILSS80.

Figure 9.

Load-displacement curves under different interlaminar shear strengths (ILSS).

According to Figure 10, it is observed that the Mises stress at the neutral plane (Coh-2) is the highest among all interfaces. The figure only provides Mises stress under ILSS10/20 MPa, but the stress distribution under ILSS40/50/80 MPa is generally consistent with that under ILSS20 MPa. From this, it can be inferred that the neutral plane is prone to damage.

Figure 10.

Interlaminar interface Mises stress: ILSS10/20 MPa.

Delamination is a typical failure mode of composite laminates. High interlaminar shear forces or low ILSS can lead to delamination, resulting in the material losing its load-bearing capacity.⁵⁰ Figure 11(a) illustrates the damage evolution process under ILSS10 MPa. At the moment δ₀, Coh-2 is damaged first. As the displacement progresses to δ₁, Coh-1 and Coh-3 also experience damage, but due to the greater extent of damage in Coh-2, it leads to the failure of Coh-2 at δ₂. At this point, Ply-1 and Ply-4 have not reached damage initiation. Therefore, the simulation suggests that the laminate experiences delamination failure under this condition. Additionally, observing Figure 9, the slope of the rising section of the ILSS10 curve does not change, and no stiffness attenuation occurs. This is because the Ply-4 fiber layer has not been damaged, which means that the part where the loading nose contacts the upper surface layer remains intact.

Figure 11.

Damage evolution cloud diagram and fracture morphology of C₄L₃₆ laminate: (a) ILSS10 MPa; (b) ILSS20 MPa; (c) ILSS40 MPa; (d) ILSS80 MPa; (e) Fracture morphology under different ILSS.

In Figure 11(b), the ILSS is increased to 20 MPa. It is observed that although Coh-2 experiences the most damage at moments δ₀ and δ₁, Coh-1 cracks first at moment δ₂. This is because at moment δ₂, the fibers in Ply-1 fracture, leading to interface cracking. Additionally, the location of interface damage under ILSS20 MPa is different compared to that under ILSS10 MPa. The damage location of each layer interface under ILSS10 MPa is in the middle area between the loading nose and the support, while the damage location under ILSS20 MPa is directly below the loading nose. Observing Figure 10, the maximum stress distribution under ILSS10 MPa is mainly between the loading nose and the support, while under ILSS20 MPa, it is more concentrated under the loading nose. Additionally, combining with Figure 9, it is evident that under ILSS20 MPa, there is stiffness degradation in the peak load region of the curve. This is attributed to the compression damage (FC) initiation criterion in the fibers of Ply-4 reaching 1 at moment δ₀.

In Figure 11(c)-(d), it is observed that when the ILSS rises to 40 MPa or 80 MPa, it becomes challenging to detect the failure of the interface, and only minor damage occurs when the displacement is loaded before the peak load. Beyond ILSS40 MPa, the peak load and fracture displacement no longer increase, as at this point, the factor controlling the material strength becomes the strength of the fiber layer.

Figure 11(e) illustrates the final fracture morphology of the laminate under different ILSS. The failure mode of ILSS10 MPa laminates is interface delamination cracking, ILSS20 and 40 MPa result in interface cracking when fiber fracture occurs, and ILSS80 MPa primarily exhibits fiber fracture. Consequently, the simulation results suggest that when the interface strength is weak, the fracture mode of the laminate tends to be delamination failure. Conversely, if the interface strength is high, the laminate is more susceptible to fiber fracture.

Figure 12 presents a comparison between the simulation and experimental failure modes of C₄L₃₆ laminates. The experimental laminates predominantly exhibit fiber fracture, which is qualitatively consistent with the simulation result of ILSS80. This suggests good interlaminar bonding quality, although direct ILSS measurement is not performed.

Figure 12.

Comparison of simulation and experimental fracture morphology of C₄L₃₆ laminates.

Under higher ILSS conditions (40–80 MPa), the interface damage rate is extremely fast, making it difficult to capture the progressive damage process numerically. Therefore, an ILSS of 20 MPa is selected for the parametric study, as the relatively slower damage rate allows for clearer observation and quantitative characterization. To further investigate the influence of stacking thickness on interface damage, three different stacking thicknesses are established by FEA, namely, C₄L₃₆, C₆L₅₄, and C₈L₇₂, under ILSS20 MPa. The following analysis and conclusions are specific to this ILSS condition.

Interlaminar interface catastrophic failure model

The failure mode of CFRP laminates is typically considered as brittle fracture, and the damage occurring at the interface of each layer can be seen as catastrophic failure. Therefore, a catastrophic damage model is employed to explore the damage evolution rules of the interface of laminates under different stacking thickness.^51,52 The purpose of the “interlaminar interface catastrophic failure model” is to analyze and predict the damage evolution and failure behavior of the interlaminar interface in CFRP laminates under bending loads. By characterizing the brittle fracture mechanism, the model identifies the power-law relationship governing damage progression and its rapid increase near the fracture point. It aids in understanding the sensitivity of damage rate to displacement, providing insights for optimizing laminate design to enhance mechanical performance and durability under varying stacking thicknesses and interlaminar shear strengths.

Damage definition

The initial damage N₀ is defined as the number of elements with an interface damage degree SDEG exceeding a given threshold at the point of deviation from linearity during the loading process, and the complete damage N_f is defined as the number of elements with an interface damage degree SDEG exceeding the threshold at the peak load. To examine the sensitivity of the damage statistics to the threshold selection, three SDEG values (0.6, 0.7, and 0.8) are analyzed and compared. The results for the SDEG threshold of 0.7 are shown in Figure 13(a). The initial damage variable is represented by $D_{0} = N_{0} / N_{f}$ , and the complete damage variable is represented by $D_{f} = 1$ . In the process of three-point bending loading, displacement δ affects the interface damage, so displacement δ is selected as the control variable, and the number of damage units corresponding to displacement δ is used to calculate the damage variable $D = N / N_{f}$ , where the displacement δ_f under peak load is called the fracture displacement corresponding to the complete damage N_f, and the initial displacement δ₀ corresponds to the initial damage N₀.

Figure 13.

(a) Damage evolution cloud diagram at the neutral plane of laminates under different stacking thickness; (b) Relationship between normalized damage rate and displacement; (c) Relationship between damage rate and normalized displacement; (d) Logarithmic relationship between damage rate and displacement.

The normalized damage variable $D / D_{f}$ increases as the normalized displacement $δ / δ_{f}$ increases, as demonstrated in Figure 13(b). The various symbols in the figure represent the simulation outcomes for different stacking thickness. It can be seen from Figure 13(b) that damage begins to appear with the increase of displacement and increases dramatically as the displacement approaches the peak value. Considering that the load of laminates drops instantaneously when the interface cracks, catastrophic failure forms will naturally occur. As a result, a catastrophic failure model is utilized to analyze the progression of damage in the interlaminar interface.

Damage model: If the control variable of a system, such as the laminate displacement δ, is continuous and differentiable before catastrophic rupture, then damage D can be expanded using Taylor’s formula as

δ = δ_{f} + {[\frac{d δ}{d D}]}_{D_{f}} (D_{f} - D) + \frac{1}{2} {[\frac{d^{2} δ}{d^{2} D}]}_{D_{f}} {(D_{f} - D)}^{2} + o {(D_{f} - D)}^{2}

(8)

where δ represents the displacement at the interface directly below the loading nose, and the damage variable D is calculated based on the number of damage elements (SDEG > threshold). When

δ = δ_{f}

and

D = D_{f}

, the above equation (8) is satisfied. When the control variable approaches the catastrophic point, the damage D increases rapidly and the damage rate R grows sharply, that is,

\lim_{δ \to δ_{f}} R = \lim_{δ \to δ_{f}} \frac{d D}{d δ} = \infty

, contrarily,

\lim_{D \to D_{f}} \frac{d δ}{d D} = 0

{[\frac{d δ}{d D}]}_{D_{f}} = 0

, substituting the relationships into equation (8) and ignoring the terms higher than the second order, the damage as a function of the control variable can be expressed as

D = D_{f} - λ {(δ_{f} - δ)}^{0.5}

(9)

where

λ = {- \frac{1}{2} {[\frac{d^{2} δ}{d^{2} D}]}_{D_{f}}}^{- 0.5}

. When

δ = δ_{0}

D = D_{0}

, substituting them into equation (9), the damage coefficient

λ = \frac{D_{f} - D_{0}}{\sqrt{δ_{f} - δ_{0}}}

. Equation (9) shows that damage D increases with displacement δ and follows a power-law relationship with an exponent of 0.5. Importantly, this exponent of 0.5 is derived theoretically from the second-order Taylor expansion of the damage variable near the catastrophic failure point (Eqs. (8) and (9)), not obtained by curve fitting. The square-root singularity is a characteristic of critical phenomena near instability. Notably, the same power-law behavior with an exponent of 0.5 has been reported for damage evolution in ceramic coating systems under three-point bending⁵¹ and under thermal loading,⁵² where damage was measured using completely different physical quantities (crack length and temperature). This consistency across different materials, loading conditions, and measurement methods strongly supports that the 0.5 exponent is a generic physical law of catastrophic failure in brittle materials, rather than an arbitrary fitting parameter.

In addition, the interface damage rate of the laminates can be obtained from equation (9).

R = \frac{d D}{d δ} = \frac{λ}{2} {(δ_{f} - δ)}^{- 0.5}

(10)

Equation (10) shows that when δ approaches δ_f, the damage rate R increases sharply.

Comparison of simulation results and damage model

Equation (11) is derived from equation (9). It can be seen from Figure 13(b) that the simulation results (symbols) agree with the theoretical power-law function curve of equation (11). It can be observed that although the initial damage points of each stacking thickness are different and the dispersion is large, the dispersion of all damage points becomes smaller when the interface tends to be completely damaged, and the damage model is more effective. For catastrophic damage, it is difficult to obtain the failure point near the fracture displacement. Therefore, by collecting partial failure points using this model, it is possible to infer the damage failure pattern.

\frac{D}{D_{f}} = 1 - [\frac{1 - \frac{D_{0}}{D_{f}}}{\sqrt{1 - \frac{δ_{0}}{δ_{f}}}}] {(1 - \frac{δ}{δ_{f}})}^{0.5}

(11)

As shown in Table 8, three SDEG thresholds (0.6, 0.7, and 0.8) are analyzed and compared. The coefficients of variation (CV) for the damage coefficient λ across the three thresholds are below 5% for all laminates. These low CV values indicate that the selection of the SDEG threshold has a negligible effect on the damage evolution equation parameters, confirming that the reported results are not threshold-dependent. This can also be visually observed in Figure 13(b)-(d), where the damage evolution curves for different SDEG thresholds show small differences. Any threshold value between 0.6 and 0.8 yields the same power-law exponent and damage evolution characteristics; the value of 0.7 is selected empirically as a representative threshold. Crucially, the power-law exponent of 0.5 is not a fitting parameter but is derived theoretically from the second-order Taylor expansion (equations (8) and (9)). The simulation data in Figure 13(b) agree with this theoretical curve, validating the derivation rather than indicating a curve-fitting exercise. Furthermore, correlation analysis is performed on the results obtained with the SDEG threshold of 0.7 across the three stacking thicknesses. The absolute values of the correlation coefficients are all greater than 0.97. These strong correlations suggest that the damage parameters are closely associated with stacking thickness in the simulation. Despite the different absolute span lengths (36 mm, 54 mm, and 72 mm) resulting from the constant span-to-thickness ratio of 60:1, the damage evolution for all three stacking thicknesses follows a consistent trend. It can be observed that at ILSS20 MPa, as the stacking thickness increases, the damage coefficient λ becomes larger, indicating a faster interface damage rate. Similarly, with a thicker thickness, the appearance of initial damage points is closer to the fracture displacement, suggesting that as the material becomes thicker, the brittle failure mode becomes more pronounced in the simulation.

Table 8.

Parameters of the damage evolution equation for three stacking thicknesses with different SDEG thresholds.

SDEG	C₄L₃₆				C₆L₅₄				C₈L₇₂				Correlation (SDEG0.7)
SDEG	0.6	0.7	0.8	CV	0.6	0.7	0.8	CV	0.6	0.7	0.8	CV	Correlation (SDEG0.7)
$D_{0} / D_{f}$	0.422	0.406	0.385	4.59%	0.127	0.118	0.115	5.20%	0.0047	0.0030	0.0023	37.03%	−0.9706
$δ_{0} / δ_{f}$	0.905	0.905	0.905	0.00%	0.939	0.939	0.939	0.00%	0.955	0.955	0.955	0.00%	0.9791
λ	1.83	1.88	1.95	3.19%	3.53	3.57	3.581	0.75%	4.69	4.70	4.71	0.21%	0.9935
α	−0.063	−0.06	2E-5	−86.7%	0.57	0.58	0.58	1.00%	0.853	0.85	0.855	0.30%	0.9735
b	0.5	0.5	0.5	0.00%	0.5	0.5	0.5	0.00%	0.5	0.5	0.5	0.00%	-

Figure 13(c) shows the normalized damage rate curve obtained with the SDEG threshold of 0.7, and the expression is shown in equation (12). It can be seen that when approaching the fracture displacement, the interface damage rate increases dramatically and theoretically becomes unbounded. Figure 13(c) solely illustrates the power-law singularity of the damage sensitivity without providing a comparison between the simulation results and the dispersion of the damage model. In addition, by fitting the simulation results, the logarithmic relationship between damage rate and displacement is

\ln R = a - b \ln (1 - δ / δ_{f})

, where a is the intercept and b is the power exponent. From Table 9, the power exponent b of each simulation result is 0.5, but the intercepts corresponding to different simulation results are scattered.

R = [\frac{1 - \frac{D_{0}}{D_{f}}}{2 \sqrt{1 - \frac{δ_{0}}{δ_{f}}}}] {(1 - \frac{δ}{δ_{f}})}^{- 0.5}

(12)

Table 9.

Simulated and experimental flexural properties.

Properties	δ (mm)			P_t (N)			σ_f (MPa)
Properties	ILSS40	ILSS20	Exp	ILSS40	ILSS20	Exp	ILSS40	ILSS20	Exp
C₄L₃₆	5.8	5.6	7.8 ± 0.51	160	157	182 ± 2.9	1841	1807	2099 ± 34
C₄L₅₄	13.3	13.5	/	106	105	/	1835	1817	/
C₆L₅₄	8.0	7.9	9.3 ± 0.23	229	225	254 ± 3.4	1761	1734	1951 ± 26
C₈L₅₄	6.1	4.6	/	406	345	/	1757	1493	/
C₈L₇₂	10.9	10.3	14.2 ± 0.50	287	283	310 ± 8.1	1654	1632	1790 ± 47

Figure 13(d) shows the logarithmic relationship between the simulation results and the damage model. The figure suggests that the interface damage rate follows a power-law singularity of −0.5 as indicated by the model, for various stacking thicknesses.

Influence of stacking thickness and ILSS on flexural properties

Understanding the influence of stacking laminate thickness and ILSS on flexural properties contributes to the design of composite materials. Based on the simulation results, when ILSS is increased to 40 MPa, the load-displacement curve of the material remains unchanged. Therefore, ILSS20 and ILSS40 are selected for analysis. Table 9 presents the average mechanical properties of each specimen under simulation and experimental conditions, where “Exp” represents “experimental.”

In combination with Figure 14(a), for the constant span-to-thickness ratio group (C₄L₃₆, C₆L₅₄, C₈L₇₂), which follows ASTM D7264 for simulation-experiment validation, the fracture displacement and peak load increase with stacking thickness. Changing ILSS from 20 MPa to 40 MPa leads to limited improvements in peak load (1.4–1.9%) and fracture displacement (1.3–5.8%). As depicted in Figure 14(b), the flexural strength decreases with increasing stacking thickness. To statistically validate this trend, a one-way analysis of variance (ANOVA) is conducted on the experimental flexural strength across the three thicknesses (C₄L₃₆, C₆L₅₄, C₈L₇₂). The results show a significant effect of stacking thickness on flexural strength (F (2,9) = 3735, p < .001), confirming that the decrease is statistically significant.

Figure 14.

Curves of different stacking thickness and ILSS: (a) Load-displacement curves; (b) Flexural strength-stacking thickness curves.

To isolate the pure thickness effect from structural scale effects, the constant span group (C₄L₅₄, C₆L₅₄, C₈L₅₄, span fixed at 54 mm) is examined. In this group, as stacking thickness increases from 0.6 mm to 1.2 mm, the peak load increases from 106 N (C₄L₅₄) to 406 N (C₈L₅₄) at ILSS40, while the flexural strength decreases from 1835 MPa to 1757 MPa (see Table 9), consistent with the trend observed in the constant span-to-thickness ratio group. This confirms that the thickness-dominated reduction in flexural strength is not an artifact of changing span length.

Combining Figure 15(a)-(b), it can be observed that with increasing stacking thickness, compression damage to the upper layer material gradually increases, leading to a decrease in the overall flexural strength. According to the weakest-link theory of Weibull,⁵³ a thicker specimen has a higher probability of encountering the weakest link during failure, resulting in lower overall strength. This phenomenon is consistent with Qiu et al.’s⁵⁴ study on the mechanical properties of nacre in Hyriopsis cumingii. In addition, increasing ILSS can also enhance flexural strength, though the effect is relatively small within the studied range.

Figure 15.

(a) Fiber compression damage cloud diagram on the upper layer before fracture; (b) Fracture morphology in the experiment.

Upon comparing simulation results with experimental data (C₄L₃₆, C₆L₅₄, C₈L₇₂), the overall trends are consistent. Although some discrepancies exist, they do not significantly affect the representation of the material’s overall performance. As discussed in Section 4.1, these discrepancies are primarily attributed to the difference in fiber volume fraction between simulation (V_f = 60%) and experiment (V_f = 67%).

Conclusion

In this study, the three-point bending damage evolution process of unidirectional carbon fiber reinforced polymer (CFRP) laminates is characterized through finite element analysis, experiment, and a catastrophic failure model. The effects of stacking thickness and interlaminar shear strength (ILSS) on the damage evolution and flexural properties of unidirectional CFRP laminates are investigated. The key findings of this research can be summarized as follows:

(1) As the stacking thickness increases (0.6 mm–1.2 mm), the experimental flexural strength decreases from 2099 MPa to 1790 MPa (see Table 9). The compressive stress on the upper layers increases, leading to more compression damage to the fibers and thereby causing a decrease in the flexural strength of the unidirectional CFRP laminates. Simulations with constant span length confirm that this thickness-induced trend is not caused by concurrent changes in span length.

(2) Simulation results indicate that for laminates with high ILSS (20, 40, 50, and 80 MPa), the predominant failure mode is fiber fracture. In contrast, laminates with low ILSS (10 MPa) primarily exhibit delamination failure and their neutral plane is most susceptible to damage (see Figure 11).

(3) At ILSS20 MPa, the catastrophic failure model reveals that the damage and failure of the interlaminar interface of unidirectional CFRP laminates obey power-law characteristics with a power index of 0.5 (see equation (11) and Table 8), indicative of brittle fracture behavior. As the displacement approaches the fracture point, damage becomes complete, and the damage rate increases dramatically without bound. The damage coefficient λ increases from 1.88 for C₄L₃₆ to 4.70 for C₈L₇₂ (see Table 8), indicating that under this specific interlaminar strength condition, thicker stacking thickness leads to faster damage progression, with initial damage points closer to the fracture point. This finding is drawn from ILSS20 MPa, where the slower damage rate enables numerical tracking; at higher ILSS levels (e.g., 80 MPa), the damage evolves too rapidly for such quantitative comparison.

The insights gained from this paper provide useful guidance for optimizing the design of carbon fiber structural components, enhancing their safety, and ensuring stability in practical engineering applications. While this study focuses on unidirectional CFRP laminates, further investigation for multidirectional or quasi-isotropic laminates is already underway. In addition, it should be noted that the conclusions regarding damage evolution mechanisms are predominantly simulation-based and should be interpreted as such. Direct experimental measurement (e.g., in situ observation of damage progression, short-beam shear testing for ILSS) is therefore recommended for future work.

Footnotes

ORCID iDs

Chengyu Guan

Lihong Liang

Author contributions

Chengyu Guan: Conceptualization, Methodology, Software, Writing - Review & Editing, Funding acquisition.

Huanyu Li: Conceptualization, Methodology, Software, Investigation, Writing - Original Draft.

Shuai Wang: Writing - Review & Editing, Funding acquisition.

Lihong Liang: Supervision, Writing - Review & Editing, Funding acquisition.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the National Natural Science Foundation of China through Grants #91860102, #12402149, #12172035, #92160203, the Fundamental Research Funds for the Central Universities of China (buctrc201930, BH202422), and Postdoctoral Fellowship Program of CPSF (GZC20230204).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article

Data Availability Statement

Data will be made available on request.*

References

Bai

Hao

Dai

, et al. A review on integrated design and manufacturing technologies for continuous fiber-reinforced polymer composite lightweight structures in aerospace. Thin-Walled Struct 2025; 217: 113750. https://doi.org/10.1016/j.tws.2025.113750

Lin

Qian

Lei

, et al. Structural integrity issues of composite materials and structures in future transportation equipment. Compos Struct 2025; 358: 118943. https://doi.org/10.1016/j.compstruct.2025.118943

Cao

, et al. Formation mechanism and detection and evaluation methods as well as repair technology of crack damage in fiber-reinforced composite wind turbine blade: a review. Int J Adv Manuf Technol 2022; 120(9): 5649–5672. https://doi.org/10.1007/s00170-022-09230-z

Wei

, et al. Characterization of progressive damage behaviour and failure mechanism of carbon fibre reinforced DP590 laminates. Thin-Walled Struct 2021; 168: 108293. https://doi.org/10.1016/j.tws.2021.108293

Zhang

Zou

Wei

, et al. Modifying stacking sequences to leverage the effects of shear thickening gel (STG) on the impact resistance of the STG applied carbon fibre-reinforced polymer (SACFRP) composite laminates. Compos Sci Technol 2026; 273: 111414. https://doi.org/10.1016/j.compscitech.2025.111414

Lei

Pan

Sun

, et al. Fatigue damage mechanisms and evolution of residual tensile strength in CFRP composites: stacking sequence effect. Compos Struct 2024; 330: 117818. https://doi.org/10.1016/j.compstruct.2023.117818

Chen

Wang

Lin

, et al. Theoretical modeling and experimental verification of the bending deformation of fiber metal laminates. Materials 2023; 16: 3486. https://doi.org/10.3390/ma16093486

Balasubramani

Nagarajan

. Experimental and optimization of interlaminar shear strength and flexural modulus of laminated composite materials under flexural loading. Mater Today Proc 2023. https://doi.org/10.1016/j.matpr.2023.09.224

Carrillo

Cantwell

. Mechanical properties of a novel fiber–metal laminate based on a polypropylene composite. Mech Mater 2009; 41: 828–838. https://doi.org/10.1016/j.mechmat.2009.03.002

10.

Husain

Khan

Mourad

AHI

. Stacking sequence effects on flexural stiffness, failure progression, and energy absorption in asymmetric CFRP laminates. Results Eng 2025; 27: 105863. https://doi.org/10.1016/j.rineng.2025.105863

11.

Gai

Yao

, et al. Mechanisms of interlaminar strength enhancement in CFRP bending beams via double-double laminate optimization: effects of ply sequence and orientation. Compos B Eng 2026; 314: 113439. https://doi.org/10.1016/j.compositesb.2026.113439

12.

Chen

Sun

Xiao

, et al. Investigation on flexural properties of intra-layer hybrid composite laminates reinforced with carbon and glass fibers. Fibers Polym 2023; 24: 1119–1130. https://doi.org/10.1007/s12221-023-00083-1

13.

Mohanavel

Suresh Kumar

Vairamuthu

, et al. Influence of stacking sequence and fiber content on the mechanical properties of natural and synthetic fibers reinforced penta-layered hybrid composites. J Nat Fibers 2022; 19: 5258–5270. https://doi.org/10.1080/15440478.2021.1875368

14.

Sun

Tong

Chen

, et al. Mechanical properties of hybrid composites reinforced by carbon and basalt fibers. Int J Mech Sci 2018; 148: 636–651. https://doi.org/10.1016/j.ijmecsci.2018.08.007

15.

Dhanaraj

Venkateshwaran

Vigneshwaran

. Unveiling synergistic and mechanical performance of hybrid fiber metal laminates by optimized aluminum/basalt/pineapple fiber stacking sequence. J Reinforc Plast Compos 2026; 45(3-4): 513–526. https://doi.org/10.1177/07316844251332222

16.

Saka Dinç

Eren

Öz

, et al. Influence of surface modifications on interfacial and mechanical properties in fibre-reinforced composites. J Reinforc Plast Compos 2025; 44(23-24): 2555–2570. https://doi.org/10.1177/07316844241259870

17.

Saurabh

Das Chakladar

Deb

. Influence of weave patterns on tensile, flexural, and short beam shear performance of FRP composites. J Reinforc Plast Compos 2025: 07316844251377629.

18.

Chen

. Residual flexural properties of surface scratched CFRP laminates with diverse span lengths and thickness profiles. J Mater Res Technol 2025; 36: 4453–4463. https://doi.org/10.1016/j.jmrt.2025.04.124

19.

Wang

, et al. Cryogenic damage mechanism of CFRP laminates under bending load via in-situ fiber-optic acoustic emission and mode decomposition. Compos B Eng 2026; 308: 112994. https://doi.org/10.1016/j.compositesb.2025.112994

20.

Ning

, et al. Toughening effect of CB-epoxy interleaf on the interlaminar mechanical properties of CFRP laminates. Compos Appl Sci Manuf 2015; 68: 226–234. https://doi.org/10.1016/j.compositesa.2014.09.030

21.

Zhang

Huan

, et al. Study on toughening modification of epoxy resin and its properties of T700 carbon fibre composites. J Reinforc Plast Compos 2026; 45(5-6): 1503–1517. https://doi.org/10.1177/07316844251331995

22.

Song

Ding

Fan

, et al. Effect of carbon nanotube on the mechanical properties of short carbon fiber reinforced polyether-ether-ketone. J Reinforc Plast Compos 2024: 07316844241295871.

23.

Yokozeki

Ueda

, et al. Simulation on the mechanical performance and fracture behavior of unidirectional carbon fiber-reinforced composites. J Compos Mater 2021; 55: 3639–3649. https://doi.org/10.1177/00219983211020118

24.

Singh

Yadav

Dodla

, et al. The influence of nanoparticle surface coatings on fiber-reinforced polymer composites: a review of mechanical, tribological, and thermal properties. Compos Interfaces 2026: 1–29. https://doi.org/10.1080/09276440.2026.2613181

25.

Mahato

Lomov

Shiverskii

, et al. A review of electrospun nanofiber interleaves for interlaminar toughening of composite laminates. Polymers 2023; 15(6): 1380. https://doi.org/10.3390/polym15061380

26.

, et al. Effect of adhesive quantity on failure behavior and mechanical properties of fiber metal laminates based on the aluminum–lithium alloy. Compos Struct 2016; 152: 687–692. https://doi.org/10.1016/j.compstruct.2016.05.098

27.

Yan

Zhu

Tong

, et al. Experimental study on the mechanical properties of laminates made of thin carbon fiber plies. Compos Struct 2020; 245: 112336. https://doi.org/10.1016/j.compstruct.2020.112336

28.

Yang

Sugahara

, et al. A study on the failure behavior and mechanical properties of unidirectional fiber reinforced thermosetting and thermoplastic composites. Compos B Eng 2016; 99: 162–172. https://doi.org/10.1016/j.compositesb.2016.06.005

29.

Tan

Naya

Yang

, et al. The role of interfacial properties on the intralaminar and interlaminar damage behaviour of unidirectional composite laminates: experimental characterization and multiscale modelling. Compos B Eng 2018; 138: 206–221. https://doi.org/10.1016/j.compositesb.2017.11.043

30.

Bhat

Mohan

Sharma

, et al. Mechanical testing and microstructure characterization of glass fiber reinforced isophthalic polyester composites. J Mater Res Technol 2019; 8: 3653–3661. https://doi.org/10.1016/j.jmrt.2019.06.003

31.

Laux

Gan

Dulieu-Barton

, et al. Ply thickness and fibre orientation effects in multidirectional composite laminates subjected to combined tension/compression and shear. Compos Appl Sci Manuf 2020; 133: 105864. https://doi.org/10.1016/j.compositesa.2020.105864

32.

Liu

Jiang

, et al. Effect of ply thickness on tensile and bending performances of carbon fiber reinforced thermoplastic unidirectional laminate. Polym Compos 2023; 44: 1889–1901. https://doi.org/10.1002/pc.27212

33.

Dong

Ranaweera-Jayawardena

Davies

. Flexural properties of hybrid composites reinforced by S-2 glass and T700S carbon fibres. Compos B Eng 2012; 43: 573–581. https://doi.org/10.1016/j.compositesb.2011.09.001

34.

Ullah

Harland

Silberschmidt

. Damage modelling in woven-fabric CFRP laminates under large-deflection bending. Comput Mater Sci 2012; 64: 130–135. https://doi.org/10.1016/j.commatsci.2012.05.036

35.

Ridha

Wang

Chen

, et al. Modelling complex progressive failure in notched composite laminates with varying sizes and stacking sequences. Compos Appl Sci Manuf 2014; 58: 16–23. https://doi.org/10.1016/j.compositesa.2013.11.012

36.

ASTM D7264/D7264M-07 . Standard test method for flexural properties of polymer matrix composite materials, 2020.

37.

Zhu

Sun

, et al. Modeling for CFRP structures subjected to quasi-static crushing. Compos Struct 2018; 184: 41–55. https://doi.org/10.1016/j.compstruct.2017.09.001

38.

Hashin

. Failure criteria for unidirectional fiber composites. J Appl Mech 1980; 47: 329–334. https://doi.org/10.1115/1.3153664

39.

Koloor

SSR

Abdullah

Tamin

, et al. Fatigue damage of cohesive interfaces in fiber-reinforced polymer composite laminates. Compos Sci Technol 2019; 183: 107779. https://doi.org/10.1016/j.compscitech.2019.107779

40.

Zou

Gong

Tian

, et al. Effect of interface angle on mode I delamination damage behavior of multidirectional fully isotropic laminates with the same global stiffness. Thin-Walled Struct 2023; 182: 110211. https://doi.org/10.1016/j.tws.2022.110211

41.

Benzeggagh

Kenane

. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos Sci Technol 1996; 56: 439–449. https://doi.org/10.1016/0266-3538(96)00005-x

42.

Xue

Kirane

. Effect of the shape of the softening damage law on the predicted tensile fracturing and energy dissipation in textile composites. Int J Damage Mech 2021; 30(5): 786–822. https://doi.org/10.1177/1056789520986849

43.

Pinho

Iannucci

Robinson

. Formulation and implementation of decohesion elements in an explicit finite element code. Compos Appl Sci Manuf 2006; 37(5): 778–789. https://doi.org/10.1016/j.compositesa.2005.06.007

44.

Ridha

Chen

, et al. On cohesive element parameters and delamination modelling. Eng Fract Mech 2019; 206: 278–296. https://doi.org/10.1016/j.engfracmech.2018.12.009

45.

Zhao

Gong

Zhang

, et al. Simulation of delamination growth in multidirectional laminates under mode I and mixed mode I/II loadings using cohesive elements. Compos Struct 2014; 116: 509–522. https://doi.org/10.1016/j.compstruct.2014.05.042

46.

Ren

Zhang

. Damage visualization and enhanced interlaminar shear strength of carbon fiber reinforced composites with Ni/PUF/PU@ER microcapsules. Compos Appl Sci Manuf 2026; 200: 109338. https://doi.org/10.1016/j.compositesa.2025.109338

47.

Lee

Seong

, et al. Effect of phenoxy-based coating resin for reinforcing pitch carbon fibers on the interlaminar shear strength of PA6 composites. Compos Appl Sci Manuf 2016; 87: 212–219. https://doi.org/10.1016/j.compositesa.2016.04.028

48.

Perogamvros

Lampeas

. Experimental and numerical investigation of AS4/8552 interlaminar shear strength under impact loading conditions. J Compos Mater 2016; 50(19): 2669–2685. https://doi.org/10.1177/0021998315610899

49.

Barbero

. Introduction to composite materials design. CRC Press, 2010.

50.

Prusty

Rathore

Singh

, et al. Experimental optimization of flexural behaviour through inter-ply fibre hybridization in FRP composite. Constr Build Mater 2016; 118: 327–336. https://doi.org/10.1016/j.conbuildmat.2016.05.054

51.

Liang

Liu

, et al. Power-law characteristics of damage and failure of ceramic coating systems under three-point bending. Surf Coating Technol 2016; 285: 113–119. https://doi.org/10.1016/j.surfcoat.2015.11.044

52.

Wang

Liang

. Size effects on damage evolution of ceramic coatings under thermal loading. Ceram Int 2024; 50(8): 12975–12984. https://doi.org/10.1016/j.ceramint.2024.01.207

53.

Weibull

. A statistical theory of the strength of material. Ingeniors Vetenskapa Acadamiens Handligar 1939; 151: 1–45.

54.

Qiu

Liang

Hao

. Thickness effect on the mechanical properties of nacre in Hyriopsis cumingii under three-point bending. Eng Fract Mech 2022; 276: 108869. https://doi.org/10.1016/j.engfracmech.2022.108869