# Progressive Failure Analysis of Laminated CFRP Composites under Three-Point Bending Load

Authors: Bing Shen; HongJun Liu; Shengli Lv; Zheng Li; Wen Cheng
Journal: Advances in Materials Science and Engineering
 (2022)

---

## Abstract

In this paper, the mechanical behavior and damage evolution of CFRP laminates under three-point bending loads were investigated. Based on the three-dimensional Hashin failure criterion and cohesive zone method, a progressive damage analysis model was constructed to predict seven different failure modes, and an improved stiffness matrix degradation law was proposed. Numerical simulations were performed using ABAQUS finite element software and UMAT user subroutines and verified by three-point bending tests. The experimental results show that the relative errors of the flexural strength and flexural modulus of the current model simulation results and the experimental data were −3.10% and −2.70%, respectively, and the damage evolution process after the ultimate load was accurately predicted. Meanwhile, the predicted final failure and damage modes were consistent with the micrographs of the failed specimens.


---

## Body

## 1. Introduction

Fibre-reinforced polymer (FRP) composite materials are widely used in aerospace, transportation, vehicle engineering, infrastructure, and construction due to their excellent specific stiffness and strength [1–4]. Generally, FRP composites are anisotropic, and the strength and stiffness of a single layer in the longitudinal direction are much higher than those in the transverse direction. Therefore, multiple layers are usually stacked in an angled orientation to manufacture a laminated structure to achieve high strength and high stiffness in different directions [5]. The anisotropy of FRP and its multiple failure modes have made the design of composite structures extremely challenging. In the design process of composite structures, designers must clearly understand the response of FRP under various loads, usually using experimental measurements. However, there is a growing realization that the “build and test” approach is too slow and too expensive, and designers of composite structures need robust and accurate tools or predictive methods. Therefore, the use of reliable failure analysis methods is essential to understand the damage mechanism of FRP materials and guide the design of composite materials [6, 7].The failure of composite laminates is a gradual deterioration process. Local failure will cause load redistribution, and different failure modes will continue to accumulate, which degrades the stiffness of the composite material, reduces the load-bearing capacity, and eventually causes the overall failure. For theoretical and numerical studies of composite laminates, progressive damage analysis (PDA) is often used to predict damage initiation and damage evolution in composite structures [8–11]. The progressive damage analysis model consists of three parts: the material constitutive model, failure criterion, and material degradation model. The failure criteria commonly used by researchers are the Puck criterion [12], Hashin criterion [13, 14], Chang-Chang criterion [15], and Hou criterion [16, 17]. A worldwide failure exercise (WWFE) was organized by Kaddour et al. [18, 20] to investigate the ability of various failure criteria to predict the damage behavior of composites. The study shows that the Hashin criterion requires fewer input parameters to perform the analysis than other failure criteria. Although the physically based Hashin criterion is conservative for strength prediction results in some cases, it provides a solid foundation for a phenomenological understanding of failure. Hashin criterion is widely used in engineering analysis because of its accurate identification of different failure modes. Gu and Chen [21] examined and evaluated the assumptions proposed by Hashin criteria, and modified the criteria of tensile fibre and tensile matrix mode based on physical considerations. Furthermore, a stress exposure factor was introduced to describe the interaction between fibre and matrix failure. Falkowicz et al. [22] used the Hashin criterion to analyze the stability, bearing capacity, and failure of thin-walled plates subjected to axial compression. The results of the numerical analysis are highly compatible with the experimental data. Linde and de Boer [23] proposed a damage model that combined the two-dimensional Hashin criterion and the cohesive zone method by using an exponential degradation law to define the material performance degradation after failure. The model can accurately predict the inter rivet buckling behavior of the composites.With the improvement of computer performance, the computational efficiency of finite element analysis software has been effectively improved, and it has become an important tool for progressive failure analysis. Developers set up various user subroutines in the ABAQUS finite element analysis software, which can be used to define complex constitutive models and failure criteria. Huang [24] studied the general constitutive relationship that described the inelastic and strength behavior of fibre-reinforced laminated composites. The user subroutine UGENS in the ABAQUS software was used to program the constitutive relations and perform the failure analysis. Massarwa [25] developed a multiscale model using ABAQUS/Implicit and user subroutine UMAT to study the mechanical behavior of glass/carbon fibre hybrid composites under different loading conditions and proved that the proposed model could capture the experimental strain distributions of the composite structure under in-plane shear loading. Hu et al. [26] used the user subroutines USDFLD and UMAT in the ABAQUS program to implement a progressive failure analysis based on elastic damage mechanics, and studied the influence of wrinkles under bending conditions on the performance of L-type laminates.The purpose of this study is to understand the mechanical response and failure behavior of asymmetric unidirectional carbon fibre-reinforced composite laminates under three-point bending loading. In this paper, the two-dimensional Hashin criterion and Hou criterion are combined into a new failure criterion. With the cohesive zone method, a progressive damage model is established that includes six intralaminar failure modes and one interlaminar failure mode. The model adopts a new stiffness degradation matrix to simulate the damage evolution process after failure. The results of the numerical analysis are verified by a three-point bending test, and the obtained load-displacement response is consistent with the test, especially in the damage evolution stage after failure has occurred. The proposed model accurately predicts the failure of laminates, which is beneficial to the study of the damage mechanism of composite laminates.


