# Computational methods for predicting impact damage in composite structures

Alastair F. Johnson, Anthony K. Pickett, Patrick Rozycki

► To cite this version:

Alastair F. Johnson, Anthony K. Pickett, Patrick Rozycki. Computational methods for predicting impact damage in composite structures. Composites Science and Technology, 2001, 61 (15), pp.2183-2192. 10.1016/S0266-3538(01)00111-7. hal-01005986

HAL Id: hal-01005986

https://hal.science/hal-01005986

Submitted on 12 Mar 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Computational methods for predicting impact damage
in composite structures

A.F. Johnson<sup>a</sup>, A.K. Pickett<sup>b</sup>, P. Rozycki<sup>c</sup>

<sup>a</sup>German Aerospace Center (DLR), Institute of Structures and Design, D-70569 Stuttgart, Germany

<sup>b</sup>Engineering Systems International GmbH, Frankfurter Str. 13-15, D-65760 Eschborn, Germany

<sup>c</sup>LAMIH, Département Génie Mécanique, Université de Valenciennes, Le Mont Houy, F-59300 Valenciennes Cedex 9, France

This paper describes recent progress in materials modelling and numerical simulation of the impact response of fibre-reinforced
composite structures. A continuum damage-mechanics (CDM) model for fabric-reinforced composites is developed as a framework
within which both in-ply and delamination failure may be modelled during impact loading. Damage-development equations are
derived and appropriate materials parameters determined from experiments. The CDM model for in-plane failure has been imple-
mented in a commercial explicit finite element (FE) code, and new techniques are used to model the laminate as a stack of shell
elements tied by contact interface conditions. This approach allows the interlaminar layers to be modelled and strength reduction
due to delamination to be represented; it also provides a computationally efficient method for the analysis of large-scale structural
parts. The code is applied to predict the response of carbon-fabric-reinforced epoxy plates impacted at different velocities by a steel
impactor. A comparison of structural response and failure modes from numerical simulations and impact tests is given which shows
a good agreement for the prediction of delamination damage at low impact energies and fracture and penetration at higher impact
energies.

Keywords: A. Polymer matrix composites; B. Impact behaviour; C. Computational simulation; C. Damage mechanics; C. Delamination

# 1. Introduction

Composite materials are now being used in primary
aircraft structures, particularly in helicopters, light air-
craft, commuter planes and sailplanes, because of numer-
ous advantages including low weight, high monotonic and
fatigue strengths and the possibility of manufacturing
large integral shell structures. Materials such as carbon-
fibre/epoxy are inherently brittle and composite struc-
tures are vulnerable to impact damage. Consequently
their use in commercial aircraft requires satisfaction of
certification procedures that cover, for example, the
high-velocity impact from runway debris or bird strike.
Impact of composite structures has been studied within
a recently completed CEC-funded research project on
'High-velocity impact of composite aircraft structures'
HICAS [1]. This paper is concerned with one aspect of

the project; namely, the development and validation of
FE codes for modelling the response of composite struc-
tures under impact loads. In polymer composite materials
there are several different failure modes such as matrix
cracking, ply delamination, fibre fracture. Furthermore
the material properties may be strain-rate dependent. Key
issues are the development of suitable constitutive laws for
the mechanical behaviour of composites to failure, the
measurement of appropriate materials parameters and
the implementation of the materials models in to com-
mercial FE codes.

Interest is centred here on *fabric-reinforced* composite
laminates, which are of considerable practical interest
for impact-resistant aircraft structures. During the high-
velocity impact loading of composite structures, failure
may occur by delamination, which is important in lower
energy impacts and in failure initiation, and by in-plane ply
failure which controls ultimate fracture and penetration
in the structure. The paper makes use of continuum-
damage mechanics (CDM) for composites as developed
by Ladevèze and his co-workers [2,3] as a framework
within which in-ply and delamination failure may be modelled. A CDM model for fabric reinforced composite plies under in-plane loads is presented, based on extending methods originally developed for unidirectional (UD) ply materials in [2]. Damage-development equations for the fabric ply are introduced which relate the damage parameters to strain energy release rates in the ply. Tension, compression and cyclic shear tests have been carried out on carbon- and glass-fabric-reinforced epoxy materials and the required mechanical and damage parameters measured for particular forms of the damage evolution equations. Delamination models for interply failure are obtained by applying the CDM framework to the ply interface, as described in [3]. Failure at the interface is modelled as in [4] by degrading stresses using an interface damage law once a critical value of strain is reached. The form of the damage law is chosen such that the total energy absorbed during delamination is equal to the interface fracture energy $G_C$.

