Modelling low velocity impact induced damage in composite laminates

Abstract

The paper presents recent progress on modelling low velocity impact induced damage in fibre reinforced composite laminates. It is important to understand the mechanisms of barely visible impact damage (BVID) and how it affects structural performance. To reduce labour intensive testing, the development of finite element (FE) techniques for simulating impact damage becomes essential and recent effort by the composites research community is reviewed in this work. The FE predicted damage initiation and propagation can be validated by Non Destructive Techniques (NDT) that gives confidence to the developed numerical damage models. A reliable damage simulation can assist the design process to optimise laminate configurations, reduce weight and improve performance of components and structures used in aircraft construction.

Introduction

Composite laminates, which are made from the continuous fibres and a polymeric resin, are used for load-carrying structures in aerospace industry due to their high specific stiffness and strength, and good fatigue and corrosion resistance when exposed to harsh environments. However, laminated constructions are vulnerable to foreign object impact, such as from debris or bird strike during landing or take-off and hail storm when in flight. Tools dropped accidently during assembly or maintenance can also introduce barely visible impact damage (BVID) that happens inside the composite (Abrate 2011, 1998). This could lead to sudden failure due to stiffness and strength degradation, if the introduced damage remains undetected. The impact behaviour of materials becomes very crucial for the safety assessment of aircraft (Shi 2014).

Fibre reinforced composite structures usually demonstrate intra- (such as fibre fracture, matrix cracking and splitting, fibre/matrix debonding) and inter-laminar (delamination) damage modes that are formed, developed and interacted in a complicated pattern internally, which becomes difficult to detect non-destructively (Diaz Valdes and Soutis 2000, 2002). In general, transverse matrix cracking/splitting is the first damage mode occurring within the laminate due to relatively low resin strength/stiffness (Nairn 2000). But this won’t result in ultimate failure of the composite, although degradation of its stiffness may occur (Tong et al. 1997). With cracks propagating in the matrix, delaminations can be initiated between plies, which could eventually lead to fibre breakage and loss of functionality depending on load or energy level available (Shi 2014). This can be a critical design issue and limitation for structural applications those require high damage resistance. Therefore, these modes of damage highlight the importance of investigating and understanding their initiation and evolution in composite laminates with the aim to select lay-up configurations for better damage resistance and tolerance.

Although the impact behaviour of a composite structure can be assessed by experimental testing, a huge consumption of time and cost will be taken because of high skilled labours required for operations and material costs. Therefore, the numerical simulation by means of finite element method (FEM) has been extensively developed to predict the complicated damage modes within composite structures when subjected to dynamic impact loadings, especially at the early design stage when such simulation can minimise the risks prior to implementation of experiments and avoid waste of mechanical tests and manufacturing of components (Shi et al. 2016). A large amount of research on developing predictive damage models have been reported in literature (Abrate 1998; Shi et al. 2012, 2016; Riccio et al. 2014; González et al. 2012; Li et al. 2014a; Schwab et al. 2016; Donadon et al. 2008; Faggiani and Falzon 2010; Feng and Aymerich 2014; Shi and Soutis 2012). Damage in composites can be usually modelled using failure criteria methods and damage mechanics approaches to predict damage evolution under external impact loading. In this paper, some recent FE methods to predict low velocity impact induced damage under high impact mass are reviewed and discussed for fibre reinforced laminates and hybrid composite structures. The accuracy of numerical modes are compared and validated by observing damage extents using non-destructive techniques (NDT).

Review

Damage criteria methods

Failure criteria method

Tow most popular approaches to estimate damages in composites are stress-based failure criteria and damage mechanics methods. The maximum stress/strain criterion (Eqs. 2.1, 2.2, 2.3, 2.4) is the easiest to numerically perform and straightforward to understand for each individual damage mode, but unable to account for mode interaction.

$$ f=\mathit{\max}\ \left(\left|\frac{\sigma_{11}}{X}\right|,\left|\frac{\sigma_{22}}{Y}\right|,\left|\frac{\sigma_{12}}{S_{12}}\right|\right) $$

(2.1)

where

$$ \begin{array}{l}{\sigma}_{11}\ge 0\ \Rightarrow\ X={X}^T;\kern1.5em {\sigma}_{11}<0\ \Rightarrow\ X={X}^C\\ {}\ {\sigma}_{22}\ge 0\ \Rightarrow\ Y={Y}^T;\kern1.5em {\sigma}_{22}<0\ \Rightarrow\ Y={Y}^C\end{array} $$

(2.2)

or

$$ f=\mathit{\max}\ \left(\left|\frac{\varepsilon_{11}}{\varepsilon_{11}^0}\right|,\left|\frac{\varepsilon_{22}}{\varepsilon_{22}^0}\right|,\left|\frac{\gamma_{12}}{\gamma_{12}^0}\right|\right) $$