## 2. Methodology

### 2.1. Intralaminar Damage Initiation

Hashin proposed a method to simulate various fibre and matrix damage modes in a single layer [13, 14]. This method has been widely used and developed by many researchers [27, 28]. Santiuste et al. [29] compared the Hou and Hashin criteria under dynamic conditions and analyzed the failure of a beam subjected to low-velocity impact in a three-point configuration to predict the failure modes (fibre tensile failure, fibre compressive failure, matrix cracking, matrix crushing, and delamination). According to the three-dimensional Hashin criterion and modified Hou’s criterion, Yang et al. [30] proposed a new progressive damage model based on the Hashin-type failure criterion and cohesive zone method, by which the degradation stage of the load-deflection curve after the ultimate bending load was well predicted.In this study, the two-dimensional Hashin criterion was used to determine the in-plane damage, and the Hou criterion was used for out-of-plane damage. The six failure modes are fibre tensile (FT) failure, fibre compressive (FC) failure, matrix tensile failure (MT), matrix compressive failure (MC), through-thickness tensile failure (TT), and through-thickness compressive failure (TC). The failure criteria are shown in equations (1)–(6).


#### 2.1.1. Fibre Tensile Failure (FT)

(1)eFT2=σ11XT2+τ12S122+τ13S132≥1,σ11≥0.


#### 2.1.2. Fibre Compressive Failure (FC)

(2)eFC2=σ11XC2≥1,σ11≤0.


#### 2.1.3. Matrix Tensile Failure (MT)

(3)eMT2=σ22YT2+τ12S122+τ23S232≥1,σ22≥0.


#### 2.1.4. Matrix Compressive Failure (MC)

(4)eMC2=σ22YCYC2S232−1+τ12S122+σ222S232≥1,σ22≤0.


#### 2.1.5. Through-Thickness Tensile Failure (TT)

(5)eTT2=σ33ZT2+τ13S132+τ23S232≥1,σ33≥0.


#### 2.1.6. Through-Thickness Compressive Failure (TC)

(6)eTC2=σ33ZC2+τ13S132+τ23S232≥1,σ33≥0,where σii (i = 1,2,3) is the normal stress; τij (i,j = 1,2,3; i≠j) is the shear stress; XP, YP, and ZP are the longitudinal strength, transverse strength, and through-thickness strength, respectively; while P = T and P=C are the tension and compression, respectively; Sij is the shear strength. In any situation, when criterion e > 1 is satisfied, the element fails as the corresponding mode.


### 2.2. Intralaminar Damage Evolution

