Open Archive Toulouse Archive Ouverte (OATAO)

OATAO is an open access repository that collects the work of Toulouse researchers
and makes it freely available over the web where possible.

This is an author-deposited version published in: http://oatao.univ-toulouse.fr/

Eprints ID: 3324

To link to this article: DOI: 10.1016/j.ijsolstr.2009.03.010

URL: http://dx.doi.org/10.1016/j.ijsolstr.2009.03.010

To cite this version: BOUVET, Christophe, CASTANI, Bruno, BIZEUL, Mathieu,

BARRAU, Jean-Jacques. Low velocity impact modelling in laminate composite panels
with discrete interface elements. International journal of solids and structures,
2009, vol. 46, no14-15, pp. 2809-2821. ISSN 0020-7683

Any correspondence concerning this service should be sent to the repository administrator:
[email protected]

Low velocity impact modelling in laminate composite panels with discrete

interface elements

Christophe Bouvet a,*, Bruno Castani a, Matthieu Bizeul b, Jean-Jacques Barrau a

a
b

Universit de Toulouse/UPS/LGMT, 118 route de Narbonne, 31062 Toulouse cedex 04, France
Universit de Toulouse/ISAE/DMSM, 10 av. E. Belin, 31055 Toulouse cedex 04, France

a b s t r a c t

Impact
Interface
Matrix crack
Delamination
Modelling

A model enabling the detection of damages developing during a low velocity/low energy impact test on
laminate composite panels has been elaborated. The ply model is composed of interface type elements to
describe matrix cracks and volumic nite elements. This mesh device allows to respect the material
orthotropy of the ply and accounts for the discontinuity experimentally observed. Afterwards delaminations are described with interfaces similar to the ones observed with matrix cracks and the coupling
between these two damages are established. In the rst step, simple stress criteria are used to drive these
interface type elements in order to assess the relevance of model principle. Nevertheless, the well known
problem of mesh sensitivity of these criteria prevents the use of this model for now as a predictive tool
but rather as a qualitative tool. An experimental validation is carried out thanks to impact experimental
tests performed by Aboissiere (2003) and a very good match has been found. However, this model could
predictivelly be used and would allow to foresee an original method to detect delaminations during an
experimental test. This modelling has been successfully tested experimentally and compared to a C-Scan
ultrasonic investigation.

1. Introduction
Composite materials are being increasingly used in airframe and
spatial applications thanks to their interesting mechanical characteristics and low specic weight. Nevertheless, for structures submitted to low energy impacts or minor objects drop, like tools
during assembly or maintenance operation, composite laminates
reveal a brittle behaviour and can undergo signicant damages in
terms of matrix cracks, bres breakages or delamination. These
damages are particularly dangerous because they drastically reduce
the residual mechanical characteristics of the structure, and at the
same time can leave very little visible mark onto the impacted surface. Many authors have consequently studied the impact behaviour of composite structures and their effects on residual
strength, both experimentally (Abrate, 1998; Aoki et al., 2006; Petit
et al., 2007; Davies and Olsson, 2004; Kwon and Sankar, 1993 . . .), as
well as numerically (Allix and Blanchard, 2006; Choi and Chang,
1992; Finn and Springer, 1993; Li et al., 2006; De Moura and Gonalves, 2004; Guinard et al., 2002; Hou et al., 2001 . . .), but a lot of
work is still necessary to improve the modelling of the damage
developing during impact on composite laminates to better assess

author. Address: Universit Paul Sabatier, Filire Gnie Mcanique, Bt. 3PN 118, route de Narbonne, 31062 Toulouse cedex 04, France. Tel.: +33
(0) 5 61 55 84 26; fax: +33 (0) 5 61 55 81 78.
E-mail address: [email protected] (C. Bouvet).

numerically their residual mechanical characteristics in order to

optimise their design.
The rst step to achieve in order to numerically design composite structures is the modelling of the impact and in particular the
impact damages like matrix cracks, bre/matrix debonding, delamination or bres breakage. These damages are classically divided in
two parts:
 The intralaminar damages, i.e. the damages developing inside
the ply like matrix cracking, bre/matrix debonding or bres
breakages.
 The interlaminar damages, i.e. the damages developing at the
interface between two consecutive plies, namely delamination.
When these two damage types are considered separately, they
are globally well understood and well simulated. The rst ones,
the intralaminar damages are often modelled thanks to damage
internal variables (Allix and Blanchard, 2006; Li et al., 2006) or more
rarely thanks to fracture mechanics (De Moura and Gonalves,
2004). These different models seem to be in correlation with experimental observations and properly explain the physical phenomena
of the different damages appearing inside ply. The second ones, the
interlaminar damages are principally modelled thanks to interface
type elements driven by damage internal variables (Allix and Blanchard, 2006) or thanks to fracture mechanics (Li et al., 2006; De
Moura and Gonalves, 2004; Collombet et al., 1996). These different

models seem to be in correlation with experimental observations

and explain the physical phenomenon of the delamination creation.
But the main problem of the physical understanding of the impact
damages formation is the interaction which exists between the two
damage types inside and between plies. The key point of impact
damage formation modelling is these interaction, consequently this
article will focus on this specic issue.
An interesting physical explanation of this interaction is proposed by Renault (1994). In this work, he proposes an experimental scenario of the different damages developing in a composite
during a low velocity/low energy impact test. This scenario
proposes to give matrix cracks a precursor role regarding the
development of specic delamination. This idea is globally admitted in the literature (Abrate, 1998; Davies and Olsson, 2004; Allix
and Blanchard, 2006). However this scenario begins with the
development of matrix cracks in the impact zone below the impactor (Fig. 1). These cracks, of transverse direction (in (l, t) plane,
where l is the longitudinal or bres direction and t the transverse
one) grow up during the loading following the bre direction.
Therefore, in each ply, a strip of bres and resin disjoints and slides
in the normal direction of this ply (z) (Fig. 1). This disjointed strip
creates an interlaminar zone of tension stress between two consecutive plies and is susceptible to induce in this zone the formation
of a delamination. In Fig. 1a, a schematic application of this scenario is proposed with a [45, 0, 45] stacking sequence, where
45 is the lower ply and 45 the upper ply. This stacking sequence
is not a real one but represents only a part of a real stacking. In
Fig. 1a, the disjointed strips of the rst two plies are drawn and
illustrate the interlaminar zone of tension stress. This zone, limited
by the disjointed strips of the two adjacent plies, has a triangular
shape with a size which grows from the impacted side to the
non-impacted side. Fig. 1b illustrates the interlaminar zones of tension stress between the 45/0 and 0/45 plies.
Another explanation of this interaction has been proposed by
Choi and Chang (1992) who explained the formation of a delamination by two main phenomena (Fig. 2):
 The rst one is the delamination induced by inner shear cracks.