(2.3)

where

$$ \begin{array}{l}{\varepsilon}_{11}\ge 0\Rightarrow {\varepsilon}_{11}={\varepsilon}_{11}^{0\mathrm{T}};\kern2em {\varepsilon}_{11}<0\Rightarrow {\varepsilon}_{11}={\varepsilon}_{11}^{0\mathrm{C}}\\ {}{\varepsilon}_{22}\ge 0\Rightarrow {\varepsilon}_{22}={\varepsilon}_{22}^{0\mathrm{T}};\kern1.3em {\varepsilon}_{22}<0\Rightarrow {\varepsilon}_{22}={\varepsilon}_{22}^{0\mathrm{C}}\end{array} $$

(2.4)

where, σ ₁₁, σ ₂₂, ε ₁₁ and ε ₂₂ are the tensile or compressive stress and strain in the axial (1) and transverse (2) directions, respectively. σ ₁₂ and γ ₁₂ are the in-plane shear stress and strain. X ^T, X ^C, Y ^T, and Y ^C represent tensile and compressive strength for failure prediction in their respective directions. S ₁₂ and ${\gamma}_{12}^0$ denote the in-plane shear strength and failure strain, respectively.

To predict the interactive effect on various damage modes, the stress/strain based failure envelop was developed by Tsai-Hill criterion (Tsai 1965) (Eqs. 2.5–2.6) and Tsai-Wu (Tsai and Wu 1971) criterion (Eqs. 2.7–2.8).

$$ \begin{array}{l}f={F}_{11}{\sigma}_{11}^2+{F}_{22}{\sigma}_{22}^2+{F}_{66}{\sigma}_{12}^2+2{F}_{12}{\sigma}_{11}{\sigma}_{22}+{F}_1{\sigma}_{11}+{F}_2{\sigma}_{22}+{F}_6{\sigma}_{12}\\ {}{F}_{11}=\frac{1}{X^2}, {F}_{22}=\frac{1}{Y^2}, {F}_{12}=-\frac{1}{2{X}^2}, \kern0.5em {F}_{66}=\frac{1}{S_{12}^2}\end{array} $$

(2.5)

$$ {F}_1=0,\kern1.25em {F}_2=0,\kern0.75em {F}_6=0 $$

(2.6)

where X, Y and S ₁₂ in Eq. (2.6) have the same definition as in Eq. (2.1).

$$ {F}_{11}=\frac{1}{X^T{X}^C},\kern2.5em {F}_{22}=\frac{1}{Y^T{Y}^C},\kern2em {F}_{66}=\frac{1}{S_{12}^2} $$

$$ {F}_1=\frac{1}{X^T}-\frac{1}{X^C},\kern1.25em {F}_2=\frac{1}{Y^T}-\frac{1}{Y^C},\kern0.75em {F}_6=0 $$

(2.7)

with the coefficient F₁₂ defined as:

$$ {F}_{12}=\frac{F_{12}^{\ast }}{\sqrt{X^T{X}^C{Y}^T{Y}^C}} $$

(2.8)

${F}_{12}^{\ast }$ in the above equation ranges from −1 to 1, and can be obtained by fitting the equibiaxial experimental data. Hoffman (Hoffman 1967) derived the same coefficients as in Eq. (2.7) but defined F ₁₂ as

$$ {F}_{12}=-\frac{1}{2{X}^T{X}^C} $$

(2.9)

Farooq et al. (Farooq and Myler 2016) attempted to use the Tsai-Hill and Tsai-Wu criteria to predict the ply level failure in their research and compared to the experimental results. However, the Tsai-Hill criterion underestimated the maximum loads compared to other failure theories. This might be due to its transverse compressive strength that is far larger than transverse tensile strength and therefore the Tsai-Hill criterion failed to predict ply based level failure. The Tsai-Wu criterion showed similar scenario that ply-by-ply failure index quantities never exceeded one, which indicated no failure occurred. The most likely reason might be due to these polynomial based criteria could not use through-thickness stresses to predict ply failure. Therefore, it can be seen the polynomial failure criteria are not ideal to model composite failure and a ply-by-ply progressive damage method is thus more popular in developing the predictive damage model.

Ply based progressive damage method

Damage initiation

Ply based individual damage mode can be generally modelled for two stages: initiation and evolution. Hashin criteria have been widely used in academic research and industry (Eqs. (2.10)–(2.13)).

Fibre tension (σ ₁₁ ≥ 0):

$$ {F}_{ft}={\left(\frac{\sigma_{11}}{X^T}\right)}^2+\kappa {\left(\frac{\sigma_{12}}{S_{12}}\right)}^2=1 $$

(2.10)

Fibre compression (σ ₁₁ < 0):