When the stress state in the composite material satisfies the failure criterion, the damage is considered to have occurred. At this time, a set of internal state variables are used to describe the discrete microstructure damage in the material. When the damage continues to increase, the stiffness coefficient of the composite material also degrades. To simulate nonlinear behavior, El-Sisi et al. [31] established a simple progressive damage model by introducing a small degradation factor and found that the degradation law (linear, quadratic, and exponential) hardly affected the load-displacement response. The exponential progressive damage evolution to describe the degradation process of the stiffness coefficient [33, 35] is shown in equation (7):(7)dk=1−1ekexpσkεk1−ekLcGk,where dkk=FT,FC,MT,MC,TT,TC represent the corresponding damage variables of different damage modes, LC is the characteristic length of the elements, and Gkis the fracture energy of the fibre or matrix. The damage-independent definition is shown in (8):(8)df=dft,σ1>0dfc,σ1<0,dm=dmt,σ2>0dmc,σ2<0,dt=dtt,σ3>0dtc,σ3<0.Therefore, the effective stiffness matrix associated with stiffness degradation after the damage initiation point can be expressed by (9) and (10):(9)Cd=C11dC12dC13d000C22dC23d000C33d000C44d00symmetricC55d0C66d.(10)C11d=1−df2C11,C22d=1−dm2C22,C33d=1−dt2C33C12d=1−df2C12,C13d=1−dm2C13,C23d=1−dt2C23C44d=1−df1−dmC44,C55d=1−df1−dtC55,C66d=1−dm1−dtC66.According to the calculation method of the stiffness degradation matrix in (10), when the element is damaged, damage coefficient d is infinitely close to 1, and the degradation factor before the stiffness coefficient in the corresponding direction instantly decays to zero. However, in the actual mechanical test, although the mechanical properties of the element with compression failure sharply drop, the fracture may be closed under continuous load, and the damaged element can maintain a low-level bearing capacity in other directions. Therefore, in (10), the calculation method of some degradation coefficients obviously has certain defects. Hu et al. [26] believed that different modes of failure had different degrees of reduction of material elastic parameters, and constant material degradation factors were used. In this paper, exponential degradation and constant parameter degradation are combined. By considering the effects of different failure modes on the composite stiffness matrix, a semiempirical form of stiffness matrix degradation was proposed, as shown in (11).(11)C11d=1−df2C11,C22d=1−0.95dm2C22,C33d,=1−0.95dt2C33C12d=1−df1−0.95dmC12,C13d=1−df1−0.95dtC13,C23d=1−0.95dm1−0.95dtC23,C44d=1−0.95df1−0.95dm1−0.1dtC44,C55d=1−0.95df1−0.1dm1−0.95dtC55,C66d=1−0.1df1−0.95dm1−0.95dtC66,where coefficients 0.1 and 0.95 of the damage variable were assumed by analyzing the influence of different failure modes on the elastic parameters.


### 2.3. Cohesive Zone Method

The cohesive zone method (CZM) originated from the research of Dugdale [33] and Barenblatt [34] in the 1960s. It has been widely used in finite element analysis. It combines damage mechanics and fracture mechanics to simulate crack growth, debonding, and layering. It is an important tool to layer composite materials [35, 36].In this study, a bilinear CZM with a quadratic damage initiation criterion is employed. The quadratic nominal stress criterion for mixed-mode loading can be expressed as(12)tntn0+tsts0+tttt0≥1,where tn, ts, and tt are the normal stress and shear stress components, respectively, and tn0, ts0, and tt0 are the interlaminar strengths. Delamination initiates when this equation is equal to or greater than 1. is the Macaulay bracket which indicates that pure compressive stress does not cause delamination.Delamination propagation of composite structures under mixed-mode loading is usually determined by the energy release rate and fracture toughness. Two different propagation criteria can be selected in ABAQUS, the power law proposed by Camanho et al. [37], and the B-K criterion proposed by Benzeggagh and Kenane [38], as shown in (13) and (14):(13)GIGICα+GIIGIICα+GIIIGIIICα=1,(14)GC=GIC+GIIC−GICGII+GIIIGI+GII+GIIIβ,where GIC and GIIC are the fracture toughness values in Mode-I and Mode-II, respectively. GI, GII, and GIII are the strain energy release rates for each related fracture mode.


## 2.1. Intralaminar Damage Initiation

Hashin proposed a method to simulate various fibre and matrix damage modes in a single layer [13, 14]. This method has been widely used and developed by many researchers [27, 28]. Santiuste et al. [29] compared the Hou and Hashin criteria under dynamic conditions and analyzed the failure of a beam subjected to low-velocity impact in a three-point configuration to predict the failure modes (fibre tensile failure, fibre compressive failure, matrix cracking, matrix crushing, and delamination). According to the three-dimensional Hashin criterion and modified Hou’s criterion, Yang et al. [30] proposed a new progressive damage model based on the Hashin-type failure criterion and cohesive zone method, by which the degradation stage of the load-deflection curve after the ultimate bending load was well predicted.In this study, the two-dimensional Hashin criterion was used to determine the in-plane damage, and the Hou criterion was used for out-of-plane damage. The six failure modes are fibre tensile (FT) failure, fibre compressive (FC) failure, matrix tensile failure (MT), matrix compressive failure (MC), through-thickness tensile failure (TT), and through-thickness compressive failure (TC). The failure criteria are shown in equations (1)–(6).


### 2.1.1. Fibre Tensile Failure (FT)