An important requirement for the aircraft designer is to develop numerical methods for the simulation of impact in composite structures, thus the composite materials models need to be implemented and validated in FE codes. Emphasis is placed in the paper on the implementation of the fabric ply CDM and delamination models into the commercial explicit FE crash and impact code PAM-CRASH [5]. The ply damage model has been implemented in shell and layered composite shell elements and a novel numerical approach for delamination modelling is developed using stacked shell elements with a contact interface condition.

The paper describes the validation of the fabric model with single elements and test specimens using the measured materials parameters. The delamination model has been calibrated against double cantilever beam (DCB) and mixed mode delamination tests. Finally, a component test case is presented for the impact loading of a simply supported square carbon fabric/epoxy flat plate by a steel spherical impactor. Test conditions and impact energies have been chosen to generate both delamination failure (low impact energies) and significant fibre damage and penetration (high impact energies). An important feature of the composites model is that it correctly distinguishes between the different failure modes in the structure, and gives good quantitative agreement with measured peak loads and plate deformations.

## 2. Fabric reinforced composites damage model

### 2.1. Elastic ply damage mechanics model

The fabric reinforced composite ply is modelled as a homogeneous orthotropic elastic or elastic-plastic damaging material whose properties are degraded on loading by microcracking prior to ultimate failure. A

CDM formulation is used in which ply degradation parameters are internal state variables which are governed by damage evolution equations. Constitutive laws for orthotropic elastic materials with internal damage parameters are described in [2] and [6], and take the general form

$$ \varepsilon^e = S\sigma \quad (1) $$

where $\sigma$ and $\varepsilon^e$ are vectors of stress and elastic strain and S the elastic compliance matrix. For shell elements a plane stress formulation with orthotropic symmetry axes ($x_1, x_2$) is required. The in-plane stress and strain components are

$$ \sigma = (\sigma_{11}, \sigma_{22}, \sigma_{12})^T \quad \varepsilon^e = (\varepsilon_{11}^e, \varepsilon_{22}^e, 2\varepsilon_{12}^e)^T. \quad (2) $$

Using a strain equivalent damage mechanics formulation, the elastic compliance matrix S may then be written:

$$ S = \begin{pmatrix} 1/E_1(1-d_1) & -v_{12}/E_1 & 0 \\ -v_{12}/E_1 & 1/E_2(1-d_2) & 0 \\ 0 & 0 & 1/G_{12}(1-d_{12}) \end{pmatrix} \quad (3) $$

where $v_{12}$ is the principal Poisson's ratio, which for simplicity is assumed here not to be degraded. The ply model introduces three scalar damage parameters $d_1, d_2, d_{12}$ which have values $0 \le d_i < 1$ and represent modulus reductions under different loading conditions due to microdamage in the material. For fabric plies $d_1$ and $d_2$ are associated with damage or failure in the principal fibre directions, and $d_{12}$ controls in-plane shear failure. In the general damage mechanics formulation [6] 'conjugate forces' or damage energy release rates $Y_1, Y_2, Y_{12}$ are introduced corresponding to 'driving' mechanisms for materials damage, and it is shown that with the compliance matrix (3) they take the form:

$$ Y_1 = \sigma_{11}^2/(2E_1(1-d_1)^2) \quad Y_2 = \sigma_{22}^2/(2E_2(1-d_2)^2), \\ Y_{12} = \sigma_{12}^2/(2G_{12}(1-d_{12})^2) \quad (4) $$

The damage parameters are defined in terms of the damage evolution functions $f_1, f_2, f_{12}$ and have the general form: $d_1=f_1(Y_1, Y_2, Y_{12})$, $d_2=f_2(Y_1, Y_2, Y_{12})$ and $d_{12}=f_{12}(Y_1, Y_2, Y_{12})$.

Specific forms for the evolution equations are required which should be consistent with test data. The elastic damage mechanics ply fabric model is based on the following assumptions:

(a) Fibre and shear damage modes are decoupled, with fibre damage determined by $Y_1$ and $Y_2$, and shear failure by $Y_{12}$.
(b) Fibre damage development may be different in tension and compression.

(c) For balanced fabrics ($E_1 = E_2$) damage development in the two fibre directions may be different, thus $d_1 = d_2$. However, it is assumed that $f_1$ and $f_2$ will have the same functional form ($f_1 = f_2$).

(d) The ply material is 'non-healing'; therefore damage during unloading is held constant until positive loading is applied which causes further damage accumulation.

(e) Damage development does not necessarily lead to ultimate failure of the ply and a global failure criterion is also necessary.

Due to condition (d) above, the damage evolution equations are based on the maximum value of the damage forces reached during the previous loading history. We thus introduce the quantities $\underline{Y}_1, \underline{Y}_2, \underline{Y}_{12}$ which are defined in terms of the maxima of $\sqrt{Y_i}$. Test data on UD composites [2] has shown that the square root of the damage forces is the quantity which arises more naturally, therefore:

$$ \begin{aligned} \underline{Y}_1(t) &= \max\{\sqrt{Y_1(\tau)}\}, & \underline{Y}_2(t) &= \max\{\sqrt{Y_2(\tau)}\}, \\ \underline{Y}_{12}(t) &= \max\{\sqrt{Y_{12}(\tau)}\}, & \tau &\le t \end{aligned} \qquad (5) $$

Taking into account (a) and (c) above, and assuming an elastic region without damage at the onset of loading, leads to the following expressions for the lower and upper thresholds of damage:

$$ \begin{aligned} d_1 &= 0, & \underline{Y}_1 < Y_{10} &\quad d_1 &= \alpha_1(\underline{Y}_1 - Y_{10}) \text{ for } Y_{10} < \underline{Y}_1 < Y_{1f} \\ d_2 &= 0, & \underline{Y}_2 < Y_{10} &\quad d_2 &= \alpha_1(\underline{Y}_2 - Y_{10}) \text{ for } Y_{10} < \underline{Y}_2 < Y_{1f} \\ d_{12} &= 0, & \underline{Y}_{12} < Y_{120} &\quad d_{12} &= \alpha_{12}(\ln \underline{Y}_{12} - Y_{120}) \text{ for } Y_{120} < \underline{Y}_{12} < Y_{12f} \end{aligned} \qquad (6) $$

Linear forms for $d_1$ and $d_2$, were found to be good approximations for fabric plies, and an equation linear in $\ln(Y_{12})$ was found to be required for modelling the shear behaviour at larger strains. Thus the evolution equations for a balanced fabric ply require the determination of two slope parameters $\alpha_1, \alpha_{12}$ and four damage threshold parameters $Y_{10}, Y_{120}, Y_{1f}, Y_{12f}$. Further refinements to allow different fibre damage behaviour in tension and compression, with an ultimate failure envelope are discussed in [7].

## 2.2. Elastic-plastic model for a fabric composite ply

For in-plane shear, deformations are controlled by matrix behaviour which may be inelastic, or irreversible, due to the presence of extensive matrix cracking or

plasticity. On unloading this can lead to permanent deformations in the ply. The extension of the fabric model to include these irreversible damage effects is now considered, based on the following main assumptions:

(f) The total strain in the ply is split into the sum of elastic and plastic (or inelastic) parts.

(g) Plastic strains are associated only with the matrix dominated in-plane shear response.

(h) A classical plasticity model is used with an elastic domain function and hardening law applied to the 'effective' stresses in the damaged material.

(i) Inelastic or plastic strain increments are assumed to be normal to the elastic domain function.

From (f) above the total strain $\varepsilon$ can be written as the sum of elastic $\varepsilon^e$ and plastic strains $\varepsilon^p$ ($\varepsilon = \varepsilon^e + \varepsilon^p$). The elastic strain component is given by (2). A plane stress model for a thin ply is assumed and only shear strain contribute to plasticity ($\varepsilon_{p1}^p = \varepsilon_{p2}^p = 0, \varepsilon_{p2}^p \neq 0$). Following [6], an elastic domain function is introduced $F(\tilde{\sigma}_{12}, R)$, where $\tilde{\sigma}_{12}$ is the 'effective' shear stress $\tilde{\sigma}_{12} = \sigma_{12}/(1 - d_{12})$ and R is an isotropic hardening function. R(p) is a function of an inelastic strain variable p. The elastic domain function has a simple form here since only the effective shear stress leads to plastic deformation:

$$ F = |\sigma_{12}|/(1 - d_{12}) - R(p) - R_0 \qquad (7) $$

where it is assumed that $R(0)=0$ and that $R_0$ is the initial threshold value for inelastic strain behaviour. The condition $F<0$ corresponds to a stress state inside the elastic domain where the material may be elastic damaging. It follows from the normality requirement (i) that $F=0, \dot{F}=0$ hence from (7) it can be shown that the plastic strain p is defined by

$$ \dot{\varepsilon}_{12}^p = \dot{p}/(1 - d_{12}) \quad \text{or} \quad p = \int_0^{\varepsilon_{12}^p} (1 - d_{12})d\varepsilon_{12}^p \qquad (8) $$

showing that p is the accumulated effective plastic strain over the complete loading cycle. The model is completed by specifying the hardening function $R(p)$. This is determined from cyclic loading tests in which both the elastic and irreversible plastic strains are measured. A typical form assumed for the hardening function [2] is an index function, which leads here to the general equation:

$$ R(p) = \beta p^m \qquad (9) $$

so that the shear plasticity model depends on the parameters $\beta$, the power index m and the yield stress $R_0$.

## 2.3. Delamination model

Delamination failures occur in composite structures under impact loads due to local contact forces in critical
regions of load introduction and at free edges. They are caused by the low, resin dominated, through-thickness shear and tensile properties found in laminated structures. In composites delamination models [3,4] the thin solid interface is modelled as a sheet of zero thickness, across which there is continuity of surface tractions but jumps in displacements. The equations of the model are given here for the case of mode I tensile failure at an interface. Let $\sigma_{33}$ be the tensile stress applied at the interface, $u_3$ the displacement across the interface, and $k_3$ the tensile stiffness. Following [4] an elastic damaging interface stress-displacement model is assumed:

$$ \sigma_{33} = k_3(1 - d_3)u_3, \\ d_3 = c_1(1 - u_{30}/u_3), \quad \text{for } u_{30} \le u_3 \le u_{3m}, \qquad (10) $$

with tensile damage parameter $d_3$, and $c_1=u_{3m}/(u_{3m}-u_{30})$. It can be verified that with this particular choice of damage function $d_3$, the stress-displacement function has the triangular form shown in Fig. 1, and $u_{30}$, $u_{3m}$ correspond to the displacement at the peak stress $\sigma_{33m}$ and at ultimate failure. The damage evolution constants are defined in terms of $\sigma_{33m}$ and $G_{IC}$, the critical fracture energy under mode I interface fracture, by $u_{30} = \sigma_{33m}/k_3$ and $u_{3m} = 2G_{IC}/\sigma_{33m}$. From these expressions it can be shown that the area under the curve in Fig. 1 is equal to the fracture energy $G_{IC}$. This interface model therefore represents an initially elastic interface, which is progressively degraded after reaching a maximum tensile failure stress $\sigma_{33m}$ so that the mode I fracture energy is fully absorbed at separation. For mode I interplay failure the interface energy $G_I$, defined as

$$ G_I = \int_0^{u_3} \sigma_{33} du_3 \qquad (11) $$

is monitored and, if this is found to exceed the critical fracture energy value $G_{IC}$, then the crack is advanced. For mode II interface shear fracture a similar damage interface law to (10) is assumed, with equivalent set of damage constants, $u_{130}$, $u_{13m}$ and critical fracture energy $G_{IIC}$.

In general there will be some form of mixed mode delamination failure involving both shear and tensile

<figure><img src="image_1.png" /><figcaption>Fig. 1. Mode I interface stress-displacement function.</figcaption></figure>

failure. This is incorporated in the model by assuming a mixed mode failure condition, which for mode I/mode II coupling could be represented by an interface failure envelope such as [4]

$$ \left(\frac{G_I}{G_{IC}}\right)^n + \left(\frac{G_{II}}{G_{IIC}}\right)^n = e_D \le 1 \qquad (12) $$

