A modified progressive damage model for simulating low-velocity impact of composite laminates

Abstract

A modified progressive damage model is constructed to analyze the damage mechanics and damage development of composite laminates induced by low-velocity impact in this article. The damage modes are judged by the 3D Hashin failure criterion. A modified damage evolution model, which was constructed based on through-thickness normal strain component $ε_{33}$ , is implemented to describe progressive damage of composites. Cohesive elements with quadratic failure criterion and B-K criterion are applied to simulate the development of delamination. The 3D Hashin criterion and modified damage evolution model are coded in VUMAT and called in the ABAQUS/Explicit package. The damage distribution and mechanical behavior, including impactor energy, react force, displacement, predicted by numerical simulation are compared with the experimental data of different impact energies (7.35, 11.03, and 14.70 J). The numerical results and experimental data are in good agreement, which suggests that the modified damage evolution model is beneficial for studying the dynamic mechanical behavior during impact. Moreover, the influence of cohesive element thickness on numerical results is discussed. It is concluded that the cohesive element thickness should be adopted between 0.001 and 0.0075 mm.

Keywords

Composite laminates dynamic impact behavior damage evolution model damage distribution finite element method

Introduction

Compared with traditional aerospace materials, composite materials are conducive to improving structural strength and reducing weight.¹ However, their poor impact resistance may raise barely visible impact damage (BVID)^2,3 subjected to low-velocity impact, which will seriously endanger the safety and overall strength of aircraft. Numerical simulation is considered as an efficient tool to describe the complicated dynamic impact behavior of composite laminates.^4,5 Therefore, a finite element model that can simulate the dynamic mechanical behavior and damage development of composite laminates induced by low-velocity impact is needed to be established.⁶

The complex mechanical properties of composite materials can be simulated by the constitutive models.⁷ Among these developed constitutive models, progressive damage model (PDM)^8,9 including failure criteria and damage evolution model has been introduced in many studies to describe the damaged mechanics of composite materials.^10,11 Considering that composite materials are usually composed of fiber and matrix, it is necessary to apply failure criteria to determine the specific damage mode of fiber and matrix. In recent years, the 3D Hashin criterion,^12,13 Chang-Chang criterion,^14,15 and Puck criterion^16,17 have been adopted in many studies. Among the three failure criteria, except the 3D Hashin criterion, the matrix damage initial conditions of the other two criteria ignore the through thickness stress component $σ_{33}$ , which makes the 3D Hashin criterion more advantageous in predicting the damage initiation of composite materials in three-dimensional cases.^2,7,12,13

Once one of fiber or matrix damage modes appears during impact, it is necessary to adopt appropriate damage evolution model to simulate the stiffness degradation of composite laminates. In previous studies, the predefined constants^18,19 or exponential degradation models²⁰ can successfully simulate the reduction of the stiffness. Lisbôa et al.²¹ determined the damage variables of composite cylinders under radial considering the winding pattern through experimental test. The experimental results are in good agreement with the finite element results. Almeida et al.²² studied the open-hole tension characteristics with variable-axial composite laminates by using a predefined constant degradation model. However, numerous of empirical parameters or empirical formulas limit the application range and prediction accuracy of these two degradation models. Sun et al.,²³ Liu et al.,²⁴ Long et al.²⁵ developed the energy evolution model in their works to study the progressive damage process of composite materials. Ferreira et al.²⁶ proposed a UEL (User Element subroutine) based on a continous damage model to evaluate the damage of composite materials. Almeida et al.^27,28 studied the buckling and post-buckling properties of filament wound composite tube under axial compression by adopting progressive damage model based on energy evolution model. Although the damage mechanical behavior of composites in 3D was studied in their works,^29,30 only the transverse strain $ε_{22}$ was considered in the matrix damage evolution model, which suggests that what they actually study is the progressive damage process of composites in two dimensions. Considering that the initial condition of the matrix damage criterion of 3D Hashin is “ $σ_{22} + σ_{33}$ ,” the influence of $ε_{22}$ and $ε_{33}$ should be comprehensively considered when calculating the matrix damage variables. In present study, an improved three-dimensional damage evolution model is constructed by calculating the equivalent strain of $ε_{22}$ and $ε_{33}$ . Compared with original energy evolution model, the proposed method is more suitable for simulating the damage mechanical behavior of composite materials under low-velocity impact in three-dimensional cases and can effectively improve the accuracy of FEM.

The inter-laminar delamination, which will seriously reduce the overall strength of composite laminates, is a typical damage mode induced by low-velocity impact.³¹ Therefore, it is necessary to apply an appropriate tool to simulate the development of inter-laminar delamination. There are two efficient tool that can be adopted to predict delamination, including virtual crack closure technique (VCCT)³² and cohesive zone model (CZM).³³ However, the crack initiation and propagation path need to be predefined when the VCCT is implemented, which cannot satisfy the requirements of damage prediction. In contrast, CZM with damage initiation criterion and damage evolution model will be more suitable for simulating the development of delamination.^2,7

In this article, a finite element model that can simulate the dynamic mechanical behavior and damage development of composite laminates is implemented in ABAQUS/Explicit package subjected to low-velocity impact. The stiffness degradation of composite materials is simulated by a proposed three-dimensional damage evolution model including $ε_{33}$ . Cohesive elements based on CZM are applied to describe the distribution of delamination. The damage distribution and the mechanical behavior of three impact energies (7.35, 11.03, and 14.70 J) predicted by finite element method (FEM) are compared with the experimental data provided by Shi et al.² to verify the advantages of proposed damage evolution model. The good agreement between FEM and experiment suggests that the modified progressive damage model is beneficial to providing reliable advices for engineering applications. Moreover, the sensitivity of cohesive element thickness is analyzed to study its influence on energy absorption and peak force.

Methodology

Three aspects should be considered in finite element model of composite materials subjected to low-velocity impact: (1) Appropriate composite failure criteria to determine each damage modes; (2) Appropriate damage evolution model is applied according to the damage modes of composite materials; (3) Appropriate numerical technology to simulate the initiation and development of inter-laminar delamination.

Failure initiation criterion