$$ {F}_{fc}={\left(\frac{\sigma_{11}}{X^C}\right)}^2=1 $$

(2.11)

Matrix tension (σ ₂₂ ≥ 0):

$$ {F}_{mt}={\left(\frac{\sigma_{22}}{Y^T}\right)}^2+{\left(\frac{\sigma_{12}}{S_{12}}\right)}^2=1 $$

(2.12)

Matrix compression (σ ₂₂ < 0):

$$ {F}_{mc}={\left(\frac{\sigma_{22}}{2{S}_{23}}\right)}^2+\left[{\left(\frac{Y^C}{2{S}_{23}}\right)}^2-1\right]\frac{\sigma_{22}}{Y^C}+{\left(\frac{\sigma_{12}}{2{S}_{12}}\right)}^2=1 $$

(2.13)

In the above equations, σ _ij (i, j = 1, 2, 3) are the stress components defined in the material coordinate system. X ^T and X ^C denote the fibre tensile and compressive strengths, Y ^T and Y ^C are the transverse tensile and compressive strength, S _{i, j} (i, j = 1, 2, 3) denote the longitudinal and transverse shear strengths of the composite, respectively. The coefficient κ in Eq. (2.10) accounts for the contribution of shear stress to fibre tensile failure and generally ranges between 0 and 1. Hashin criteria comprehensively considers the various damage modes such as fibre ruptures in tension, matrix cracking due to transverse tension and shear, fibre compressive damage (buckling or kinking) and matrix crushing under transverse compression and shear effect (Hoffman 1967). Hashin criteria have been proved that it is effective method to predict ply based damage initiation (Shi et al. 2012; Feng and Aymerich 2014; Farooq and Myler 2016; Johnson et al. 2006; Liu et al. 2016; Li et al. 2014b), however, based on the experimental observation, the matrix compressive crack is always formed with an angle through the ply thickness whereas Hashin criteria are less able to simulate this process. Puck and Schurmann (1998) developed a damage model for transverse compression. They proposed to use the failure criteria of Mohr (Salencon 2001) instead of the yield criterion of von Mises which is normally applied. Puck’s damage criterion for compression damage mode can be expressed as Eq. (2.14).

$$ {F}_{mc}={\left(\frac{\sigma_{NT}}{S_{23}^A+{\mu}_{NT}{\sigma}_{NN}}\right)}^2+{\left(\frac{\sigma_{NL}}{S_{12}+{\mu}_{NL}{\sigma}_{NN}}\right)}^2=1 $$

(2.14)

In Eq. (2.14) σ _ij (i, j = L, T, N) are the stresses σ _ij (i, j = 1, 2, 3) rotated to the fracture plane, by reference to the axes shown in Fig. 1 (Shi 2014):

$$ {\sigma}_{NN}={\sigma}_2{m}^2+{\sigma}_3\left(1-{m}^2\right)+2{\tau}_{23}mn $$

(2.15)

$$ {\sigma}_{NT}=-{\sigma}_2mn+{\sigma}_3mn+{\tau}_{23}\left(2{m}^2-1\right) $$

(2.16)

$$ {\sigma}_{TT}={\sigma}_2\left(1-{m}^2\right)+{\sigma}_3{m}^2-2{\tau}_{23}mn $$

(2.17)

$$ {\sigma}_{NL}={\tau}_{12}m+{\tau}_{13}n $$

(2.18)

$$ {\sigma}_{LT}=-{\tau}_{12}n+{\tau}_{13}m $$

(2.19)

$$ {S}_{23}^A=\frac{Y_C}{2}\left(\frac{1- \sin \varphi }{cos\varphi}\right) $$

(2.20)

$$ \varphi =2\alpha -{90}^{\mathrm{o}} $$

(2.21)

where m = cos(α) and n = sin(α) in Eqs. (2.15)–(2.19). ${S}_{23}^A$ is the transverse shear strength along the fracture plane, which can be determined by the transverse compression strength Y _C and the angle of the fracture plane as shown in Eqs. (2.20) and (2.21).

The friction coefficients μ _NT and μ _NL in Eq. (2.14) can be defined based on the material friction angle (see Eq. (2.21)), φ, and material properties by reference to the Mohr failure criteria.

$$ {\mu}_{NT}= tan\varphi = \tan \left(2\alpha -{90}^{\mathrm{o}}\right) $$

(2.23)

$$ {\mu}_{NL}=\frac{\mu_{NT}}{S_{23}^A}{S}_{12} $$

(2.24)

In general, the fracture plane is oriented at α = 53° through the thickness direction under a uniaxial compressive load (Puck and Schurmann 1998). However, impact loading usually leads to various values of the fracture angle which can be numerically determined by executing the damage initiation index in a swept angle range. This method was performed by Shi et al. (2012), Faggiani and Falzon (2010) and Feng and Aymerich (2014) where the fracture plane and increase of shear strength due to normal compressive stress acted on fracture plane have been successfully simulated.

