# A Macro-Meso Correlation Model for Numerical Simulation of CFRP Tensile Notched Strength

Authors: Fang Chen; Weixing Yao; Wen Jiang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, the tensile mechanical behavior of Carbon Fiber Reinforced Plates (CFRP) with central open-holes was studied by experimental and simulation ways. A correlation model was built by macromechanical analysis and mesodamage analysis to form a progressive analysis architecture. Composite laminates were disassembled into two kinds: Representative Volume Element (RVE), which were 3D intralayer orthotropic, and the 2D interlayer cohesive element. The macromechanical analysis built connections between external loading cases and RVE stress distribution filed. The mesodamage analysis aimed to determine multimode damage initiation and evolution inside RVE. By comparing simulation results with experimental data, the prediction accuracy on failure mode, ultimate load, and fracture morphology were good enough to show the effectiveness and rationality of this new model. In addition, this model’s applicability to different material and geometrical parameters was also verified by simulating the experiment results in the literature.


---

## Body

## 1. Introduction

The application of composite materials in aerocrafts has expanded from secondary load bearing components to primary load carrying structures. With the growing demand for composite materials in engineering, the study of its notched strength has drawn more and more attention.In recent decades, many experimental or theoretical investigations have been conducted. Auerbuch and Madhukar [1] classified and summarized eleven semiempirical models and compared the difference among linear elastic fracture mechanics, point stress criteria, and average stress criteria. Yao [2] proposed a stress field intensity model to theoretically predict the notched strength. In the late 1990s, with the development of finite element analysis technology, researchers began to use simulation methods to conduct the mechanical analysis of notched plates. For example, Tan [3] used Tsai-Wu criterion [4] to analyze the progressive damage of composite perforated plates, but this criterion had difficulty in predicting failure modes. Also, Sleight [5] chose Hasin criterion [6] and obtained good simulation results. Then, Chang et al. [7] studied the influences of in-situ strength and shear nonlinearity on notched strength numerical prediction. In the 21st century, with the improvement of computational science, researchers gradually extended their investigations from two-dimensional numerical analysis to three-dimensional ones. Hallet and Wisnom [8] introduced the damage analysis of cohesive elements into the simulation system. Falzon and Apruzzese [9] selected Puck criterion [10], a promising phenomenological model, as failure criterion. Although many kinds of numerical methods were proposed, there was little synthetic simulation architecture that could predict the notched strength of specimens with different material and geometrical parameters. In recent years, based on multiscale strategy, some valuable investigations have been conducted. For instance, Nerilli and Vairo [11] developed a nonlinear computational approach, which was proved to be a very effective method for the progressive damage analysis of composite bolted joints. Marfia and Sacco [12] proposed an efficient procedure for analyzing the mechanical responses of elastoplastic and viscoplastic composite materials.In this paper, a macro-meso correlation model was proposed to form a general simulation architecture and its applicability to different material system was verified. Four groups of NHT specimens with different stacking sequences were manufactured and tested according to ASTM-D5766 standard [13]. The numerical progressive analysis was conducted through a macro-meso correlation model. The corresponding prediction results were compared with experimental data, and the effectiveness of this model was verified. Moreover, in order to show the applicability of this model, another numerical validation example was conducted according to the parameters provided by Tian et al. [14]. Simulation results of this model showed higher prediction accuracy than Tian’s method.


## 2. Macro-Meso Correlation Model

The Macro-Meso Correlation Model is a two-phase progressive damage analysis architecture based on the Representative Volume Element (RVE).The external phase is macrolevel mechanical analysis:(I)
During the external phase, composite laminates are equivalent to a homogeneous orthotropic plate, which is combined by two types of representative volume elements, i.e., the 3D orthotropic plate element (noted as RVEA) and the 2D adhesive layer element (noted as RVEB).(II)
The purpose of the external phase is to establish the connection between external loading and the stress and strain distributions on RVEA and RVEB.The internal phase is mesolevel damage analysis:(I)
During the internal phase, RVEA and RVEB are deconstructed into nonhomogenous fiber reinforced matrix and cohesive layers, respectively.(II)
The purpose of the internal phase is to determine damage initiation and damage evolution within RVEA and RVEB.The simplified logical framework of macro-meso correlation analysis model is shown in Figure1.Figure 1
The correlation between macrolevel and mesolevel analysis.


### 2.1. Macrolevel Mechanical Analysis

The main purpose of the macrolevel analysis is to determine the continuous internal displacement field and stress distribution under external loading cases. Figure2 shows the framework of macrolevel analysis. RVEA are eight-node elements under 3D stress state, which is described by three independent normal stresses σii=1,2,3 and three shear stresses τjj=12,23,31. RVEB are four-node elements under 2D stress state, which is described by interfacial normal stress σn and two interfacial shear stresses τkk=s,t.Figure 2
Macrolevel mechanical analysis framework.The displacement filed determination is the first step for transmitting external load into RVE stress distribution. For progressive damage analysis, Newton–Rapson iteration method or explicit two-order central timing difference method can be used to determine RVE displacement distribution. After then, the stress information of RVEA and RVEB can be obtained by combining the displacement field with the physical constitutive equation. Equations (1) and (2a) show two different constitutive equations for RVEB and RVEA, respectively.(1)σnτsτt=KnKsKtδnδsδt,(2a)σ1σ2σ3τ23τ13τ12=C11C12C13C12C22C23C13C23C33C44C55C66ε11ε22ε33γ23γ13γ12,where, δii=n,s,t are interfacial relative displacements, εjj=11,22,33 are normal strains, γkk=12,13,23 are shear strains, and Kii=n,s,t are interfacial stiffness parameters. In addition, Cjj=11,22,33, Ckk=12,13,23, and Cll=44,55,66 are intralaminar stiffness parameters, which could be determined by equation (2b).(2b)C11=E11−ν23ν32ν,C12=E2ν12+ν32ν13ν,C13=E3ν13+ν12ν23ν,C22=E21−ν13ν31ν,C23=E3ν23+ν21ν13ν,C33=E31−ν12ν21ν,C44=G23,C55=G13,C66=G12,ν=1−ν12ν21−ν23ν32−ν13ν31−2ν21ν32ν31,where Eii=1,2,3 are elastic modulus parameters, Gjj=12,13,23 are shear modulus parameters. In addition, νjj=12,13,23 and νkk=21,31,32 are primary and secondary Poisson ratio, respectively. These parameters satisfy Maxwell–Betti rules (equation (2c)).(2c)ν21E2=ν12E1,ν31E3=ν13E1,ν32E3=ν23E2.


### 2.2. Mesolevel Damage Analysis

According to stress field distribution obtained in macrolevel analysis, damage behaviors within RVEA and RVEB are analyzed. Figure3 shows the framework of mesolevel analysis. The purpose of this internal level analysis is to determine the damage field and update material properties of RVE. The damage field distribution can be divided into damage initiation and damage evolution. Within RVEA, damage analysis shall be subdivided into fiber failure and inter fiber failure. For RVEB, the initiation and evolution of delamination of adhesive layer shall be discussed in detail.Figure 3
Mesolevel damage analysis framework.


#### 2.2.1. Damage Initiation