A shear matrix crack located in the inner plies of the laminate
generates a substantial delamination along the bottom interface
and a small conned delamination along the upper interface of
the cracked plies.
Impacted Side

 The second one is the delamination induced by a surface bending crack. A bending matrix crack located at the surface ply of
the laminate generates a delamination along the rst interface
of the cracked ply.
These two damages scenarios underline the fundamental role
played by the interaction between these two damages which exists
during the impact damage progression. Consequently this interaction must be taken into account in modelling to correctly simulate
the experimental observations. To account for this interaction,
many solutions have been suggested in the literature.
For instance one put forward by Choi and Chang, with reference
to the impact damage scenario mentioned above, has used a very
interesting delamination criterion evaluated only with mean stresses of the upper (noted n + 1) and lower (noted n) plies (z is the
interface normal direction):

sn1
tz
Da  4 n1
Si

!2

snlz

2

Sni

rnt 2 5
1

where sn1
is the shear stress in the upper ply with t and z,
tz
respectively, the transverse and normal direction of the considered
ply, snlz and rnt , respectively, the shear and normal stress of the lower ply expressed in the considered ply reference, Sni and Sin1 the
in situ interlaminar shear strength, respectively, in the lower and
upper plies and Y n the in situ transverse tensile or compressive
strength in the lower ply. The two rst terms in the square brackets
represent the effects of a shear matrix crack on the creation of a
delamination. These terms allow to account for the phenomena
named by Choi and Chang the delamination induced by inner
n1
stress charactershear cracks (top of Fig. 2). The rst term in stz
izes more specially the opening in fracture mode I of the delamination due to the shear matrix cracks of the upper ply and the second
term in snlz stress characterizes the propagation in fracture mode II
of the delamination due to the important stiffness in the longitudinal direction (due to bres) of the lower ply. This last phenomenon
is meant to explain, according to the authors, the propagation of
delamination in the bres direction of the lower ply.
The third term in the square brackets represents the effects of a
matrix transverse crack on the creation of a delamination in the
upper interface. This term allows to account for the phenomenon

Ply n2 (0)

Ply n3 (45)

Propagation
direction
Impact
zone

Impact
zone

Ply n3 (45)
Ply n2 (0)
Ply n1 (-45)

45

0
-45

Impact
zone

Intralaminar Matrix Cracks :

Creation of Disjointed Strips

Interlaminar Tension Stress

Zone of Triangular Shape

Ply n2 (0)

Section B-B

Ply n1 (-45)

Impact
zone
45

Propagation
direction

Impact
zone

Disjointed
strip

45
0
-45

Section A-A

45
Zones of interlaminar
0
-45
tension stress

0
-45

Non-Impacted Side

-a-

-b-

Fig. 1. Formation mechanism of delaminations (a) and interface tension stress zones (b).

C. Bouvet et al. / International Journal of Solids and Structures 46 (2009) 28092821

2811

Fig. 2. Impact damage mechanisms. Top: delamination induced by inner shear cracks. Bottom: delamination induced by surface bending crack.

named by Choi and Chang the delamination induced by surface

bending crack indicated at the bottom of Fig. 2.
And the Da parameter is an empirical coupling coefcient to
account for the precursor effect of the matrix cracking on the
delamination which is practically taken equal to 1.8. Then the quadratic delamination criterion written inside the square brackets is
considered reached if it equals 0.56 1=Da in presence of matrix
cracks and if it equals 1 in their absence. However this delamination criterion is only assessed if a matrix cracking has been detected in the lower or upper ply, the matrix cracking being
classically evaluated by the help of a stress quadratic function.
Then for this criterion, the presence of matrix crack is necessary
before the delamination formation. The authors have concluded
that model prediction is in correlation with test data on different
stacking sequences.
One drawback of this criterion is the different role played by the
lower and upper ply. We want to point out that this designation of
a lower or upper ply is purely articial and in the model proposed
in this paper, which is only partially inuenced by this criterion, an
equivalent role is played by lower or upper ply. But in addition we
want to demonstrate that it is the bending effect induced by an
impact which gives the lower plies a more important role in
delamination progression.
Another problem of the model is that it detects only the delamination onset and does not account for the degradation due to it
which can be a limitation for other bigger delamination area. In
this article, the modelled delamination areas are practically less
than 1000 mm2. Then to represent the degradations due to delamination and to completely simulate the impact damages, an interface element seems to be necessary. Consequently, in the present
model, an interface element with a total degradation possibility
is elaborated.
This idea is equally used by Allix and Blanchard (2006) which
are looked at an interface element driven by three damage internal
variables; the rst one, representing the damage created by the
normal stress, is associated with the rst mode of fracture, and
the other two, representing the damage created by the shear stresses, are associated with the second and third modes of fracture.
Then a non-associated damage evolution law is introduced to substitute these three damage internal variables with only one damage variable. The intraply damage is simulated with a similar
way thanks to three damage internal variables, one for the matrix

cracks in transverse direction, one for the in-plane shear matrix

cracks and one for the bres breakages. Then experimental evolution laws are introduced to drive these three damage internal variables. Finally, in this model, many couplings have been introduced
between the variables of intraply damages and the ones for delamination but no coupling is explicitly written between these two
damage types.
Afterwards this model is used to simulate the failure of holed
sandwich structures and shows good match with experimental results monitored but it shows a problem to simulate the creation of
delamination during an impact test. The authors conclude that due
to the complexity of the phenomena involved in what is called
delamination, the identication of the interface is still in its early
stage.
Li et al. (2006) have used failure mechanics to simulate the
delamination in a laminated panel and a pipe during an impact test.
The energy release rates are evaluated thanks to VCCT (Virtual Crack
Closure Technique) methods and the transverse matrix cracking
damage is taken into account thanks to a degradation of the strength
modulus. The simulations are in correlation with experiment even if
the stacking sequences are particularly simple. Afterwards, the
authors have underlined the difculty to implement this model in
a commercial FE (Finite Element) code because it requires accessing
and transmitting information between elements in the neighbourhood of the delamination front which is a problem equally met in
the implementation of the proposed model. Finally they have proposed the use of decohesive interfaces between sub laminates,
whose properties are degraded as described by a damage variable
in the context of continuum damage mechanics. This model has
been applied to a few slightly different problems, like DCB (Double
Cantilever Beam) tests, but its applicability to more sophisticated
scenarios, such as modelling the damage in impact problem is yet
to be explored which is seen as next phase of development.
Moura and Gonalves (2004) have equally built an interface element with a softening behaviour between the three stresses; one
normal and two shear stresses; and the three corresponding displacements. The coupling between these three damage variables
is made simply thanks to linear interpolation which allows to simulate the fracture mixed modes and the area under the softening
curve is classically linked to the energy release rate (Mi et al.,
1998). These interface elements are equally used to account for
the matrix cracking in the transverse direction. This idea, rarely