where $G_I$ and $G_{II}$ are the monitored interface strain energy in modes 1 and 2 respectively, $G_{IC}$ and $G_{IIC}$ are the corresponding critical fracture energies and the constant $n$ is chosen to fit the mixed mode fracture test data. Typically $n$ is found to between 1 and 2. Failure at the interface is imposed by degrading stresses when $e_D < 1$ using (10) and the corresponding shear relation. When $e_D \ge 1$ there is delamination and the interface separates. Fig. 2 depicts the energy interaction curves for the general case of mixed mode loading. Here a linear interaction of fracture energies for modes I and II is assumed so that $n=1$. Some experimental testing work has been undertaken [8], using a mixed mode test [9], which does indicate the interaction of energies is approximately linear in typical composites materials. The energy absorption for pure modes I and II are shown together with the response under typical mixed mode loading. For the mixed mode loading points (A) and (B) represent the limits of elastic response and full fracture respectively. The interface traction law must now use modified curves shown in Fig. 2 for mixed mode loading.

### 3. Code implementation and validation of composites model

#### 3.1. The explicit FE method

In general two finite element solution techniques are used in commercial FE software packages; namely, the 'implicit' and 'explicit' methods [11]. The implicit variation is more common and used in most general purpose FE codes for the solution of linear and moderate non-linear problems. In recent years, however, the explicit method has proved successful for the analysis of dynamic, highly non-linear problems, particularly where contact plays an important role. Automotive crashworthiness and metal stamping analysis are two notable examples. Recent experience has also found that this technique provides a robust and reliable method for composites impact analysis. In this case geometric instabilities due to ply buckling and material softening due to damage are both easily treated using the dynamic formulation. Furthermore, potential contact between delaminated plies, or contact with impacting bodies, can be efficiently and reliably treated. The ply and delamination models described in Section 2 have both been implemented in
<figure><img src="image_2.png" /><figcaption>Fig. 2. Fracture interaction curves for modes I and II, and the corresponding traction stress-displacement (crack opening) curves.</figcaption></figure>

the commercial explicit impact code PAM-CRASH [5].
The following two section give some brief details of this
implementation and validation.

3.2. Ply in-plane fabric model and validation

The proposed ply model has been implemented in the multi-layered ‘Mindlin–Reissner’ shell element [12]. As such it can either be used to model the complete laminate or, in the case of specifying only one layer, it can be used to represent an individual ply of the laminate. Using the latter approach a stack of individual plies may be constructed which are then joined together using the delamination model to create the improved laminate model. In order to reduce CPU costs, at the expense of accuracy, alternative modelling schemes could be envisaged. For example the laminate could be constructed using a reduced number of multi-layered shells, with each shell now representing a sub-group of plies; these are then joined using fewer delamination interfaces. This simplification would still allow the principle ply and delamination failure modes to be captured.

A predictor-corrector algorithm is used for the
implementation of the ply elastic-plastic damaging law.
At calculation time *t* for each finite element and for each
time step, the known variables are the total strains, the
associated strain increments and the set of variables
(strains, stresses and damage parameters) calculated at
the previous time step *t*−1. Since the elastic orthotropic
behaviour (1) requires the knowledge of the current
elastic strains, these are computed from the plastic
strains at the previous time *t*−1. It is then possible to
determine the three damage parameters corresponding
to this state of elastic strain, using the compliance
matrix (3), the definition of the damage energy release
rates (4) and the damage evolution functions (5)–(6).

<figure><img src="image_3.png" /><figcaption>Fig. 3. Cyclic shear curves for CF/epoxy (test versus simulation).</figcaption></figure>