In this paper, Puck criterion was used to determine RVEA damage initiation. Puck criterion gives the expression of the damage risk factorfE according to the relationship between local stress and local strength at the material point. Damage initiation occurs when fE is greater than or equal to the critical value. Equation (3) is divided into equations (3a) and (3b) according to the location of damage initiation, respectively, corresponding to fiber failure (FF) and inter fiber failure (IFF).(3a)fE,FF=1XTσ1−ν12−ν12fmσfE1E1fσ2+σ3,σ1≥0,1−XCσ1−ν12−ν12fmσfE1E1fσ2+σ3,σ1≥0.(3b)fE,IFF=1YT−p⊥ψTR⊥ψAσθn2+σθsS122+σθtR⊥⊥A2+p⊥ψTR⊥ψAσθn,IFFσn≥0,p⊥ψCR⊥ψAσθn2+σθsS122+σθtR⊥⊥A2+p⊥ψCR⊥ψAσθn,IFFσn<0,where σθn is the normal stress on the potential IFF plane, and τθjj=s,t are the shear stresses on the IFF plane. Their determination method is shown in equation (4). Parameters R⊥ψA, R⊥⊥A, and p⊥ψTC can be acquired by using equation (5). In addition, the value of mean magnification factor mσf is recommended in Table 1.(4)σθn=σ2cos2θ+σ3cos2θ+2τ23sinθcosθ,σθs=τ31sinθ+τ12cosθ,σθt=σ3−σ2sinθcosθ+τ23cos2θ−sin2θ,(5)R⊥⊥A=YC21+p⊥⊥C,p⊥ψTCR⊥ψA=p⊥⊥TCR⊥⊥Aσθt2σθs2+σθt2+p⊥∥TCS12σθs2σθs2+σθt2.Table 1
Recommendation values of magnification and inclination factors [10, 15].
Materialmσfp⊥⊥Tp⊥⊥Cp⊥∥Tp⊥∥CCarbon-fiber1.10.25 to 0.300.25 to 0.300.350.30Glass-fiber1.30.20 to 0.250.20 to 0.250.300.25In the above formulae,p⊥⊥TC and p⊥∥TC are inclination factors, and their values can be determined according to the recommendation given by Puck in his paper (i.e., refs [10, 15]).In addition, for RVEB, the Quadratic Stress (Quads) criterion is used to determine the initiation of Delamination (DEA) damage, and its specific formula is given in equation (6):(6)σnRn2+tsRs2+ttRt2=1,where Rn is the interfacial normal stress and Rkk=s,t are the interfacial shear stresses.


#### 2.2.2. Damage Evolution

Except for FF damage, both IFF and DEA damage have a gradual propagation process in RVEA and RVEB. In this paper, the damage evolution behavior is supposed to satisfy the linear softening rule. Linear priori-damage and posterior-damage behaviors constitute the bilinear constitutive relationship curves (Figure4). The initial point A (or A′) and ending point B (or B′) of softening evolution law are two characteristic points for RVEA (or RVEB).(I)
PointAσm0,εm0 and A′t0,δ0 are determined by Puck criterion and Quads criterion, respectively. They correspond to the initiation moment of IFF and DEA damage.(II)
PointB0,εmf and B′0,δf are determined by the area of triangle OAB and OA′B′, respectively. The values of these two triangle areas are equal to the critical strain energy density gc (for RVEA) and critical strain energy release rate Gc (for RVEB).Figure 4
Bilinear constitutive relationship of RVE.When the strain energy density reaches to its critical valuegc or the strain energy release rate reaches Gc, the damage propagation is completed, which means final failure occurs within RVEA or RVEB. The damage variable under different evolution states can be synthetically represented by equation (7):(7)dIFF=0,εm<εm0,εmfεm−εm0εmεmf−εm0,εm0≤εm<εmf,1,εm≥εmf,dDEA=0,δ<δ0,δfδ−δ0δδf−δ0,δ0≤δ<δf,1,δ≥δf,where dIFF and dDEA are damage state variables for IFF and DEA damage, respectively.According to equation (7), the damage state variables are correlated to the strain and displacement field of RVEA and RVEB, respectively. These two fields are determined by corresponding stress distributions that depend on the external loading cases. If the strain and displacement of reference point (i.e., integration point) within representative volume elements are lower than εm0 or δ0, the corresponding damage variables equal to zero, which represents undamaged condition of the representative volume elements. If the strain and displacement of reference point are larger than εmf or δf, the corresponding damage variables equal to one, which represents fully damaged condition.


#### 2.2.3. Damage Characterization

The damage characterization is the data transition module from mesolevel damage analysis to macrolevel analysis (Figure5). Through the multiplication with damage state variable and initial stiffness parameters, the influences caused by damage initiation or evolution of FF, IFF, and DEA can be put on RVEA and RVEB macromechanical behaviors during macrolevel analysis. Considering variant damage modes have different effects on stiffness parameters, the specific damage characterization rules are represented as shown in equation (8):(8)Kid=1−dDEAKi,i=n,s,t,E1d=1−dFFE1,Ejd=1−dIFFEj,j=2,3,G23d=1−dIFFG23,Gkd=1−dFF1−dIFFGk,k=12,13.Figure 5
Damage characterization under different damage states.According to equation (8), if representative volume elements are in undamaged conditions, the elements stiffness parameters maintain their initial value. If the elements are in partially damaged conditions, the elements’ stiffness parameters are reduced by a linear coefficient. Specially, if the elements are in fully damaged conditions, the elements’ stiffness parameters are reduced to zero, which means the elements lose the bearing capacity of external loading.


### 2.3. Progressive Damage Analysis Based on Macro-Meso Model

By correlating macrolevel and mesolevel analyses, a progressive damage analysis flow chart can be proposed as shown in Figure6. As mentioned in above sections, stress distribution is the kernel module of macrolevel analysis. Damage initiation, damage evolution, and damage characterization constitute the mesolevel damage analysis.Figure 6
Progressive damage analysis flow chart based on macro-meso correlation model.


## 2.1. Macrolevel Mechanical Analysis

The main purpose of the macrolevel analysis is to determine the continuous internal displacement field and stress distribution under external loading cases. Figure2 shows the framework of macrolevel analysis. RVEA are eight-node elements under 3D stress state, which is described by three independent normal stresses σii=1,2,3 and three shear stresses τjj=12,23,31. RVEB are four-node elements under 2D stress state, which is described by interfacial normal stress σn and two interfacial shear stresses τkk=s,t.Figure 2
Macrolevel mechanical analysis framework.The displacement filed determination is the first step for transmitting external load into RVE stress distribution. For progressive damage analysis, Newton–Rapson iteration method or explicit two-order central timing difference method can be used to determine RVE displacement distribution. After then, the stress information of RVEA and RVEB can be obtained by combining the displacement field with the physical constitutive equation. Equations (1) and (2a) show two different constitutive equations for RVEB and RVEA, respectively.(1)σnτsτt=KnKsKtδnδsδt,(2a)σ1σ2σ3τ23τ13τ12=C11C12C13C12C22C23C13C23C33C44C55C66ε11ε22ε33γ23γ13γ12,where, δii=n,s,t are interfacial relative displacements, εjj=11,22,33 are normal strains, γkk=12,13,23 are shear strains, and Kii=n,s,t are interfacial stiffness parameters. In addition, Cjj=11,22,33, Ckk=12,13,23, and Cll=44,55,66 are intralaminar stiffness parameters, which could be determined by equation (2b).(2b)C11=E11−ν23ν32ν,C12=E2ν12+ν32ν13ν,C13=E3ν13+ν12ν23ν,C22=E21−ν13ν31ν,C23=E3ν23+ν21ν13ν,C33=E31−ν12ν21ν,C44=G23,C55=G13,C66=G12,ν=1−ν12ν21−ν23ν32−ν13ν31−2ν21ν32ν31,where Eii=1,2,3 are elastic modulus parameters, Gjj=12,13,23 are shear modulus parameters. In addition, νjj=12,13,23 and νkk=21,31,32 are primary and secondary Poisson ratio, respectively. These parameters satisfy Maxwell–Betti rules (equation (2c)).(2c)ν21E2=ν12E1,ν31E3=ν13E1,ν32E3=ν23E2.


