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

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 ﬁled. 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 eﬀectiveness and rationality of this new model. In addition, this model’s applicability to diﬀerent material and geometrical parameters was also veriﬁed by simulating the experiment results in the literature.


Introduction
e 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. en, Chang et al. [7] studied the influences of in-situ strength and shear nonlinearity on notched strength numerical prediction. In the 21 st 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]. e numerical progressive analysis was conducted through a macro-meso correlation model. e 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.

Macro-Meso Correlation Model
e Macro-Meso Correlation Model is a two-phase progressive damage analysis architecture based on the Representative Volume Element (RVE). e 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) e purpose of the external phase is to establish the connection between external loading and the stress and strain distributions on RVEA and RVEB. e 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) e purpose of the internal phase is to determine damage initiation and damage evolution within RVEA and RVEB.
e simplified logical framework of macro-meso correlation analysis model is shown in Figure 1.

Macrolevel Mechanical Analysis.
e main purpose of the macrolevel analysis is to determine the continuous internal displacement field and stress distribution under external loading cases. Figure 2 shows the framework of macrolevel analysis. RVEA are eight-node elements under 3D stress state, which is described by three independent normal stresses σ i (i � 1, 2, 3) and three shear stresses τ j (j � 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 τ k (k � s, t).
e 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) where, δ i (i � n, s, t) are interfacial relative displacements, ε j (j � 11, 22, 33) are normal strains, c k (k � 12, 13, 23) are shear strains, and K i (i � n, s, t) are interfacial stiffness parameters. In addition, C j (j � 11, 22, 33), C k (k � 12, 13, 23), and C l (l � 44, 55, 66) are intralaminar stiffness parameters, which could be determined by equation (2b). (2c)

Mesolevel Damage Analysis.
According to stress field distribution obtained in macrolevel analysis, damage behaviors within RVEA and RVEB are analyzed. Figure 3 shows the framework of mesolevel analysis. e purpose of this internal level analysis is to determine the damage field and update material properties of RVE. e 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.

Damage Initiation.
In this paper, Puck criterion was used to determine RVEA damage initiation. Puck criterion gives the expression of the damage risk factor f E according to the relationship between local stress and local strength at the material point. Damage initiation occurs when f E 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).
where σ θn is the normal stress on the potential IFF plane, and τ θj (j � s, t) are the shear stresses on the IFF plane. eir determination method is shown in equation (4). Parameters R A ⊥ψ , R A ⊥⊥ , and p T(C) ⊥ψ can be acquired by using equation (5). In addition, the value of mean magnification factor m σf is recommended in Table 1. σ θn � σ 2 cos 2 θ + σ 3 cos 2 θ + 2τ 23 sin θ cos θ, σ θs � τ 31 sin θ + τ 12 cos θ, In the above formulae, p T(C) ⊥⊥ and p T(C) ⊥‖ 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): where R n is the interfacial normal stress and R k (k � s, t) are the interfacial shear stresses.

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 When the strain energy density reaches to its critical value g c or the strain energy release rate reaches G c , the damage propagation is completed, which means final failure occurs within RVEA or RVEB. e damage variable under different evolution states can be synthetically represented by equation (7): where d IFF and d DEA 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. ese 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 ε 0 m 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 ε f m or δ f , the corresponding damage variables equal to one, which represents fully damaged condition.

Damage Characterization.
e damage characterization is the data transition module from mesolevel damage analysis to macrolevel analysis ( Figure 5).
rough 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): 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.  Damage characterization Mathematical Problems in Engineering

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 Figure 6. 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.

Numerical Implementation.
e Macro-Meso Correlation Model was numerically implemented as shown in Figure 7. e 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. e mesolevel damage analysis, including damage initiation, evolution, and characterization is realized by operating user-defined material behavior subroutine.
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. e corresponding results are shown in Section 3.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.

Specimen and Equipment.
e geometric parameters of U3160/5284 carbon fiber-reinforced laminates are shown in Figure 8(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°p ly 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. e 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. e 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.

Numerical
Model. e finite element model of U3160/ 5284 notched plate is shown in Figure 9. e 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 twodimensional 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.

Results and Discussions.
rough 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.
ese processes show several similar phenomena which can be illustrated by Figures 10 and 11, where the blue cells represent undamaged materials and the red ones represent complete damage: (I) e 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 ( Figure 10), IFF on 90°p ly 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 (Figure 10), 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. en, the delamination crack grew forward and led to the interfiber failure of 0°p ly. 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 (Figure 10), 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. 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. is kind of fracture morphology shown in Figure 12 shall be denoted as MGM failure, according to the failure mode definition listed in ASTM-D5766 [13]. is 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 13 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. e simulation model gives similar ultimate fracture load by comparing with the mean experimental value for each group. e 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.

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.

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. e stacking sequence of laminates was [45/0/ −45/90] 3S , and the nominal thickness of lamina was 0.1 mm. e corresponding material parameters are listed in Table 5. Figures 14 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.

Numerical Results and Discussions.
As shown in Figure 14, 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.

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. e macro-meso damage model proposed in this paper was applied to simulate the tensile responses of specimens. e 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.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
ere are no conflicts of interest in the research work presented in this manuscript.