Besides, Sun et al. (1996) also appropriately modified Hashin’s criteria to improve the accuracy of modelling matrix compressive damage mode:

$$ {F}_{mc}={\left(\frac{\sigma_{22}}{Y^C}\right)}^2+{\left(\frac{\sigma_{12}}{S_{12}-{\varsigma \sigma}_{22}}\right)}^2=1 $$

(2.25)

where ς is a constant determined experimentally and generally regarded as an internal material friction parameter.

Camanho et al. (2005) developed a failure criterion called LaRC03 based on continuum damage mechanics (CDM) method. In order to demonstrate their damage model, the failure envelope plotted between the transverse stress σ ₂₂ and in-plane shear stress τ ₁₂ was generated and compared with other criteria. All predicted results have been validated by World Wide Failure Exercise (WWFE) test results. For the tensile damage mode, it rarely find big difference for simulation accuracy by proposed criteria, except the maximum stress criterion because the maximum stress criterion only defines the single failure without interaction modelled. Therefore, it cannot give a satisfactory prediction of the failure of composites, especially when damage is matrix dominated.

For matrix compression damage mode, it was illustrated as similar conclusion from experimental observation. Due to the fracture plane existed, Puck’s envelope offered the most accurate prediction while both Sun’s criterion (Sun et al. 1996) and LaRC03 (Davila et al. 2005) also demonstrate an appropriately predictive capacity. Hashin’s criteria showed an underestimation due to fracture plane was ignored. Topac et al. (2017) employed the damage criteria LaRC04 developed by Pinho et al. (2005) to further improve the accuracy of modelling matrix cracking using such continuum damage mechanics based failure criteria where the transverse cracks through the thickness have been accurately simulated compared to experimental measurement under the relatively fine meshing strategy.

Damage evolution

To predict the damage evolution within composites, a material degradation strategy should be generally defined for each individual damage mode. The most common way is to apply a degradation parameter corresponding to the damage mode for simulating the softening effect. The parameter for degradation of material might be from experimental measurement. Tita et al. (2008) reduced the stiffness by using appropriate factors with respect to the various failure modes observed experimentally. For instance, the transverse Young’s modulus E₂₂ and the in-plane Poisson’s ratio ν₁₂ were reduced to zero directly to represent the complete damage in their work. Farooq et al. also defined the degraded modulus for composite IM7/8552 with complete failure of each single damage mode such as the tensile moduli of E₁, E₂ and E₃ as well as in-plane and out-of-plane shear moduli to simulate the impact induced damage within the laminate (Farooq and Myler, 2015, 2016). In order to find the accurate degradation parameters for tensile and compressive moduli of the composite, Louca et al. performed cyclic tests and eventually concluded the average values for stiffness loss once the composite had failed (Johnson et al. 2006). It can be seen that this method is straightforward and always accurate for predicting the composite moduli degradation, however, experimental measurements have to be necessary for different damage modes and different kinds of fibre and resin systems. Moreover, this method cannot numerically demonstrate a progressive damage evolution, which could be a constraint to composite engineers, since it relies heavily on experiments.

A progressive damage evolution based on fracture mechanics was reported by Chang and Chang (1987), Chang et al. (1991). Considering the matrix failure, they proposed a degradation law to reduce the moduli E ₁₁ and G ₁₂ based on an exponential decaying, but other moduli were reduced abruptly to zero once the damage initiation happened.

$$ \begin{array}{l}\ {E}_{11}^d={E}_{11}\ \mathit{\exp}\left[-\left(\frac{A}{A_o}\right)H\right]\\ {}{G}_{12}^d={G}_{12}\ \mathit{\exp}\left[-\left(\frac{A}{A_o}\right)H\right]\end{array} $$

(2.26)

where A wass the area of the damage zone; A _o was the area of the interaction zone of the fibre failure from Chang and Chang (1987), Chang et al. (1991). H was a factor to control the degradation of the material stiffness. Following this work, Matzenmiller et al. (1995) developed a damage model called ‘MLT model’ for the non-linear analysis of the composite laminates. They constructed the model using damage variables with respect to the individual failure modes in the material principal directions. The model assumes that each unidirectional lamina in the composite acted as a continuum irrespective of the damage state. The damage growth was controlled based on a Weibull distribution. The post-damage softening behaviour of the composite can be predicted by an exponential function:

$$ d=1-\mathit{\exp}\left(-\frac{1}{me}{\left(\frac{E^0\varepsilon }{X}\right)}^m\right) $$

(2.27)