## 2.2. Mesolevel Damage Analysis

According to stress field distribution obtained in macrolevel analysis, damage behaviors within RVEA and RVEB are analyzed. Figure3 shows the framework of mesolevel analysis. The purpose of this internal level analysis is to determine the damage field and update material properties of RVE. The damage field distribution can be divided into damage initiation and damage evolution. Within RVEA, damage analysis shall be subdivided into fiber failure and inter fiber failure. For RVEB, the initiation and evolution of delamination of adhesive layer shall be discussed in detail.Figure 3
Mesolevel damage analysis framework.


### 2.2.1. Damage Initiation

In this paper, Puck criterion was used to determine RVEA damage initiation. Puck criterion gives the expression of the damage risk factorfE according to the relationship between local stress and local strength at the material point. Damage initiation occurs when fE is greater than or equal to the critical value. Equation (3) is divided into equations (3a) and (3b) according to the location of damage initiation, respectively, corresponding to fiber failure (FF) and inter fiber failure (IFF).(3a)fE,FF=1XTσ1−ν12−ν12fmσfE1E1fσ2+σ3,σ1≥0,1−XCσ1−ν12−ν12fmσfE1E1fσ2+σ3,σ1≥0.(3b)fE,IFF=1YT−p⊥ψTR⊥ψAσθn2+σθsS122+σθtR⊥⊥A2+p⊥ψTR⊥ψAσθn,IFFσn≥0,p⊥ψCR⊥ψAσθn2+σθsS122+σθtR⊥⊥A2+p⊥ψCR⊥ψAσθn,IFFσn<0,where σθn is the normal stress on the potential IFF plane, and τθjj=s,t are the shear stresses on the IFF plane. Their determination method is shown in equation (4). Parameters R⊥ψA, R⊥⊥A, and p⊥ψTC can be acquired by using equation (5). In addition, the value of mean magnification factor mσf is recommended in Table 1.(4)σθn=σ2cos2θ+σ3cos2θ+2τ23sinθcosθ,σθs=τ31sinθ+τ12cosθ,σθt=σ3−σ2sinθcosθ+τ23cos2θ−sin2θ,(5)R⊥⊥A=YC21+p⊥⊥C,p⊥ψTCR⊥ψA=p⊥⊥TCR⊥⊥Aσθt2σθs2+σθt2+p⊥∥TCS12σθs2σθs2+σθt2.Table 1
Recommendation values of magnification and inclination factors [10, 15].
Materialmσfp⊥⊥Tp⊥⊥Cp⊥∥Tp⊥∥CCarbon-fiber1.10.25 to 0.300.25 to 0.300.350.30Glass-fiber1.30.20 to 0.250.20 to 0.250.300.25In the above formulae,p⊥⊥TC and p⊥∥TC are inclination factors, and their values can be determined according to the recommendation given by Puck in his paper (i.e., refs [10, 15]).In addition, for RVEB, the Quadratic Stress (Quads) criterion is used to determine the initiation of Delamination (DEA) damage, and its specific formula is given in equation (6):(6)σnRn2+tsRs2+ttRt2=1,where Rn is the interfacial normal stress and Rkk=s,t are the interfacial shear stresses.


### 2.2.2. Damage Evolution

Except for FF damage, both IFF and DEA damage have a gradual propagation process in RVEA and RVEB. In this paper, the damage evolution behavior is supposed to satisfy the linear softening rule. Linear priori-damage and posterior-damage behaviors constitute the bilinear constitutive relationship curves (Figure4). The initial point A (or A′) and ending point B (or B′) of softening evolution law are two characteristic points for RVEA (or RVEB).(I)
PointAσm0,εm0 and A′t0,δ0 are determined by Puck criterion and Quads criterion, respectively. They correspond to the initiation moment of IFF and DEA damage.(II)
PointB0,εmf and B′0,δf are determined by the area of triangle OAB and OA′B′, respectively. The values of these two triangle areas are equal to the critical strain energy density gc (for RVEA) and critical strain energy release rate Gc (for RVEB).Figure 4
Bilinear constitutive relationship of RVE.When the strain energy density reaches to its critical valuegc or the strain energy release rate reaches Gc, the damage propagation is completed, which means final failure occurs within RVEA or RVEB. The damage variable under different evolution states can be synthetically represented by equation (7):(7)dIFF=0,εm<εm0,εmfεm−εm0εmεmf−εm0,εm0≤εm<εmf,1,εm≥εmf,dDEA=0,δ<δ0,δfδ−δ0δδf−δ0,δ0≤δ<δf,1,δ≥δf,where dIFF and dDEA are damage state variables for IFF and DEA damage, respectively.According to equation (7), the damage state variables are correlated to the strain and displacement field of RVEA and RVEB, respectively. These two fields are determined by corresponding stress distributions that depend on the external loading cases. If the strain and displacement of reference point (i.e., integration point) within representative volume elements are lower than εm0 or δ0, the corresponding damage variables equal to zero, which represents undamaged condition of the representative volume elements. If the strain and displacement of reference point are larger than εmf or δf, the corresponding damage variables equal to one, which represents fully damaged condition.


### 2.2.3. Damage Characterization

The damage characterization is the data transition module from mesolevel damage analysis to macrolevel analysis (Figure5). Through the multiplication with damage state variable and initial stiffness parameters, the influences caused by damage initiation or evolution of FF, IFF, and DEA can be put on RVEA and RVEB macromechanical behaviors during macrolevel analysis. Considering variant damage modes have different effects on stiffness parameters, the specific damage characterization rules are represented as shown in equation (8):(8)Kid=1−dDEAKi,i=n,s,t,E1d=1−dFFE1,Ejd=1−dIFFEj,j=2,3,G23d=1−dIFFG23,Gkd=1−dFF1−dIFFGk,k=12,13.Figure 5
Damage characterization under different damage states.According to equation (8), if representative volume elements are in undamaged conditions, the elements stiffness parameters maintain their initial value. If the elements are in partially damaged conditions, the elements’ stiffness parameters are reduced by a linear coefficient. Specially, if the elements are in fully damaged conditions, the elements’ stiffness parameters are reduced to zero, which means the elements lose the bearing capacity of external loading.


## 2.2.1. Damage Initiation