(1)eFT2=σ11XT2+τ12S122+τ13S132≥1,σ11≥0.


### 2.1.2. Fibre Compressive Failure (FC)

(2)eFC2=σ11XC2≥1,σ11≤0.


### 2.1.3. Matrix Tensile Failure (MT)

(3)eMT2=σ22YT2+τ12S122+τ23S232≥1,σ22≥0.


### 2.1.4. Matrix Compressive Failure (MC)

(4)eMC2=σ22YCYC2S232−1+τ12S122+σ222S232≥1,σ22≤0.


### 2.1.5. Through-Thickness Tensile Failure (TT)

(5)eTT2=σ33ZT2+τ13S132+τ23S232≥1,σ33≥0.


### 2.1.6. Through-Thickness Compressive Failure (TC)

(6)eTC2=σ33ZC2+τ13S132+τ23S232≥1,σ33≥0,where σii (i = 1,2,3) is the normal stress; τij (i,j = 1,2,3; i≠j) is the shear stress; XP, YP, and ZP are the longitudinal strength, transverse strength, and through-thickness strength, respectively; while P = T and P=C are the tension and compression, respectively; Sij is the shear strength. In any situation, when criterion e > 1 is satisfied, the element fails as the corresponding mode.


## 2.1.1. Fibre Tensile Failure (FT)

(1)eFT2=σ11XT2+τ12S122+τ13S132≥1,σ11≥0.


## 2.1.2. Fibre Compressive Failure (FC)

(2)eFC2=σ11XC2≥1,σ11≤0.


## 2.1.3. Matrix Tensile Failure (MT)

(3)eMT2=σ22YT2+τ12S122+τ23S232≥1,σ22≥0.


## 2.1.4. Matrix Compressive Failure (MC)

(4)eMC2=σ22YCYC2S232−1+τ12S122+σ222S232≥1,σ22≤0.


## 2.1.5. Through-Thickness Tensile Failure (TT)

(5)eTT2=σ33ZT2+τ13S132+τ23S232≥1,σ33≥0.


## 2.1.6. Through-Thickness Compressive Failure (TC)

(6)eTC2=σ33ZC2+τ13S132+τ23S232≥1,σ33≥0,where σii (i = 1,2,3) is the normal stress; τij (i,j = 1,2,3; i≠j) is the shear stress; XP, YP, and ZP are the longitudinal strength, transverse strength, and through-thickness strength, respectively; while P = T and P=C are the tension and compression, respectively; Sij is the shear strength. In any situation, when criterion e > 1 is satisfied, the element fails as the corresponding mode.


## 2.2. Intralaminar Damage Evolution

When the stress state in the composite material satisfies the failure criterion, the damage is considered to have occurred. At this time, a set of internal state variables are used to describe the discrete microstructure damage in the material. When the damage continues to increase, the stiffness coefficient of the composite material also degrades. To simulate nonlinear behavior, El-Sisi et al. [31] established a simple progressive damage model by introducing a small degradation factor and found that the degradation law (linear, quadratic, and exponential) hardly affected the load-displacement response. The exponential progressive damage evolution to describe the degradation process of the stiffness coefficient [33, 35] is shown in equation (7):(7)dk=1−1ekexpσkεk1−ekLcGk,where dkk=FT,FC,MT,MC,TT,TC represent the corresponding damage variables of different damage modes, LC is the characteristic length of the elements, and Gkis the fracture energy of the fibre or matrix. The damage-independent definition is shown in (8):(8)df=dft,σ1>0dfc,σ1<0,dm=dmt,σ2>0dmc,σ2<0,dt=dtt,σ3>0dtc,σ3<0.Therefore, the effective stiffness matrix associated with stiffness degradation after the damage initiation point can be expressed by (9) and (10):(9)Cd=C11dC12dC13d000C22dC23d000C33d000C44d00symmetricC55d0C66d.(10)C11d=1−df2C11,C22d=1−dm2C22,C33d=1−dt2C33C12d=1−df2C12,C13d=1−dm2C13,C23d=1−dt2C23C44d=1−df1−dmC44,C55d=1−df1−dtC55,C66d=1−dm1−dtC66.According to the calculation method of the stiffness degradation matrix in (10), when the element is damaged, damage coefficient d is infinitely close to 1, and the degradation factor before the stiffness coefficient in the corresponding direction instantly decays to zero. However, in the actual mechanical test, although the mechanical properties of the element with compression failure sharply drop, the fracture may be closed under continuous load, and the damaged element can maintain a low-level bearing capacity in other directions. Therefore, in (10), the calculation method of some degradation coefficients obviously has certain defects. Hu et al. [26] believed that different modes of failure had different degrees of reduction of material elastic parameters, and constant material degradation factors were used. In this paper, exponential degradation and constant parameter degradation are combined. By considering the effects of different failure modes on the composite stiffness matrix, a semiempirical form of stiffness matrix degradation was proposed, as shown in (11).(11)C11d=1−df2C11,C22d=1−0.95dm2C22,C33d,=1−0.95dt2C33C12d=1−df1−0.95dmC12,C13d=1−df1−0.95dtC13,C23d=1−0.95dm1−0.95dtC23,C44d=1−0.95df1−0.95dm1−0.1dtC44,C55d=1−0.95df1−0.1dm1−0.95dtC55,C66d=1−0.1df1−0.95dm1−0.95dtC66,where coefficients 0.1 and 0.95 of the damage variable were assumed by analyzing the influence of different failure modes on the elastic parameters.