The traditional failure criteria, which can only judge the failure of composite materials,² are mainly based on equivalent stress and equivalent strain, such as Tsai-Wu and Tsai-Hill criterion. However, the manufacturing characteristics of composite materials require that the damage modes of matrix and fiber should be clearly observed through failure criteria. The 3D Hashin criterion, as the most popular failure criterion, can accurately predict four main damage modes of composite materials, and has been widely used in numerous studies because the influence of $σ_{33}$ on matrix damage is taken into consideration.^2,7,29 The 3D Hashin criterion including tensile and compressive damage of fiber and matrix is shown in following four equations, respectively:

Fiber tensile damage ( $σ_{11} \geq 0$ ):

F_{ft} = {(\frac{σ_{11}}{X_{T}})}^{2} + {(\frac{σ_{12}}{S_{12}})}^{2} + {(\frac{σ_{13}}{S_{13}})}^{2} \geq 1

(1)

Fiber compressive damage $(σ_{11} < 0)$ :

F_{fc} = {(\frac{σ_{11}}{X_{C}})}^{2} \geq 1

(2)

Matrix tensile damage $(σ_{22} + σ_{33} \geq 0)$ :

\begin{matrix} F_{mt} = {(\frac{σ_{22} + σ_{33}}{Y_{T}})}^{2} + \frac{1}{S_{23}^{2}} (σ_{23}^{2} - σ_{22} σ_{33}) + {(\frac{σ_{12}}{S_{12}})}^{2} \\ + {(\frac{σ_{13}}{S_{13}})}^{2} \geq 1 \end{matrix}

(3)

Matrix compressive damage $(σ_{22} + σ_{33}) < 0$ :

\begin{matrix} F_{mc} = {(\frac{σ_{22} + σ_{33}}{2 S_{23}})}^{2} + \frac{σ_{22} + σ_{33}}{Y_{C}} [{(\frac{Y_{C}}{2 S_{23}})}^{2} - 1] \\ + \frac{1}{S_{23}^{2}} (σ_{23}^{2} - σ_{22} σ_{33}) + {(\frac{σ_{12}}{S_{12}})}^{2} + {(\frac{σ_{13}}{S_{13}})}^{2} \geq 1 \end{matrix}

(4)

Where $σ_{11}, σ_{22}, σ_{33}$ are the stress component in normal and transverse directions, respectively; $σ_{12}, σ_{13}, σ_{23}$ are the corresponding out of plane shear modulus; $X_{T}, X_{C}$ are fiber tensile and compressive strength, respectively; $Y_{T} . Y_{C}$ matrix tensile and compressive strength, respectively; $S_{12}, S_{13}, S_{23}$ are shear strengths of corresponding out of plane.

Damage evolution model

When the value of damage criterion is equal to 1, it is necessary to apply suitable damage evolution model to describe the continuous degradation process of composite stiffness. The damage evolution model including fracture energy of fiber and matrix has been introduced in many studies.^34,35 The strength and damage variables of each element nodes are closely related to their strains $δ_{I}$ . The relationships between the strain of element node and damage variables and composite strength are illustrated in Figure 1. The computation of damage variables for each damage modes is shown in equation (5). When the value of strain $δ_{I}$ increases to the initial failure strain $δ_{I}^{0}$ , the element has reached the bearing limit and entered the stage of damage evolution. With the increase of strain $δ_{I}$ , the element stiffness decreases continuously, and the value of damage variable increases monotonically. The sign of the end of the damage evolution stage is that the strain $δ_{I}$ increases to the final failure strain $δ_{I}^{f}$ , which indicates that the element has completely lost its load-bearing capacity.

d_{I} = \frac{δ_{I}^{f} (δ_{I} - δ_{I}^{0})}{δ_{I} (δ_{I}^{f} - δ_{I}^{0})} (I = f t, f c, m t, m c)

(5)

Figure 1.

The constitutive model of damage evolution model based on fracture energy of composite materials.

Fiber tensile and compressive damage evolution model

According to the initial conditions of fiber tensile damage ( $σ_{11} \geq 0$ ) and fiber compression damage $(σ_{11} < 0)$ of the 3D Hashin failure criterion, it is found that the fiber damage modes are only related to the $σ_{11}$ . Therefore, the value of strain $δ_{I}$ in equation (5) is the strain component $ε_{11}$ along the fiber direction of each element node. The damage evolution model corresponding to composite fiber tensile damage and fiber compression damage is given as follows:

d_{I} = \frac{δ_{I}^{f} (ε_{11} - δ_{I}^{0})}{ε_{11} (δ_{I}^{f} - δ_{I}^{0})} (d_{I} \in [0, 1], I = ft, fc)

(6)

δ_{ft (fc)}^{0} = \frac{X_{T (C)}}{E_{11}}

(7)

δ_{ft (fc)}^{f} = \frac{2 G_{ft (fc)}}{l_{c} X_{T (C)}}

(8)

δ_{I} = m a x (ε_{11}, δ_{f t (f c)}^{0})

(9)

Where $X_{T}$ represents the fiber tensile strength and $X_{C}$ denotes fiber compressive strength. $G_{ft}$ is fracture energy of fiber tensile damage. $G_{fc}$ is fracture energy of fiber compressive damage. $l_{c}$ represents the characteristic length of elements, which is equal to the cube root of the element volume.

Matrix tensile and compressive damage evolution model

The original energy evolution model only introduced the transverse strain $ε_{22}$ in the matrix damage evolution model.^29,30 Considering the obvious deformation along the thickness direction of composite laminates under impact load and the matrix damage initial condition of 3D Hashin criterion, ignoring the influence of $ε_{33}$ on matrix damage variables makes the original model unable to describe the progressive damage process of composites under low-velocity impact in three dimensions. In this paper, an improved three-dimensional damage evolution model is constructed by computing the equivalent strain of $ε_{22}$ and $ε_{33}$ when the matrix damage variables are calculated. Compared with the original model, the proposed method can more accurately describe the dynamic mechanical response of composite laminates in three-dimensional cases. Table 1 shows the difference between the original damage evolution model and the modified damage evolution model in the computation of matrix equivalent strain. The proposed damage evolution model corresponding to matrix tensile damage and matrix compression damage is given as follows:

d_{I} = \frac{δ_{I}^{f} (δ_{I} - δ_{I}^{0})}{δ_{I} (δ_{I}^{f} - δ_{I}^{0})} (d_{I} \in [0, 1], I = mt, mc)

(10)

δ_{mt (mc)}^{0} = \frac{Y_{T (C)}}{E_{22}}

(11)

δ_{mt (mc)}^{f} = \frac{2 G_{mt (mc)}}{l_{c} Y_{T (C)}}

(12)

δ_{mt (mc)} = \max (δ_{mt (mc)}, δ_{mt (mc)}^{0})

(13)

Where $ε_{22}$ and $ε_{33}$ are the strain of each element node in transverse direction. $G_{mt}$ represents fracture energy of matrix tensile damage and $G_{mc}$ represents fracture energy of matrix compressive damage. The symbol <> denotes the Macaulay operator.

Table 1.

Computation of equivalent strain for the original and modified damage evolution model.

	$δ_{mt}$	$δ_{mc}$
Original	$\sqrt{{ε_{22}}^{2}}$	$\sqrt{- {ε_{22}}^{2}}$
Modified	$\sqrt{{ε_{22}}^{2} + {ε_{33}}^{2}}$	$\sqrt{- {ε_{22}}^{2} + - {ε_{33}}^{2}}$

Damage matrix of composite materials

The damage variables of fiber and matrix should be introduced in the damage matrix of composite materials to capture the changes of the mechanical behavior. The degraded compliance matrix $S_{d}$ is shown in equation (14) and the degraded stiffness matrix $C_{d}$ is illustrated in equation (15):

S_{d} = [\begin{matrix} \begin{matrix} \frac{1}{d_{f} E_{11}} & - \frac{v_{21}}{E_{22}} & - \frac{v_{31}}{E_{33}} \\ - \frac{v_{12}}{E_{11}} & \frac{1}{d_{m} E_{22}} & - \frac{v_{32}}{E_{33}} \\ - \frac{v_{13}}{E_{11}} & - \frac{v_{23}}{E_{22}} & \frac{1}{E_{33}} \end{matrix} & \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \\ \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} & \begin{matrix} \frac{1}{d_{f} d_{m} G_{12}} & 0 & 0 \\ 0 & \frac{1}{d_{f} d_{m} G_{23}} & 0 \\ 0 & 0 & \frac{1}{d_{f} d_{m} G_{31}} \end{matrix} \end{matrix}]

(14)

C_{d} = \frac{1}{Δ} [\begin{matrix} \begin{matrix} d_{f} E_{11} (1 - d_{m} v_{23} v_{32}) & d_{f} d_{m} E_{11} (v_{12} + v_{32} v_{13}) & d_{f} E_{11} (v_{13} + d_{m} v_{12} v_{23}) \\ 0 & d_{m} E_{22} (1 - d_{f} v_{13} v_{31}) & d_{m} E_{22} (v_{23} + d_{f} v_{21} v_{13}) \\ 0 & 0 & E_{33} (1 - d_{f} d_{m} v_{12} v_{21}) \end{matrix} & \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \\ \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} & \begin{matrix} Δ d_{f} d_{m} G_{12} & 0 & 0 \\ 0 & Δ d_{f} d_{m} G_{23} & 0 \\ 0 & 0 & Δ d_{f} d_{m} G_{13} \end{matrix} \end{matrix}]

(15)