In this paper, Puck criterion was used to determine RVEA damage initiation. Puck criterion gives the expression of the damage risk factorfE according to the relationship between local stress and local strength at the material point. Damage initiation occurs when fE is greater than or equal to the critical value. Equation (3) is divided into equations (3a) and (3b) according to the location of damage initiation, respectively, corresponding to fiber failure (FF) and inter fiber failure (IFF).(3a)fE,FF=1XTσ1−ν12−ν12fmσfE1E1fσ2+σ3,σ1≥0,1−XCσ1−ν12−ν12fmσfE1E1fσ2+σ3,σ1≥0.(3b)fE,IFF=1YT−p⊥ψTR⊥ψAσθn2+σθsS122+σθtR⊥⊥A2+p⊥ψTR⊥ψAσθn,IFFσn≥0,p⊥ψCR⊥ψAσθn2+σθsS122+σθtR⊥⊥A2+p⊥ψCR⊥ψAσθn,IFFσn<0,where σθn is the normal stress on the potential IFF plane, and τθjj=s,t are the shear stresses on the IFF plane. Their determination method is shown in equation (4). Parameters R⊥ψA, R⊥⊥A, and p⊥ψTC can be acquired by using equation (5). In addition, the value of mean magnification factor mσf is recommended in Table 1.(4)σθn=σ2cos2θ+σ3cos2θ+2τ23sinθcosθ,σθs=τ31sinθ+τ12cosθ,σθt=σ3−σ2sinθcosθ+τ23cos2θ−sin2θ,(5)R⊥⊥A=YC21+p⊥⊥C,p⊥ψTCR⊥ψA=p⊥⊥TCR⊥⊥Aσθt2σθs2+σθt2+p⊥∥TCS12σθs2σθs2+σθt2.Table 1
Recommendation values of magnification and inclination factors [10, 15].
Materialmσfp⊥⊥Tp⊥⊥Cp⊥∥Tp⊥∥CCarbon-fiber1.10.25 to 0.300.25 to 0.300.350.30Glass-fiber1.30.20 to 0.250.20 to 0.250.300.25In the above formulae,p⊥⊥TC and p⊥∥TC are inclination factors, and their values can be determined according to the recommendation given by Puck in his paper (i.e., refs [10, 15]).In addition, for RVEB, the Quadratic Stress (Quads) criterion is used to determine the initiation of Delamination (DEA) damage, and its specific formula is given in equation (6):(6)σnRn2+tsRs2+ttRt2=1,where Rn is the interfacial normal stress and Rkk=s,t are the interfacial shear stresses.


## 2.2.2. Damage Evolution

Except for FF damage, both IFF and DEA damage have a gradual propagation process in RVEA and RVEB. In this paper, the damage evolution behavior is supposed to satisfy the linear softening rule. Linear priori-damage and posterior-damage behaviors constitute the bilinear constitutive relationship curves (Figure4). The initial point A (or A′) and ending point B (or B′) of softening evolution law are two characteristic points for RVEA (or RVEB).(I)
PointAσm0,εm0 and A′t0,δ0 are determined by Puck criterion and Quads criterion, respectively. They correspond to the initiation moment of IFF and DEA damage.(II)
PointB0,εmf and B′0,δf are determined by the area of triangle OAB and OA′B′, respectively. The values of these two triangle areas are equal to the critical strain energy density gc (for RVEA) and critical strain energy release rate Gc (for RVEB).Figure 4
Bilinear constitutive relationship of RVE.When the strain energy density reaches to its critical valuegc or the strain energy release rate reaches Gc, the damage propagation is completed, which means final failure occurs within RVEA or RVEB. The damage variable under different evolution states can be synthetically represented by equation (7):(7)dIFF=0,εm<εm0,εmfεm−εm0εmεmf−εm0,εm0≤εm<εmf,1,εm≥εmf,dDEA=0,δ<δ0,δfδ−δ0δδf−δ0,δ0≤δ<δf,1,δ≥δf,where dIFF and dDEA are damage state variables for IFF and DEA damage, respectively.According to equation (7), the damage state variables are correlated to the strain and displacement field of RVEA and RVEB, respectively. These two fields are determined by corresponding stress distributions that depend on the external loading cases. If the strain and displacement of reference point (i.e., integration point) within representative volume elements are lower than εm0 or δ0, the corresponding damage variables equal to zero, which represents undamaged condition of the representative volume elements. If the strain and displacement of reference point are larger than εmf or δf, the corresponding damage variables equal to one, which represents fully damaged condition.


## 2.2.3. Damage Characterization

The damage characterization is the data transition module from mesolevel damage analysis to macrolevel analysis (Figure5). Through the multiplication with damage state variable and initial stiffness parameters, the influences caused by damage initiation or evolution of FF, IFF, and DEA can be put on RVEA and RVEB macromechanical behaviors during macrolevel analysis. Considering variant damage modes have different effects on stiffness parameters, the specific damage characterization rules are represented as shown in equation (8):(8)Kid=1−dDEAKi,i=n,s,t,E1d=1−dFFE1,Ejd=1−dIFFEj,j=2,3,G23d=1−dIFFG23,Gkd=1−dFF1−dIFFGk,k=12,13.Figure 5
Damage characterization under different damage states.According to equation (8), if representative volume elements are in undamaged conditions, the elements stiffness parameters maintain their initial value. If the elements are in partially damaged conditions, the elements’ stiffness parameters are reduced by a linear coefficient. Specially, if the elements are in fully damaged conditions, the elements’ stiffness parameters are reduced to zero, which means the elements lose the bearing capacity of external loading.


## 2.3. Progressive Damage Analysis Based on Macro-Meso Model

By correlating macrolevel and mesolevel analyses, a progressive damage analysis flow chart can be proposed as shown in Figure6. As mentioned in above sections, stress distribution is the kernel module of macrolevel analysis. Damage initiation, damage evolution, and damage characterization constitute the mesolevel damage analysis.Figure 6
Progressive damage analysis flow chart based on macro-meso correlation model.


## 3. Validation

### 3.1. Numerical Implementation

The Macro-Meso Correlation Model was numerically implemented as shown in Figure7. The macrolevel mechanical analysis can be conducted by common finite element method, where RVEA and RVEB are represented by 3D stress element and 2D-cohesive element, respectively. The mesolevel damage analysis, including damage initiation, evolution, and characterization is realized by operating user-defined material behavior subroutine.Figure 7
Numerical implementation of macro-meso correlation model.In order to verify the rationality of the macro-meso correlation model, we conducted tensile tests and numerical simulations on U3160/5284 specimens with four types of stacking sequences. The corresponding results are shown in Section3.2. Moreover, in order to verify the applicability of the proposed method to other issues with different material properties and geometric parameters, another simulation was also carried out in Section 3.3 to predict the experimental results provided by investigator Tian et al. [14]. All five numerical finite element models were built in ABAQUS commercial software.


### 3.2. Experiment and Simulation for U3160/5284 Plates

#### 3.2.1. Specimen and Equipment

The geometric parameters of U3160/5284 carbon fiber-reinforced laminates are shown in Figure8(a). According to different ply ratios (Table 2), four groups of specimens were manufactured and tested. In this paper, the ply ratio is defined as the quotient obtained when the quantity of plies with a certain fiber orientation angle is divided by the total quantity of whole laminated layers. For instance, specimens from group E totally own twenty plies, including six 0° plies, twelve ±45° plies, and two 90° plies. Accordingly, for specimens from group E, the 0° ply ratio is 30%, the ±45° ply ratio is 60%, and the 90° ply ratio is 10%. All four groups in Table 2 own the same number of total plies, but the 0° ply ratio is gradually increased in order to study the influences of ply ratio on specimens’ notched tensile strength. In addition, each group contains six specimens to reduce the dispersion effects of material properties on experimental results.Figure 8
Image of notched specimen and tension test machine: (a) U3160/5284 composite notched plates; (b) MTS370 hydraulic experiment system.
(a)(b)Table 2
Ply orientations of different groups’ U3160/5284 notched plates.
Group numberPly ratio [0/±45/90]Stacking sequenceLamina thickness (mm)E30/60/10[45/−45/0/45/90/0/45/0/−45/45]S0.167F40/50/10[45/−45/0/0/90/45/0/0/−45/45]SG50/40/10[45/0/0/−45/90/0/−45/0/0/45]SH60/30/10[45/0/0/90/0/0/−45/0/0/45]SThe tensile tests were carried out on the MTS370.10 hydraulic experiment system according to ASTM-D5766 standard [13]. As shown in Figure 8(b), the experiment system has both a stationary head and a movable head. During the tensile tests, specimens were clamped between movable head and stationary head, and the movable head was controlled at a stable velocity (i.e., 2 mm/min) with respect to the stationary head. The tension load was applied align the long axis of each specimen. Above the stationary head, there is a force indicator that owns a maximum measure range of ±100 kN, and it can indicate the total force being carried by the specimen with ±1% precision.