## 2.3. Cohesive Zone Method

The cohesive zone method (CZM) originated from the research of Dugdale [33] and Barenblatt [34] in the 1960s. It has been widely used in finite element analysis. It combines damage mechanics and fracture mechanics to simulate crack growth, debonding, and layering. It is an important tool to layer composite materials [35, 36].In this study, a bilinear CZM with a quadratic damage initiation criterion is employed. The quadratic nominal stress criterion for mixed-mode loading can be expressed as(12)tntn0+tsts0+tttt0≥1,where tn, ts, and tt are the normal stress and shear stress components, respectively, and tn0, ts0, and tt0 are the interlaminar strengths. Delamination initiates when this equation is equal to or greater than 1. is the Macaulay bracket which indicates that pure compressive stress does not cause delamination.Delamination propagation of composite structures under mixed-mode loading is usually determined by the energy release rate and fracture toughness. Two different propagation criteria can be selected in ABAQUS, the power law proposed by Camanho et al. [37], and the B-K criterion proposed by Benzeggagh and Kenane [38], as shown in (13) and (14):(13)GIGICα+GIIGIICα+GIIIGIIICα=1,(14)GC=GIC+GIIC−GICGII+GIIIGI+GII+GIIIβ,where GIC and GIIC are the fracture toughness values in Mode-I and Mode-II, respectively. GI, GII, and GIII are the strain energy release rates for each related fracture mode.


## 3. Finite Element Modelling

Figure1 shows the finite element model established in the ABAQUS 2020 software. A UMAT subroutine and cohesion model were integrated into ABAQUS/Standard to simulate the mechanical behavior of composite laminates under three-point bending loads. The loading nose and supports were modeled as analytical rigid bodies, the six degrees of freedom of the supports were fixed, and the loading nose was allowed to move only in the direction of the laminate thickness. The contact among the loading nose, supports, and laminate adopted a general contact algorithm, and the friction coefficient between the contact parts was set as 0.2 [39]. The specimens were asymmetrical laminates of 20 layers. The ply angles are [0/90]10. C3D8R solid elements with enhanced hourglass control were used for each layer of the laminate, and a 0.001-mm-thick cohesive solid layer composed of COH3D8 elements was inserted between the layers to simulate interlaminar damage. During the simulation, the properties of the cohesive layer were defined by the traction-separation response. The onset of delamination damage in the mixed mode follows a quadratic nominal stress criterion [40]; once the delamination damage onset condition is satisfied, it evolves according to the BK criterion. The state of delamination is represented by an ABAQUS scalar stiffness degradation variable (SDEG), which ranges from 0 (no delamination) to 1 (delamination). The mechanical properties of the unidirectional CF/epoxy composite for numerical simulation input are shown in Table 1. Material properties are mainly provided by suppliers and estimated by homogenization methods. To ensure the calculation accuracy and efficiency, the global size of the mesh is 1 mm × 1 mm, and three elements are set in the thickness direction of every single layer. The entire model contains 133,510 elements in total. The numerical simulations were performed on a 64-bit Intel i7-10700k processor with 64 GB RAM.Figure 1
Finite element model.
(a)(b)Table 1
Input material properties for carbon fibre/epoxy composite ply.
PropertiesValuesPhysical propertiesρf = 1.78 g/cm3; ρm = 1.2 g/cm3; Vf = 50%Composite modulus (GPa)E11 = 115; E22 = E33 = 7.7; G12 = G13 = 4.8; G23 = 3.8;Poisson’s ratiosν12 = ν13 = 0.32; ν23 = 0.35;Tensile and compressive strength (MPa)XT = 1680; XC = 1250; YT = ZT = 60; YC = ZC = 165;Shear strength (MPa)S12 = S13 = 98; S23 = 55;Intralaminar fracture toughness (kJ/m2)Gf = 28.2; Gm = 1.4;Cohesive stiffness (N/mm)Kn = 1 × 104; Ks = Kt = 4 × 103;Cohesive strength (MPa)tn = ts = tt = 17.3;Interlaminar fracture toughness (J/m2)GIC = 0.325; GIIC = GIIIC = 0.68;B-K coefficientβ = 2.07The flowchart of the UMAT subroutine is schematically shown in Figure2. The damage variables di were stored as solution-dependent state variables (SDVs), which can display the damage status of the elements in the visualization module of the ABAQUS software.Figure 2
The implementation of the finite element model.