where E ⁰ was the individually fibre or transverse modulus, ε was the strain related to the progressive damage at different time steps, “e” was a mathematical constant. X was the tensile or compressive strength with regard to the different damage modes in the different loading directions. “m” was the strain softening parameter during damage progression, which is a key factor to determine the accuracy of prediction eventually. In general, a high value of “m” could result in brittle failure of the material while a low value of “m” indicated a ductile failure response that could lead to a high absorbed energy due to damage.

It is an effective approach to model the damage growth for a composite laminate using strain softening parameter “m” in MLT model (Williams and Vaziri 2001; Gama and Gillespie 2011; Tabiei and Aminjikarai 2009; Jung et al. 2017). The appropriate value of “m” was usually related to the mesh size and load conditions. The value of “m” for various damage modes can be determined using uniaxial tensile or compressive tests. A set of values of “m” can be applied for the strain softening to model the complex failure process of composites for different damage modes. Obviously, the “m” value has a strong effect on determining the accuracy of prediction of damage progression. An inappropriate value of “m” could give rise to numerical difficulties to simulate the process of damage growth and eventually lead to the unrealistic results (Williams and Vaziri 2001). For example, if a relatively small value of “m” was applied in the damage model for composite laminate, it will exhibit a ductile behaviour which is obviously contradictory to its brittle characteristics. Therefore, in order to avoid experiment dependent “m”, the damage variable was proposed by an exponential function that only referred to characteristic length and composite properties (Dassault Systemes Simulia Corp 2008).

$$ {d}_{\alpha }=1-\frac{1}{r_{\alpha }}\mathit{\exp}\left[-{A}_{\alpha}\left({r}_{\alpha }-1\right)\right]\kern2em and\kern1.25em {d}_{\alpha}\ge 0 $$

where the coefficient A_α was defined as

$$ {A}_{\alpha }=\frac{2{g}_0^{\alpha }{L}_c}{G_{IC}^{\alpha }-{g}_0^{\alpha }{L}_c}\kern1.25em with\kern2em {g}_0^{\alpha }=\frac{X_{\alpha}^2}{2{E}_{\alpha }} $$

(2.29)

All parameters used in this exponential function correspond to material properties and therefore the experiment for evaluating “m” can be effectively avoided. Schwab et al. employed this method to successfully predict the impact induced damage with perforation for woven composites under intermediate velocity impact (Schwab et al., 2015, 2016; Schwab and Pettermann 2016). More specifically, using this model, they simulated the impact behaviour of a large composite fan containment of a jet engine impacted by deformable bodies (Schwab and Pettermann 2016).

In addition, the energy based damage mechanics approach was also extensively developed and used to model the progressive failure in composite laminates (Iannucci and Willows 2006, 2007; Iannucci and Ankersen 2006). The damage model used the strains at damage onset and at complete damage, effectively defined the damage variable for degradation and accurately captured the damage progression for composite laminates (Shi et al. 2012; Donadon et al. 2008; Faggiani and Falzon 2010; Feng and Aymerich 2014; Liu et al. 2016; Topac et al. 2017), which a typical damage variable defined for tensile damage mode is shown in Eq. (2.30).

$$ {d}_{1,2}^t=\frac{\varepsilon_{1,2}^{ft}}{\varepsilon_{1,2}^{ft}-{\varepsilon}_{1,2}^{0t}}\left(1-\frac{\varepsilon_{1,2}^{0t}}{\varepsilon_{1,2}}\right) $$

(2.30)

where the subscript 1 and 2 denoted the fibre and transverse directions, respectively; ${\upvarepsilon}_{1,2}^{0\mathrm{t}}$ was the strain when the damage initiation condition was fulfilled. Due to the irreversibility of the damage variable, the strain calculated at each time step could be updated in comparison with the strain at damage initiation ${\varepsilon}_{1,2}=\mathit{\max}\left({\varepsilon}_{1,2},{\varepsilon}_{1,2}^{0t}\right)$ in Eq. (2.30). In order to avoid a zero or even negative energy absorption, the complete failure strain could be also defined to be greater than the initial failure strain ${\varepsilon}_{1,2}^{ft}>{\varepsilon}_{1,2}^{0t}$.

The failure initiation strain is given by the following equation:

$$ {\varepsilon}_{1,2}^{0t}=\frac{\sigma_{1,2}^T}{E_{1,2}} $$

(2.31)

For tensile failure in fibres, σ ^T was the tensile strength X ^T while Y ^T was used for the matrix tensile failure mode.

ε ^ft can be derived from the fracture energy ${G}_{1,2}^T$ for the individual failure mode, the failure strength of the material and the characteristic length:

$$ {\varepsilon}_{1,2}^{ft}=\frac{2{G}_{1,2}^T}{\sigma^T{l}^{\ast }} $$