used in the literature, seems to be very interesting because it allows to simply introduce the interaction between intra and inter
plies damages and consequently is suggested in the proposed
model.
Afterwards, these interfaces are included in a laminate panel
mesh but only where damages were detected during an experimental test which is of course a limitation for a predicting model.
The simulations are quite in accordance with experiments even if
the delamination is overestimated and the matrix cracks underestimated. According to the authors, these differences can be attributed to the mixed-mode damage model used and to the inter/
intralaminar values of fracture mechanics characteristics which
are considered identical. However, the apparent simplicity of the
stacking sequences does not allow to test the predictive capacity
of this model.
Collombet et al. (1996) have equally used an interface element
driven by a simple criterion in normal stress which is implemented in an explicit FE code. The simplicity of the delamination
criterion driven only with normal stress is explained by the big
importance given to the rst opening mode of fracture compared
to the second and third fracture modes. This hypothesis is difcult to be experimentally tested and is largely discussed in the literature but it is used in the proposed model because it enables a
very good delamination prediction comparing to experimental results. More experimental investigations are yet necessary to be
done to evaluate the inuence of each fracture mode on delamination formation.
Afterwards, in the modelling of Collombet et al., the matrix
cracking is taken into account by a simple transverse stress criterion and the coupling between inter and intra laminar damage is
imposed by a precursor role of the matrix cracks: a delamination
is possible only if the lower ply is saturated in the matrix cracking.
The comparison between experiment and simulation is quite in
accordance even if the apparent simplicity of the stacking sequence does not allow to test the predictive capacity of this model.
Before concluding on this brief review of a few impact models,
one is going to focus on two works dealing specially with the interaction between the inter and intraply damages.
Ladevze et al. (2006) have studied the bridge between the
mesomodel of impact damage mentioned above (Allix and Blanchard, 2006) and a micro model in microscale. Their main conclusion is that the mesomodel can be interpreted as the homogenized
result of micro model. Thanks to this homogenisation, they have
assessed the inuence of intraply micro damage on the inter ply
meso damage, i.e. the interaction between the intra and inter plies
damages. They have concluded in particular that intralaminar
damages have a negligible inuence on interface damage under
normal stress loading but a major role under shear stress loading.
This conclusion is in correlation with experimental scenario of
Renault (1994) who has given to shear matrix cracks and consequently to shear stresses a major role on the delamination formation. This idea has been largely used for the model building.
Then the authors have concluded that in their opinion a more
specic link between delamination and transverse cracking must
be derived.
A work with similar conclusions has equally been performed by
Lammerant and Verpoest (1996). In their approach, they have
modelized a laminate panel with volumic elements and damage
springs for matrix cracks, as well as for delamination cracks. These
damage spring elements are driven by critical energy release rate
in mixed mode with linear interaction. The energy release rates
are identied for 0/0, 0/90 and 45/45 interfaces and for
intraply matrix cracks. They have studied only the propagation of
these damages, i.e. initial damages are rst present in the structure. They have shown in particular that the delamination shape
depends strongly on the rst matrix cracks and have concluded

that the existence of matrix cracks cannot be neglected when calculating the delamination development. This idea is present in the
proposed model even if the cases of initiation and propagation of
delamination are distinguished. In fact the existence of precursor
matrix cracks seems to be necessary for a delamination initiation
but not for its propagation.
Before concluding from this bibliography, which is of course not
exhaustive, the authors invite the interest readers to consult review papers of Abrate (1998) or Davies and Olsson (2004).
The conclusions which are made from this bibliography are the
following:
 An interface element is necessary to correctly simulate a delamination and in particular its degradation.
 An interface element is necessary to correctly simulate matrix
cracks. Indeed we think that a continuum degradation model
of the ply cannot take into account the effect of the matrix crack
on the delamination. This important remark, which can be
proved, is in our opinion the reason of an appropriate model
derived from experiment.
 A coupling between the intra and inter ply damages is necessary
and informations must be exchanged between the interfaces
elements of the matrix cracks and delamination.
Then the aim of this work is to build a FE model which allows to
account for these different remarks:
 Interface elements are used to simulate delamination.
 Interface elements are used to simulate matrix cracks in each
ply and consequently the FE mesh must respect the bres direction of each ply.
 The FE proposed allows the discussing between these two interface element types in order to account for the interaction
between delamination and matrix cracks.
Finally this model is used to simulate experimental impact
tests performed by Aboissiere (2002, 2003) on laminate composite 100  150 mm2 plates simply supported by a 75  125 mm2
shadow IGC04:26:383N 4 Airbus) with a spherical impactor of
16 mm diameter. The material used is a prepreg with carbon unidirectional bres and epoxy matrix HTA/EH24 manufactured by
HEXCEL of around 0.25 mm thickness ply. The material characteristics evaluated by test are summarized on Table 1. Where El and
Et are the Young Modulus in longitudinal and transverse direction, respectively, mlt the Poisson ratio, efl the failure strain in longitudinal direction, rft the failure stress in transverse direction, Glt
the shear modulus, sflt the failure shear stress and GI the critical
energy release rate obtained in propagation with a 0/0
interface.
Among the numerous experimental investigations in the literature (Aoki et al., 2006; Davies and Olsson, 2004; Petit et al., 2007;
Kwon and Sankar, 1993 . . .), Aboissieres work (2002, 2003) was
used to set this model because of the following interests:
Firstly, the main stacking sequence used in these tests
0 2 ; 45 2 ; 90 2 ; 45 2 S respects the classical stacking laws, as the
mirror symmetry and no angle exceeds 45 between two consecutive plies. This stacking is very simple, but does not allow to have a
too complicated model; other model verications are in progress in
more specic industrial sequences.

Table 1
Material parameters.
El (GPa)

Et (GPa)

mlt

efl (%)

rft (MPa)

Glt (GPa)

sflt (MPa)

GI (N/m)

143

100

0.29

1.357

80

5.1

77

280 50

Secondly, the damage observations, in particular delaminations,