## 4. Experimental Validation

The composite laminate raw material in the experiment is unidirectional T300 carbon fibre/6509 epoxy prepreg with a fibre content of 50%, which is provided by Weihai Guangwei Composite Materials Co., Ltd. The composite material performance parameters used in the finite element analysis are provided by the manufacturer. The single-layer prepreg is 0.15 mm, the layup angle is [0°/90°]10, and the thickness of the laminate is approximately 3 mm thick. Laminates are manufactured using vacuum bag molding with a curing pressure of 0.8 bar. The curing temperature was raised from room temperature to 120°C at a heating-up rate of 1–2°C/min, maintained at 120°C for 90 minutes, and naturally cooled to below 60°C for demolding. After demolding, the manufactured composite material plate is cut into corresponding dimensions with diamond cutters according to the requirements in the test standard.A three-point bending test was performed according to the procedures in standard ASTM D7264. The obtained experimental data were used to calculate the ultimate strength and flexural modulus of the specimens according to (15).(15)σ=3PD2bh2,ε=6δhD2,E=∆σ∆ε,where σ and ε are the stress and strain at the outer surface in the load region, respectively; P is the allied force; b and h are the width and thickness, respectively; δ is the midspan deflection; D is the span length. After the strict screening of the manufactured specimens by geometric size and molding quality, five specimens were tested according to the testing configuration in Figure 3. The collected ultimate loads and corresponding deflection are listed in Table 2, which were used to calculate the flexural strength and modulus.Figure 3
Schematic diagram of the specimen size and test configuration.Table 2
Three-point bending test data of carbon fibre/epoxy composites laminate.
Ultimate loadPDeflectionδFlexural strengthσFlexural modulusEUnitNmmMPaGPaX1483.0015.182594.46358.735X2507.7105.107624.87462.646X3515.5075.222634.47062.208X4510.0125.449627.70758.981X5490.2935.209603.43859.313Average501.3055.234616.99060.377


## 5. Results and Discussion

### 5.1. Load-Deflection Response

According to the finite element model established by the theory introduced in Methodology section, the load-deflection response curve of the laminate under a three-point bending load was obtained and compared with the mechanical test results. The average ultimate load of the eight mechanical specimens is 501.305 N, the average failure deflection is 5.234 mm, and the average ultimate strength is 616.99 MPa. The ultimate loads of the simulation results of the two finite element models are 483.54 N and 485.776 N, respectively, and the ultimate strengths are 595.127 MPa and 597.878 MPa, respectively, as shown in Table3. The curves of specimens X2, X4, and X5 are strongly similar to the simulation results, as shown in Figure 4.Table 3
Comparison of mechanical properties between simulation and experiments.
Ultimate loadErrorDeflectionErrorFlexural strengthErrorFlexural modulusErrorUnitN%mm%MPa%GPa%Average501.305/5.234/616.990/60.377/FEM1483.540−3.545.202−0.61595.127−3.5458.573−2.99FEM2485.776−3.105.211−0.44597.878−3.1058.746−2.70Figure 4
Comparison of the load-deflection curves of the numerical model and experiments.From the results in Figure4, both finite element models can accurately predict the ultimate strength. Compared with the FEM2 model, the FEM1 model without stiffness degradation factor correction cannot effectively predict the damage evolution process after ultimate loading. For the data of characteristic points B ∼ E, the comparison between simulation results and average data of three typical specimens is shown in Table 4. From the data, the prediction of point B is the most accurate. The load and deflection prediction errors of the model in the damage evolution stage are -5.02% and 1.84% at point C and 9.93% and 1.3% at point D, respectively. Due to the difference between the termination judgment in the simulation program and the mechanical tests, the load and deflection errors at point E are 6.58% and -7.63%, respectively. Therefore, the model can be used to predict the load-deflection curves of composite laminates under a three-point bending load.Table 4
Comparison of characteristic points between simulation and experiments.
PointDataX2X4X5AverageStandard deviationFEM2Error (%)BLoad (N)325.292354.886376.019352.06620.805349.896−0.62Deflection (mm)5.1125.4305.2665.2690.1305.265−0.08CLoad (N)379.181394.069422.259398.50317.864378.485−5.02Deflection (mm)5.9996.4616.0766.1790.2026.2931.84DLoad (N)275.632285.186298.690286.5039.459314.9679.93Deflection (mm)6.0146.4646.2686.2480.1846.3291.30ELoad (N)305.305318.540312.393312.0805.408332.6176.58Deflection (mm)7.6087.0398.0687.5710.4216.993−7.63