(2.32)

In Eq. (2.32), l ^* is the characteristic length which can maintain an energy release rate per unit area of rack constant and also keep the predicted results independent of the mesh size in a FE model. The definition of characteristic length can be found for more details in (Shi 2014; Bazǎnt and Oh 1983; Olivier 1989; Pinho 2005).

Shi et al. identified the most effective damage evolution laws for individual damage mode and developed the damage model including non-linear shear with damage for prediction of low velocity impact induced damage of reinforced carbon composite laminates (Shi et al. 2012). The damage model was implemented into commercial FE code ABAQUS/EXPLICIT by user defined subroutine VUMAT. The various impact energies were tested to explore the multi-damage modes of the composite, which was experimentally recorded by X-ray radiography. The damage of the laminate compared well with X-ray images taking at different impact loads (see Fig. 2). However, in Fig. 3 it can be seen that there is a difference between experimental measurement and numerical prediction, since splitting was not simulated but observed in the X-ray radiographs. The solution to this issue was proposed by Shi et al. (2014a, 2014b), Shi and Soutis (2015) using cohesive zone elements, which will be discussed in the following section.

Fig. 1

Fracture plane for matrix compressive failure relative to the material coordinate system

Fig. 2

X-Ray radiographs and numerically predicted delamination area for cross-ply [0/90]_2S: a-b Impact energy of 7.35 J; c-d 11.03 J; e-f 14.7 J (Shi 2014)

Fig. 3

Experimental and numerical impact energy-time histories for cross-ply [0/90]_2S: a Impact energy of 7.35 J; b 11.03 J; c 14.7 J (Shi 2014)

Cohesive zone elements (CZE) method

The cohesive zone elements are found to be an effective way to capture delamination initiation and propagation at the interface of adjacent plies with definition of initiation and propagation laws for damage (Shi et al. 2012; Camanho and Matthews 1999; Dávila and Johnson 1993; Camanho and Dávila 2002). Delamination can be formed under the effect of mixed-mode including all of the damage modes (see Fig. 4) when subjected to impact loading.

Fig. 4

Three damage modes a mode I (opening) b mode II (sliding) c mode III (tearing) (Shi 2014)

A quadratic stress failure criterion can be used for the estimation of delamination onset:

$$ {\left(\frac{\sigma_n}{N}\right)}^2+{\left(\frac{\sigma_s}{S}\right)}^2+{\left(\frac{\sigma_t}{T}\right)}^2=1 $$

(3.1)

where σ _i (i = n, s, t) denoted the traction stress vectors in the normal n and shear directions, s and t, respectively, while N, S and T were defined as the corresponding inter-laminar normal and two shear strengths.

The damage evolution can be defined with a similar way represented previously in Eq. (3.2) as:

$$ d=\frac{\delta_m^f\left({\delta}_m-{\delta}_m^0\right)}{\delta_m\left({\delta}_m^f-{\delta}_m^0\right)} $$

(3.2)

where δ _m refers to the value of the mixed-mode displacement attained during the loading history, as illustrated in Fig. 5. The δ _m parameter corresponds to the total mixed-mode displacement (normal, sliding, tearing) given by:

$$ {\delta}_m=\sqrt{\delta_n^2+{\delta}_{shear}^2}=\sqrt{\delta_n^2+{\delta}_s^2+{\delta}_t^2} $$

(3.3)

Fig. 5

Typical traction-displacement relationship (Shi 2014)

In Eq. (3.2) ${\delta}_m^f$ is the mixed-mode displacement at complete failure and ${\delta}_m^0$ is the effective displacement at damage initiation. The displacement ${\delta}_m^f$ can be defined based on fracture energy represented by several reported approaches such as power law (Eq. 3.4) or the Benzeggagh–Kenane (BK) criterion (Eq. 3.5) which can offer the specific high efficiency for composite materials with the use of identical critical fracture energies along the two shear directions (ABAQUS 2010).