#### 3.2.2. Numerical Model

The finite element model of U3160/5284 notched plate is shown in Figure9. The whole model was 300 mm in length, 36 mm in width, 3.34 mm in thickness, and the diameter of the opening hole was 6 mm. Within numerical models, RVEA and RVEB were simulated by a three-dimensional C3D8R mesh and a two-dimensional COH3D8 mesh, respectively. In the thickness direction, the numerical model was discretized into twenty plies of C3D8R meshes. Each three-dimensional C3D8R ply was 0.167 mm in thickness. In addition, between two adjacent C3D8R plies which owned different fiber orientations, there was a two-dimensional cohesive layer discretized by COH3D8 meshes. Displacement constraints were applied at the left end and uniform tension loads were applied at the right end. Table 3 gives the basic material mechanical properties used in numerical models.Figure 9
Numerical mesh of U3160/5284 notched plates.Table 3
Material properties’ symbols and values of U3160/5284 layer.
ModulusSymbol/unitValueLongitudinal elastic modulusE1/MPa116320Transverse elastic modulusE2, E3/MPa8400In-plane shear modulusG12, G13/MPa4570Out-plane shear modulusG23/MPa3075In-plane Poisson ratioν12,ν130.15Out-plane Poisson ratioν230.30Mass densityρkgmm−31.6E − 06Fiber longitudinal modulusE1f/MPa239000Fiber in-plane shear modulusν12f0.2StrengthSymbol/unitValueLongitudinal tension strengthXT/MPa1413Longitudinal compression strengthXC/MPa993Transverse tension strengthYT, ZT/MPa43Transverse compression strengthYC, ZC/MPa184In-plane shear strengthS12, S13/MPa172Out-plane shear strengthS23/MPa102Fracture toughnessSymbol/UnitValueModel-I critical energy release rateGIC/N (mm)−10.52Model-II critical energy release rateGIIC/N (mm)−10.92Model-III critical energy release rateGIIIC/N (mm)−10.92InterfaceSymbol/unitValueNormal stiffnessKN/MPa (mm)−1120000Shear stiffnessKS, KT/MPa (mm)−143000Normal strengthRN/MPa43Shear strengthRS, RT/MPa60


#### 3.2.3. Results and Discussions

Through conducting the simulation analysis of four numerical models, the initiation and evolution process of multiple damage modes, including fiber failure, interfiber failure, and delamination were predicted. These processes show several similar phenomena which can be illustrated by Figures10 and 11, where the blue cells represent undamaged materials and the red ones represent complete damage:(I)
The interfiber failure (IFF) on 90° ply was always the first damage mode that initiated in all four numerical models. Taking E-Groups’ numerical analysis results for instance (Figure10), IFF on 90° ply initiated at 53% fracture load, and it gradually evolved with increasing load until the final fracture occurred (Figure 11).(II)
In addition to the evolution of monodamage mode that mentioned above, the coupling evolution processes between different damage modes were also found. Taking E-Groups’ numerical analysis results for instance (Figure10), as the matrix crack caused by the IFF on 90° ply also propagated alongside the thickness direction, it promoted the initiation of delamination damage on the adjacent cohesive layer 90/0. Then, the delamination crack grew forward and led to the interfiber failure of 0° ply. Finally, the IFF on 0° ply evolved and caused localized stress concentration, which resulted in the fiber failure.(III)
Fiber failure on Layer 0° obviously reduced plates’ axial tensile loading capability. Taking E-Groups’ numerical analysis results for instance (Figure10), there was only 1% load increment between fiber failure of 0° ply and final fracture of whole plates.Figure 10
Numerical prediction on multiple-modes damage initiation and evolution process.Figure 11
Numerical prediction on initiation and evolution process of inter fiber failure on 90° ply.Figure12 shows four numerical models’ prediction results on final damage distribution of all plies and cohesive layers. The projection image of all damaged elements within each ply and cohesive layer forms corresponding numerical models’ fracture morphology. The experimental records on fractured specimen of four groups are also shown in Figure 12. By comparing these figures, two points can be briefly concluded:(I)
For all four groups, the simulation predictions on fracture morphology are similar with corresponding experimental results.(II)
From both experimental and simulation results, we can find that specimens of four groups own the same type of fracture morphology. All composite specimens fail in tension at the hole and exhibit multiple modes of failure in various plies. This kind of fracture morphology shown in Figure12 shall be denoted as MGM failure, according to the failure mode definition listed in ASTM-D5766 [13]. This notation given by ASTM-D5766 uses the first place to describe damage type (i.e., Multimode damage), the second to describe damage area (i.e., Gage of the specimen), and the last to describe damage location (i.e., Middle of the gage).Figure 12
Simulation predictions and experimental records on notched plates’ fracture morphology: (a) group E; (b) group F; (c) group G; (d) group H.
(a)(b)(c)(d)Figure13 and Table 4 show the experimental data and numerical prediction on ultimate tensile fracture load. For all four groups, the standard deviation of experimental ultimate load values is small, which shows the reliability of the whole test system. The simulation model gives similar ultimate fracture load by comparing with the mean experimental value for each group. The biggest prediction error is less than 10%, which is a satisfactory performance for the mechanical numerical analysis of composites. Synthetically judging from the numerical prediction on damage modes, fracture morphology, and ultimate tensile load, the effectiveness of this proposed model can be verified.Figure 13
Comparison of simulation and experimental data on fracture load.Table 4
Relative error calculation of numerical prediction on fracture load.
Group numberSimulation result (kN)Experimental data (kN)Relative error (%)E33.2331.46 ± 2.535.62F34.2337.60 ± 1.398.95G40.3542.44 ± 2.344.94H43.5643.07 ± 3.681.15


### 3.3. Simulation and Comparison with Tian’s Model

Tian et al. [14] conducted CFRP tensile notched tests and proposed a numerical analysis model. According to geometrical and material parameters given in reference [14], we conducted another numerical validation model by adopting macro-meso correlation model. Comparing with experimental data, this model gives better prediction results than Tian’s model.


#### 3.3.1. Geometrical and Material Parameters

According to reference [14], the length L of notched plate was 220 mm, the diameter D of central hole was 6.3 mm, and the width W was 38 mm. The stacking sequence of laminates was [45/0/−45/90]3S, and the nominal thickness of lamina was 0.1 mm. The corresponding material parameters are listed in Table 5.Table 5
Material Properties used in the calculation of Tian’s Model [14].
Fiber reinforced plateCohesiveXT/MPa2806ν12,ν130.34RN/MPa43XC/MPa1400ν230.43RS, RT/MPa60YT/MPa60G12, G13/MPa5170KN/MPa (mm)−1120000YC/MPa185G23/MPa3980KS, KT/MPa (mm)−143000S12, S13/MPa90E1/MPa161000GIC/J (mm)−20.293S23/MPa90E2, E3/MPa11400GIIC/J (mm)−20.631ρkgmm−31.6E − 06GIIIC/J (mm)−20.631


#### 3.3.2. Numerical Results and Discussions