{\begin{matrix} d_{f} = (1 - d_{ft}) (1 - d_{fc}) \\ d_{m} = (1 - S_{mt} d_{mt}) (1 - S_{mc} d_{mc}) \\ Δ = 1 - d_{f} d_{m} v_{12} v_{21} - d_{m} v_{23} v_{32} - d_{f} v_{13} v_{31} \\ - 2 d_{f} d_{m} v_{21} v_{32} v_{13} \end{matrix}

(16)

Where the total fiber damage variable $d_{f}$ is calculation by fiber tensile damage variables $d_{ft}$ and fiber compression damage variable $d_{fc}$ . $d_{mt}, d_{mc}, S_{mt}, S_{mc}$ represent matrix tensile damage variables, matrix compression damage variable, matrix tensile damage shear stiffness degradation coefficient, and matrix compressive damage shear stiffness degradation coefficient, respectively and these four parameters can determine the total matrix damage variables $d_{m}$ . In this study, the value of $S_{mt}, S_{mc}$ are equal to 0.9 and 0.5, respectively.

Delamination

In this paper, the cohesive elements with bilinear traction-separation relationship are inserted between layers to predict the initiation and development of inter-laminar delamination. The Quadratic failure criterion is used to judge the initiation of delamination, as shown in equation (17). The Benezeggaph and Kenane (B-K) criterion^36,37 is used to describe the evolution of delamination, as shown in Eq. 18.

\frac{{t_{n}}^{2}}{T_{n}^{2}} + \frac{t_{s}^{2}}{T_{s}^{2}} + \frac{t_{t}^{2}}{T_{t}^{2}} = 1

(17)

Where $t_{n}, t_{s}, t_{t}$ represent the normal and shear tractions; $T_{n}, T_{s}, T_{t}$ represent the strength in the normal and shear direction.

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

(18)

Where $G^{C}, G_{n}^{C}, G_{s}^{C}$ are the total, normal, and shear critical fracture energy, respectively. $G_{S}$ is the dissipated energy in the out-of-plane direction; $G_{T}$ is the total dissipated energy; $η$ is material coefficient which is relate to the type of composite materials and set to 1.45.

Finite element model for low-velocity impact

Dimension parameters and boundary conditions

In this study, the finite element model composed of composite laminate and impactor is established in ABAQUS/Explicit package to capture the changes of mechanical behavior. The geometric parameters are illustrated in Figure 2. The geometric dimensions of the finite element model are provided by Shi et al.² The composite laminate with a diameter of 75 mm is composed of eight layers with a thickness of 0.25 mm and seven interfaces with a thickness of 0.0075 mm and the stacking sequence is [0/90]_2S. Considering that the semicircular impactor with a diameter of 15 mm will not be significantly deformed during the low-velocity impact, it is regarded as an analytical rigid impactor to reduce the computation cost. Constrain the freedom of impactor in the X and Y directions to ensure that the direction of movement of the impactor is always perpendicular to the composite laminate. In order to obtain the three impact energies of 7.35, 11.03, and 14.7 J in the experiment test, the mass of the impactor with a speed of 3.83 m/s is set to 1, 1.5 and 2 kg respectively.

Figure 2.

The dimensional parameters and boundary conditions of finite element model.

Mesh types and contact algorithm

In this model, seven interfaces with a thickness of 0.0075 mm, which are composed of 13,860 eight-node three-dimensional cohesive elements (COH3D8), are inserted into layers to simulate the initiation and evolution of delamination. In order to avoid excessive deformation of composite laminates, the maximum degradation of cohesive element is controlled at 0.99. The material properties of cohesive element provided by Zhou et al.⁷ are shown in Table 2. Eight layers, which are composed of 15,840 eight-node elements with reduced integration (C3D8R), are used to simulate the intra-laminar damage of composite materials. The distortion control is activated to avoid excessive element deformation. The properties of composite materials provided by Shi et al.² is shown in Table 3. Appropriate mesh density can ensure the accuracy and efficiency of the computation. Therefore, the size of 1 mm × 1 mm elements are applied to simulate the experimental damage zone.

Table 2.

Mechanical properties of interface.^2,7

Elastic modulus (MPa)	$E_{n} = E_{s} = E_{t} = 5000$
Normal and shear strength (MPa)	$T_{n} = T_{s} = T_{t} = 30$
Fracture energy (N/mm)	$G_{n}^{C} = 0.6; G_{s}^{C} = 2.1$

Table 3.

Mechanical properties of the carbon fiber/epoxy composite materials.^2,7

Density	1600/(kg/m³)
Engineering constant	Young’s modulus (GPa)	$E_{11} = 153; E_{22} = E_{33} = 10.3$
	Poisson’s ratio	$v_{12} = v_{13} = 0.3; v_{23} = 0.4$
	Shear modulus (GPa)	$G_{12} = G_{13} = 6; G_{23} = 3.70$
Strength properties	Fiber tensile and compressive strength (MPa)	$X_{T} = 2537; X_{C} = 1580$
	Matrix tensile and compressive strength (MPa)	$Y_{T} = 82; Y_{C} = 236$
	Shear strength (MPa)	$S_{12} = S_{13} = 90; S_{23} = 40$
Fracture Energy	Fiber tensile and compressive fracture energy (N/mm)	$G_{ft} = 91.6; G_{fc} = 79.9$
	Matrix tensile and compressive fracture energy (N/mm)	$G_{mt} = 0.22; G_{mc} = 1.1$

A suitable contact algorithm shoule be implemented to simulate the contact between the impactor and the composite laminate, which will be beneficial to obtain more accurate numerical results. The general contact algorithm in ABAQUS/Explicit package is adopted in this article, and the friction coefficient between the impactor and the composite laminate is 0.5.^38,39

Numerical results and discussion

In this section, the results of FEM and experiment results of Shi et al.² under three impact energies (7.35, 11.03, and 14.70 J) are compared to verify the effectiveness of the modified progressive damage model. Moreover, the sensitivity of cohesive element thickness is studied.

Model validation

In this part, the mechanical behavior changes of composite laminates predicted by FEM, including impactor energy, react force, displacement, and damage development, are compared to experimental data under three impact energies (7.35, 11.03, and 14.70 J). By comparing the differences between original and modified damage evolution model, the advantages of $ε_{33}$ on the numerical simulation is discussed.

Mechanical behavior changes of composite laminates

Figure 3 illustrates the energy-times curves of impactor predicted by original and modified damage evolution model. The original damage evolution model in method A only considers $ε_{22}$ when calculating the matrix damage variables. The modified damage evolution model in method B includes $ε_{22}$ and $ε_{33}$ . It is found that the trend of absorbed energy shows a “slow-fast-slow” law. In the initial impact stage, the lower energy absorption rate is due to the insufficient contact of semicircular impactor and composite laminate, which is determined by the shape of impactor. When the impactor is in full contact with composite laminates, the kinetic energy of the impactor will be quickly diminished, resulting in the area of each damage mode began to appear and spread. As the kinetic energy of the impactor decreases and the damage area increases, the overall stiffness of the composite laminate decreases, which causes the energy absorption rate to decrease. The curve no longer rises indicating that the kinetic energy of the impactor has been completely converted into friction energy, internal energy, and other energy modes. Meantime, the stored internal energy in composite laminate will be released, which will accelerate the impactor to rebound. Finally, the absorbed energy in composite laminates will tend to a stable value.

Figure 3.

Energy-time curves of FEM and experiment under three impact energies: (a) 7.35 J, (b)11.03 J, and (c) 14.7 J.

Comparing Figure 3(a) to (c), the FEM results of three different impact energies have a good agreement with experimental data. The final absorbed energy of two damage evolution models is always lower than experiment, but the difference between FEM and experiment will gradually decrease with the increase of impact energy. The reasons for the phenomenon may come from two aspects: (1) It is difficult for FEM to accurately describe the complex contact behavior induced by low-velocity impact, such as the interaction between the matrix and the fiber in the composite material, and the friction force between the layers. (2) There is still a certain residual bonding strength between the layers because the maximum degradation of cohesive elements is controlled at 0.99. Method B leads to a higher energy absorption rate of composite laminates than method A, which is reflected in the energy-time curve of method B always reaches the peak point before method A, which makes the result of method B closer to the experimental results. The energy-curves at the rebound stage predicted by method B have good relationship with the experimental results, especially under the impact energy of 14.7 J. The reason is that the computation of matrix damage variables in modified damage evolution model is controlled by $ε_{22}$ and $ε_{33}$ , which leads to the trend of damage stiffness matrix $C_{d}$ is closer to experiment.

It is found that considering $ε_{33}$ in damage evolution model is beneficial for numerical simulation. To explain the specific influence of $ε_{33}$ on the computation of matrix damage variables, Table 4 shows the strain components of the eighth layer center element without any damage mode appears on composite laminate at 0.5 ms, including $ε_{11}$ , $ε_{22}$ , and $ε_{33}$ . It is found that although $ε_{33}$ must always be much smaller than $ε_{22}$ , the influence of $ε_{33}$ on the calculation of damage variable $d_{I}$ becomes larger and larger when the impact energy increases from 7.35 to 14.7 J, which indicates that the matrix damage variable of the modified damage evolution model reaches 1 earlier than original damage evolution model. Therefore, the influence of $ε_{33}$ on the computation of matrix damage variables cannot be ignored.

Table 4.

The strain component of central element on the back side.

Impact energy (J)	$ε_{11}$ /mm	$ε_{22}$ /mm	$ε_{33}$ /mm
7.35	0.1711	0.0210	0.0008
11.03	0.1964	0.0317	0.0012
14.70	0.2058	0.0321	0.0015

Figure 4 shows force-time curves of FEM results and experiments of different impact energies (7.35, 11.03, and 14.7 J). It is observed that there will be a more severe vibration amplitude in impact phase when the method A is implemented. In contrast, the overall vibration amplitude predicted by method B is relatively stable. The phenomenon indicates that the modified progressive damage model can make the damage stiffness matrix $C_{d}$ degenerate more continuously, which makes the stiffness degradation process of composite laminates predicted by the FEM closer to the experimental data. It is found that the vibration amplitude is more serious around the peak of force-time curves, while the curve fluctuation in the rebound phase is smaller. The reason for this phenomenon could be that the development of damage modes leads to a decrease in the overall strength of composite laminate. Therefore, the impact load is constantly redistributed, resulting in different vibration amplitude in the force-time curves. Due to large-area damage, the overall strength of composite laminates reaches a minimum near the peak of force-time, and smaller loads can cause larger fluctuations. The area of each damage mode hardly increases at the rebound stage, which makes the overall mechanical behavior and energy release rate relatively stable, and the force-time curves at the rebound stage are smoother.

Figure 4.

Force-time curves of FEM and experiment under three impact energies: (a) 7.35 J, (b) 11.03 J, and (c) 14.7 J.

Comparing Figure 4(a) to (c), the FEM results of three impact energies have a good agreement with experiment force-time curves. It is found that when the impact energy increases from 7.35 to 14.7 J, the differences of the peak force between method A and experimental result becomes larger, while the maximum force predicted by method B has the opposite trend, which agrees well with the experiment. In the rebound phase, the difference between the force-time curves predicted by the two methods under low level impact energy is small. However, under impact energy 14.7 J, there is a significant difference between method A and method B. The reason is that the $ε_{33}$ is introduced in damage evolution model in method B. According to previous discussion, the influence of $ε_{33}$ on the computation of the damage variable $d_{I}$ becomes greater when the level of impact energy increases. The damage variable $d_{I}$ of method B will reach 1 earlier, which makes the stiffness matrix $C_{d}$ enter the degradation phase earlier, and the strength decline of composite laminates is more realistic. Therefore, the vibration of force-time curves predicted by method B becomes smaller and the stability of the analysis is improved.

Figure 5 illustrates the force-displacement curves of FEM and experiment under three impact energies. It is found that the force-displacement curves of method A and method B under impact energy 7.35 and 11.03 J are greater than experimental curves in rebound phase. The force-displacement curve predicted by method B are in good agreement with the experimental data under the impact energy of 14.7 J. The maximum displacement obtained by method B which represents the maximum deformation of composite laminate is closer to the experimental results, and the maximum displacement is almost the same as the experimental results when the impact energies are 7.35 and 11.03 J. This verifies the effectiveness of the proposed damage evolution model which was constructed based on through-thickness normal strain component $ε_{22}$ and $ε_{33}$ again. Under the impact energy 14.7 J, the maximum displacement predicted by method A will have larger error when predicting the mechanical behavior of higher impact energy. The introduction of $ε_{33}$ can improve the mechanics of composite laminates under higher impact energy, especially in rebound stage.

Figure 5.

Force-displacement curves of FEM and experiment under three impact energies: (a) 7.35 J, (b) 11.03 J, and (c)14.7 J.

Damage modes and damage distribution

One advantage of FEM is the convenience to observe dynamic impact process and damage modes. In this part, the dynamic process of impact and failure, matrix tensile damage, matrix compressive damage, and inter-laminar delamination of composite laminate are studied in detailed.

Figure 6 illustrate the process of impact and failure under the impact energy of 11.03 J. At 2.12 ms, the kinetic energy of the impactor has been completely absorbed. It is found that during the entire impact process, multiple failure modes appeared on the composite laminate. From the point of damage propagation rate, the entire impact process could be divided into three stages. (1) The damage develops slowly. Due to the excellent design mechanical properties of composite laminates, the damage propagation rate is low. However, the duration of this stage is short, because the occurrence of fiber damage makes the load-bearing capacity of the composite laminates drop sharply; (2) The damage develops rapidly. The larger impact energy makes each damage mode spread rapidly at this stage. Among them, the propagation rate of matrix tensile damage is the largest, which propagates from the bottom of the composite laminate to the impact surface. The development direction of matrix compression damage is opposite to that of matrix tensile damage. The end of this stage is marked by the kinetic energy of the impactor being 0. (3) Rebound stage. Due to the absence of external impact loads at this stage, the damage modes on the composite laminate hardly propagate.

Figure 6.

The process of impact and failure under the impact energy of 11.03 J.

The distribution of matrix tensile damage on each layer predicted by two damage evolution model under the impact energy of 14.7 J are shown in Figure 7. The elements in the red region have been completely damaged, while the blue denotes undamaged region. It is found that the area of matrix tensile damage is the largest on the fifth, sixth, and seventh layers. Furthermore, the area of matrix tensile damage increases from first layer to seventh layer. The reason for the phenomenon could come from two aspects: (1) The bend deformation of the composite laminate gradually increases from the first layer to the eighth layer, resulting in a larger tensile load; (2) The impact load transmitted by layer aggravates the matrix tensile damage. Also, more serious matrix tensile damage areas appeared on each layer predicted by method B. This can be explained by the introduction of $ε_{33}$ in method B, which makes the matrix tensile damage variables reach 1 earlier than original damage model.

Figure 7.

The distribution of matrix tensile damage on composite laminate under the impact energy of 14.7 J: (a) Method A and(b) Method B.

The matrix tensile damage on each layer shows different shapes. The matrix tensile damage shape of the first to seventh layers is similar to rectangle, but not all the long sides of the rectangle are parallel to the fiber orientation of the corresponding layer, such as the first and sixth layers. The tensile damage shape of the matrix on the eighth layer is similar to butterfly shape. These phenomena suggest that the propagation direction of matrix tensile damage is influenced by fiber direction and impact load. There are tiny hollows shapes with different areas on each layer predicted by method A, but this phenomenon has been effectively alleviated after applying method B, especially in the top four layers. The tiny hollows shape can be attributed to the Hashin criterion which uses the summation of transverse and out-of-plane normal stress to evaluate the tension and compression failure for matrix.⁷ However, the modified progressive damage model including $ε_{33}$ can suppress the appearance of tiny hollowed shape, which denotes that the proposed damage evolution model in this article can overcome the defects of 3D Hashin and improve the accuracy of numerical simulation.

The distribution of matrix tensile damage of composite laminates predicted by two progressive damage model under the impact energy of 14.7 J are shown in Figure 8. It is found that compared with the matrix compression damage, the matrix tensile damage shows a more serious damage area, which indicates that the main reason for the decrease of the overall strength of composite laminates is the excessive area of matrix tensile damage. The area of matrix compressive damage gradually decreases from the first layer to last layer. The reason could be that the C3D8R element and COH3D8 element whose damage variables reaches 1 is not deleted in this paper so that the first layer keeps contact with impactor, which makes the propagation orientation of matrix compressive damage is from first layer to eighth layer. A phenomenon needs to be noted is that although the method A and method B show the similar matrix compression damage area, there are scattered damage near the impact center of each layer predicted by method B. The reason could be that the through thickness strain component $ε_{33}$ is considered in the progressive damage model of method B, which makes the matrix compression damage variable $d_{mc}$ easier to reach 1. On the other hand, method B can increase the safety factor of composite structures. Considering the damage distribution shape and comparing the fiber direction of each layer, it is found that the development orientation of the matrix compression damage is not parallel to the fiber direction of corresponding layer, but spreads from center to edge. The matrix compression damage on the first and second layer is distributed in a butterfly shape, while the damage on the fourth, fifth, and sixth layers is affected by the fiber direction of adjacent layers. The above phenomena suggest that the impact load plays a major role in the development of matrix compression damage.

Figure 8.

The distribution of matrix compressive damage on composite laminate under the impact energy of 14.7 J: (a) Method A and (b) Method B.

The distribution of inter-laminar delamination predicted by FEM and experiment under 14.7 J impact are shown in Figure 9. By comparing Figure 9(a) to (c), it is found that delamination area predicted by FEM is composed of complete damaged region (red region), undamaged region (blue region), and partially damaged region (the rest color region). Among them, the partially damage region surrounds the complete failed region, which indicates that delamination damage extends outward from the center of each layer. The delamination damage area which includes complete damage region and partially damage region predicted by two methods is larger than the X-ray radiograph. Compared to method A, method B shows a more severe area of complete delamination damage. The reason could be that method B shows lager area of matrix tensile damage, which leads to higher interface force.

Figure 9.

The distribution of inter-laminar delamination of 14.7 J impact energy predicted by FEM and experiment: (a) X-ray radiograph, (b) overall delamination distribution of method A, (c) overall delamination distribution of method B, (d) delamination on seven interfaces of method A, and (e) delamination on seven interfaces of method B.

The distribution of delamination on each interface predicted by two damage evolution models are illustrated in Figure 9(d) and (e), respectively. It is found that both two methods capture that the development orientation of delamination on the last five interfaces is parallel to the fiber direction of the layer below the corresponding interface, but the first interface has the opposite trend. The reason for this phenomenon could be that the development orientation of delamination on interface 1 is controlled by the high-level friction force, while the propagation orientation of delamination on last five interfaces is influenced by bending deformation. Both methods have observed that the delamination on the interface 2 is mainly partial damaged region, which indicates that the interface 2 still maintains certain load-bearing capacity. Furthermore, the area displayed on the interface 6 is mainly the complete delamination damage, with few partial failures. The main difference between the two methods for predicting delamination is the area of complete damaged region. Method A shows larger delamination area than method B because the maximum deformation of the composite laminates predicted by two methods shows the similar trend. On the remaining six interfaces, the delamination area predicted by method B is larger than that of method A. It should be noted that in the numerical analysis, more accurate mechanical properties and refined meshes could make the computation results closer to the experiment. However, the certain error between FEM and experiment is beneficial to improving the design allowable value of composite structure and improve the safety factor.

In summary, the modified three-dimensional progressive damage model considering $ε_{33}$ can accurately predict the intra-laminar damage and inter-laminar delamination under low-velocity impact, making it more suitable for simulating the damage mechanical behavior of composite materials in three-dimensional cases and can effectively improve the accuracy of FEM.

Research on degradation model versatility

Another impact event is taken to verify the versatility of the modified damage evolution model in this article.⁴⁰ The dimension parameters of the Graphite/epoxy composite laminate with the stacking sequence of ${[0_{2}^{o} / 90_{2}^{o} / 0_{2}^{o} / 90_{2}^{o} / 0_{2}^{o}]}_{S}$ – ${[0_{2} / 90_{2} / 0_{2} / 90_{2} / 0_{2}]}_{S}$ is 40 mm × 40 mm × 2.1 mm. Table 5 gives the mechanical properties of composite materials and cohesive elements. The semicircular impactor is simplified into a rigid body, whose mass and diameter are 6.5 kg and 12.7 mm, respectively, as shown in Figure 10. The force-time and displacement-time curves were compared with experimental data provided by Kim et al.⁴⁰ and given in Figure 11. Similar to the previous description, method A and method B respectively represent the original and the modified damage evolution model. It is found that the force-time curves and the deformation of composite laminate predicted by method B are more consistent with the experimental results, which indicates that the modified progressive damage evolution model has good versatility and can be applied to the simulate the mechanical behavior changes of different types composite laminates.

Table 5.

Material properties of the graphite/epoxy composites.^40–42

Compositelamina properties	Density	1600 $kg / m^{3}$
	Engineering constant (GPa)	$E_{11} = 143.4; E_{22} = E_{33} = 9.27; v_{12} = v_{13} = 0.31; v_{23} = 0.52; G_{12} = G_{13} = 3.8; G_{23} = 3.2$
	Strength properties (MPa)	$X_{T} = 2945; X_{C} = 1650; Y_{T} = 54; Y_{C} = 240; S_{12} = S_{13} = S_{23} = 89.9$
	Fracture energy (N/mm)	$G_{ft} = 91.6; G_{fc} = 79.9; G_{mt} = 0.22; G_{mc} = 0.76$
Interface properties	Elastic modulus (GPa)	$E_{n} = E_{s} = E_{t} = 5$
	Strength properties (MPa)	$N = S = T = 30$
	Fracture energy (N/mm)	$G_{n}^{C} = 0.6; G_{S}^{C} = 2.1$
	Material coefficient $η$	1.45

Figure 10.

Finite element model subjected to low velocity impact in Kim et al.⁴⁰

Figure 11.

Mechanical behavior of FEM and experiment under different impact energies: (a) RF-T curves of 5 J impact energy,(b) RF-T curves of 15 J impact energy, (c) RF-U curves of 5 J impact energy, and (d) RF-U curves of 15 J impact energy.

Sensitivity study on cohesive element thickness

In this part, the sensitivity of cohesive element thickness is analyzed to study its influence on energy absorption and peak force. Four groups of values (0.001, 0.004, 0.0075, 0.01 mm) under the 14.7 J impact are applied to obtain the changes of energy absorption and peak force, shown in Figure 12. It is found that when the cohesive element thickness is 0.004 mm, the values of absorbed energy is the largest, but the peak force is the smallest. As the cohesive element thickness increases, the energy absorption value and peak force of the composite laminate begin to stabilize. However, the excessively large thickness of the cohesive element not only causes the element to warpage during the calculation process, but also increases the size of the finite element model, making the analysis results of no reference significance. Therefore, according to the analysis results of this paper, the thickness of the cohesive element in the finite element model should be adopted between 0.001 and 0.0075 mm.

Figure 12.

The absorbed energy and peak force of cohesive elements with four thicknesses under the impact energy of 14.7 J.

Conclusion

This paper establishes a finite element model which can simulate the dynamic impact mechanics and damage development of composite laminates induced by low-velocity impact. The damage modes and damage initiation are judged by the 3D Hashin criterion. The stiffness loss of composite materials is described by the proposed modified damage evolution model which includes the $ε_{33}$ . Cohesive elements with Quadratic failure criterion and B-K damage evolution model are applied to simulate the development and distribution of inter-laminar delamination. The 3D Hashin criterion and modified damage evolution model are coded in VUMAT and implemented in ABAQUS/Explicit. Compared with original damage model ignored $ε_{33}$ , the damage distribution and the mechanical behavior of three impact energies (7.35, 11.03, and 14.70 J) predicted by proposed modified damage evolution model are in good agreement with experimental data provided by Shi et al.² The conclusions are followed:

Compared with experimental results, the errors of the peak force predicted by proposed method under three impact energies are 4.7%, 6.9%, and 3.2%, while the errors of the original model are 5.1%, 10.4%, and 7.7%, respectively. Therefore, the proposed method is more suitable for simulating the mechanical behavior of composite laminates under low-velocity impact. The positive influence of considering $ε_{33}$ will be increases as the impact energy increases.

The limitations of 3D Hashin criterion to compute the matrix damage has been well overcome by the modified damage evolution model, which reduces the area of the tiny hollowed region and improves the application and prediction accuracy of the 3D Hashin criterion.

The delamination area predicted by the modified progress damage model is larger than the experiment results, which indicates that the finite element model based on modified progressive damage model with $ε_{33}$ can provide conservative design suggestions for engineering design.

A sensitivity study on cohesive element thickness shows that the absorbed energy and peak force will tend to be stable as the thickness increases, but considering the actual size of structure, the recommended cohesive element thickness should be adopted between 0.001 and 0.0075 mm.

In summary, the modified progressive damage model has good versatility and can be directly applied to analyze the damage mechanics and mechanical behavior of composite laminates with other stacking sequences induced by impact load.

Footnotes

Handling Editor: Chenhui Liang

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work has been supported by the National Natural Science Foundation of China (11702046), Major Theme Project of Chongqing Technology Innovation and Application Development (cstc2019jscx-zdztzx0028, cstc2019j scx-zdztzx0020).

ORCID iD

Niu Liu

References

Chen

Lin

. Evaluation on crashworthiness and energy absorption of composite light airplane. Adv Mech Eng 2018; 10: 1687814018794080.

Shi

Swait

Soutis

. Modelling damage evolution in composite laminates subjected to low velocity impact. Compos Struct 2012; 94: 2902–2913.

Zhang

Duodu

. Finite element modeling of damage development in cross-ply composite laminates subjected to low velocity impact. Compos Struct 2017; 173: 219–227.

Tuo

, et al. An experimental and numerical investigation on low-velocity impact damage and compression-after-impact behavior of composite laminates. Compos Part B Eng 2019; 167: 329–341.

El Moumen

Tarfaoui

Hassoon

, et al. Experimental study and numerical modelling of low velocity impact on laminated composite reinforced with thin film made of carbon nanotubes. Appl Compos Mater 2018; 25: 309–320.

Liao

Jia

Zhou

, et al. An explicit–implicit combined model for predicting residual strength of composite cylinders subjected to low velocity impact. Compos Struct 2020; 247: 112450.

Zhou

Wen

Wang

. Finite element analysis of a modified progressive damage model for composite laminates under low-velocity impact. Compos Struct 2019; 225: 111113.

Kolanu

Raju

Ramji

. A unified numerical approach for the simulation of intra and inter laminar damage evolution in stiffened CFRP panels under compression. Compos B Eng 2020; 190: 107931.

Rozylo

. Stability and failure of compressed thin-walled composite columns using experimental tests and advanced numerical damage models. Int J Numer Methods Eng 2021; 122: 5076–5099.

10.

Donadon

Iannucci

Falzon

, et al. A progressive failure model for composite laminates subjected to low velocity impact damage. Comput Struct 2008; 86: 1232–1252.

11.

Lee

C-S

Kim

J-H

Kim

, et al. Initial and progressive failure analyses for composite laminates using puck failure criterion and damage-coupled finite element method. Compos Struct 2015; 121: 406–419.

12.

Hashin

Rotem

. A fatigue failure criterion for fiber reinforced materials. Compos Mater 1973; 7: 448–464.

13.

Hashin

. Failure criteria for unidirectional fiber composites. J Appl Mech 1980; 47: 329–334.

14.

Chang

. A progressive damage model for laminated composites containing stress concentrations. Compos Mater 1987; 21: 834–855.

15.

Huang

Wang

. Numerical and experimental investigations on the axial crushing response of composite tubes. Compos Struct 2009; 91: 222–228.

16.

Puck

Schürmann

. Failure analysis of FRP laminates by means of physically based phenomenological models. Compos Sci Technol 2002; 62: 1633–1662.

17.

Chen

. Micro–macro FE modeling of damage evolution in laminated composite plates subjected to low velocity impact. Compos Struct 2016; 147: 111–121.

18.

Zhang

. An efficient approach for predicting low-velocity impact force and damage in composite laminates. Compos Struct 2015; 130: 85–94.

19.

Zhang

. Simulating low-velocity impact induced delamination in composites by a quasi-static load model with surface-based cohesive contact. Compos Struct 2015; 125: 51–57.

20.

Singh

Namala

Mahajan

. A damage evolution study of E-glass/epoxy composite under low velocity impact. Compos Part B Eng 2015; 76: 235–248.

21.

Lisbôa

Almeida

JHS

Jr Spickenheuer

, et al. FEM updating for damage modeling of composite cylinders under radial compression considering the winding pattern. Thin-Walled Struct 2022; 173: 108954.

22.

Almeida

JHS

Bittrich

Spickenheuer

. Improving the open-hole tension characteristics with variable-axial composite laminates: optimization, progressive damage modeling and experimental observations. Compos Sci Technol 2020; 185: 107889.

23.

Sun

Tie

, et al. Experimental and numerical investigations of adhesively bonded CFRP single-lap joints subjected to tensile loads. Int J Adhes Adhes 2019; 95: 102402.

24.

Liu

Ding

, et al. The behaviour of thermoplastic and thermoset carbon fibre composites subjected to low-velocity and high-velocity impact. J Mater Sci 2020; 55: 15741–15768.

25.

Long

Yao

Zhang

. Delamination prediction in composite laminates under low-velocity impact. Compos Struct 2015; 132: 290–298.

26.

Ferreira

GFO

Almeida

JHS

Ribeiro

, et al. A finite element unified formulation for composite laminates in bending considering progressive damage. Thin-Walled Struct 2022; 172: 108864.

27.

Almeida

JHS

Jr St-Pierre

Wang

, et al. Design, modeling, optimization, manufacturing and testing of variable-angle filament-wound cylinders. Compos Part B Eng 2021; 225: 109224.

28.

Almeida

JHS

Tonatto

MLP

Ribeiro

, et al. Buckling and post-buckling of filament wound composite tubes under axial compression: linear, nonlinear, damage and experimental analyses. Compos Part B 2018; 149: 227–239.

29.

Liu

Liao

Jia

, et al. Finite element analysis of dynamic progressive failure of carbon fiber composite laminates under low velocity impact. Compos Struct 2016; 149: 408–422.

30.

Fang

Liang

Wang

B-L

. Progressive damage and nonlinear analysis of 3D four-directional braided composites under unidirectional tension. Compos Struct 2009; 89: 126–133.

31.

Feng

Aymerich

. Finite element modelling of damage induced by low-velocity impact on composite laminates. Compos Struct 2014; 108: 161–171.

32.

Krueger

. Virtual crack closure technique: history, approach, and applications. Appl Mech Rev 2004; 57: 109–143.

33.

Camanho

Davila

de Moura

. Numerical simulation of mixed-mode progressive delamination in composite materials. J Compos Mater 2003; 37: 1415–1438.

34.

Riccio

Di Felice

LaManna

, et al. A global–local numerical model for the prediction of impact induced damage in composite laminates. Appl Compos Mater 2014; 21: 457–466.

35.

Lapczyk

Hurtado

. Progressive damage modeling in fiber-reinforced materials. Compos Part A Appl Sci Manuf 2007; 38: 2333–2341.

36.

Benzeggagh

Kenane

. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos Sci Technol 1996; 56: 439–449.

37.

Yan

, et al. 3D explicit finite element analysis of tensile failure behavior in adhesive-bonded composite single-lap joints. Compos Struct 2018; 201: 261–275.

38.

Schön

. Coefficient of friction of composite delamination surfaces. Wear 2000; 237: 77–89.

39.

Bing

Sun

. Effect of compressive transverse normal stress on mode II fracture toughness in polymeric composites. Int J Fract 2007; 145: 89–97.

40.

Kim

Rim

Lee

, et al. Composite damage model based on continuum damage mechanics and low velocity impact analysis of composite plates. Compos Struct 2013; 95: 123–134.

41.

Guan

Yang

. Low-velocity impact and damage process of composite laminates. J Compos Mater 2002; 36: 851–871.

42.

Singh

Mahajan

. Modeling damage induced plasticity for low velocity impact simulation of three dimensional fiber reinforced composite. Compos Struct 2015; 131: 290–303.