$$ {\delta}_m^f=\left\{\begin{array}{c}\hfill \frac{2\left(1+{\beta}^2\right)}{K{\delta}_m^0}{\left[{\left(\frac{1}{G_{IC}}\right)}^{\alpha_p}+{\left(\frac{\beta^2}{G_{IIC}}\right)}^{\alpha_p}\right]}^{-1/{\alpha}_p}\kern1em {\delta}_n>0\hfill \\ {}\hfill \sqrt{{\left({\delta}_s^f\right)}^2+{\left({\delta}_t^f\right)}^2}\kern9.3em {\delta}_n\le 0\hfill \end{array}\right. $$

(3.4)

$$ {\delta}_m^f=\left\{\begin{array}{c}\hfill \frac{2}{K{\delta}_m^0}\left[{G}_{IC}+\left({G}_{IIC}-{G}_{IC}\right){\xi}^{\eta}\right]\kern3em {\delta}_n>0\hfill \\ {}\hfill \sqrt{{\left({\delta}_s^f\right)}^2+{\left({\delta}_t^f\right)}^2}\kern8em {\delta}_n\le 0\hfill \end{array}\right. $$

(3.5)

where α _p was power law parameter and β represented the mixity ratio. η in Eq. (3.5) was the B-K power law parameter that can be determined using a least-square fit from a set of mixed-mode bending experimental data. Obviously, the damage evolution law was strongly dependent on the fracture energy for mix-mode loading. Damage was assumed to grow when the energy release rate was equal or greater than the critical energy release rate, Eq. (3.6).

$$ G={G}_{IC}+\left({G}_{IIC}-{G}_{IC}\right){\left(\frac{G_{shear}}{G}\right)}^{\eta } $$

(3.6)

where G = G _I  + G _shear and G _shear  = G _II  + G _III in Eq. (3.6) are the energy release rate for mixed-mode and shear loading, respectively. $\frac{G_{shear}}{G}$ can be expressed as $\xi =\frac{\beta^2}{1+{\beta}^2}$, which is a function of mixed ratio, β. ξ generally took values between zero and one. When ξ = 0 the crack was mode I driven, while as ξ → 1 propagation was mode II dominated (and this could be also the case when η = 0).

Compared Eqs. (3.1) and (3.2) with the interactive quadratic failure criterion (as Eq. (3.7) (Pinho et al. 2006a, 2006b)) and Eq. (2.30), the constitutive equations applied were found to be almost consistent for predicting the damage initiation and evolution process. Therefore, Shi et al. (2014a, 2014b), Shi and Soutis (2015) proposed to simulate matrix cracking/splitting by inserting cohesive elements within each lamina, at space intervals determined from experimental observations.

$$ {F}_{mt}={\left(\frac{\sigma_{22}}{Y^T}\right)}^2+{\left(\frac{\sigma_{12}}{S_{12}}\right)}^2+{\left(\frac{\sigma_{23}}{S_{23}}\right)}^2=1 $$

(3.7)

In Shi’s work, cohesive zone elements were used to simulate matrix cracking/splitting as well as interfacial delamination while the damage criteria were used to predict the rest of failure modes. The predicted matrix crack density was well predicted by Shi, comparing to the analytical and experimental results under uniaxial tensile load (Shi et al. 2014b). To simulate the four point bending test, the crack density modelled in 90 plies of laminates showed a good agreement with experimental measurement, in particular, the partial and full matrix cracks and their propagation process have been successfully captured by cohesive zone elements (Shi and Soutis 2015). To improve the predicted accuracy for low velocity impact behaviour of composite laminates (Shi et al. 2012), the splitting in bottom 0 ply was accurately simulated by inserting cohesive zone elements within individual ply of the model, which was observed in experiment (Shi et al. 2014a; Shi and Soutis 2015), as shown in Fig. 6. In addition, the absorbed energy for each energy level was better captured when matrix cracking and splitting were included in the calculations. In particular, for the 7.35 J impact case, the difference is reduced from almost 14% to less than 5.2% (Shi et al. 2014a).

Fig. 6

Experimentally observed and numerically captured damage under impact energy of 14.7 J for cross-ply [0/90]_2S a X-ray radiograph and b proposed model (Shi 2014)

Moreover, the proposed method by Shi et al. (2014a) offered an outstanding benefit to predict composite damage that cracks may close once the load was removed, becoming undetectable by X-ray radiography or any other non-destructive detection technique and hence underestimating the severity and extent of internal damage, as shown in Fig. 7.

Fig. 7

Initial form and propagation of matrix cracking in composite laminate subjected to impact energy of 7.35 J for cross-ply [0/90]_2S. a Initial contact b Maximum deflection c Rebounding d End of impact event (Shi 2014)

In this example, the impact energy of 7.35 J was selected to simulate the damage evolution process in the cross-ply laminate since 90° matrix cracking was the main damage mode expected with crack induced delamination (the impact energy is too low to cause fibre breakage). Matrix cracking was initially formed in the bottom 90° ply with few cracks found in the middle 90° plies, in Fig. 7a. As the impactor further contacted the composite laminate the surface indentation was formed with a bending crack, appearing on the top 0° surface ply due to concentrated compressive load, while complete cracks through the thickness of the 90° layers were developed accompanied by delamination at the 0/90 interfaces, in Fig. 7b. As the impact event progressed, the growth rate of cracks in the middle and bottom 90° plies was found to decrease, while the projectile rebound occurred, in Fig. 7c. The simulation showed that cracks started to close once the complete rebound was reached, in Fig. 7d with extensive delamination that matched the detected result by X-ray radiography.

Sun et al. (2016) also employed cohesive zone elements to predict the quasi-static indentation tests for carbon fibre composites. The inserted cohesive zone elements, especially in angled ply, well captured the matrix cracking and splitting.

Modelling impact induced damage for hybrid composites

In order to further improve mechanical properties of composites such as energy absorption, damage resistance, fatigue properties and further reduced weight, the hybrid composite was developed. In this section honeycomb sandwich structures and fibre metal laminates are discussed and their impact behaviour is simulated by the FE method.

Honeycomb sandwich structure

In general, the FE process was more complicated for such hybrid composites as more damage modes were involved to simulate. The basic strategy to model the individual damage mechanisms was to define the damage model for different components, respectively. To simulate the low velocity impact response of composite sandwich structure, the continuum damage mechanics based damage criteria were defined for predicting the intra-laminar damage of composite face sheet with delamination using cohesive zone elements while the mechanical behaviour of honeycomb core can be modelled, referring to the experimentally measured data (Aktay et al. 2008; Liu et al. 2017) or damage model for Nomex (Chen et al. 2017; Liu et al. 2015; Giglio et al. 2012) or Aluminium honeycomb (Lee et al. 2010; Kay 2003).

Chen et al. successfully predicted the perforation process of a composite sandwich structure where the damage of composite face skin and Nomex honeycomb core were numerically captured compared to experimental observation (Chen et al. 2017). Liu et al. (2017) applied the composite damage model and elastic–plastic material model for filled aluminium honeycomb core to predict the crushing and bending responses where the failure process for both responses were well simulated when compared to experimental measurements.

Fibre metal laminates

Fibre metal laminates (FMLs) are hybrid composites consisting of thin metal sheets and fibre reinforced composite layers, which the mechanical properties such as impact resistance, and fatigue properties can be further improved. As developed FMLs the commercially manufactured Glass fibre aluminium reinforced (Glare) laminate has been used for the upper fuselage of the Airbus A380, which consists of thin aluminium alloy layers (2024, 7075 etc.) and reinforced glass fibre plies (S2-glass/FM94). Glare usually exhibits a high compressive strength and damage resistance, especially under impact loading. Under low velocity impact most of the energy is dissipated by plastic deformation of the metal layers, followed by delamination of neighbouring layers. Glare thus shows improved impact behaviour when compared to monolithic aluminium but is far better than carbon fibre composites (Vlot and Gunnink 2001; Wu and Yang 2005). To model the impact behaviour of FMLs, the multi-damage modes could be simulated for various material components: a) For metal sheet, the plastic damage model like the Johnson-Cook Kay (2003) or experimentally measured data can be defined; b) For the fibre reinforced composite layer we can employ the damage criteria that were discussed in section of 2 of this work; c) Interfacial damage such as the debonding between metal and composite or delamination at interfaces of the composite plies can be modelled by introducing cohesive zone elements (Moriniere et al. 2014; Mohamed et al. 2012; Bikakis and Savaidis 2016; Zarei et al. 2017). Based on the observation of experiment (Leonard et al. 2014), composite damage like splitting might promote deflection, even crack of their neighboured aluminium, and therefore the cohesive zone element discussed in section 3 of this paper could help to predict the interaction of damages between composite and metal layer in FMLs for improved accuracy.