were realised with ultrasonic investigation of high quality (Fig. 13)
which allows to simply identify each interface. On top of that, CScan was equally performed on the non-impacted side (Fig. 13b)
which gives supplementary informations on the delamination
shape and in particular the exact shape of the rst interface of
non-impacted side.
2. Numerical modelling principle
In our opinion and according to the concluding remarks mentioned above, the FE modelling of each ply of the composite laminate panel must account for the discontinuity created by matrix
cracks and in particular by shear matrix cracks. Consequently the
ply mesh must respect the exact material orthotropy in order to
allow for the creation of these cracks. These matrix cracks are globally created by three stress components (Fig. 3):
 A transverse stress rtt which creates cracks normal to transverse direction (t).
 A out-of-plane shear stress stz which creates cracks inclined of
about 45 .
 And an in-plane shear stress slt which creates cracks normal to
transverse direction (t).
In order to respect at best these cracks direction, the ply is
meshed in little strips in the longitudinal direction with one volumic element in its width and thickness of the ply (Fig. 4). In this
model, only inclined direction is not respected to avoid too complex modelling. This hypothesis seems reasonable if the ply thickness is small compared to the laminate thickness. Moreover, this
hypothesis justies equally the unique FE in the ply thickness.
Afterwards these little strips are connected together with interface type elements of null length (Fig. 4). There are practically not
interface elements but spring elements which are used for the sim-

Transverse cracks

plicity of programming in the used FE code SAMCEF . The main

difference between these two element types is that interface element can easily take into account the surface element which is
very important for a specic non-homogeneous mesh. In our case,
the proposed ply mesh imposes a homogeneous size mesh and
these two element types are very similar but it is easier to develop
spring element in a software like SAMCEF .
Then null length springs of very high stiffness (typically 107 N/
mm) are used to join together two consecutive strips and these
stiffnesses are put to zero if matrix crack exists. In order to avoid
numerical instabilities, the stiffnesses are cancelled progressively
with time, and during this delay, the global displacement is kept
constant. This classical way, to avoid instabilities, is similar to
delay effect proposed by Allix and Blanchard (2006).
The matrix crack is derived from a classical quadratic equation:

hrt i

!2

s2 s2tz

lt

f
t

sflt

tz

z
l
t

where rt is the transverse stress, slt and stz the shear stresses in the
(lt) and (tz) planes, < > the positive value and rft and sflt the failure
stresses mentioned above. A particularity of this model is to drive
this matrix cracking criterion of these springs thanks to stresses
in adjacent elements and not to stresses in the spring. Then, the criterion is evaluated in each volumic element of the ply, and a spring
is broken if this criterion is reached in at least one of its 4 neighbouring volumic elements. For example, for the spring R1 of the
Fig. 4, the criterion is evaluated thanks to the mean stress in the 4
volumic elements E1 ; E2 ; E3 ; E4 , and not to the force in the spring,
which allows to avoid stress concentration at the tip of the matrix
cracks. This criterion can be considered as an average stress criterion as the one proposed by Whitney and Nuismer (1974). In fact,
this is similar to average stresses over a distance which depends
on the mesh size. This mesh sensitivity will have to be further studied, nevertheless in the present simulated impact test, the matrix

In-plane shear cracks

Out-of-plane shear cracks

2 6 1

lt

tz
z

tz

lt

lt

Fig. 3. The three types of matrix cracks inner the ply.

Fibres
direction

Model

E4 A
E1

E3
R1
B E2

Matrix Cracks
Planes

Spring R1 : driven by
criterion in volumic
elements E1, E2, E3, E4

Model

Disjointed
Strips
Fig. 4. Model of the ply.

Broken
Springs

Safe
Springs

cracks revealed to be important especially for damage initiation and

not for damage propagation. However we want to point out that
this propagation is determinant in the nal damage morphology.
This ply model enables to correctly take into account the proposed experimental observations in respecting the orthotropic
behaviour of the ply, but it obliges to build a rened and fastidious
mesh. In fact this model imposes a constant size of the FE mesh in
the studied zone. Another particularity of this ply model is to
impose a meshing for the 45 plies with a diamond-shaped element to allow for the coincidence between two consecutive plies
(Figs. 7 and 8). Therefore only stacking sequences with 0, 90
and 45 plies can be meshed, which is not a limit because most
industrial applications are of this type.
In this model, the bres failure is equally taken into account
because this damage can cause drastic decreasing of material stiffness in the bres direction. It is classically simulated at each Gauss
point of volumic elements with a strain failure criterion:

el 6 efl

and Blanchard, 2006), similar to the one presented above for matrix crack, is implemented.
In fact, four springs are necessary (Fig. 6a), because there are
two nodes for the upper ply and two nodes for the lower ply. Then
each spring makes the bonding between one upper ply node and
one lower ply node. Physically, each spring represents one quarter
of the delamination concerned surface (Fig. 6b). For example, in
Fig. 6b, the spring R1 between the nodes A and D represents the
surface R1 . Now a criterion must be dened to drive interface crack
spring and to simulate the delamination. Two approaches are classically possible, one in fracture mechanics with a criterion in energy release rate and another one in limit stress. In the actual

0 ply

90 ply

3
f
l

where el is the strain in the bre direction and e is given in Table 1.