Figures14 and 15, respectively, show the stress-strain curves and fracture mode obtained from the experiment and two simulation models. Moreover, Table 6 gives the prediction error of ultimate tensile strength obtained by comparing numerical prediction results and experimental mean value.Figure 14
Tensile responses comparison with Tian’s model.Figure 15
Fracture mode prediction comparison with Tian’s model.Table 6
Fracture load prediction comparison with Tian’s model.
ModelSimulation result (MPa)Experimental data (MPa)Relative error (%)Tian534.14477.0211.97This paper455.084.60As shown in Figure14, a deviation band formed with two gray curves represents the experimental results of five different specimens. Due to the initial light nonlinear behaviors of experimental data, both simulation models’ prediction curves locate outside of the deviation band, but in contrast, this model’s results is closer to experimental band than Tian’s model. In addition, on the prediction accuracy of ultimate tensile strength, as shown in Table 6, the relative error of this model is below 5% while Tian-model’s error is larger than 10%. Moreover, Figure 15 shows that the fracture mode predicted by this model is also in good agreement with the experimental results. From above tables and figures, the model proposed in this paper performs better in simulating Tian’s experiments.


## 3.1. Numerical Implementation

The Macro-Meso Correlation Model was numerically implemented as shown in Figure7. The macrolevel mechanical analysis can be conducted by common finite element method, where RVEA and RVEB are represented by 3D stress element and 2D-cohesive element, respectively. The mesolevel damage analysis, including damage initiation, evolution, and characterization is realized by operating user-defined material behavior subroutine.Figure 7
Numerical implementation of macro-meso correlation model.In order to verify the rationality of the macro-meso correlation model, we conducted tensile tests and numerical simulations on U3160/5284 specimens with four types of stacking sequences. The corresponding results are shown in Section3.2. Moreover, in order to verify the applicability of the proposed method to other issues with different material properties and geometric parameters, another simulation was also carried out in Section 3.3 to predict the experimental results provided by investigator Tian et al. [14]. All five numerical finite element models were built in ABAQUS commercial software.


## 3.2. Experiment and Simulation for U3160/5284 Plates

### 3.2.1. Specimen and Equipment

The geometric parameters of U3160/5284 carbon fiber-reinforced laminates are shown in Figure8(a). According to different ply ratios (Table 2), four groups of specimens were manufactured and tested. In this paper, the ply ratio is defined as the quotient obtained when the quantity of plies with a certain fiber orientation angle is divided by the total quantity of whole laminated layers. For instance, specimens from group E totally own twenty plies, including six 0° plies, twelve ±45° plies, and two 90° plies. Accordingly, for specimens from group E, the 0° ply ratio is 30%, the ±45° ply ratio is 60%, and the 90° ply ratio is 10%. All four groups in Table 2 own the same number of total plies, but the 0° ply ratio is gradually increased in order to study the influences of ply ratio on specimens’ notched tensile strength. In addition, each group contains six specimens to reduce the dispersion effects of material properties on experimental results.Figure 8
Image of notched specimen and tension test machine: (a) U3160/5284 composite notched plates; (b) MTS370 hydraulic experiment system.
(a)(b)Table 2
Ply orientations of different groups’ U3160/5284 notched plates.
Group numberPly ratio [0/±45/90]Stacking sequenceLamina thickness (mm)E30/60/10[45/−45/0/45/90/0/45/0/−45/45]S0.167F40/50/10[45/−45/0/0/90/45/0/0/−45/45]SG50/40/10[45/0/0/−45/90/0/−45/0/0/45]SH60/30/10[45/0/0/90/0/0/−45/0/0/45]SThe tensile tests were carried out on the MTS370.10 hydraulic experiment system according to ASTM-D5766 standard [13]. As shown in Figure 8(b), the experiment system has both a stationary head and a movable head. During the tensile tests, specimens were clamped between movable head and stationary head, and the movable head was controlled at a stable velocity (i.e., 2 mm/min) with respect to the stationary head. The tension load was applied align the long axis of each specimen. Above the stationary head, there is a force indicator that owns a maximum measure range of ±100 kN, and it can indicate the total force being carried by the specimen with ±1% precision.


### 3.2.2. Numerical Model

The finite element model of U3160/5284 notched plate is shown in Figure9. The whole model was 300 mm in length, 36 mm in width, 3.34 mm in thickness, and the diameter of the opening hole was 6 mm. Within numerical models, RVEA and RVEB were simulated by a three-dimensional C3D8R mesh and a two-dimensional COH3D8 mesh, respectively. In the thickness direction, the numerical model was discretized into twenty plies of C3D8R meshes. Each three-dimensional C3D8R ply was 0.167 mm in thickness. In addition, between two adjacent C3D8R plies which owned different fiber orientations, there was a two-dimensional cohesive layer discretized by COH3D8 meshes. Displacement constraints were applied at the left end and uniform tension loads were applied at the right end. Table 3 gives the basic material mechanical properties used in numerical models.Figure 9
Numerical mesh of U3160/5284 notched plates.Table 3
Material properties’ symbols and values of U3160/5284 layer.
ModulusSymbol/unitValueLongitudinal elastic modulusE1/MPa116320Transverse elastic modulusE2, E3/MPa8400In-plane shear modulusG12, G13/MPa4570Out-plane shear modulusG23/MPa3075In-plane Poisson ratioν12,ν130.15Out-plane Poisson ratioν230.30Mass densityρkgmm−31.6E − 06Fiber longitudinal modulusE1f/MPa239000Fiber in-plane shear modulusν12f0.2StrengthSymbol/unitValueLongitudinal tension strengthXT/MPa1413Longitudinal compression strengthXC/MPa993Transverse tension strengthYT, ZT/MPa43Transverse compression strengthYC, ZC/MPa184In-plane shear strengthS12, S13/MPa172Out-plane shear strengthS23/MPa102Fracture toughnessSymbol/UnitValueModel-I critical energy release rateGIC/N (mm)−10.52Model-II critical energy release rateGIIC/N (mm)−10.92Model-III critical energy release rateGIIIC/N (mm)−10.92InterfaceSymbol/unitValueNormal stiffnessKN/MPa (mm)−1120000Shear stiffnessKS, KT/MPa (mm)−143000Normal strengthRN/MPa43Shear strengthRS, RT/MPa60


### 3.2.3. Results and Discussions