The error made with this approximation for the deter-
mination of the updated damage parameters d₁, d₂, d₁₂
was found to be less than 3%, which is considered a
reasonable approximation for dynamic studies. Having
determined the elastic strains and the damages which
are associated to them, the state of stress can then be
determined. At this point, the elastic domain function
(7) is used to check if the predicted elastic strains are
correct or not. In the affirmative case, calculations
finish, and if not a correction of the stress state is then
applied. Using an iterative algorithm of Newton–
Raphson type (with a maximum of fifty iterations), the
plastic strains are re-calculated (8). The schema previously
proposed for the updating of elastic strain is then used
again. This process is repeated until a convergence
criteria is satisfied. Generally, it is found that the stress
state converges after a few iterations, usually less than
ten. These calculation steps are applied for each time
step and for each layered composite shell element in the
structure.
In the HICAS project [1] an extensive materials test programme has been carried out on carbon and glass fabric reinforced epoxy materials including in-plane and through-thickness tension, compression and shear tests. Cyclic shear tests to determine the plastic strain contribution have also been conducted. Test data were used first to justify the chosen forms for the damage evolution equations (6) and shear plasticity hardening law (9), then to obtain the required materials parameters for the model. The analysis procedures for determining the ply damage and plasticity parameters are described in detail in [7]. The materials dataset obtained has been used as the basis for code validation on single elements and materials test specimens. As an example Fig. 3 shows results of a single element simulation to validate the shear damage plasticity model. Using measured properties of CF/epoxy material a 4-node rectangular element was simulated under cyclic shear loads using the load method applied in a tension-shear test. That is the rectangular element had fibre directions at ±45° to the element axes and was loaded in cyclic tension. The figure shows significant material nonlinearity due to elastic damage and the extent of plastic shear strain measured in this material on unloading. It is not possible to synchronise exactly the test cycle with the simulation cycle and the load hysteresis effect is not included in the ply model. Nevertheless there is good agreement between test data and numerical simulation, which verifies the plasticity model and choice of hardening function.

In Fig. 4 results from simulations of actual test specimens of GF/epoxy are shown under 0°- and 45°-tension loads. The test specimen geometry dimensions were 220×15×2 mm for the 0°-tension specimens and 100×25.8×2.3 mm for the 45°-tension specimens, where in each case the length is the gauge length in the test. The specimens were modelled with about 100 laminated shell elements and the specimens were loaded at constant velocity. Fig. 4 compares computed load- deflection curves with actual test data. The displacements shown are clip

<figure><img src="image_4.png" /><figcaption>Fig. 4. Test specimen simulations under 0 and 45° tension.</figcaption></figure>

gauge displacements based on a 50 mm gauge length. The figure shows good agreement in the 0°-specimen up to brittle failure, and in the 45°-specimen up to about 5 mm displacement. The failure load was predicted well, but the computation broke off at about 8 mm displacement below the final failure deflection in the specimen. The deviation here at larger displacements is due to a number of reasons. The 45°-tension test is not considered valid for determining shear properties at larger strains, and there are different shear stress and strain measures being used in the analysis of the test to determine the damage and plasticity parameters and in the PAM-CRASH code. The results are however considered to be encouraging for the new ply damage model and confirm that the model and implementation are stable for structural calculations.

### 3.3. Interply delamination model

A novel approach has been developed to implement the delamination model of § 2.3 into the PAM-CRASH code. One approach for modelling delamination in composites would be to construct a detailed solid FE model of the laminate including discrete layers for both the plies and resin rich interfaces. This, however, leads to a CPU intensive computational model due to the enormous number of elements needed and, in the case of using an explicit integration scheme for the FE solution, an unfavourable timestep is associated with the thin solid elements in the interface layers. Consequently, the approach is not practical for large scale structural analysis. An alternative compromise approach is used here in which the laminate is treated as a stack of shell elements. Each ply or sublaminate ply group is now represented by one layer of shell elements and the individual ply layers are tied together using a 'sliding interface' with an interface traction-displacement law. This approach gives a good approximation for delamination stresses and failure, and has the added advantage that the critical integration timestep is much larger since it depends on the area size of the shell elements. An advantage of this numerical modelling approach is that large composite structures may be modelled efficiently with shells, or stacked shells, requiring fewer elements than solid models, and computationally expensive interface solid elements are eliminated. In this stacked shell laminate model the interface failure law (10)–(12) is applied to determine tractions and displacement discontinuities at the interface under mixed mode tension/shear loading.

The delamination model requires interface fracture energy under mode I ($G_{IC}$), mode II ($G_{IIc}$) and the fracture response under mixed mode loading. Pure modes I and II data may be readily obtained from Double Canti-lever Beam (DCB) and End Notched Flexure (ENF) tests, respectively. The delamination model with stacked
shell elements under pure mode I loading based on the interface law (10) was validated with DCB fracture mechanics test data for a composite laminate, as discussed in more detail in [10]. Mixed mode data requires