Conclusions

This paper reviewed the techniques developed recently for modelling the impact induced damage of composite laminates. Damage criteria were discussed for initiation and evolution of individual damage modes. The developed damage models such as Hashin and Puck were compared and discussed for the accuracy of prediction while the damage evolution laws were reviewed based on continuum damage mechanics approach. The experimentally based parameter degradation of damage composite was easily performed whereas the progressive damage evolution was recommended to appropriately capture the damage process and eventual damage extent under low velocity impact loadings.

The cohesive zone elements method proposed by Shi et al. was shown to be an effective way to predict matrix cracking/splitting which was commonly ignored by damage criteria method that could underestimate energy absorbed. Moreover in order to understand the crack evolution process in composites the impact test with low velocity impact energy level was numerically studied in more detail, since matrix cracking can lead to delamination. The model showed regions in the ply where resin cracks initiated, propagated and then closed during projectile rebounding, something that is difficult to detect by X-ray radiography or any other non-destructive detection technique, hence underestimating the severity and extent of internal damage. This study can be extended to estimate the residual strength properties and fatigue life of such laminates, something that could be investigated in future work.

For hybrid composites such as a honeycomb sandwich structure or fibre metal laminates, the strategy to predict the multi-damage modes is to define the individual damage model for composite, metal or interfacial damage (adhesive), respectively. It will be a great challenge to capture accurately the whole damage process accounting for mode interaction, while the cohesive zone elements approach could potentially achieve more accurate damage predictions.