Through conducting the simulation analysis of four numerical models, the initiation and evolution process of multiple damage modes, including fiber failure, interfiber failure, and delamination were predicted. These processes show several similar phenomena which can be illustrated by Figures10 and 11, where the blue cells represent undamaged materials and the red ones represent complete damage:(I)
The interfiber failure (IFF) on 90° ply was always the first damage mode that initiated in all four numerical models. Taking E-Groups’ numerical analysis results for instance (Figure10), IFF on 90° ply initiated at 53% fracture load, and it gradually evolved with increasing load until the final fracture occurred (Figure 11).(II)
In addition to the evolution of monodamage mode that mentioned above, the coupling evolution processes between different damage modes were also found. Taking E-Groups’ numerical analysis results for instance (Figure10), as the matrix crack caused by the IFF on 90° ply also propagated alongside the thickness direction, it promoted the initiation of delamination damage on the adjacent cohesive layer 90/0. Then, the delamination crack grew forward and led to the interfiber failure of 0° ply. Finally, the IFF on 0° ply evolved and caused localized stress concentration, which resulted in the fiber failure.(III)
Fiber failure on Layer 0° obviously reduced plates’ axial tensile loading capability. Taking E-Groups’ numerical analysis results for instance (Figure10), there was only 1% load increment between fiber failure of 0° ply and final fracture of whole plates.Figure 10
Numerical prediction on multiple-modes damage initiation and evolution process.Figure 11
Numerical prediction on initiation and evolution process of inter fiber failure on 90° ply.Figure12 shows four numerical models’ prediction results on final damage distribution of all plies and cohesive layers. The projection image of all damaged elements within each ply and cohesive layer forms corresponding numerical models’ fracture morphology. The experimental records on fractured specimen of four groups are also shown in Figure 12. By comparing these figures, two points can be briefly concluded:(I)
For all four groups, the simulation predictions on fracture morphology are similar with corresponding experimental results.(II)
From both experimental and simulation results, we can find that specimens of four groups own the same type of fracture morphology. All composite specimens fail in tension at the hole and exhibit multiple modes of failure in various plies. This kind of fracture morphology shown in Figure12 shall be denoted as MGM failure, according to the failure mode definition listed in ASTM-D5766 [13]. This notation given by ASTM-D5766 uses the first place to describe damage type (i.e., Multimode damage), the second to describe damage area (i.e., Gage of the specimen), and the last to describe damage location (i.e., Middle of the gage).Figure 12
Simulation predictions and experimental records on notched plates’ fracture morphology: (a) group E; (b) group F; (c) group G; (d) group H.
(a)(b)(c)(d)Figure13 and Table 4 show the experimental data and numerical prediction on ultimate tensile fracture load. For all four groups, the standard deviation of experimental ultimate load values is small, which shows the reliability of the whole test system. The simulation model gives similar ultimate fracture load by comparing with the mean experimental value for each group. The biggest prediction error is less than 10%, which is a satisfactory performance for the mechanical numerical analysis of composites. Synthetically judging from the numerical prediction on damage modes, fracture morphology, and ultimate tensile load, the effectiveness of this proposed model can be verified.Figure 13
Comparison of simulation and experimental data on fracture load.Table 4
Relative error calculation of numerical prediction on fracture load.
Group numberSimulation result (kN)Experimental data (kN)Relative error (%)E33.2331.46 ± 2.535.62F34.2337.60 ± 1.398.95G40.3542.44 ± 2.344.94H43.5643.07 ± 3.681.15


## 3.2.1. Specimen and Equipment

The geometric parameters of U3160/5284 carbon fiber-reinforced laminates are shown in Figure8(a). According to different ply ratios (Table 2), four groups of specimens were manufactured and tested. In this paper, the ply ratio is defined as the quotient obtained when the quantity of plies with a certain fiber orientation angle is divided by the total quantity of whole laminated layers. For instance, specimens from group E totally own twenty plies, including six 0° plies, twelve ±45° plies, and two 90° plies. Accordingly, for specimens from group E, the 0° ply ratio is 30%, the ±45° ply ratio is 60%, and the 90° ply ratio is 10%. All four groups in Table 2 own the same number of total plies, but the 0° ply ratio is gradually increased in order to study the influences of ply ratio on specimens’ notched tensile strength. In addition, each group contains six specimens to reduce the dispersion effects of material properties on experimental results.Figure 8
Image of notched specimen and tension test machine: (a) U3160/5284 composite notched plates; (b) MTS370 hydraulic experiment system.
(a)(b)Table 2
Ply orientations of different groups’ U3160/5284 notched plates.
Group numberPly ratio [0/±45/90]Stacking sequenceLamina thickness (mm)E30/60/10[45/−45/0/45/90/0/45/0/−45/45]S0.167F40/50/10[45/−45/0/0/90/45/0/0/−45/45]SG50/40/10[45/0/0/−45/90/0/−45/0/0/45]SH60/30/10[45/0/0/90/0/0/−45/0/0/45]SThe tensile tests were carried out on the MTS370.10 hydraulic experiment system according to ASTM-D5766 standard [13]. As shown in Figure 8(b), the experiment system has both a stationary head and a movable head. During the tensile tests, specimens were clamped between movable head and stationary head, and the movable head was controlled at a stable velocity (i.e., 2 mm/min) with respect to the stationary head. The tension load was applied align the long axis of each specimen. Above the stationary head, there is a force indicator that owns a maximum measure range of ±100 kN, and it can indicate the total force being carried by the specimen with ±1% precision.


## 3.2.2. Numerical Model

The finite element model of U3160/5284 notched plate is shown in Figure9. The whole model was 300 mm in length, 36 mm in width, 3.34 mm in thickness, and the diameter of the opening hole was 6 mm. Within numerical models, RVEA and RVEB were simulated by a three-dimensional C3D8R mesh and a two-dimensional COH3D8 mesh, respectively. In the thickness direction, the numerical model was discretized into twenty plies of C3D8R meshes. Each three-dimensional C3D8R ply was 0.167 mm in thickness. In addition, between two adjacent C3D8R plies which owned different fiber orientations, there was a two-dimensional cohesive layer discretized by COH3D8 meshes. Displacement constraints were applied at the left end and uniform tension loads were applied at the right end. Table 3 gives the basic material mechanical properties used in numerical models.Figure 9
Numerical mesh of U3160/5284 notched plates.Table 3
Material properties’ symbols and values of U3160/5284 layer.
ModulusSymbol/unitValueLongitudinal elastic modulusE1/MPa116320Transverse elastic modulusE2, E3/MPa8400In-plane shear modulusG12, G13/MPa4570Out-plane shear modulusG23/MPa3075In-plane Poisson ratioν12,ν130.15Out-plane Poisson ratioν230.30Mass densityρkgmm−31.6E − 06Fiber longitudinal modulusE1f/MPa239000Fiber in-plane shear modulusν12f0.2StrengthSymbol/unitValueLongitudinal tension strengthXT/MPa1413Longitudinal compression strengthXC/MPa993Transverse tension strengthYT, ZT/MPa43Transverse compression strengthYC, ZC/MPa184In-plane shear strengthS12, S13/MPa172Out-plane shear strengthS23/MPa102Fracture toughnessSymbol/UnitValueModel-I critical energy release rateGIC/N (mm)−10.52Model-II critical energy release rateGIIC/N (mm)−10.92Model-III critical energy release rateGIIIC/N (mm)−10.92InterfaceSymbol/unitValueNormal stiffnessKN/MPa (mm)−1120000Shear stiffnessKS, KT/MPa (mm)−143000Normal strengthRN/MPa43Shear strengthRS, RT/MPa60


## 3.2.3. Results and Discussions