a specialised test [9], which can propagate an interface crack under a predetermined ratio of $G_I$ and $G_{II}$, as shown schematically in Fig. 5. This test can provide pure mode I data if loaded upward at point A (≈ DCB test),

<figure><img src="image_5.png" /><figcaption>Fig. 5. Mixed mode delamination test: setup, simulation results at two deformation states and comparison of test and simulation force time histories during crack propagation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6. Drop tower impact tests on CF/epoxy plates.</figcaption></figure>
pure mode II data if loaded downward at B (≈ ENF test) and mixed mode data if loaded downward at C. The ratio of $G_I$ to $G_{II}$ is varied by selecting a different lever arm length. Calibration of the numerical model is necessary to get the failure stress $\sigma$ since this is not directly available from the test. Fig. 5 also shows example simulation results for this test for a quasi-isotropic UD carbon/epoxy laminate [8] and compares test and simulation force histories at the load point for the special case $G_I = G_{II}$ for the case of a linear interaction with $n=1$ in (12). The good agreement between test and simulation in this mixed mode test case supports the code development of mixed mode interface failure and validates the materials parameters used in the simulations.

4. Simulation of composite plate impact tests

Encouraged by the code validations on single elements and test specimens, trial impact simulations with the damage/delamination model in PAM-CRASH are now considered on idealised composite structures in the form of CF/epoxy plates. The plates had nominal dimensions 300×300×4.5 mm and were made up of 16 plies of carbon fabric/epoxy with quasi-isotropic layup. The impact test setup in the DLR drop tower is shown in Fig. 6a, with the plates simply supported on a 250×250 mm square steel frame. The impactor head was a 50mm diameter steel sphere with an added mass of 21 kg, and various impact velocities were used in the range

<figure><img src="image_7.png" /><figcaption>b) Delamination damage at centre ply</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>a) Plate damage contours at failure</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>a) Delamination damage at centre ply</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>b) Plate damage during rebound</figcaption></figure>

<caption>Fig. 8. Simulation of CF/epoxy plate under impact by steel impactor at 6.28 m/s.</caption>
2.33–6.28 m/s, to give different impact energies 57–414 J and failure modes from rebound to full penetration.

For the simulation a stacked shell laminate model is used with 16 plies. However, in order to reduce CPU time the laminate has been simplified by using four sub-laminates of shell elements tied with the delamination sliding interface. For each shell a sub-group of 4 plies is used so that the correct quasi-isotropic laminate layup is represented. This model allows delamination at three interfaces in the plate, which is in agreement with typical observed failure behaviour. The composite plate is simply supported over the frame and impacted at its centre with the sphere, which is modelled here as a rigid impactor. Numerical results are presented here for two impact test cases, corresponding to impact velocities of 6.28 m/s when the impactor penetrated the plate, and 2.33 m/s when the impactor rebounded causing some back face cracking and extensive delamination.

### 4.1. Impact velocity 6.28 m/s

Fig. 7 shows typical simulated deformation plots during the impact event, with contours of ply damage and delamination damage parameters. Post processing allows information on both the ply damage and ply interface delamination to be assessed. The simulations predict well the penetration of the plate and the fibre fracture at the impactor, seen in Fig. 6b.

### 4.2. Impact velocity 2.33 m/s

The second simulation refers to a lower impact velocity, which led in the test to a rebound failure mode. Fig. 8 shows typical simulated deformation plots during the impact event. At this lower impact velocity the impactor rebounds and the plate has some fibre cracking on the back face and C-scanning indicated extensive delamination in the impact region. The rebound was successfully modelled in the FE simulation. The delamination contours and size of delamination predicted at the plate middle plane are also found to be in good agreement with the results of the C-scan test.

### 4.3. Comparison of measured and predicted loads

Figs. 7 and 8 demonstrate that PAM-CRASH with the new ply damage and delamination model can successfully predict the failure modes in the impacted plates, by penetration or rebound, with convincing contour plots of ply damage and delamination damage. However, to fully validate the code it is necessary to compare quantitative predictions of impact loads, plate deflections and energy absorbed with measured test data. These data have been measured in the low velocity plate impact test programme. The load cell on the impactor head records the contact load-time pulse during impact.