## 5.1. Load-Deflection Response

According to the finite element model established by the theory introduced in Methodology section, the load-deflection response curve of the laminate under a three-point bending load was obtained and compared with the mechanical test results. The average ultimate load of the eight mechanical specimens is 501.305 N, the average failure deflection is 5.234 mm, and the average ultimate strength is 616.99 MPa. The ultimate loads of the simulation results of the two finite element models are 483.54 N and 485.776 N, respectively, and the ultimate strengths are 595.127 MPa and 597.878 MPa, respectively, as shown in Table3. The curves of specimens X2, X4, and X5 are strongly similar to the simulation results, as shown in Figure 4.Table 3
Comparison of mechanical properties between simulation and experiments.
Ultimate loadErrorDeflectionErrorFlexural strengthErrorFlexural modulusErrorUnitN%mm%MPa%GPa%Average501.305/5.234/616.990/60.377/FEM1483.540−3.545.202−0.61595.127−3.5458.573−2.99FEM2485.776−3.105.211−0.44597.878−3.1058.746−2.70Figure 4
Comparison of the load-deflection curves of the numerical model and experiments.From the results in Figure4, both finite element models can accurately predict the ultimate strength. Compared with the FEM2 model, the FEM1 model without stiffness degradation factor correction cannot effectively predict the damage evolution process after ultimate loading. For the data of characteristic points B ∼ E, the comparison between simulation results and average data of three typical specimens is shown in Table 4. From the data, the prediction of point B is the most accurate. The load and deflection prediction errors of the model in the damage evolution stage are -5.02% and 1.84% at point C and 9.93% and 1.3% at point D, respectively. Due to the difference between the termination judgment in the simulation program and the mechanical tests, the load and deflection errors at point E are 6.58% and -7.63%, respectively. Therefore, the model can be used to predict the load-deflection curves of composite laminates under a three-point bending load.Table 4
Comparison of characteristic points between simulation and experiments.
PointDataX2X4X5AverageStandard deviationFEM2Error (%)BLoad (N)325.292354.886376.019352.06620.805349.896−0.62Deflection (mm)5.1125.4305.2665.2690.1305.265−0.08CLoad (N)379.181394.069422.259398.50317.864378.485−5.02Deflection (mm)5.9996.4616.0766.1790.2026.2931.84DLoad (N)275.632285.186298.690286.5039.459314.9679.93Deflection (mm)6.0146.4646.2686.2480.1846.3291.30ELoad (N)305.305318.540312.393312.0805.408332.6176.58Deflection (mm)7.6087.0398.0687.5710.4216.993−7.63


## 6. Damage Initiation and Evolution