Through conducting the simulation analysis of four numerical models, the initiation and evolution process of multiple damage modes, including fiber failure, interfiber failure, and delamination were predicted. These processes show several similar phenomena which can be illustrated by Figures10 and 11, where the blue cells represent undamaged materials and the red ones represent complete damage:(I)
The interfiber failure (IFF) on 90° ply was always the first damage mode that initiated in all four numerical models. Taking E-Groups’ numerical analysis results for instance (Figure10), IFF on 90° ply initiated at 53% fracture load, and it gradually evolved with increasing load until the final fracture occurred (Figure 11).(II)
In addition to the evolution of monodamage mode that mentioned above, the coupling evolution processes between different damage modes were also found. Taking E-Groups’ numerical analysis results for instance (Figure10), as the matrix crack caused by the IFF on 90° ply also propagated alongside the thickness direction, it promoted the initiation of delamination damage on the adjacent cohesive layer 90/0. Then, the delamination crack grew forward and led to the interfiber failure of 0° ply. Finally, the IFF on 0° ply evolved and caused localized stress concentration, which resulted in the fiber failure.(III)
Fiber failure on Layer 0° obviously reduced plates’ axial tensile loading capability. Taking E-Groups’ numerical analysis results for instance (Figure10), there was only 1% load increment between fiber failure of 0° ply and final fracture of whole plates.Figure 10
Numerical prediction on multiple-modes damage initiation and evolution process.Figure 11
Numerical prediction on initiation and evolution process of inter fiber failure on 90° ply.Figure12 shows four numerical models’ prediction results on final damage distribution of all plies and cohesive layers. The projection image of all damaged elements within each ply and cohesive layer forms corresponding numerical models’ fracture morphology. The experimental records on fractured specimen of four groups are also shown in Figure 12. By comparing these figures, two points can be briefly concluded:(I)
For all four groups, the simulation predictions on fracture morphology are similar with corresponding experimental results.(II)
From both experimental and simulation results, we can find that specimens of four groups own the same type of fracture morphology. All composite specimens fail in tension at the hole and exhibit multiple modes of failure in various plies. This kind of fracture morphology shown in Figure12 shall be denoted as MGM failure, according to the failure mode definition listed in ASTM-D5766 [13]. This notation given by ASTM-D5766 uses the first place to describe damage type (i.e., Multimode damage), the second to describe damage area (i.e., Gage of the specimen), and the last to describe damage location (i.e., Middle of the gage).Figure 12
Simulation predictions and experimental records on notched plates’ fracture morphology: (a) group E; (b) group F; (c) group G; (d) group H.
(a)(b)(c)(d)Figure13 and Table 4 show the experimental data and numerical prediction on ultimate tensile fracture load. For all four groups, the standard deviation of experimental ultimate load values is small, which shows the reliability of the whole test system. The simulation model gives similar ultimate fracture load by comparing with the mean experimental value for each group. The biggest prediction error is less than 10%, which is a satisfactory performance for the mechanical numerical analysis of composites. Synthetically judging from the numerical prediction on damage modes, fracture morphology, and ultimate tensile load, the effectiveness of this proposed model can be verified.Figure 13
Comparison of simulation and experimental data on fracture load.Table 4
Relative error calculation of numerical prediction on fracture load.
Group numberSimulation result (kN)Experimental data (kN)Relative error (%)E33.2331.46 ± 2.535.62F34.2337.60 ± 1.398.95G40.3542.44 ± 2.344.94H43.5643.07 ± 3.681.15


## 3.3. Simulation and Comparison with Tian’s Model

Tian et al. [14] conducted CFRP tensile notched tests and proposed a numerical analysis model. According to geometrical and material parameters given in reference [14], we conducted another numerical validation model by adopting macro-meso correlation model. Comparing with experimental data, this model gives better prediction results than Tian’s model.


### 3.3.1. Geometrical and Material Parameters

According to reference [14], the length L of notched plate was 220 mm, the diameter D of central hole was 6.3 mm, and the width W was 38 mm. The stacking sequence of laminates was [45/0/−45/90]3S, and the nominal thickness of lamina was 0.1 mm. The corresponding material parameters are listed in Table 5.Table 5
Material Properties used in the calculation of Tian’s Model [14].
Fiber reinforced plateCohesiveXT/MPa2806ν12,ν130.34RN/MPa43XC/MPa1400ν230.43RS, RT/MPa60YT/MPa60G12, G13/MPa5170KN/MPa (mm)−1120000YC/MPa185G23/MPa3980KS, KT/MPa (mm)−143000S12, S13/MPa90E1/MPa161000GIC/J (mm)−20.293S23/MPa90E2, E3/MPa11400GIIC/J (mm)−20.631ρkgmm−31.6E − 06GIIIC/J (mm)−20.631


### 3.3.2. Numerical Results and Discussions

Figures14 and 15, respectively, show the stress-strain curves and fracture mode obtained from the experiment and two simulation models. Moreover, Table 6 gives the prediction error of ultimate tensile strength obtained by comparing numerical prediction results and experimental mean value.Figure 14
Tensile responses comparison with Tian’s model.Figure 15
Fracture mode prediction comparison with Tian’s model.Table 6
Fracture load prediction comparison with Tian’s model.
ModelSimulation result (MPa)Experimental data (MPa)Relative error (%)Tian534.14477.0211.97This paper455.084.60As shown in Figure14, a deviation band formed with two gray curves represents the experimental results of five different specimens. Due to the initial light nonlinear behaviors of experimental data, both simulation models’ prediction curves locate outside of the deviation band, but in contrast, this model’s results is closer to experimental band than Tian’s model. In addition, on the prediction accuracy of ultimate tensile strength, as shown in Table 6, the relative error of this model is below 5% while Tian-model’s error is larger than 10%. Moreover, Figure 15 shows that the fracture mode predicted by this model is also in good agreement with the experimental results. From above tables and figures, the model proposed in this paper performs better in simulating Tian’s experiments.


## 3.3.1. Geometrical and Material Parameters

According to reference [14], the length L of notched plate was 220 mm, the diameter D of central hole was 6.3 mm, and the width W was 38 mm. The stacking sequence of laminates was [45/0/−45/90]3S, and the nominal thickness of lamina was 0.1 mm. The corresponding material parameters are listed in Table 5.Table 5
Material Properties used in the calculation of Tian’s Model [14].
Fiber reinforced plateCohesiveXT/MPa2806ν12,ν130.34RN/MPa43XC/MPa1400ν230.43RS, RT/MPa60YT/MPa60G12, G13/MPa5170KN/MPa (mm)−1120000YC/MPa185G23/MPa3980KS, KT/MPa (mm)−143000S12, S13/MPa90E1/MPa161000GIC/J (mm)−20.293S23/MPa90E2, E3/MPa11400GIIC/J (mm)−20.631ρkgmm−31.6E − 06GIIIC/J (mm)−20.631


## 3.3.2. Numerical Results and Discussions

Figures14 and 15, respectively, show the stress-strain curves and fracture mode obtained from the experiment and two simulation models. Moreover, Table 6 gives the prediction error of ultimate tensile strength obtained by comparing numerical prediction results and experimental mean value.Figure 14
Tensile responses comparison with Tian’s model.Figure 15
Fracture mode prediction comparison with Tian’s model.Table 6
Fracture load prediction comparison with Tian’s model.
ModelSimulation result (MPa)Experimental data (MPa)Relative error (%)Tian534.14477.0211.97This paper455.084.60As shown in Figure14, a deviation band formed with two gray curves represents the experimental results of five different specimens. Due to the initial light nonlinear behaviors of experimental data, both simulation models’ prediction curves locate outside of the deviation band, but in contrast, this model’s results is closer to experimental band than Tian’s model. In addition, on the prediction accuracy of ultimate tensile strength, as shown in Table 6, the relative error of this model is below 5% while Tian-model’s error is larger than 10%. Moreover, Figure 15 shows that the fracture mode predicted by this model is also in good agreement with the experimental results. From above tables and figures, the model proposed in this paper performs better in simulating Tian’s experiments.


## 4. Conclusion

In this paper, the model proposed builds the correlation between mesoscopic damage initiation and macroscopic damage propagation. Four kinds of carbon fiber reinforced laminates with different stacking sequences were prepared, and their ultimate notched strength was measured. The macro-meso damage model proposed in this paper was applied to simulate the tensile responses of specimens. The ultimate failure load and fracture morphology prediction results are in good agreement with experimental data, which shows the effectiveness and rationality of this new model. Moreover, in order to verify the applicability to other material and geometrical parameters, another simulation was conducted and acceptable predictions were also acquired.



---

*Source: 5812912-2020-05-18.xml*