When this criterion is reached, the longitudinal Young modulus El
and the shear modulus Glt ; Gtz and Glz are drastically decreased,
practically divided by 20 to avoid numerical instability. This bre
criterion, which can seem secondary compared to delamination criterion, strongly inuences the delamination propagation in particular in the 90 direction during the impact test mentioned below
(Fig. 17).
When the different plies are meshed with volumic and matrix
cracks spring elements, springs are added to join them together
and allow delamination (Fig. 5).
Therefore these delamination spring elements have the following properties; if no interface cracks exist, two consecutive plies
are attached together with four null length springs of very high
stiffness (typically 107 N/mm in the three directions) otherwise
the stiffness is put to zero (Fig. 5). And to avoid numerical instabilities during the cancellation of stiffness, the delay effect (Allix

45 ply

x
y

-45 ply

x
Fig. 7. The four mesh types of 0, 90, 45 and 45 plies.

Disjointed
Strips
Model

Delamination
Fig. 5. Model of the interfaces.

Fig. 6. The four springs of delamination (a) and the covered surface (b).

Fig. 8. Finite element mesh.

model, a limit stress criterion is implemented for its simplicity and

an energy release rate criterion is actually studying. This choice of
very simple criterion was made in the rst step in order to assess
the relevance of model principle. Nevertheless the well known
problem of mesh sensitivity (Mi et al., 1998) of this criterion prevents the use of this model for now as a predictive tool but rather
as a qualitative tool. The main objective of this model is by the way
the physical understanding of the impact damage creation rather
than the predictive model building.
On top of that, with regard to the bibliography remarks and in
particular the experimental delamination scenario proposed by
Renault (1994), only the interlaminar normal stress is considered
in the criterion. This hypothesis neglects the second and third fracture modes as regard to the rst opening mode of fracture. With
this hypothesis, good results on delamination morphology compared to the experiment are obtained.
This result, which is conrmed by other authors (Collombet
et al., 1996), is of course available only for the considered impact
test and is probably not applicable to other damage scenarios like
impact on thick plates or tension of holed composite plates. So this
delamination criterion is used:

For i 1 to 4 :

rzi

F zi
6 rlim
Surf =4

where rzi is the interlaminar normal stress in the z direction for the
spring i, F zi the force in the spring i, Surf the surface concerned with
the group of four springs and r lim a limit stress of interface crack.
This limit stress, fundamental for this model, is characteristic of
the delamination initiation, with or without preliminary matrix
cracking as well as its propagation. Then three delamination cases
can be distinguished:

 The rst case is the initiation which is dened by a delamination

creation without preliminary damage. In this case the limit
stress rlim is chosen equal to the matrix failure stress
rft 80 MPa. But practically in every simulated test this criterion
is never reached before the propagation one, which is coherent
with the work performed by Choi and Chang.
 The second case is the propagation which is dened if a neighbouring delamination spring element is broken. Practically, each
delamination spring has four neighbours and is considered in
propagation if at least one of its four neighbours is broken. In
this case, the rlim value has been identied thanks to the delaminated area measured experimentally: rlim 10:5 MPa (Fig. 11b).
 The last case is a pseudo-propagation which is dened if a
matrix cracking exists preliminary. Practically, each delamination spring is confounded with two matrix cracking springs,
one for the upper ply and one for the lower ply. Then a delamination spring is considered in pseudo-propagation if at least one
of these two matrix cracking springs is broken. In this case, the
rlim value has been chosen equal to the propagation one:
rlim 10:5 MPa.
This value of rlim could seem very low compared to the matrix
failure stress rft 80 MPa but represents a different phenomenon: the failure propagation. So this limit stress should rather be
related with the critical energy released rate GI 28050N=m
measured experimentally. In order to achieve this relation, a DCB
test was simulated using this criterion. However to be representative of the impact test, two plies at 0 and 45 of 0.5 mm thickness
were meshed and jointed together with delamination interface elements with the limit stress of 10.5 MPa above mentioned (Fig. 9).
The force versus displacement curve can be obtained (Fig. 10a)

Fig. 9. DCB test.

450
100

G (N/m)

F (N)

75
50

300

150

25
0

0
0

15

30

45

60

12.5

25

37.5

Y (mm)
-b-

d (mm)
-a-

Fig. 10. Force versus displacement curve (a) and critical energy release rate at the crack tip (b).

and a local calculation of the critical energy release rate with a

VCCT method was made at the crack tip (Fig. 10b). This calculation
was performed for each spring element, thanks to stress obtained
just before failure and displacement just after, along the crack tip
for a crack of about 60 mm length. Outside boundary effects
(Fig. 10b) due to particular mesh, a constant critical energy release
rate of about 300 N/m is obtained. This value is in correlation with
the experimental value of 280 50 N/m measured by Aboissiere
and conrms the relevance of the used limit stress value. Due to
the VCCT method, the fracture modes II and III were also evaluated,
and this calculation showed that these two modes are less than 5%
of the fracture mode I and can be considered negligible for this test.
However we want to underline that the critical energy release
rate evaluated thanks to this numerical DCB test is of course mesh
sensitive but is also sensitive to the stiffness of the structure. In
fact, the energy release rate is dissipated in the delay effect of
the springs as mentioned above, and depends totally on the local
stiffness around the spring and therefore also on the local mesh.
This calculation is only valid with the presented mesh and for this
type of solicitation, therefore the same mesh is used for both the
DCB test and the impact test. In particular, if this numerical DCB
test is performed with another stacking sequence or with other
plies thicknesses, the obtained critical energy release rate will be
different, which is due to the stress monitoring of the criterion.
To avoid that, softening springs should be used to drive the delamination propagation thanks to the critical energy release rate. In
fact, the area under the stress/ displacement curves of these softening spring allows to indirectly introduce the energy release rate of
the interface (Mi et al., 1998). In the present model, this factor is
not taken into account and to avoid mesh size problem, all FE have
the same size (about 2  2  0:5 mm3 ) and the limit stress rlim was
identied with reference to this mesh size. On top of that, as mentioned above, this mesh sensitivity prevents the use of this model
for now as a predictive tool but rather as an understanding tool.

3. Experimental validation

Finally, this model was set up in the FE software Samcef and

Aboissieres impact tests were simulated. A half composite plate
is meshed with one volumic element by plies sequence of the same
orientation,
then
seven
elements
in
the
thickness
0 2 ; 45 2 ; 90 2 ; 45 4 ; 90 2 ; 45 2 ; 0 2  with a double thickness element in the middle are presented. An axial symmetry condition
around the z-axis is imposed (Fig. 8). Later in this paper, the term
ply will be used for a plies sequence of the same orientation.
The mesh size in x and y directions is constant and equals to
2  2 mm2 . The volumic elements are afterwards attached together
by the discrete interface element mentioned above. In all, the FE
model counts 135000 degrees of freedom and the calculation lasts
for about 6 h on a PC dual core 3.5 GHz with 2 Gb of RAM but this
size and this time should certainly be decreasing by optimising the
programmation. Then non-linear static numerical calculation with
an imposed displacement of the impactor is performed, since a lot
of authors (Abrate, 1998; Kwon and Sankar, 1993) showed the static/dynamic equivalence for this type of low velocity/low energy
impacts. The displacement of the impactor, which is meshed with
innite stiffness element, is imposed between 0 and 6 mm and to
compare this result with experimental test, the equivalent impact
energy is evaluated thanks to the integration of the displacement/
force curve. Then the maximum impact force versus equivalent impact energy curve is drawn in Fig. 11a and is compared to the
experimental results. We want to point out that this comparison
should be considered with caution because of the criterion mesh
sensitivity mentioned above. However there is good correlation
up to 15 J but deteriorates after, it is certainly due to the impactor
perforation phenomenon which is not simulated in the present
model. In particular the impactors perforation creates bres failures by out-of-plane shear stresses slz which are not taken into
account by the bres failure criterion. The simulation should over-

2500

12

Model with fibres failure

Sdelam (mm2)

Fmax (kN)

2000

Model without fibres failure

Experiment

Tests detailed in next figures

1500

Model with fibres failure

Model without fibres failure

1000

Experiment

500
0

0
0

10

20

Energy (J)
-a-

30

10

Fmax (kN)
-b-

Fig. 11. Maximum impact force versus energy curve (a) and delaminated area versus maximum impact force curve (b).

stress rlim mentioned above. Moreover the limit stress used for
delamination springs was identied due to this curve and the correlation allows only to verify the good identication of this
parameter.
The different damage types given by the model are summarized
in Fig. 12. At the left of this gure, the matrix cracking is drawn
from the rst ply, non-impacted side, to the seventh ply, impacted
side. These matrix cracking damages show an axial symmetry imposed by the model, then the impact point is the central point of

estimate the impact force and could explain a part of the difference
between the experiments and modelling.
Afterwards, the evolution of the delaminated area versus maximum impact force is drawn in Fig. 11b and is compared to the
experiments. Like the previous case and for the same reason, the
experimental and model curves are in accordance up to 6 kN but
there is no proper match afterwards. We want equally to point
out that this comparison should be considered with caution because it is this curve which has been used to evaluate the limit

Stacking sequence :
[0 2 , 45 2 , 90 2 , -45 4 , 90 2 , 45 2 , 0 2 ]

Ply n7
impacted side

Impactor

= 16 mm

90

Ply n1 non
impacted side

0
7 plies (with one double
ply in the middle) :

90

45

-45

Ply 7 : 0

Fib res90
failure

45

90
45

M atr ix
cracks

-45

-45
90
45

Fib res
failure

90

Interface 6:
Ply 6 : 45

(45/0 )
-45
90

45

45

M atr ix
cracks

-45

-45

Interface 5: (90/45)
90

45

Scale :
20 mm

45

Fib res
failure

Ply 5 : 90

-45

90
45

M atr ix
cracks

0
90
-45

45

-45
0

Interface 4: (-45/90)
90
45

Fib res
failure

Ply 4 : -45

Scale :
20 mm

90

-45

45

M atr ix
cracks

90

0
-45

45

-45

Interface 3: (90/-45)
Scale :
20 mm

90
45

Ply 3 : 90

Fib res
failure

90
45

-45

M atr ix
cracks

Scale :
20 mm

-45

45

90

-45

90

45

Fib res
failure

90

Interface 2: (45/90)
Ply 2 : 45

M atr ix
cracks

0
-45

45
0
90
-45
-45

45

-45
0

Interface 1: (0/45)
M atr ix
cracks

Scale :
20 mm

90

45

Fib res
failure
-45

-45

90

Ply 1 : 0
45
0

0.25
0.5
0.75
1
Fibres failure dama
ge

Delamination

Modelling

: 28J / disp. 6 mm / 1800 mm2

Fig. 12. Numerical modelling of matrix cracking, bres failures and delamination for 6 mm displacement.

damage maps. These damage maps represent the matrix cracking

indicator which equals 0 in safe material and 1 in damaged material. This indicator is evaluated for each matrix cracking spring and
a constant value in thickness of the ply is imposed due to the
adopted matrix cracking model. This damaged zone by matrix
cracking is bigger in the non-impacted side than the impacted size
due to exion effects.
In the middle of this gure, the bres failure is drawn from the
rst to the seventh ply. These damage maps represent the bres
failure indicator which equals 0 in safe material and 1 in broken
material. This indicator is evaluated at each Gauss point of volumic
elements which counts 8 Gauss points due to the used FE. The indicators value is between 0 and 1 in mean per element. This damage
is equally bigger in non-impacted side than impacted size due to
exion effects.
Finally at the right of this gure, the delamination damage is
drawn for each interface. These damage maps represent the delamination indicator which equals 0 in safe interface and 1 in delaminated interface. This indicator is evaluated for each delamination
spring and is taken equal for each group of four delamination
springs mentioned above. Afterwards, in order to compare these
delamination maps with experimental ones obtained by C-Scan,
these delamination surfaces are brought together and lled with
the colour obtained by C-Scan (Fig. 13).
In fact, to test for the relevance of a modelling, the global
delamination area or the force/impact energy curve are insufcient, the morphology of the delaminated area of each interface
gives much more information on the impact damage scenario. So
the simulated delaminations are compared to the experimental
one obtained by C-Scan at 25 J for experiment and 28 J for modelling (Fig. 13). The greatest value for the impact energy of the model
is explained by the overestimation of the impact force mentioned
above. Even if the mesh sensitivity of the delamination criterion
mentioned above prevents to wholly appreciate the quantitative
prediction of this model, in a denite way the qualitative comparison can be made and in particular the comparison of the delamination shape of each interface. The comparison between these two

delamination pictures (Fig. 13) is very good, the proposed model

allows to account for a lot of experimental observations:
Firstly the propagation of a delamination is always driven by
the direction of the lower ply (Fig. 12). In fact, it is principally
the tension in the bre direction of the lower ply which creates
the opening of the delamination fracture. This phenomenon is very
well illustrated in Fig. 14. This gure represents the displacement
eld obtained by FE calculation at 6 mm displacement for two
cross sections in 0 and 45 directions. It can be noted in these sections that only elements with bres in the cross section direction
are well represented. In fact, the other ones have no parallel sides
to the section plan which gives a wrong impression of penetration
between elements. Nevertheless, it can be seen on this picture the
large opening of the last delamination between the 0 and 45 obtained numerically thanks to a cross section in the lower ply direction (here in 0 direction) and the relatively large opening of the
previous delamination between the 45 and 90 thanks to a cross
section in the lower ply direction (here in 45 direction).
To study more particularly this phenomenon, the displacement
in z direction of the last ply is drawn versus the x-position (Fig. 15)
for a cross section along 0 and 45 directions. It clearly appears on
these curves a knee point at the location where the last (or the before last) delamination ends.
To conrm this result, an experimental test was performed in
the laboratory with a carbon/epoxy prepreg material (T700/M21)
and with the same stacking sequence. The main difference
between this material and Aboissieres lies in the resin which is a
pure epoxy resin for Aboissiere and a mix epoxy resin with thermoplastic constituents for the M21. This test is similar to Aboissieres but in static conditions to allow for the observations of
the 3D shape of the last ply thanks to an image correlation system
with two CCD cameras (Fig. 16a). That is we compare the experimental displacement eld with the model at the same imposed
displacement of 5 mm.
Afterwards the curves of z-displacement along 0 and 45 lines
were drawn (Fig. 15) and compared with the model. These curves
are in good t for the 0 cross section as for the 45 cross section

Fig. 13. Experimental and modelling delamination in the impacted (a) and non-impacted side (b).

Fig. 14. Displacement response of the model.

and in particular the knee point foresaw by modelling is clearly

observed experimentally which conrms the large opening of the
delamination and in particular of the last one simulated by the
model. So this test validates the model because its prediction
was exactly foreseeable from the experimental test.
Afterwards, other experimental curves of z-displacement along
different direction were drawn and enabled us to observe each
time two knee points which are plotted on the z-displacement eld
(Fig. 16a). A criterion on the derivate of the z-displacement curve is
used to evaluate automatically these points. Then this knee point
curve is compared to a C-Scan performed on the non-impacted side
of the laminate panel (Fig. 16b). The curve, evaluated thanks to
knee points, is in very good agreement with C-Scan and allows to
conrm the following method of delamination. Moreover it allows
to observe the delamination evolution during the experimental
test, contrary to C-Scan observation, which should normally give
interesting supplementary informations on the damage formation.
Other experiments must be carried out to conrm this method.
Moreover, even if for instance this observation is performed during
a static test, the use of two rapid cameras could allow us to perform this observation during a dynamic test. This procedure is
actually in progress.
Secondly, the delaminations do not propagate (or propagate
only slightly) in the upper ply direction of the impacted zone due
to a compression zone (Fig. 14). This phenomenon is properly illustrated by the experimental C-Scan obtained on the non-impacted
side (Fig. 13) which clearly shows the clamping of the delamination in the 45 of the impacted zone direction and it is well illustrated by the model.
Thirdly, the global direction of the projected delaminated area is
strongly directed in the 0 direction (Fig. 13). This experimental
observation is globally well taken into account by the model even
if it gives a big delamination in the 90 direction. This is perhaps
due to bres failures which are imperfectly modelled by this model. It is not easy to experimentally conrm this hypothesis because
bres failure is very difcult to observe experimentally. Nevertheless micrographic cuts or deplying (Sztefek and Olsson, 2008) could
allow to situate the bres failure and to compare them with the
modelling (Fig. 12).
Another possibility to indirectly conrm this hypothesis is to
remove the bre failures of the model. In Fig. 17, the delamina-

tion obtained without bre failure is drawn, and the curves

force/impact energy and delaminated area/maximum impact
force are equally reported on Fig. 11. The lack of bre failure increases the impact force and energy for the same imposed displacement, and for the same displacement of 6 mm the
delaminated area increases strongly from 1800 to 2850 mm2
(Fig. 17). However the delamination morphology is very different,
the projected delaminated area is not yet directed by 0 direction
but quasi-circular without the bres failure (Fig. 17). The delamination propagation in the 90 direction is fast and stopped close
to the boundary conditions. This phenomenon equally allows to
better understand why the projected delaminated area is in the
0 direction. This is due to the rectangular shape shadow which
preferentially induces bre failures in the 45 and 90 directions,
in the non-impacted side. Then, in the studied stacking sequence,
the 45 and 90 plies (plies no 2 and 3), in the non-impacted side,
break rst (the 45 ply breaks before the 90 ply due to exion
effect), and decrease stresses in these directions which stops
delamination propagation. At the same time, the last 0 ply (ply
no 1) breaks very late during the impact test modelling. Practically, this ply is still safe until 5 mm displacement, which tends
to propagate delamination in its direction. Now tests with different boundary conditions to check this result are under
development.
Finally the proper resemblance between the experimental and
numerical damage morphology allows us to conclude that the
hypothesis used in this model are globally representative of the
physical phenomena. In particular, this model can help to better
understand the phenomenon of delamination formation during
an impact test. In fact, three cases of delaminations are considered:
Pure initiation, propagation and pseudo-propagation (cf. Section
2). As already mentioned above, the rst case of pure initiation,
i.e. a delamination without preliminary damage is never found in
the different simulated impact tests. And the third case of pseudo-propagation, i.e. a delamination with preliminary matrix cracking is found only for one element (2 with the axial symmetry) of
each interface. This case is necessary in the model to initiate the
delamination but once it begins, the propagation phenomenon is
always reached before the pseudo-propagation.
This suggests two different types of matrix delamination
formation:

Cross section in 45 direction

Cross section in 0 direction

Impactor

45 Knee
point

4 mm

4 mm

Impactor

-45

0 45
90
45 90

0 Knee point

0
Fibres direction

Cross
section 45

Z Y
X

Cross
section 0

Z-displacement (mm)

0
-1
-2
-3
-4
-5
-62.5

Model :
Cross section 0
Cross section 45

Experiment :
Cross section 0
Cross section 45
-50

-37.5

-25

-12.5

12.5

25

37.5

50

62.5

X (mm)

-a-

-b-

Fig. 15. Z-displacement of the non-impacted side composite panel in the 0 and 45 directions.

Fig. 16. Knee point curve with z-displacement eld (a) and with non-impacted side C-Scan (b).

Modelling without fibres failure: 33 J / 6 mm / 2850 mm2

90

50 mm

45

Boundary conditions

0
-45

6 : 45/0
5 : 90/45
1 : 0/45
4 : -45/90
2 : 45/90

3 : 90/-45

Impacted side

 A delamination formation which corresponds to the delamination propagation: when a disjointed strip is created by preliminary matrix cracking, a zone of tension stress is created and
involves a delamination propagation in opening mode I fracture
before the matrix cracking. This scenario corresponds to the one
proposed by Renault (Fig. 1). It is particularly visible in Fig. 12
where the delaminated areas are bigger than the matrix cracked
ones.
In a real laminate panel damaged by an impact test, the delamination initiation will be present in the zone directly under the impact point and the delamination propagates far from the impact
point.
Consequently, it seems necessary to take into account these two
types of delamination formation to correctly simulate the impact
test damage.

Fig. 17. Modelling delamination without bre failure impacted side.

4. Conclusion
 A delamination formation which corresponds to the delamination initiation: when a ply develops an important matrix cracking damage with safe interfaces and when these matrix cracks
reach this interface, a delamination is initiated at this specic
point. These matrix cracks can be created by normal stresses,
which correspond to the surface bending cracks mentioned by
Choi and Chang or by shear stresses, which correspond to the
inner shear cracks mentioned by Choi and Chang (Fig. 2).

An impact damage model has been set up to simulate the different damage types forming during an impact test on laminate composite panel. The main ideas used to build this model are:
 The matrix crack damage is modelized thanks to localized damage to take into account discontinuity created by this
phenomenon.

 The delamination damage is modelized with interface type elements, and the failure criteria depends on the stresses of the
interface but also on the stresses of the adjacent layers.
 The bre failure is modelized thanks to continuum variable.
The main objective of this model is to better phycically understand the impact damage creation. In the rst step, very simple criteria written with stresses, was used in order to test the relevance
of fundamental ideas of this model. The well known problem of
mesh sensitivity of these criteria limits for now the use of this
model as a predictive tool.
Finally this model has been used to simulate experimental tests
performed by Aboissiere (2002, 2003) and shows very good agreement with experiment for the global response of the structure as
well as for the delamination morphology in each interface obtained
by C-Scan. This good result allows to justify the used hypothesis
and shows the relevance of this model. Then, it has been used to
underline different points:
 The bres failure is fundamental in the impact damage development and in particular the bres failure in a direction locks the
delamination propagation in this direction and induces it in the
perpendicular direction. Experimental conrmation must be
done to conrm this effect, for example the use of bres with
very different failure strain.
 The fracture mode I is fundamental on the delamination propagation. Indeed, the proposed model allows us to take into
account the delamination with only an interlaminar normal
stress criterion. Even if additional experiments must be done
to study the effect of the fracture mode II on delamination propagation, the present model allows to clearly conrm the predominance of the fracture mode I on this phenomenon. Of
course this result is true only for impact on this type of panel
and may be inappropriate on other congurations like for example impact on thick panel or tension on holed panel.
 The displacement eld of the last ply is inuenced by the existence of the delamination and it is possible to evaluate the delaminated area during an indentation test thanks to an image
correlation system with two CCD cameras. A correlation of this
measure with a C-Scan allowed us to validate this method and
other tests are in progress to conrm this result.
Now the fundamental ideas of this model were conrmed by
the good relevance of the damage morphology, the different criteria should be improved to obtain a really predictive tool. The rst
point to study is the use of softening springs for delamination elements to eliminate the mesh sensitivity. This mesh sensitivity will
be tested on different meshes, for delamination elements, as well
as for matrix cracking elements. Another point to focus is the
development of a new criterion for damage of bres failure to better simulate the panel perforation. Finally the fracture mode II will
be taken into account to modelize other damages types like impact
on thick plate.
Afterwards this model will be used as an initial condition to
simulate a test of compression after impact, indeed the loss of
strength in compression is directly due to impact damage. But
the initial shape of the panel after impact inuence the residual

compression strength and this model must be improved to simulate the permanent indentation after impact.
However this permanent indentation is a dominating parameter
to certify a composite structure in the eld of aeronautics. Since it
is the damage tolerance concept: The structure must withstand
ultimate loads with a permanent indentation smaller than the
BVID (Barely Visible Impact Damage).
A lot of work is still necessary to wholly simulate the damage
tolerance of a composite panel and to take into account, at the
same time, the damage during impact and the permanent indentation to evaluate the residual strength and to optimise the design of
composite structures in damage tolerance.
References
Aboissiere, J., 2003. Propagation de dommages dimpact dans un matriau
composite strati bres de carbone et rsine poxyde. Thesis of the
University Paul Sabatier of Toulouse, France.
Aboissiere, J., Michel, L., Eve, O., Barrau, J.J., 2002. Matrix cracking and delamination
under fatigue loading. ECCM10, Brugge, Belgium.
Abrate, S., 1998. Impact on Composites Structures. Cambridge University Press.
Allix, O., Blanchard, 2006. Mesomodeling of delamination: towards industrial
applications. Composites Science and Technology 66, 731744.
Aoki, Y., Iwahori, Y., Ishikawa, T., Kondo, H., Hiraoka K., 2006. Dent depth and CAI
property of CFRP laminates subjected to low velocity impact. ECCM12, Biarritz,
France.
Choi, H.Y., Chang, F.K., 1992. A model for predicting damage in graphite/expoxy
laminated composites resulting from low-velocity point impact. Journal of
Composite Materials 26 (14), 21342169.
Collombet, F., Bonini, J., Lataillade, J.L., 1996. A three dimensional modelling of low
velocity impact damage in composite laminates. International Journal for
Numerical Methods in Engineering 39, 14911516.
Davies, G.A.O., Olsson, R., 2004. Impact on composite structures. The Aeronautics
Journal 108, 541563.
De Moura, M.F.S.F., Gonalves, 2004. Modelling the interaction between matrix
cracking and delamination in carbon-epoxy laminates under low velocity
impact. Composites Science and Technology 64, 10211027.
Finn, S.R., Springer, G.S., 1993. Delaminations in composite plates under transverse
static or impact loads. Composite Structures 23, 177204.
Guinard, S., Allix, O., Gudra-Degeorges, D., Vinet, A., 2002. A 3D damage analysis of
low-velocity impacts on laminated composites. Composite Science and
Technology 62, 585589.
Hou, J.P., Petrinic, N., Ruiz, C., 2001. A delamination criterion for laminated
composites under low-velocity impact. Composite Science and Technology 61,
20692074.
Kwon, Y.S., Sankar, B.V., 1993. Indentation exure and low velocity impact damage
in graphite epoxy laminate. Journal of Composite Technology and Research 15
(2), 101111.
Ladevze, P., Lubineau, G., Marsal, D., 2006. Towards a bridge between the microand mesomechanics of delamination for laminated composites. Composite
Science and Technology 66, 698712.
Lammerant, L., Verpoest, I., 1996. Modelling of the interaction between matrix
cracks and delaminations during impact of composite plates. Composites
Science and Technology 56, 11711178.
Li, S., Reid, S.R., Zou, Z., 2006. Modelling damage of multiple delaminations and
transverse matrix cracking in laminated composites due to low velocity lateral
impact. Composites Science and Technology 66, 827836.
Mi, Y., Criseld, M.A., Davies, G.A.O., 1998. Progressive delamination using interface
elements. Journal of Composite Materials 32 (14), 12461272.
Petit, S., Bouvet, C., Bergerot, A., Barrau, J.J., 2007. Impact and compression after
impact of a composite laminate with a cork thermal shield. Composites Science
and Technology 67, 32863299.
Renault, M., 1994. Compression aprs impact dune plaque stratie carbone
poxyde Etude exprimentale et modlisation lments nis associe.
Rapport interne EADS CCR.
Sztefek, P., Olsson, R., 2008. Tensile stiffness distribution in impacted composite
laminates determined by an inverse method. Composites Part A 39 (8), 1283
1293.
Whitney, J.M., Nuismer, R.J., 1974. Stress fracture criteria for laminated composites
containing stress concentrations. Journal of Composite Materials 8, 253265.