Test data for these impactor loads for the 6.28 and 2.33 m/s impact tests are plotted in Fig. 9a & b. These are compared with simulations using the ply damage laminate model with a single laminated shell (without delamination), and with the 4 sublaminate stacked shells (with delamination). The figures show that the single laminated shell model which has ply damage and plasticity overpredicts the measured peak failure loads by a factor of two or higher. The simulation with the 4-sublaminate delamination model has significantly lower peak loads than the single shell result and is much closer to the test data, particularly for the 6.28 m/s penetration test. The simulations show that delamination in the plate reduces the plate bending stiffness near the impactor which leads to lower peak loads, greater deflections and reduced energy absorption. It follows that the simulation with only laminate ply failure and no delamination gives poor agreement with test data, with predicted peak loads and energy

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Comparison of impactor load-time pulses from test, with simulations using the single shell and 4 ply delamination models.</figcaption></figure>
absorption too high. The results for the delamination model are encouraging and show that delamination effects are significant in plate impact and that these can be successfully modelled with the improved code. Modifications to the delamination failure initiation stresses and interface fracture energies in the model leads to further changes in the peak loads, showing that the impact simulations are sensitive to the interface parameters, which are not well documented for these materials. More experience and relevant test data are now needed to further improve the simulation results.

## 5. Conclusion

This paper has presented a materials failure model for composites with fibre fabric reinforcement which includes both intraply damage and plasticity, and interply delamination. This is implemented in the dynamic FE code PAM-CRASH and code validations have been successfully carried out on single elements and both in-plane and fracture mechanics materials test specimens. Trial structural simulations of low velocity drop tower impact tests on carbon fabric/epoxy composite plates is described, which gives good correlation with test data. An important feature of the model is that it distinguishes clearly between different failure modes in the structure. In post-processing it is possible to follow the progression during impact of the fibre and shear damage parameters in the shell, the irreversible plastic shear strains, the fibre strains and the extent of delamination. Thus it is possible to simulate both the impact failure modes and failure progression during impact loading in composite structures. Ongoing work is concerned with studying methods for bringing composites rate dependent properties into the materials models and the FE code, and with the application of the model to impact in larger composite shell structures.

## Acknowledgements

This work was developed in the EU project HICAS [1]. The authors wish to acknowledge financial support from the CEC and the HICAS partners for their scientific contribution. Mr P. Rozycki thanks ESI for funding the activities at LAMIH and Dr Haug (ESI) and Professor Coutellier (LAMIH) for their support.

## References

[1] HICAS. High velocity impact of composite aircraft structures, CEC DG XII BRITE-EURAM Project BE 96-4238, 1998.

[2] Ladevèze P, Le Dantec E. Damage modelling of the elementary ply for laminated composites. Composites Science and Technology 1992;43:257-67.

[3] Allix O, Ladevèze P. Interlaminar interface modelling for the prediction of delamination. Composites Structures 1992;22:235-42.

[4] Crisfield MA, Mi Y, Davies GAO, Hellweg HB. Finite element methods and the progressive failure modelling of composite structures. In: Owen DRJ, et al., editors. Computational plasticity—fundamentals and applications' Part 1. Barcelona: CIMNE, 1997. p. 239-54.

[5] PAM-CRASH™ FE Code, Engineering Systems International, F-94578 Rungis Cedex, France.

[6] Ladevèze P. Inelastic strains and damage. In: Talreja R, editor. Damage mechanics of composite materials. Composite Materials Series, 9. Amsterdam: Elsevier, 1994, Chapter 4.

[7] Johnson AF. Modelling fabric reinforced composites under impact loads. Composites A, 2001;32(9):1207-1215.

[8] Gruber J. Private communication, IPF Dresden e. V. TU Dresden, Germany. July 2000.

[9] Reeder JR, Crews JH. Mixed-mode bending method for delamination testing. AIAA Journal 1989;28(7):1270-6.

[10] Johnson AF, Pickett AK, Impact and crash modelling of composite structures European Conf. Computational Mechanics, ECCM'99, Munich, 1999.

[11] Cook RD, Malkud DS, Plesha ME. Concepts and applications of finite element analysis, Wiley, 1989.

[12] Belytschko T, Lin JL, Tsay C. Explicit algorithms for the nonlinear dynamics of shells. Computer methods for applied mechanics and engineering 1984;42:225-51.