In the load-deflection response curve, five characteristic points can be clearly identified: points A ∼ E. Point A is the ultimate load point and initial damage point. At point A, initial damage is mainly concentrated near the loading nose, and the failure modes are mainly longitudinal fibre compressive failure on the upper surface and transverse matrix tensile failure at the bottom, as shown in Figure5. In stage AB, although the increase in strain is small, the stress dramatically changes in the damaged area. With the increase in strain, the damaged area expands to the bottom of the specimen along the thickness direction, which leads to the fibre compression failure in the 0° plies and matrix compression damage in the 90° plies, as shown in Figures 6(a) and 6(b). The damaged area forms a continuous crack on the side of the specimen, accompanied by delamination on both sides of the crack, as shown in Figure 6(c).Figure 5
Failure modes at characteristic point A. (a) FC failure and (b) MT failure.
(a)(b)Figure 6
Failure modes at characteristic point B. (a) FC failure, (b) MC failure, and (c) delamination failure.
(a)(b)(c)After point B, the load gradually increases and exhibits linear elasticity. When the load grows into stage CD, the 0° ply near the bottom surface of the specimen begins to show fibre tensile failure, as shown in Figures7(a) and 7(b). In this stage, the FC failure in the 0° plies and the MC failure in the 90° plies no longer significantly evolve, as shown in Figures 7(c) and 7(d). With the fibre tensile failure of the 0° plies, severe delamination failure appears on both sides of the damaged area and tends to expand to both ends of the specimen, as shown in Figures 7(e) and 7(f).Figure 7
Failure modes at characteristic points C and D. (a) FT failure at point C, (b) FT failure at point D, (c) FC failure at point C, (d) FC failure at point D, (e) delamination failure at point C, and (f) delamination failure at point D.
(a)(b)(c)(d)(e)(f)When approaching pointE, more than 50% of the plies in the laminate specimen are damaged, and the received load is approximately 40% lower than the ultimate load, which implies that the entire composite laminate structure has almost failed. The difference in the performance of different plies due to processing and complex coupling of different failure modes makes the appearance of point E difficult to accurately predict in numerical simulations. Compared with the mechanical tests, the E point in the simulation appears earlier and is more serious because the generated cracks on the upper surface of the specimen are closed under the action of compressive stress, so a certain degree of bearing capacity can be maintained. However, when the analysis program determines the occurrence of fibre compression failure and matrix compression failure, the mechanical properties of the failure area will be greatly reduced.The simulation and experimental results of the seven failure modes at pointE are shown in Figure 8. The side views of the five specimens subjected to the mechanical test are shown in Figure 9. Near the middle of the top surface of the specimen, the compressive stress under the loading nose causes the fibre compressive failure of ply 1 and matrix compressive failure of ply 2 which indicates the initiation of damage. With the increase in deflection, the compressive damage area expands from ply 1 down the thickness direction to ply 11 and ply 12 and forms cracks, as shown in region 1 in Figure 9(a). At the bottom of the specimen, region 2 shows a fracture that extends from ply 20 to ply 16 due to tensile stress. Delamination damage is accompanied on both sides of the compression crack and tensile fracture, as shown in regions 3, 4, 5, and 6. Similar damage is also clearly observed in the specimens in Figures 9(b)∼9(d) with the same number of damaged plies. The comparison between simulation and test results shows that the simulation model in this paper can effectively predict the damage distribution and evolution of composite laminates under three-point bending loads.Figure 8
Failure modes at characteristic point E. (a) FT failure, (b) FT failure, (c) FC failure, (d) FC failure, (e) MT failure, (f) MT failure, (g) MC failure, (h) MC failure, (i) TT failure, (j) TT failure, (k) TC failure, and (l) TC failure.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)Figure 9
Failure modes at characteristic point E. (a) Specimen X1, (b) Specimen X2, (c) Specimen X3, (d) Specimen X4, and (e) Specimen X5.
(a)(b)(c)(d)


## 7. Conclusions

In this study, a progressive failure analysis method was employed to investigate the mechanical response of unidirectional CFRP composite laminates under three-point bending loads. By using the improved 3D Hashin failure criterion and cohesive zone method, an analytical model was created to predict seven different failure modes. Considering the influence of different failure modes on the stiffness coefficient of composite materials, combining exponential and constant damage variables, the original stiffness degradation matrix was improved. Numerical simulations were performed using the ABAQUS finite element software and UMAT user subroutines and verified by three-point bending tests. The results show that the load-deflection curve of the improved model is consistent with the experimental data. Compared with the average value of the experimental data, the errors of the ultimate load and failure displacement obtained by simulation are −3.10% and −0.44%, respectively, and the errors of ultimate bending strength and bending elastic modulus are −3.10% and −2.70%, respectively. In the damage evolution stage after the ultimate load, the analytical model can accurately predict the appearance of characteristic points. Meanwhile, the predicted final failure and damage modes are consistent with the micrographs of the failed specimens. The model can be fully used to understand the damage mechanism of composite materials and provides a theoretical basis for the structural application of composite laminates. [38].



---

*Source: 3047319-2022-09-05.xml*