Donadon2009.pdf

Hindawi Publishing Corporation
International Journal of Aerospace Engineering
Volume 2009, Article ID 486063, 22 pages
doi:10.1155/2009/486063
Research Article
A Three-Dimensional Ply Failure Model for Composite Structures
Maurı́cio V. Donadon, Sérgio Frascino M. de Almeida,
Mariano A. Arbelo, and Alfredo R. de Faria
Department of Mechanical Engineering, Aeronautical Institue of Technology (ITA), CTA-ITA-IEM, Praça Mal. Eduardo Gomes 50,
12228-900 São José dos Campos-SP, Brazil
Correspondence should be addressed to Maurı́cio V. Donadon, [email protected]
Received 21 October 2008; Revised 28 January 2009; Accepted 23 June 2009
Recommended by Ever Barbero
A fully 3D failure model to predict damage in composite structures subjected to multiaxial loading is presented in this paper. The
formulation incorporates shear nonlinearities effects, irreversible strains, damage and strain rate effects by using a viscoplastic
damageable constitutive law. The proposed formulation enables the prediction of failure initiation and failure propagation by
combining stress-based, damage mechanics and fracture mechanics approaches within an unified energy based context. An
objectivity algorithm has been embedded into the formulation to avoid problems associated with strain localization and mesh
dependence. The proposed model has been implemented into ABAQUS/Explicit FE code within brick elements as a userdefined
material model. Numerical predictions for standard uniaxial tests at element and coupon levels are presented and discussed.
Copyright © 2009 Maurı́cio V. Donadon et al. This is an open access article distributed under the Creative Commons Attribution
License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly
cited.
1. Introduction
Damage in composite structures is a very complex phenomenon which can occur through a different number of
failure mechanisms such as fibre breakage, fibre buckling,
matrix cracking, fibre-matrix debonding and delamination
either combined or individually. The increasing computational resources have allowed reliable prediction of such
phenomenon with a certain degree of accuracy by using the
finite element method. However, there is still a lot of work
to be done in this field to better understand the physics
of failure in order to improve numerical failure models for
composite materials. The damage modelling in composites
can be broadly divided into four different approaches:
(i) failure criteria approach,
(ii) fracture mechanics approach,
(iii) plasticity approach,
(iv) damage mechanics approach.
Failure criteria approaches were initially developed for
unidirectional materials and restricted to the static regime.
They are divided into two categories: interactive and noninteractive criteria [1]. Noninteractive criteria assume that
the failure modes are decoupled and specific expressions
are used to identify each failure mechanism. In the stress
based criteria, for example, each and every one of the stresses
in the principal material coordinates must be less than
the respective strengths; otherwise, fracture is said to have
occurred. In a similar fashion, the maximum strain failure
criteria states that the failure occurs if one of the strains in the
principal material coordinates exceed its respective failure
strain.
On the other hand, interactive criteria assume an
interaction between two or more failure mechanisms and
they describe the failure surface in the stress or strain
space. Usually stress or strain polynomial expressions are
used to describe the boundaries for the failure surface or
envelope. Any point inside the envelope shows no failure
in the material. Several interactive failure criteria can be
found in the open literature, such as the Tsai and Wu [2]
criterion, Tsai-Hill criterion, Hoffman criteria among others
[1]. Nevertheless, the disadvantage in using such polynomial
criteria is that they do not say anything about the damage
mechanisms themselves, therefore modified versions have
been used to distinguish between failure modes. Hashin [3]
proposed a three-dimensional failure criteria for unidirectional composites. In his model four distinct failure modes
2
associated with fibre failure in tension or compression and
matrix cracking in tension or compression are modelled
separately. Engblom and Havelka [4] proposed the use of a
combination of the Hashin [3] and Lee [5] failure criteria.
They used the Hashin criterion to detect in-plane failure and
the Lee criteria to predict delaminations. The degradation
was performed by reducing the stresses associated with
each failure mechanism to zero in the composite material
constitutive law.
Shivakumar et al. [6] used the Tsai-Wu failure criterion
and maximum stress criterion to model low-velocity impact
damage in composites. In this approach the Tsai-Wu criterion was used in order to determine whether the damage
has occurred and the maximum stress criteria was used to
identify the failure mode. Good agreement was obtained
between computed and experimental damage areas.
F.-K. Chang and K.-Y. Chang [7, 8] proposed a progressive in-plane damage model for predicting the residual
strength of notched laminated composites. Three in-plane
failure modes are considered: matrix cracking, fibre-matrix
shearing and fibre breakage. Their model is currently available in the material model library of the LS-DYNA3D explicit
finite element code. In their model the shear stress-shear
strain relation is assumed to be nonlinear and an expression
proposed by Hahn and Tsain [9, 10] is used to represent
composite behaviour in shear. For fibre breakage and/or
fibre-matrix shearing the degree of property degradation
within the damaged area depends on the size of the damage
predicted by the fibre failure criterion. F.-K. Chang and
K.-Y. Chang [7] proposed a property reduction model for
fibre failure based on the micromechanics approach for fibre
bundle failure. It is postulated that for fibre failure both E y
and ν yx are reduced to zero, but Ex and the shear modulus
Gxy degenerate according to the Weibull distribution.
Choi et al. [11] investigated the low velocity impacts
on composite plates using a line-nose impactor. They used
the Chang-Chang failure criterion to detect the initiation of
matrix cracking and delamination. Good agreement between
experimental and numerical results was achieved.
Davies and Zhang [12] studied the low-velocity impact
damage in carbon/epoxy composite plates impacted by
an hemispherical impactor. The plates tested were quasiisotropic having a symmetrical lay-up. Different impact
energy levels, thickness, dimensions and boundary conditions were considered. The finite element modelling was
carried out using plate elements.
The Chang-Chang failure criteria were used for the inplane damage predictions and in general good agreement
was obtained between the computed force histories and
experiments.The Chang-Chang failure criteria have shown a
reasonable performance for in plane damage predictions in
composite structures and therefore, their applications have
been widely reported in the literature [13–18].
The prediction of the onset of delamination depends on
the interlaminar stress state and interlaminar strength of
the laminate. Kim and Soni [19] used the distribution of
interlaminar normal stress, and averaged that stress along
a ply thickness distance in all cases they studied. They
assumed that failure occured when the average of the normal
International Journal of Aerospace Engineering
interlaminar stress value over the fixed distance reached the
interlaminar tensile strength.
Jen et al. [20] developed a model based on boundary layer
theory to predict initiation and propagation of delamination
in a composite laminate containing a central circular hole.
The Hashin-Rotem failure criterion was adopted in their
model to predict the loading and location at which the
initiation of delamination occurs.
Brewer and Lagace [21] proposed a quadratic stress
criterion for initiation of delamination using the approach
suggested by Kim and Soni [19]. According to their criterion
only out-of-plane stresses contribute to delamination and
they assume that the predicted stresses should be independent of the sign of the interlaminar shear stress. A similar
criteria was also proposed by Liu et al. [22] to predict
delamination in composite laminates.
The disadvantage in using stress based criteria for
composite materials is that the scale effects relating to the
length of cracks subject to the same stress field cannot be
modelled correctly [13, 23]. In the failure criteria approach
either the position and size of the cracks are unknown. For
these reasons the fracture mechanics approach may be more
attractive. Fracture mechanics considers the strain energy at
the front of a crack of a known size and compares the energy
with critical quantities such as critical strain energy release
rate.
The fracture mechanics approach has been used to
predict compression after impact strength of composite
laminates [24–27]. In such models the damaged area is
replaced by an equivalent hole and the inelastic deformation
associated with fibre microbuckling that develops near the
hole edge is replaced with a equivalent crack loaded on
its faces by a bridging traction which is linearly reduced
with the crack closing displacement. The diameter of the
hole is obtained from X-radiographs and/or ultrasonic Cscan images. The results showed good correlation between
analytical and experimental values.
Another potential application of the fracture mechanics
approach is its indirect use to predict progressive delamination in composites [28–31]. In such models stressdisplacement constitutive laws describe the interfacial material behaviour and fracture mechanics concepts are used. The
area defined by the constitutive relationship is equal to the
fracture energy or energy release rate and once the stresses
have been reduced to zero the fracture energy has been
consumed and the crack propagates. Linear and quadratic
interaction relationships were assumed to describe the crack
propagation in mixed-mode delamination. Comparisons
were made with experimental and closed-form results and
good agreement was obtained.
However, the fracture mechanics approach cannot be
easily incorporated into a progressive failure methodology
because its application requires an initial flaw. A possible
solution is to use a hybrid approach by using a stress or
strain-based criterion for the failure initiation and a fracture
mechanics approach for the failure propagation.
The plasticity approach is suitable for composites
that exhibit ductile behaviour such as boron/aluminium,
graphite/PEEK and other thermoplastic composites.
International Journal of Aerospace Engineering
Vaziri et al. [32] proposed an orthotropic plane stress
material model that combines the classical flow theory
of plasticity with a failure criterion. In their work the
material constitutive law is assumed to be elastic-plastic
and it has two stages. The first stage is the post-yield and
pre-failure where an orthotropic plasticity model is used
to model the nonlinear material behaviour. The second
stage is the postfailure where brittle or ductile failure modes
start to occur. Favourable agreement is obtained between
experiment and the model results.
The damage mechanics approach has been investigated
by many researchers in recent years and its application to
damage modelling in composites has shown to be efficient.
The method was originally developed by Kachanov [33] and
Rabotnov [34] and it has the potential to predict different
composite failure modes such as matrix cracking, fibre
fracture and delamination.
Ladeveze and Dantec [35] proposed an in-plane model
based on damage mechanics to predict matrix microcracking
and fibre/matrix debonding in unidirectional composites.
Two internal damage variables were used to degrade the
ply material properties, one of them associated with the
transverse modulus and another with the in-plane shear
modulus. A linear elastic-damage behaviour was assumed for
tensile and compressive stresses and a plasticity model was
developed to account for the inelastic strains in shear. Damage evolution laws associated with each failure mechanism
were introduced which relate the damage variables to strain
energy release rates in the ply. Tension, compression and
cyclic shear tests were performed to determine the constants
required in the damage-development laws. Comparisons
with experiments were performed by the authors and a good
correlation between numerical and experimental results was
obtained.
Johnson [36] applied the model suggested by Ladeveze
and Dantec for the prediction of the in-plane damage
response of fibre reinforced composite structures during
crash and impact events. The damage model was implemented for shell elements into the PAM-CRASH explicit
finite element code. An experimental programme was carried
out in order to validate the model. Low velocity impact
tests were performed using a drop test rig. The plates were
simply supported on a square steel frame and they were
impacted at their centre using a hemispherical impactor
with 50 mm diameter. The mass of the impactor was 21 kg.
The plates were fabricated from 16 plies of carbon/epoxy
with a quasi-isotropic layup. Different energy levels were
considered to give a range of different failure modes
from rebound to full penetration. Based on force history
comparisons between experiments and numerical results the
model overpredicted the peak load. According to Johnson,
the discrepancy between measured and predicted peak load
is explained by the neglected delamination in the model. This
fact was demonstrated by the author in his subsequent work
[37] where the delamination was included in the damage
modelling using contact interface conditions.
Williams and Vaziri [38] implemented an in-plane
damage model based on continuum damage mechanics
(CDM) into LS-DYNA3D for impact damage simulation in
3
composite laminates. The model was originally developed
by Matzenmiller et al. [39]. The model considers three
damage parameters: fibre failure damage parameter, matrix
failure damage parameter and an extra damage parameter
to account for the effect of damage in shear response.
Individual stress based criteria for each failure mechanism
are used for the damage initiation. The failure criterion
defines certain regions in stress (or strain) space where
the damage state does not change. The damage growth
law adopted by Matzenmiller is a function of the strain
and it assumes an exponential form. The stress/strain curve
predicted by this damage function is a Weibull distribution
that can be derived from statistical analysis of the probability
of the failure of a bundle of fibres with initial defects. Impact
simulations with different energy levels were performed and
the performance of the proposed model was checked against
the Chang-Chang failure criteria and experimental results.
Some limitations of the model were pointed out by the
authors. Firstly, the response is predicted by a single equation
and, as a result the loading and postfailure responses cannot
be separated, thus restricting the versatility of the model.
Also, the dependence of the damage growth on the loading
rate as well as mesh size.
In his recent work Williams et al. [40] proposed a planestress continuum damage mechanics model for composite
materials. The model was implemented for shell elements
into LSDYNA3D explicit finite element code. Also, the
model is an extension of their previous work [41] and
it deals with some issues related to the limitation of the
model suggested by Matzenmiller et al. [39]. The approach
adopted by the authors uses the sub-structuring concept
where one integration point is used for each sub-laminate
and composite laminate theory is used to obtain its effective
material properties. The model assumes a bilinear damage
growth law and the damage process has two phases, one
associated with matrix/delamination damage and another
with fibre breakage. The driving force for damage growth is
assumed to be a strain potential function and the threshold
values for each damage phase are experimentally determined. High-velocity and low-velocity impact simulations
were performed in order to assess the performance of the
model. Also, the results were compared with the model
proposed by Matzenmiller et al. [39] and the existing ChangChang failure criteria. Pinho et al. [42] proposed a threedimensional failure model to predict failure in composites.
Their model is an extension of the plane-stress failure model
proposed by Davila and Camanho [43] and it accounts for
shear nonlinearities effects. The model was implemented
into LS-DYNA3D finite element code and the authors
obtained good correlation between numerical predictions
and experimental results for static standard in-plane tests.
However, neither strain rate effects nor strain localization
problems associated with irregular mapped meshes were
addressed by the authors.
This paper presents a formulation for a threedimensional ply failure model for composite laminates.
The proposed model is an extension of the authors previous
work [44]. The new version of the model incorporates
strain rate effects in shear by means of a semiempirical
4
International Journal of Aerospace Engineering
viscoplastic constitutive law and combines a quadratic stress
based criterion with the Mohr-Coulomb failure criterion
to predict interfibre-failure (IFF) without knowing a priori
the orientation of the fracture plane. The new version of
the model also uses the Second Piola-Kirchhoff stress tensor
which potentially avoids material deformation problems.
High-order damage evolution laws are proposed to avoid
both numerical instabilities and artificial stress waves
propagation effects commonly observed in the numerical
response of Finite Element Codes based on Explicit Time
Integration Schemes. These laws compared to the widely
used bilinear law ensure smoothness at damage initiation
and fully damaged stress onsets leading to a more stable
numerical response. Iterative energy based criteria have
been also added to the previous version [44] to better
predict damage under multiaxial stress states. The model
also incorporates a three dimensional objectivity algorithm
[45] which accounts for orthotropic and crack directionality
effects, making the model to be completely free of strain
localization and mesh dependence problems.
system. The compliance matrix C −1 is defined in terms of the
material axes as
C −1
⎡ 1
−ν21 −ν31
0
0
⎢ E11
E22 E33
⎢
⎢ ν12
1
−ν32
⎢−
0
0
⎢ E11 E22
E33
⎢
⎢ ν13 −ν23
1
⎢−
0
0
⎢ E
E22 E33
11
=⎢
⎢
1
⎢ 0
0
0
0
⎢
G12
⎢
⎢
1
⎢ 0
0
0
0
⎢
G23
⎢
⎣
0
0
0
0
νi j = νxi xj
Eii = Exi .
This section presents the failure model formulation. The
formulation is based on the Continuum Damage Mechanics
(CDM) approach and enables the control of the energy
dissipation associated with each failure mode regardless of
mesh refinement and fracture plane orientation by using a
smeared cracking formulation. Internal thermodynamically
irreversible damage variables were defined in order to
quantify damage concentration associated with each possible
failure mode and predict the gradual stiffness reduction
during the fracture process. The material model has been
implemented into ABAQUS explicit finite element code
within brick elements as an user defined material model.
2.1. 3-D Orthotropic Stress-Strain Relationship. The orthotropic material law that relates second Piola-Kirchhoff stress
S to the Green-St. Venant strain e in the local material
coordinate system is written as
S = T t CTe,
(1)
where e = [e11 , e22 , e33 , e12 , e23 , e31 ] is the Green-St. Venant
strain vector and T is the transformation matrix given by
k2
⎢
⎢ 2
⎢ l
⎢
⎢
⎢ 0
⎢
T=⎢
⎢ 0
⎢
⎢
⎢ 0
⎢
⎣
−kl
l
2
0 0 0
k2 0 0 0
0 1 0 0
0 0 k −l
0 0 l k
kl 0 0 0
2kl
⎤
⎥
⎥
−2kl ⎥
⎥
⎥
0 ⎥
⎥
⎥,
0 ⎥
⎥
⎥
0 ⎥
⎥
⎦
2
2
k −l
(2)
where k = cos(β) and l = sin(β). β is the inplane rotation
angle around the Z-direction which defines the orientation
of the material axes with respect to the reference coordinate
⎥
⎥
0 ⎥
⎥
⎥
⎥
0 ⎥
⎥
⎥,
⎥
0 ⎥
⎥
(3)
⎥
⎥
0 ⎥
⎥
⎥
1 ⎦
G31
where the subscripts denote the material axes, that is,
2. Formulation
⎡
0
⎤
0 ⎥
(4)
Since C is symmetric,
νi j
ν ji
=
.
Eii
Ej j
(5)
The Green-St. Venant strain tensor e is given by
e=
1 t
F F−I ,
2
(6)
where F is the deformation gradient and I is the identity
matrix. The Cauchy stress tensor σi j can be determined in
terms of the second Piola-Kirchhoff stress tensor S as follows:
σi j = J −1 FSFt ,
(7)
where J is defined as the Jacobian determinant which is
the determinant of the deformation gradient F. The present
formulation will predict realistic material behaviour for finite
displacements and rotation as long as the strains are small.
2.2. Degraded Stresses. The degraded or damaged stresses are
defined as stresses transmitted across the damaged part of
the cross-section in a representative volume element of the
material. Based on the isotropic damage theory as originally
proposed by Kachanov [33], an orthotropic relationship
between local stresses acting on the damaged configuration
can be written in terms of the local effective stresses in the
undamaged configuration at ply level as follows:
International Journal of Aerospace Engineering
⎧ ⎫
⎪
σd ⎪
⎪
⎪
⎪
⎪ 1⎪
⎪
⎪
⎪
⎪
d⎪
⎪
⎪
⎪
⎪
σ
⎪
⎪
2
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
d
⎨ σ3 ⎪
⎬
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
=⎢
⎪
⎪
d
⎢
⎪
⎪
τ
⎪
⎢
12 ⎪
⎪
⎪
⎪
⎪
⎢
⎪
⎪
⎪
⎢
⎪
⎪
d ⎪
⎪
⎪
τ23 ⎪
⎪
⎪ ⎢
⎪
⎪
⎪
⎪
⎪ ⎣
⎪
⎩ d⎪
⎭
τ13
(1 − d1 (1 ))
0
t t
1 − dm
m
0
5
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
⎪
⎪ dc ⎪
⎥⎪
⎪
⎪σ m ⎪
⎥⎪
⎪
⎪
2 ⎪
0
0
0⎥⎪
⎪
⎪
⎪
⎥⎪
⎪
c ⎪
⎪
⎥⎪
d
⎪
m⎪
⎪
⎬
⎨ σ3 ⎪
0
0
0⎥
⎥
⎥ dc ,
⎪
⎪
t
⎥
m
t 1 − dm
1 − d12 γ12
0
0⎥⎪
⎪
⎪ τ12 ⎪
m
⎪
⎪
⎪
⎪
⎥
⎪
⎪
⎥⎪
c ⎪
⎪
⎪
d
t
m
t
⎪
⎪
⎪
⎪
1 − dm m
0⎥
0
τ
⎪
⎪
⎦⎪
⎪
⎪ 23 ⎪
⎪
c ⎪
⎪
d
⎭
0
0
1 ⎩τ m ⎪
(i) It accounts for fibre damage effects due to tensile
and compression loadings by means of the internal
variable d1 .
(ii) It provides the coupling between shear induced
damage and matrix cracking during the degradation
t
process using the internal damage variables d12 , dm
c
and dm .
(iii) It accounts for damage coupling between the normal
stress σ2 , out-of-plane shear stress τ23 (transverse
shear cracking, which leads to delaminations) and
inplane shear stress τ12 for matrix cracking predictions.
2.3. Fibre Failure. The failure index to detect fibre failure in
tension is given by
σ1
≥ 1.
Xt
(9)
In order to detect the catastrophic failure in compression
related to the total instability of the fibres, the maximum
stress criteria is used to detect damage initiation:
F1c (σ1 ) =
| σ1 |
Xc
≥ 1,
(8)
13
t , and d
where d1 , dm
12 are internal variables introduced to
quantify the damage concentration within the Representative
Volume Element (RVE) [46, 47]. A detailed description on
the physical meaning and definition of these variables will
be given in the following sections. The stress components
c are the stresses degraded locally on
with the superscript dm
the action fracture plane. Details about the determination of
the local action plane and the local degradation procedure is
given in Section 2.4.
The proposed orthotropic relationship between
degraded and intact stresses ensures that the material
stiffness matrix is positive defined during the degradation
process. Moreover, it is physically based. By inspecting
carefully (8) one can notice that transverse compression
c degrades the through-the-thickness stresses
damage dm
σ3 , τ23 and τ13 . From the physical point of view, this
relationship incorporates the the main features observed in
the experiments, such as follows:
F1t (σ1 ) =
⎫
⎧
⎤ σ1 ⎪
0 ⎪
⎪
⎪
⎪
⎪
0
(10)
where Xt and Xc are the longitudinal strengths in tension
and compression, respectively. When one of criteria given
above is met, damage commences and grows according to the
damage evolution law proposed by Donadon et al. [44]:
d1 (1 ) = d1t 1+ + d1c 1− − d1t 1+ d1c 1− ,
(11)
where d1t (1+ ) and d1c (1− ) are the contributions of the
irreversible damage due to fibre breakage in tension and fibre
kinking in compression, respectively:
d1t 1+ = 1 −
d1c
t 1,0
2
1 + κ1,t
1+ 2κ1,t 1+ − 3 ,
1+
c −
− 1,0
2
1 = 1 − − 1 + κ1,c
1 2κ1,c 1− − 3
1
(12)
with
κ1,t 1+ =
t
1+ − 1,0
t
t
1,
f − 1,0
,
c
− − 1,0
κ1,c 1− = c1
c ,
1, f − 1,0
(13)
t
c
where 1,0
and 1,0
are the failure strains in tension and
t
) and 1− =
compression, respectively. 1+ = max(1 (t), 1,0
c
max(|1 (t)|, 1,0 ) are the maximum achieved strains in the
strain time history in tension (σ1 > 0) and compression
t
c
and 1,f
are the final strains in
(σ1 < 0), respectively. 1,f
tension and compression which are written as a function
of the tensile fibre breakage and compression fibre kinking
fracture toughnesses, respectively, as follows,
t
1,
f =
c
1,
f
2Gtfibre
,
Xt l ∗
2Gcfibre
=
Xc l ∗
(14)
where Gtfibre and Gcfibre are the intralaminar fracture toughnesses associated with fibre breakage in tension and compression, respectively and l∗ is characteristic element related
to the size of the process zone. Experimental procedures
and data reduction schemes to characterise the intralaminar
fracture toughness for composites can be found in [48].
The proposed damage evolution law d1 (1 ) results in a
thermodynamically consistent degradation procedure, regularizing damage and combining stress, damage mechanics
and fracture mechanics based approaches within an unified
way similarly to the approach proposed by Miami et al. [49]
6
International Journal of Aerospace Engineering
2.4. Interfibre Failure (IFF). The interfibre failure modes
consist of transverse matrix cracking either in tension or
compression. Based on the Worldwide Failure Exercise
experimental results [50] Pinho et al. [42, 51] found that
the failure envelope defined between the transverse stress
σ2 and in-plane shear stress τ12 is accurately described by a
quadratic interaction criteria. Thus, a failure index based on
the interactive quadratic failure criterion given in [42, 51] has
been used to predict tensile transverse matrix cracking. For
tensile matrix cracking a failure index based on an interactive
quadratic failure criterion written in terms of tensile and
shear stresses is proposed in the following form:
F2t (σ2 , τ23 , τ12 )
σ2
=
Yt
2
τ
+ 23
S23
2
τ
+ 12
S12
2
≥ 1.
tm,0 tm
2
1 + κm,t
tm
2κm,t tm − 3
σ2t = σ tm cos(θ),
where θ is the angle defined between the resultant stress
and tensile normal stress in the transverse direction, that is,
θ = a cos[max(0, σ2t )/σ tm ]. In a similar way the tensile and
resultant shear strain components are given as follows:
2t = tm cos(θ),
where
tm
tm − tm,0
tm, f − tm,0
For the damage evolution law given by (16) the specific
fracture energies associated with tensile and shear stresses are
respectively given by
(16)
gmt =
gms =
,
(17)
is the defined as resultant strain which is given by,
tm = 22 + γ2s,m
with
(18)
2
2
γs,m = γ12
+ γ23
,
tm,0
(19)
σ tm,0
and
are the damage onset resultant strain and
stress, respectively, that is,
tm,0
= tm F2t =1
,
(20)
σ tm,0 = σ tm F2t =1 .
In order to account for damage irreversibility effects
tm = max(tm (t), tm,0 ) must be used in (16), where tm is the
maximum achieved resultant strain in the strain time history.
The derivation of the resultant failure strain tm,f associated
with tensile/shear matrix cracking is based on a power law
criterion, which accounts for interactions between energies
per unit of volume of damaged material within the RVE
subjected to tensile and shear loadings. The power law energy
criterion is given in the following form,
gmt
gmt c
λ
gs
+ sm
gm c
λ
=1
(21)
with λ = 1 for UD laminates. The resultant stress in the
transverse direction (or matrix direction) due to combined
tensile and shear loadings is given by
2
2
2
σ tm = σ22 + τ12
+ τ23
= σ22 + τ tm .
(22)
(24)
γts,m = tm sin(θ).
with
κm,t tm =
(23)
τ t = σ tm sin(θ),
(15)
Once the criterion above is met, the proposed expression
for damage growth due to tensile matrix cracking is given by
t
tm = 1 −
dm
The tensile and resultant shear stress components can be
written in terms of the resultant stress as follows:
2f
0
γf
σ2t d2 =
s,m
0
σ tm,0 tm, f cos2 (θ)
,
2
σ tm,0 tm, f sin2 (θ)
τ t dγ =
,
2
(25)
where the area under the stress-strain curves defined by the
proposed polynomial damage evolution laws is identical to
the one defined by the widely used bilinear softening law
given in [42, 44]. Substituting (25) into (21) we obtain the
following expression for the final strain due to the combined
tensile and shear stress state,
⎡
tm, f
2
cos2 (θ)
= t ⎣
σ m,0
gmt c
λ
sin2 (θ)
+
gms c
λ ⎤−1/λ
⎦
(26)
where gmt c and gms c are the critical specific fracture energies.
These energies are related to the intralaminar fracture toughnesses. By using a smeared cracking formulation [52] and
assuming that for UD laminates the values of intralaminar
toughnesses associated with tensile matrix cracking and
shear matrix cracking are comparable with mode I and mode
II interlaminar fracture toughnesses, a relationship between
specific critical fracture energies and intralaminar fracture
toughnesses can be written as follows:
gmt c =
gms c
GI c
,
l∗
GII
= ∗c ,
l
(27)
where l∗ is the characteristic length associated with the
length of the process zone for each particular failure
mode. A detailed description about the characteristic length
calculation will be presented in the following section.
The failure index to detect matrix cracking in compression failure is based on the criterion proposed by Puck and
Schürmann [53, 54]. Their criterion is based on the MohrCoulomb theory and it enables the prediction of fracture
planes for any given stress state related to Interfibre-Failure
International Journal of Aerospace Engineering
7
(IFF) Modes. This criterion is currently the state of the art
to predict transverse compression response of composite
laminates. The failure index to detect matrix cracking in
compression based on the failure criterion proposed by Puck
and Shürmann [53, 54] can be written as follows,
F2c (σnt , σnl ) =
σnt
SA23 + μnt σnn
2
+
σnl
S12 + μnl σnn
2
≥ 1,
(28)
here the subscripts n, l and t refer to the normal and tangential directions in respect to the fracture plane direction. S12
is the inplane shear strength and SA23 is the transverse shear
strength in the potential fracture plane (Action Plane), which
is given by [55],
SA23 =
Yc 1 − sin φ
2
cos φ
(29)
also confirmed by Puck and Schurmann [53, 54]. However,
θf ≈ 530 is only valid for uniaxial compression loading.
This implies that θf changes for different stress states and
alternatives are needed in order to handle such a problem.
By examining (28) it is possible to see that the criterion has
also the potential of predicting transverse intralaminar shear
cracking for values of θf different from 530 . The transverse
intralaminar shear cracking is a very important failure mode
because it leads to delamination between adjacent layers. In
order tackle this problem we have used an iterative procedure
to compute the fracture plane orientation for a given stress
state. The procedure consists of incrementally varying θf
within the interval [−900 ≤ θf ≤ 900 ] for a given stress
state defined at ply level and check if the failure index for
matrix cracking in compression being reached. Once the
failure index is reached the local shear components σnl and
σnt acting on the candidate fracture plane are degraded to
zero according to the following damage evolution law:
with
φ = 2θ f − 900 ,
(30)
where Yc is the transverse compression strength. The fracture
angle can be determined either experimentally or alternatively using (30), where θ f maximises the failure criterion.
Following the Mohr-Coulomb failure theory the friction
coefficients can be determined as a function of the material
friction angle as follows:
μnt = tan φ = tan 2θ f − 900 .
(32)
The stress components acting on the potential fracture
plane are written in terms of angle θ f which defines the
orientation of the fracture plane in respect to throughthe-thickness direction (direction-3 in the local material
coordinate system):
σnt θ f = −σ2 mn + σ3 mn + τ23 2m2 − 1 ,
σtt θ f = σ2 1 − m2 + σ3 m2 − 2mnτ23 ,
(33)
σnl θ f = τ12 m + τ13 n,
where m = cos(θf ) and n = sin(θf ). θf is the rotation around
the local fibre direction (direction-1 in the local material
coordinate system). It is clear that in order to apply (28)
θf must be known. Many authors have defined θf ≈ 530
based on experimental results for standard compression tests
in UD laminates [42, 44, 51]. This is true and has been
(34)
cm − cm,0
cm, f − cm,0
,
(35)
where cm is the resultant shear strain on the action plane
which is defined as
2
2
cm = nl
+ nt
,
2
σ cm = σnl2 + σnt
(36)
with
nl θ f = γ12 m + γ13 n,
nt θ f = −2 mn + 3 mn + γ23 2m2 − 1 ,
(37)
where cm,0 and σ cm,0 are the damage onset resultant strain and
stress, respectively, that is,
(38)
The resultant final strain for transverse compression
failure is defined in terms of mode II interlaminar fracture
toughness as follows:
σlt θ f = −τ12 n + τ13 m,
2
1 + κm,c
cm 2κm,c cm − 3
cm,0 = cm F2c =1 ,
σ cm,0 = σ cm F2c =1 .
σnn θ f = σ2 m2 + σ3 1 − m2 + 2τ23 mn,
cm
κm,c cm =
(31)
μnl
μnt
=
S12 .
SA23
cm,0 with
In the absence of experimental values an orthotropic
relationship for the friction coefficients can be used [54],
c
cm = 1 −
dm
cm, f =
2Gcmatrix
2GII
= c c∗ .
σ cm,0 l∗
σ m,0 l
(39)
In order to account for damage irreversibility effects
cm = max(cm (t), cm,0 ) must be used in (34), where cm is the
maximum achieved resultant strain on the action plane in the
strain time history. After degrading the shear stresses acting
on the potential fracture angle, the stresses are rotated back
8
International Journal of Aerospace Engineering
to the local material coordinate system using the following
transformation:
dc
c
cm
σ2 m = m2 σnn − 2mnσnt 1 − dm
dc
=
c
mσnl 1 − dm
c
dm
2
τ13 = mnσnn + 2m − 1
c
σnt 1 − dm
dc
c
cm
τ23m = nσnl 1 − dm
γ12,0 = γ12 S12 ,
+ σtt m2 ,
(40)
c m − nmσtt ,
τ12 = αG12 γ12
(41)
G12 = G012 + c1 (e−c2 γ12 − 1),
(42)
with
where G012 is the initial shear modulus and c1 , c2 are material
constants obtained from static/quasistatic in-plane shear
tests. α is the strain-rate enhancement given by the following
law,
α = 1 + e(γ̇12 /c3 )
(43)
where c3 is another material constant obtained from dynamic
in-plane shear tests. By decomposing the total shear-strain
i
e
and γ12
elastic components, the inelastic
into inelastic γ12
shear-strain can be written in terms of the elastic and total
strain components as follows:
=
e
γ12 − γ12
τ12 γ12
= γ12 −
.
G012
(44)
The failure index for in-plane shear failure is based on the
maximum stress criterion and it is given by:
F12 (τ12 ) =
|τ12 |
S12
≥ 1.
(45)
The proposed damage evolution law for in-plane shear
failure is given by
d12 γ12
i
γ12,0 − γ12,0
2
=1−
1 + κ12
γ12 2κ12 γ12 − 3
i
γ12 − γ12,0
(46)
with
κ12 γ12 =
γ12, f is written in terms of the intralaminar toughness in
shear:
γ12, f =
+ mσlt .
2.5. In-Plane Shear Failure. The observed behaviour of glass
and carbon fibres laminates generally shows marked rate
dependence in matrix-dominated shear failure modes and
for this reason a rate dependent constitutive model has been
used to model the in-plane shear behaviour. The constitutive
model formulation is based on previous work carried out
by Donadon et al. [44, 56] and it accounts for shear
nonlineatities, irreversible strains and damage within the
RVE. The stress-strain behaviour for in-plane shear failure
is defined as follows,
i
γ12
(48)
i
i γ12,0
= γ12
,
S
12
c m − nσlt ,
+ σtt 1 − m2 ,
c
cm
σ3 m = σnn 1 − m2 + 2mnσnt 1 − dm
dc
τ12m
i
are the total strain and total inelastic
where γ12,0 and γ12,0
strain at failure (τ12 = S12 ), respectively, that is,
i
γ12 − γ12,0 − 2γ12,0
,
i
γ12,0 − γ12,0
− γ12, f
(47)
2Gshear
.
S12 l∗
(49)
In the absence of experimental results for Gshear , is
reasonable to assume Gshear = GIIc for UD plies, where GIIc is
the mode II interlaminar fracture toughness.
2.6. Objectivity Algorithm. The smeared formulation
described in the previous sections relates the specific
energy within a Representative Volume Element (RVE)
with the fracture energy of the material for each particular
failure mode. Since finite elements are volume based, mesh
dependency problems will arise as a result of the mesh
refinement. The correction of the postfailure softening slope
according to the finite element size, as reported by Bazant
[52] seems an attractive solution for the problem. However
the approach has some limitations. Firstly, the crack growth
direction must be parallel to one edge of the finite element,
which is not the case for multidirectional composite
laminates where layers can have arbitrary directions.
Secondly, it cannot handle nonstructured meshes required
in most of the complex finite element models with geometric
discontinuities. In order to overcome such limitations and to
ensure the objectivity of the model for generalized situations,
a methodology originally developed by Oliver [57] has been
used and extended to handle composite layers [58]. The
dependence of the characteristic length on the fracture
energy as well as its mathematical expression, were derived
based in the work proposed by Oliver [57]. The method
ensures a constant energy dissipation regardless of mesh
refinement, crack growth direction and element topology so
that, it is still applicable to nonstructured meshes.
2.6.1. Crack Modelling in the Continuous Medium. Imagine a
singular line in a two dimensional domain as a continuous
material line, across which displacements are continuous but
displacement gradients are discontinuous. The condition for
a point belonging to a singular line with unit normal n at this
point is that the determinant of the acoustic tensor in the n
direction be zero, that is [57],
det ni Ci jkl nl = 0.
(50)
Nonpositive materials bifurcate, producing singular lines and
the equation above permits their direction to be determined
at each point. In the context of standard finite elements of
C 0 continuity, a singular line can be modelled only by the
sides of the elements, these being the only points in the mesh
where displacement gradient discontinuities can be obtained.
International Journal of Aerospace Engineering
9
However, a crack produces not only displacement gradient
discontinuities but also displacement discontinuities. This
latter kind of discontinuity cannot be modelled by a C 0 finite
element mesh for finite levels of discretization. However, a
displacement discontinuity can be modeled as the limit of
two parallel singular lines Γ− and Γ+ which tend to coincide
with each other. The band delimited by these lines is known
as singular band, and h is its width.
By assuming an orthogonal curvilinear coordinate system (x , y ) in the interior of the band, where y coordinates
lines are parallel to the singular lines Γ− and Γ+ , and x are
the straight coordinates lines. Let u+ (y ) and u− (y ) be the
displacement vectors on Γ+ and Γ− the relative displacement
vector can be written as,
ω y = u+ y − u− y (51)
as a vector representing the displacement “jump” between
the two singular lines and
δ y = limω y = lim u+ y − u− y ,
h→0
h→0
(52)
If δ =
/ 0, the singular band is modelling a discontinuous
displacement field as the limit of a continuous one. This
allows a crack to be idealised as a limit (with mesh
refinement) of a band of finite elements where, by means of
some numerical mechanisms, the condition δ =
/ 0 is satisfied.
2.6.2. Displacement and Traction Vectors in the Singular Band.
Consider a singular band in the solid, with a width h
according to Figure 1.
Along a coordinate line x , the displacement vector u can
be expanded from its value in the Γ− line using Taylor’s series
as
u x , y = u− y +
∂u
∂x
−
Δx + O h2 ,
(53)
The equilibrium across the singular band will be enforced by
assuming the traction vector t acting on the plane defined by
the normal n:
t i = σi j n j ,
which is constant in the x direction that is
u+ y = u− y +
∂u
∂x
−
(54)
∞
gf =
u x , y = u− y +
u x , y = u− y +
−
∼
=u
∂u
∂x
−
Δx + O h2 ,
Δx 2 ω y +0 h
h
(55)
0
gf =
(57)
(59)
∞
0
σi j ˙ i j dτ.
(60)
=
0
∞
0
1 ∂u̇i ∂u̇ j
+
dτ
2 ∂x j ∂xi
σi j
∂u̇
σi j i dτ −
∂x j
∞
0
(61)
1 ∂u̇i ∂u̇ j
σi j
+
dτ.
2 ∂x j ∂xi
The last integrand in (61) is zero, being the product of a
symmetric and antisymmetric tensor, so that
gf =
∞
0
∂u̇
σi j i dτ =
∂x j
∞
0
∂ σi j u̇i dτ
∂x j
(62)
where the Cauchy’s equations for quasistatic processes and
negligible body forces (∂σi j /∂xi = 0) have been considered.
The total dissipated energy in the domain Ω∗ is
Ω∗
∞
∗
g f dΩ =
Ω∗
0
!
∂ σi j u̇i dτ dΩ∗ .
∂x j
(63)
By applying the Gauss’s theorem to (63) and using (58)
we obtain,
" ∞
Γ∗
0
#
σi j n j u̇i dτ dΓ∗ =
" ∞
Γ∗
0
#
ti u̇i dτ dΓ∗ (64)
Owing to the infitesimal width of the band, the curvilinear integral in (64) can be evaluated only on the lines Γ∗+ and
Γ∗− (see Figure 1):
W∗ =
" ∞
Γ∗+ ∪Γ∗−
0
#
ti u̇i dτ d y ,
(65)
and taking into account (56) and (59) we obtain,
W∗ =
" ∞
Γ∗+ ∪Γ∗−
φ = 0 =⇒ Γ− ,
∞
y + φ x , y ω y ,
φ = 1 =⇒ Γ+ .
σi j x , y , τ d i j x , y , τ =
(56)
where φ is a function to be determined, which approximates
Δx /h when h → 0.
From (51) and (56) it can be seen that
For uniaxial deformation process g f would be, for a given
point, the area under the stress-strain curve at that point. By
taking the linearized geometric equations g f can be expressed
as,
W∗ =
From the (51), (53), (54) we can write
2.6.3. Energy Dissipation Within the Band. For a generic
deformation process which takes place over a time τ(0 ≤ τ ≤
∞) the specific energy dissipation (energy per unit volume)
within a closed domain Ω∗ (see Figure 1) is given by
W =
h y + O h2 .
t x , y = t+ y = t− y .
∗
and consenquently, for a point in the Γ+ line
(58)
+
Γ∗+ ∪Γ∗−
0
#
ti− y , τ u̇−i y , τ dτ d y " ∞
0
#
ti− y , τ φ x , y ω̇ y , τ dτ d y .
(66)
10
International Journal of Aerospace Engineering
y’
y’
Γ−
B
Γ−
Γ+
Γ+
Γ∗−
h
Ω∗
n
Γ∗+
x’
x’
Ω
A
Γ∗
(a)
(b)
Figure 1: Analysis within the singular band.
The first integral vanishes because the contributions on
Γ∗+ and Γ∗− cancel each other out. Thus, the dissipated
energy on Ω∗ is
W∗ =
"
Γ∗+ ∪Γ∗−
φ x , y
∞
0
−
ti y , τ ω̇ y , τ dτ d y .
(67)
W=
Ω
g f dΩ =
yB " ∞
yA
0
#
ti− y , τ ω̇ y , τ dτ d y .
(68)
Equation (68) establishes that the energy dissipated
within the idealized band can be written as a curvilinear
integral along its length. The integrand of (68) represents the
energy dissipated per unit of area, which in terms of fracture
mechanics, is the fracture energy G f :
G f y =
∞
0
ti− y , τ ω̇ y , τ dτ.
W∗ =
Γ∗+ ∪Γ∗−
G f φ x , y d y = G f
Γ∗
W∗ = Gf
Ω∗
∂φ
dΩ∗ =
∂x
Ω∗
g f dΩ∗ .
Gf
∂φ
= ∗,
∂x
l
∂φ
∂x
−1
.
(73)
σ 0 f
,
2
(74)
where σ 0 is the material strength and the degraded stress is
given by
σ d = σ(1 − d()),
(75)
where the damage evolution law d() is given in terms of the
strain as follows:
d() = 1 −
φ x , y d y . (70)
0 1 + κ2 ()(2κ() − 3) ,
(76)
with
(71)
The local form of (71) is
gf = Gf
gf =
By applying Green’s theorem and taking into account
(60), we can write the energy dissipated within the band as,
The parameter l∗ plays the role of relating the specific
energy (per unit of volume) and the fracture energy (per unit
of area). l∗ is also identified as the characteristic length or
crack band width used in existing cracking models [52]. For
the unidimensional case and using the proposed Hermitian
stress-strain softening law, the specific energy can be written
as
(69)
If G f is assumed to be a material property independent
of the spatial position of the point from (67) and (69) we
obtain,
l∗ x , y =
#
Now, if the case where Ω∗ = Ω is considered; that is,
the whole band between points A and B in Figure 1, the total
energy dissipated within the band between points A and B is
where
(72)
κ() =
− 0
.
f − 0
(77)
Using (74) the failure strain f can be written in terms of
the fracture energy and the characteristic length as it follows:
f =
2G f
.
σ 0 l∗
(78)
International Journal of Aerospace Engineering
11
Φ1 = 1
ζ
(ϕ1 , ζ) = (1, 0)
ξ
η
Node 8
Φ2 = 1
x’, n
Node 5
Node 7
ξ
η
θj
y’
(ϕ4 , ζ) = (1, 0)
Node 6
n
Φ4 = 1
(ϕ2 , ζ) = (0, 0)
Φ3 = 0
Node 3
Figure 3: Finite element band modelling.
(ϕ3 , ζ) = (0, 0)
Node 1
Node 2
Figure 2: Determination of the characteristic length for hexahedron elements.
2.6.4. Determination of the Function φ and the Characteristic
Length l∗ . In order to apply the theory presented in the
previous section to the discretized medium, consider a mesh
of C 0 continuous hexahedron solid finite elements (see
Figure 2).
A set of cracked elements is determined by using
failure criteria for detecting the crack initiation, the crack
orientation depends on the fibre direction. The cracked plane
is defined here as a normal vector which is parallel to the
fibres for fibre failure and normal to the fibre direction for
matrix failure. The algorithm described in this section for
determining the characteristic length was proposed by Oliver
[57] and it has the advantages of calculating the characteristic
length for arbitrary crack directions and any finite element
geometry.
From (71) the function φ has to be continuous and
derivable, satisfying (57). A simple function defined in
the isoparametric coordinates ξ and η which fulfils these
requirements is
φ ξ, η =
nc
$
Ni ξ, η φi ,
Γ(e)
integration points. The following algorithm proposed by
Oliver [57] has been used for determining the characteristic
length at each integration point j, as shown in Figure 4,
(1) A set of local cartesian axes x , y is defined at the
centre of the element, this being identified by the values of
the isoparametric coordinates (ξ = 0, η = 0 and ζ = 0).
The direction of the local axis x is defined by the normal to
the fracture plane, which is the fibre angle for fibre failure
(θ j = θ f ) and (θ j = θ f + 900 ) for matrix failure.
⎡
⎤⎧ ⎫
sin θ j ⎥⎨xi ⎬
⎢ cos θ j
⎦
=
⎩ yi ⎭
⎩ y⎭ ⎣
− sin θ j cos θ j
i
⎧ ⎫
⎨ xi ⎬
(2) Values of φ at each corner node are established
according to their position with respect to the local axis x ,
y (φi = 1 if xi ≥ 0, otherwise φi = 0).
(3) The characteristic length, at the present integration
point j with isoparametric coordinates ξ j and η j and
cracking angle θ j is obtained as
l
∗
ξj, ηj =
−1
∂φ(ξ j , η j )
∂x
⎛ ⎡
nc ∂N ξ , η
$
i j
j
=⎝ ⎣
cos θ j
(79)
(81)
∂x
i=1
+
∂Ni ξ j , η j
∂y
i=1
where nc is the number of corner nodes of a virtual plane
located at the midplane of the element (nc = 4 for our case),
Ni are the standard C 0 shape functions of an element of nc
virtual nodes in its midplane and φi is the value of the φ at
corner i. If the crack location inside the element is known, φi
takes the value +1 if the corner node i is ahead the crack, and
0 otherwise. The function defined by (79) fulfils the required
condition of continuity within elements and takes the values
+1 for the nodes on the boundaries ahead the crack and 0 for
the nodes on the boundaries behind the crack (see Figure 3).
In general, however, the exact crack location is not
known, and usually only some indication of the onset
of cracking and the crack directions is obtained at the
(80)
⎤ ⎞−1
sin θ j ⎦φi ⎠
where
⎧
⎫
∂N ⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎬
∂x
∂N ⎪
⎪
⎪
⎪
⎪
⎪
⎩
⎭
∂y
⎧
⎫
∂N ⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎬
∂ξ
−1
= Jxy
,
∂N ⎪
⎪
⎪
⎪
⎪
⎪
⎩
⎭
(82)
∂η
and Jxy is defined as the Jacobian matrix given by
⎡
Jxy
∂x ∂y
⎤
⎥
⎢
⎢ ∂ξ ∂ξ ⎥
= ⎢ ∂x ∂y ⎥,
⎦
⎣
∂η ∂η
(83)
12
International Journal of Aerospace Engineering
(φi = 1 if xi ≥ 0, otherwise φi = 0) in a similar way as the
one described previously.
The characteristic length associated with transverse compression at the present integration point j with isoparametric
coordinates ξ j and ζ j and fracture angle θ f is given by
x’, n
η
Φ1 = 1
ξ
θj
y’
Φ4 = 1
Φ2 = 0
l
∗
ξj, ζj =
∂φ(ξ j , ζ j )
∂x
−1
⎡
nc ∂N ξ , ζ
$
i j j
=⎝ ⎣
cos θ f
⎛
∂x
i=1
Φ3 = 0
(87)
) −1
∂Ni (ξ j , ζ j )
+
sin(θ f ) φi
∂z
Figure 4: Computation of φ values at the virtual midplane of the
element.
where
where the partial derivatives with respect to the isoparametric coordinates are written as
⎧
⎫
∂N ⎪
⎪
⎪
⎪
⎨
⎬
∂x
1
1
1
1
=
1 + η x1 − 1 + η x2 − 1 − η x3 + 1 − η x4 ,
∂ξ
4
4
4
4
⎪
⎪
⎪
⎩ ∂N ⎪
⎭
∂x 1
1
1
1
= (1 + ξ)x1 + (1 − ξ)x2 − (1 − ξ)x3 − (1 + ξ)x4 ,
∂η 4
4
4
4
∂x
∂z
where the pairs (xi , yi ) refers to the global coordinates of the
virtual nodes defining the midplane of the element.
For transverse compression failure a set of cracked
elements is determined by using the stress based criterion
defined by (28) and the cracked plane is defined by the
fracture angle θ f . The function φ is given b
φ(ξ, ζ) =
nc
$
Ni (ξ, ζ)φi ,
(85)
i=1
where nc is the number of corner nodes of a virtual plane
located at the midplane of the element according to Figure 5,
Ni are the linear shape functions defined previously and nc
virtual nodes defining the virtual cracking midplane and φi
is the value of the φ at corner i.
x , z is an auxiliar coordinate system defined at the
centre of the element, this being identified by the values of
the isoparametric coordinates (ξ = 0, η = 0 and ζ = 0) with
the direction of the local axis x is defined by the normal to
the fracture plane,
⎤⎧ ⎫
sin θ j ⎥⎨xi ⎬
⎢ cos θ j
⎦
.
=
⎩ zi ⎭
⎩z ⎭ ⎣
− sin θ j cos θ j
i
⎧ ⎫
⎨xi ⎬
⎡
(86)
The values of φ at each corner node are established
according to their position with respect to the local axis x , z
∂ξ
−1
= Jxz
,
∂N ⎪
⎪
⎪
⎪
⎪
⎪
⎩
⎭
(88)
∂ζ
and Jxz is defined as the Jacobian matrix given by
⎡
∂y 1 1
1
1
=
1 + η y1 − 1 + η y2 − 1 − η y3 + 1 − η y4 ,
∂ξ
4
4
4
4
∂y 1
1
1
1
= (1 + ξ)y1 + (1 − ξ)y2 − (1 − ξ)y3 − (1 + ξ)y4 ,
∂η 4
4
4
4
(84)
⎧
⎫
⎪
∂N ⎪
⎪
⎪
⎪
⎪
⎨
⎬
Jxz
∂x ∂z ⎤
⎢ ∂ξ ∂ξ ⎥
⎥
=⎢
⎣ ∂x ∂z ⎦,
(89)
∂ζ ∂ζ
where the partial derivatives with respect to the isoparametric coordinates are written as
∂x
1
1
1
1
= (1 + ζ)x1 − (1 + ζ)x2 − (1 − ζ)x3 + (1 − ζ)x4 ,
∂ξ
4
4
4
4
∂x
1
1
1
1
= (1 + ξ)x1 + (1 − ξ)x2 − (1 − ξ)x3 − (1 + ξ)x4 ,
∂ζ
4
4
4
4
∂z 1
1
1
1
= (1 + ζ)z1 − (1 + ζ)z2 − (1 − ζ)z3 + (1 − ζ)z4 ,
∂ξ
4
4
4
4
∂z 1
1
1
1
= (1 + ξ)z1 + (1 − ξ)z2 − (1 − ξ)z3 − (1 + ξ)z4 ,
∂ζ
4
4
4
4
(90)
where the pairs (xi , zi ) refer to the global coordinates of
the virtual nodes defining the midplane of the element.
In-plane shear cracking is strongly dependent on the fibre
orientarion within the element therefore, the characteristic
length associated with in-plane shear failure has been
assumed to be the same as the one defined for fibre failure or
failure in the warp direction. For out-of-plane shear failure
modes, cracks are assumed to be smeared over the thickness
of the element with a crack band defined between upper
and lower faces of the element which is equivalent to assume
θ f = 900 in (87),
l∗ ξ j , ζ j =
∂φ(ξ j , ζ j )
∂x
−1
⎛ *
) ⎞−1
nc
$
∂N
(ξ
,
ζ
)
i
j
j
=⎝
φi ⎠ .
i=1
∂z
(91)
International Journal of Aerospace Engineering
13
(5) Compute failure index for fibre kinking (σ1n+1 < 0):
ζ
η
x’, n
Node 8
c
F1 σ1
(ϕ1 , η) = (1, 0)
ξ
z’
Node 5
(ϕ2 , η) = (0, 0)
(ϕ4 , η) = (1, 0)
Node 6
Node 3
Node 1
n+1 σ1 =
Xc
(95)
≥ 1.
c n+1
Store (1,0
) when F1c (σ1n+1 ) = 1.
(6) Compute IFF index for tensile/shear matrix cracking
(σ2n+1 > 0):
Node 7
θf
n+1
F2t
n+1 n+1
σ2n+1 , τ23
, τ12
=
σ2n+1
Yt
2
τ n+1
+ 23
S23
2
τ n+1
+ 12
S12
2
≥ 1.
(96)
(ϕ3 , η) = (0, 0)
n+1 n+1
Store (tm,0 )n+1 and (σ tm,0 )n+1 when F2t (σ2n+1 , τ23
, τ12 ) = 1.
(7) Compute IFF index for matrix crushing and
intralaminar shear cracking ( σ2n+1 < 0):
Node 2
Figure 5: Computation of the characteristic length for transverse
compression.
n+1 n+1
, σnl =
F2c σnt
n+1
σnt
A
n+1
S23 + μnt σnn
2
+
σnln+1
n+1
S12 + μnl σnn
2
≥ 1.
(97)
3. Numerical Implementation
This section presents details about the numerical implementation of the proposed failure model into ABAQUS
FE code. The code formulation is based on the updated
Lagrangian formulation which is used in conjunction with
the central difference time integration scheme for integrating
the resultant set of nonlinear dynamic equations. The
method assumes a linear interpolation for velocities between
two subsequent time steps and no stiffness matrix inversions
are required during the analysis. The explicit method is
conditionally stable for nonlinear dynamic problems and the
stability for its explicit operator is based on a critical value of
the smallest time increment for a dilatational wave to cross
any element in the mesh. The numerical implementation
steps are listed as follows:
(1) Stresses and strain-rates are transformed from the
element axes to material axes:
σld
n
= T σgd
n
n+1 n+1
Store (cm,0 )n+1 , (σ cm,0 )n+1 and (θ f )n+1 when F2c (σnt
, σnl ) =
1.
(8) Compute failure index for inplane shear cracking:
F12 τ12
=
l∗ ξ j , η j =
∂φ(ξ j , η j )
∂x
(4) Compute failure index for fibre failure
σ n+1
F1t σ1n+1 = 1 ≥ 1.
Xt
> 0):
(94)
(98)
≥ 1.
−1
(99)
∂x
i=1
(92)
(σ1n+1
S12
⎛ ⎡
nc ∂N ξ , η
$
i j
j
=⎝ ⎣
cos θ j
,
where the indices l and g refer to the material (local) and
element (global) coordinate systems, respectively and T is the
transformation matrix defined by (2). The superscripts n and
n + 1 refer to the previous and current time, respectively:
(2) Compute constitutive matrix C from (3) for trial
stresses.
(3) Based on the current time step compute strain
increments and update the elastic stresses and strains
σ n+1 = σ n + CΔ,
(93)
n+1 = n + Δ.
n+1 τ12 i
n+1
)n+1 when F12 (τ12
) = 1.
Store (γ12,0 )n+1 , (γ12,0
(9) Compute characteristic lengths for each failure mode
using objectivity algorithm:
In-plane failure modes
+
˙ ln+1 = T ˙ gn+1 ,
t n+1
) when F1t (σ1n+1 ) = 1.
Store (1,0
n+1
∂Ni ξ j , η j
∂y
⎤ ⎞−1
sin(θ j )⎦φi ⎠
θ j = θfibre for fibre failure in tension and compression, θ j =
θfibre + 900 for matrix cracking in tension/shear, θ j = θfibre
for in-plane shear failure.
Out-of-plane failure modes
⎛
l∗ ξ j , ζ j = ⎝
∂φ ξ j , ζ j
∂x
⎞−1
⎠
⎤ ⎞−1
⎛ ⎡
nc ∂N ξ , ζ
$
i j j
⎦φ i ⎠ ,
=⎝ ⎣
i=1
∂z
(100)
θ j = θ f , where θ f is the fracture plane orientation when
F2c (σnt , σnl ) = 1 for matrix crushing/intralaminar shear
failure modes.
(10) Compute final strains based on intralaminar fracture toughness and characteristic length associated with each
failure mode.
14
International Journal of Aerospace Engineering
Fibre failure
U1 = 0
2Gtfibre
,
Xt l ∗
t
1,
f =
2Gcfibre
c
1,
f =
Xc l ∗
U1 = 0
U1 = Up
(101)
,
U1 = 0
IFF (Matrix cracking in tension/shear)
⎡
tm, f
2
cos2 (θ)
= t ⎣
σ m,0
gmt c
λ
sin2 (θ)
+
gms c
λ ⎤−1/λ
⎦
,
(102)
U1 = Up
U1 = 0
U1 = Up
IFF (Matrix cracking in compression/shear)
2Gcmatrix
2GII
= c c∗ ,
σ cm,0 l∗
σ m,0 l
cm, f =
(103)
U1 = Up
in-plane shear failure
Undeformed shape
Deformed shape
2Gshear
.
S12 l∗
γ12, f =
(104)
Figure 6: Boundary conditions for tensile loading.
(11) Update damage parameters d1 , dtm , dcm , d12 when the
failure criteria are met with,
⎧
n+1
⎪
⎪
⎪(d1 )
⎪
⎨
if (d1 )n+1 > (d1 )n ,
F1t
d1 = ⎪
⎪
⎪
⎪
⎩(d )n
1
t
dm
=
t
dm
n
⎧ n+1
c
⎪
dm
⎪
⎪
⎪
⎨
c
dm
=
⎪
⎪
⎪
⎪
⎩
c
dm
σ1n+1
> 1 or
F1c
U1 = 0
σ1n+1
U1 = Up
> 1,
otherwise,
⎧ n+1
t
⎪
dm
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎩
U1 = 0
n
⎧
⎪
(d12 )n+1
⎪
⎪
⎪
⎨
t
if dm
F2t
n
U1 = Up
t
> dm
,
n+1 n+1
σ2n+1 , τ23
, τ12
U1 = 0
> 1,
otherwise,
c
if dm
n+1
(105)
n
U1 = 0
U1 = Up
c
> dm
,
n+1 n+1
F2c σnt
, σnl > 1,
otherwise,
U1 = Up
if (d12 )n+1 > (d12 )n ,
Undeformed shape
Deformed shape
n+1
F12 τ12
> 1,
d12 = ⎪
⎪
⎪
⎪
⎩(d )n
12
n+1
Figure 7: Boundary conditions for compression loading.
otherwise.
(12) Stress update
⎧ ⎫n+1
⎪
σd ⎪
⎪
⎪
⎪
⎪ 1⎪
⎪
⎪
⎪
⎪
d⎪
⎪
⎪
⎪
⎪
σ
⎪
⎪
2
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
d
⎨ σ3 ⎪
⎬
⎪
d ⎪
⎪
⎪
τ12
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
d ⎪
⎪
⎪
τ23 ⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩τ d ⎪
⎭
13
=
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
(1 − d1 )σ1
dc
t σ m
1 − dm
2
dc
σ3 m
⎫
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎬
c
⎪
t (1 − d )τ dm ⎪
⎪
1 − dm
⎪
⎪
12 12 ⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
dmc
⎪
⎪
⎪
⎪
t
⎪
⎪
1
−
d
τ
⎪
⎪
m 23
⎪
⎪
⎪
⎪
⎪
⎪
c
⎪
⎪
d
⎩
⎭
m
τ13
(13) Stress transformation from the material(local) to the
element coordinate system (global)
.
σgd
n+1
= (T)
−1
σld
n+1
.
(107)
(106)
(14) Compute element nodal forces, add element hourglass control nodal forces, transform element nodal forces
to global coordinates, solve for accelerations via equation of
motion and update velocities and displacements.
International Journal of Aerospace Engineering
15
1
U1 = U2 = U3 = 0
U1 = Up
U2 = 0
U3 = 0
U1 = U2 = U3 = 0
0.75
U1 = Up
U2 = 0
U3 = 0
Damage
U1 = U2 = U3 = 0
U1 = Up
U2 = 0
U3 = 0
0.5
U1 = U2 = U3 = 0
U1 = Up
U2 = 0
U3 = 0
0.25
Undeformed shape
Deformed shape
0
Figure 8: Boundary conditions for in-plane shear loading.
–0.08
–0.04
0
0.04
0.08
Strain (m/m)
Damage evolution (compression only)
Damage evolution (tension only)
Combined damage evolution (tension + compression)
1
Figure 10: Damage evolution law due to combined tension and
compression loading in the fibre direction.
0.75
Normalised stress
0.5
0.25
0
–0.25
–0.5
–0.08
–0.04
0
0.04
0.08
Strain (m/m)
Loading-unloading-reloading (tension + compression)
Loading in compression only
Loading in tension only
Figure 9: Element loaded-unloaded-reloaded in tension and
compression in the fibre direction.
4. Single Element Validation
Numerical simulations at element level were carried out
to validate the proposed damage model. The numerical
tests consisted of exciting each failure mode individually
and verify if damage irreversibility conditions and energy
concepts are always satisfied for each failure mode. The
material properties used in model were taken from [58] and
they are summarised in Tables 1, 2 and 3. The material system
consists of a quasi unidirectional carbon UD tape supplied by
EUROCARBON, which has a plain weave pattern with T700
carbon fibres in the warp direction and a small fraction of
PPG glass fibres in weft direction to hold the carbon fibres
together embedded into an infusible PRIME 20 LV epoxy
resin system.
The element was loaded under displacement control
in each direction with prescribed displacements U1 =
U p and boundary conditions shown in Figures 6, 7 and
8. The stress in the fibre direction was normalised with
respect to the fibre tensile strength whilst the stress in the
matrix (transverse) direction was normalised with respect
to the matrix compression strength. The shear stress was
normalised with respect to the shear strength. The structural
responses for one element loaded in tension, compression and combined tension-compression under loadingunloading-reloading conditions in the fibre and matrix
directions are shown in Figures 9 and 11, respectively. As
damage commences the stresses are gradually reduced to
zero. The material stiffness in each direction are also reduced
as damage cummulates and during unloading the damage
irreversibility condition is fully satisfied in order to avoid
material self-healing as shown in Figures 10, 12 and 14. The
nonlinear behaviour in shear including damage combined
with plasticity effects is shown in Figure 13.
5. Coupon Tests Validation
In this section the predictions obtained using the proposed
failure model are compared against experimental results
16
International Journal of Aerospace Engineering
Table 1: Mechanical properties of each ply.
E22 = E33 (GPa)
8.11
E11 (GPa)
100
G12 (GPa)
4.65
G13 (GPa)
4.65
G23 (GPa)
5.0
ν12 = ν13
0.3
ν23
0.4
1
0.6
0.4
0.8
0
0.6
Damage
Normalised stress
0.2
–0.2
0.4
–0.4
–0.6
0.2
–0.8
–1
–0.04
–0.02
0
0.02
Strain (m/m)
0.04
0
0.06
–0.04
Loading in compression only
Loading-unloading-reloading (tension + compression)
Loading in tension only
Figure 11: Element loaded-unloaded-reloaded in tension and
compression in the matrix direction.
–0.02
0
0.02
Strain (m/m)
0.04
0.06
Damage evolution (compression only)
Damage evolution (tension only)
Damage evolution due to compression (tension + compression)
Damage evolution due to tension (tension + compression)
Figure 12: Damage evolution law due to combined tension and
compression loading in the matrix direction.
Table 2: Ply strengths.
Xt (MPa)
2000
Xc (MPa)
1000
Yt (MPa)
100
Yc (MPa)
160
S12 (MPa)
140
1
Table 3: Intralaminar fracture toughnesses.
Gtm (kJ/m2 )
2.0
Gcm (kJ/m2 )
2.0
Gs (kJ/m2 )
2.0
at coupon level in tension compression and shear. The
experimental results were taken from [58] in which the first
author tested composite laminates manufactured using the
Resin Infusion under Flexible Tooling (RIFT) process. The
mechanical properties of the material are listed in Tables 1, 2
and 3. The dimensions of the coupons used by the author for
the tensile tests in the fibre and matrix directions were 200 ×
20 × 1.35 mm3 and 200 × 20 × 2.25 mm3 , respectively, with a
gauge length of 100 mm for both cases. For the in-plane shear
tests the specimen had dimensions of 200 × 20 × 2.7 mm3
with a gauge length of 100 mm and lay-up of [+450 / − 450 ]3 .
The material constants c1 and c2 required by the model to
predict the inplane shear behaviour were obtained by finding
the best-fit between the model and the experimental shear
Normalised stress
Gtf (kJ/m2 ) Gcf (kJ/m2 )
100.0
25.0
0.8
0.6
0.4
0.2
0
0
0.02
0.04
0.06
Absolute strain
0.08
Normalised stress
Figure 13: Element loaded-unloaded-reloaded in shear.
International Journal of Aerospace Engineering
17
1
0.8
Damage
0.6
Figure 16: Finite element mesh used in the tensile test simulation.
0.4
0.2
0
0.02
0.04
0.06
0.08
Absolute strain
Damage evolution
Figure 14: Damage evolution law in shear.
Figure 17: Finite element mesh used in the compression test
simulation.
×10 6
150
Shear stress (Pa)
120
90
C1 = 3.71 GPa
C2 = 24
60
30
0
0
0.02
0.04
0.06
0.08
Strain (m/m)
0.1
0.12
Experimental model
Figure 15: Comparison between numerical and experimental
curves for the in-plane shear behaviour.
stress and shear strain curve. A comparison between the
experimental and numerical shear stress-shear strain curves
together with the material model parameters c1 and c2 are
presented in Figure 15. As the experiments were carried out
quasistatically, a high value in the order of 109 was assigned
to c3 in order to neglect the strain rate effects in the numerical
response in shear.
The lay-ups used for the specimens loaded in tension
in both fibre and matrix directions were [00 ]3 and [00 ]5 ,
respectively. The mechanical properties of each layer are
listed in Tables 1, 2, 3. Each ply with a nominal thickness
of 0.45 mm was individually modelled using a single solid
hexahedral linear element through the thickness direction.
The tabs behaviour was assumed to be orthotropic linear
elastic. The virtual specimens were assumed to be fixed in one
end and loaded quasistatically under displacement control
using the dynamic relaxation method and the meshes used
in the tensile and compression test simulations are shown in
Figures 16 and 17.
Figures 18–22 show the correlation between numerical
predictions and experimental results for the tension and
compression in both fibre and matrix directions and in-plane
shear tests, respectively.
In general, the numerical predictions agree very well
with the experimental results for all simulations carried
out at coupon level. It can be clearly seen from Figure 18
the stiffening effect in tension in the longitudinal direction.
This effect has not been taken into account in the model
development and it is due to the fabric arquitecture which
leads to the development of high interlaminar shear stresses
between fibres located in the fill and warp direction during
the loading process which results in a reduction of the initial
crimp angle of the fabric therefore increasing the laminate
stiffness along the loading history.
Comparisons between the experimental and predicted
failure locations are shown in Figures 23–26. The numerical
predictions indicate that the use of sharp corners end tabs
leads to failure near to the tab due to the high stress
concentrations in these regions. This effect was also observed
18
International Journal of Aerospace Engineering
×104
5
4.5
4
Load (N)
3.5
3
2.5
2
1.5
1
0.5
0
0
0.5
1
1.5
2
×10−3
Displacement (m)
Numerical
Experimental (specimen-1)
Experimental (specimen-2)
Experimental (specimen-3)
Experimental (specimen-4)
Experimental (specimen-5)
Figure 18: Comparison between numerical predictions and experimental results for the fibre tensile tests.
3000
2500
Load (N)
2000
1500
1000
500
0
0
1
2
3
4
5
Displacement (m)
6
7
8
×10−4
Numerical
Experimental (specimen-1)
Experimental (specimen-2)
Experimental (specimen-3)
Experimental (specimen-4)
Experimental (specimen-5)
Figure 19: Comparison between numerical predictions and experimental results for the matrix tensile tests.
in the experiments. The prediction of the failure location
for transverse compression, Figure 25, correlate very well
with the experiments indicating that the failure in transverse
direction is a shear dominated failure mode with a fracture
plane orientated about 500 in respect to the through the
thickness direction. In order to avoid element distortion
International Journal of Aerospace Engineering
19
12000
10000
Load (N)
8000
6000
4000
2000
0
0
0.2
0.4
0.6
0.8
1
1.2
×10−4
Displacement (m)
Numerical
Experimental (specimen-1)
Experimental (specimen-2)
Experimental (specimen-3)
Experimental (specimen-4)
Figure 20: Comparison between numerical predictions and experimental results for the fibre compression tests.
3500
3000
Load (N)
2500
2000
1500
1000
500
0
0
0.5
1
1.5
2
Displacement (m)
2.5
3
×10−4
Numerical
Experimental (specimen-1)
Experimental (specimen-2)
Experimental (specimen-3)
Experimental (specimen-4)
Experimental (specimen-5)
Figure 21: Comparison between numerical predictions and experimental results for the matrix compression tests.
problems, the failed elements in transverse compression have
been deleted from the mesh using an erosion algorithm
available in the ABAQUS Finite Element Code.
been presented and discussed in this paper. The relevant
contributions of the present formulation are the following:
6. Conclusions
(i) incorporation of strain rate effects on in-plane shear
behaviour by means of a semiempirical viscoplastic
constitutive law,
A detailed formulation for a three dimensional failure
model for composite laminates under multiaxial loading has
(ii) implementation of thermodynamically consistent
degradation laws coupled with an advanced
20
International Journal of Aerospace Engineering
6000
5000
Load (N)
4000
3000
2000
1000
0
0
0.2
0.4
0.6
0.8
1
1.2 1.4
Displacement (m)
1.6
1.8
2
×10−3
Numerical
Experimental (specimen-1)
Experimental (specimen-2)
Experimental (specimen-3)
Experimental (specimen-4)
Experimental (specimen-5)
Figure 22: Comparison between numerical predictions and experimental results for the in-plane shear tests.
(a) Numerical prediction
(b) Experimental
Figure 23: Comparison between predicted and experimental failure
locations for the tensile tests.
(a) Numerical prediction
objectivity algorithm which avoids strain localization
and mesh dependence problems regardless of mesh
refinement and element topology,
(iii) high order damage evolution laws have been incorporated into the model in order to avoid generation
of artifitial stress waves during damage growth,
(iv) the model formulation combines a quadratic stress
based criterion with the matrix compression criterion
based on the action plane to predict interfibre-failure
without knowing a priori the orientation of the
fracture plane
(v) usage of the Second Piola Kirchhoff stresses which
potentially avoids material deformation problems,
such as trellising,
(vi) implementation of the new failure model into
ABAQUS/Explicit FE code within brick elements as
an user-defined material model.
(b) Experimental
Figure 24: Comparison between predicted and experimental failure
locations for the compression test in the fibre direction.
Single element simulations were carried out to validade
the model at element level. Comparisons between numerical
predictions and experimental results at coupon level under
different loading conditions were carried out to validade the
proposed failure model. A very good correlation between
numerical predictions and experimental results was obtained
using the proposed failure model.
International Journal of Aerospace Engineering
21
[8]
(a) Numerical prediction
[9]
[10]
[11]
(b) Experimental
Figure 25: Comparison between predicted and experimental failure
locations for the compression test in the transverse direction.
[12]
[13]
(a) Numerical prediction
[14]
[15]
(b) Experimental
Figure 26: Comparison between predicted and experimental failure
locations for the in-plane shear test.
[16]
Acknowledgments
[17]
The authors acknowledge the financial support received for
this work from Fundação de Amparo Pesquisa do Estado de
São Paulo (Fapesp), contract no. 2006/06808-6, 2007/027104 and CNPq Grant 305601/2007-5.
References
[1] R. M. Jones, Mechanics of Composite Materials, Taylor and
Francis, Philadelphia, Pa, USA, 2nd edition, 1999.
[2] S. W. Tsai and E. Wu, “A general theory of strength for
anisotropic materials,” Journal of Composite Materials, vol. 5,
pp. 58–72, 1971.
[3] Z. Hashin, “Failure criteria for uni-directional fibre composites,” Journal of Applied Mechanics, vol. 47, no. 1, pp. 329–334,
1980.
[4] J. Engblom and J. Havelka, “Transient reponse predictions for
transversely loaded laminated composite plates,” AIAA paper
no.89-1302-CP , 1989.
[5] J. Lee, “Three dimensional finite element analysis of layered
fiber-reinforced composite materials,” Computers and Structures, vol. 12, no. 3, pp. 319–333, 1980.
[6] K. Shivakumar, W. Elber, and W. Illg, “Prediction of lowvelocity impact damage in thin circular laminates,” AIAA
Journal, vol. 23, no. 3, pp. 442–449, 1985.
[7] F.-K. Chang and K.-Y. Chang, “Post-failure analysis of bolted
composite joints in tension or shear-out mode failure mode,”
[18]
[19]
[20]
[21]
[22]
[23]
[24]
Journal of Composite Materials, vol. 21, no. 9, pp. 809–833,
1987.
F.-K. Chang and K.-Y. Chang, “A progressive damage model
for laminated composites containing stress concentrations,”
Journal of Composite Materials, vol. 21, no. 9, pp. 834–855,
1987.
H. Hahn and S. W. Tsai, “Nonlinear elastic behavior of
unidirectional composite laminate,” Journal of Composite
Materials, vol. 7, pp. 102–110, 1973.
H. Hahn and S. W. Tsai, “Nonlinear behavior of laminated
composite,” Journal of Composite Materials, vol. 7, pp. 257–
271, 1973.
H. Choi, H. Wu, and F. Chang, “New approach toward understanding damage mechanisms and mechanics of laminated
composites due to low-velocity impact—part II: analysis,”
Journal of Composite Materials, vol. 25, no. 8, pp. 1012–1038,
1991.
G. Davies and X. Zhang, “Impact damage prediction in
carbon composite structures,” International Journal of Impact
Engineering, vol. 16, no. 1, pp. 149–170, 1995.
S. Watson, The modelling of impact damage in Kevlar-reinforced
epoxy composite structures, Ph.D. thesis, Department of Aeronautics, Imperial College of Science Technology and Medicine,
London, UK , 1993.
G. Davies, D. Hitchings, and G. Zhou, “Impact damage and
residual strengths of woven fabric glass/polyester laminates,”
Composites Part A, vol. 27, no. 12, pp. 1147–1156, 1996.
J. Wiggenraad, X. Zhang, and G. Davies, “Impact damage
prediction and failure analysis of heavily loaded, bladestiffened composite wing panels,” Composite Structures, vol.
45, no. 2, pp. 81–103, 1999.
X. Zhang, G. Davies, and D. Hitchings, “Impact damage with
compressive preload and post-impact compression of carbon
composite plates,” International Journal of Impact Engineering,
vol. 22, no. 5, pp. 485–509, 1999.
J. Hou, N. Petrinic, C. Ruiz, and S. Hallett, “Prediction of
impact damage in composite plates,” Composites Science and
Technology, vol. 60, no. 2, pp. 273–281, 2000.
J. Hou, N. Petrinic, and C. Ruiz, “A delamination criterion for
laminated composites under low-velocity impact,” Composites
Science and Technology, vol. 61, no. 14, pp. 2069–2074, 2001.
R. Kim and S. R. Soni, “Experimental and analytical studies on
the onset of delamination in laminated composites,” Journal of
Composite Materials, vol. 18, no. 4, pp. 70–80, 1984.
H. Jen, Y. Kau, and J. Hsu, “Initiation and propagation of
delamination in a centrally notched composite laminate,”
Journal of Composite Materials, vol. 27, no. 3, pp. 272–285,
1993.
J. Brewer and P. Lagace, “Quadratic stress criterion for
initiation of delamination,” Journal of Composite Materials,
vol. 22, no. 12, pp. 1141–1155, 1988.
S. Liu, Z. Kutlu, and F.-K. Chang, “Matrix cracking-induced
delamination propagation in graphite/epoxy laminated composites due to a transverse concentrated load,” in Comp Materials: Fatigue and Fracture, Fourth Volume, W. W. Stinchcomb
and N. E. Ashbaugh, Eds., vol. 1156 of ASTM STP, pp. 86–
101, American Society for Testing and Materials, Philadelphia,
USA, 1993.
L. Iannucci, “Progressive failure modelling of woven carbon
composite under impact,” International Journal of Impact
Engineering, vol. 32, no. 6, pp. 1013–1043, 2006.
C. Soutis and P. Curtis, “A method for predicting the fracture
toughness of CFRP laminates failing by fibre microbuckling,”
Composites Part A, vol. 31, no. 7, pp. 733–740, 2000.
22
[25] C. Soutis, F. Smith, and F. Matthews, “Predicting the compressive engineering performance of carbon fibrereinforced
plastics,” Composites Part A, vol. 31, pp. 531–536, 2000.
[26] Y. Zhuk, I. Guz, and C. Soutis, “Compressive behaviour of
thin-skin stiffened composite panels with a stress raiser,”
Composites Part B, vol. 32, no. 8, pp. 696–709, 2001.
[27] V. J. Hawyes, P. Curtis, and C. Soutis, “Effect of impact
damage on the compressive response of composite laminates,”
Composites Part B, vol. 32, no. 9, pp. 1263–1270, 2001.
[28] Y. Mi, M. Crisfield, and G. Davies, “Progressive delamination
using interface elements,” Journal of Composite Materials, vol.
32, no. 14, pp. 1247–1271, 1998.
[29] J. Chen, M. Crisfield, A. Kinloch, E. Busso, F. Matthews, and
Y. Qiu, “Predicting progressive delamination of composite
material specimens via interface elements,” Mechanics of
Composite Materials and Structures, vol. 6, no. 4, pp. 301–317,
1999.
[30] M. R. Wisnom and F.-K. Chang, “Modelling of splitting
and delamination in notched cross-ply laminates,” Composites
Science and Technology, vol. 60, no. 15, pp. 2849–2856, 2000.
[31] P. P. Camanho and C. G. Davila, “Mixed-mode decohesion
finite elements for the simulation of delamination in composite materials,” NASA/TM-2002-211737, 2002.
[32] R. Vaziri, M. Olson, and D. Anderson, “Damage in composites:
a plasticity approach,” Computer and Structures, vol. 44, pp.
103–116, 1992.
[33] L. Kachanov, “Time of rupture process under creep conditions,” Izvestiya Akademii Navk SSSR. Odtelemie Tekhniheskikh
Nauk, vol. 8, pp. 26–31, 1958.
[34] Y. Rabotnov, “Creep rupture,” in Proceedings of the12th
International Congress of Applied Mechanics, pp. 342–349,
1968.
[35] P. Ladeveze and E. L. Dantec, “Damage modelling of the
elementary ply for laminated composites,” Composites Science
and Technology, vol. 43, no. 3, pp. 257–267, 1992.
[36] A. F. Johnson, “Modelling fabric reiforced composites under
impact loads,” Composites, vol. 32, pp. 1197–1206, 2001.
[37] A. Johnson, A. Pickett, and P. Rozycki, “Computational methods for predicting impact damage in composite structures,”
Composites Science and Technology, vol. 61, no. 15, pp. 2183–
2192, 2001.
[38] K. V. Williams and R. Vaziri, “Application of a damage
mechanics model for predicting the impact reponse of composite materials,” Computers and Structures, vol. 79, pp. 997–
1011, 2001.
[39] A. Matzenmiller, J. Lubliner, and R. Taylor, “A constitutive
model for anisotropic damage in fiber composites,” Mechanics
of Materials, vol. 20, no. 2, pp. 125–152, 1995.
[40] K. V. Williams, R. Vaziri, and A. Poursartip, “A physically based
continuum damage mechanics model for thin laminated
composite structures,” International Journal of Solids and
Structures, vol. 40, no. 9, pp. 2267–2300, 2003.
[41] K. V. Williams, “Simulation of damage progression in laminated composite plates using LS-DYNA,” in Proceedings of
the 5th International LS-DYNA Conference, Southfield, Mich,
USA, September 1998.
[42] S. Pinho, L. Iannucci, and P. Robinson, “Physically-based
failure models and criteria for laminated fibrereinforced composites with emphasis on fibre kinking—part I: development,”
Composites Part A, vol. 37, no. 1, pp. 63–73, 2006.
[43] C. G. Davila and P. P. Camanho, “Physically based failure
criteria for FRP laminates in plane stress,” NASA-TM, 2003.
[44] M. V. Donadon, L. Iannucci, B. G. Falzon, J. M. Hodgkinson,
and S. F. M. de Almeida, “A progressive failure model
International Journal of Aerospace Engineering
[45]
[46]
[47]
[48]
[49]
[50]
[51]
[52]
[53]
[54]
[55]
[56]
[57]
[58]
for composite laminates subjected to low velocity impact
damage,” Computers and Structures, vol. 86, no. 11-12, pp.
1232–1252, 2008.
M. V. Donadon and L. Iannucci, “An objectivity algorithm
for strain softening material models,” in Proceedings of the
9th International LS-DYNA Users Conference, Dearborn, Mich,
USA, June 2006.
G. Voyiadjis and P. Kattan, Advances in Damage Mechanics:
Metals and Metals Matrix Composites, Elsvier, Oxford, UK,
1999.
P. Kattan and G. Voyiadjis, Damage Mechanics with Finite
Elements: Pratical Applications with Composite Tools, Springer,
New York, NY, USA, 2001.
M. V. Donadon, B. G. Falzon, L. Iannucci, and J. M. Hodgkinson, “Intralaminar toughness characterisation of unbalanced
hybrid plain weave laminates,” Composites Part A, vol. 38, no.
6, pp. 1597–1611, 2007.
P. Maimi, P. Camanho, J. Mayugo, and C. Davila, “A continuum damage model for composite laminates: pat II—
computational implementation and validation,” Mechanics of
Materials, vol. 39, pp. 909–919, 2007.
P. Soden, M. J. Hinton, and A. S. Kaddour, “Biaxial test
results for strength and deformation of a range of E-glass and
carbon fibre reinforced composite laminates: failure exercise
benchmark data,” Composites Science and Technology, vol. 62,
no. 12-13, pp. 1489–1514, 2002.
S. Pinho, P. Robinson, and L. Iannucci, “Physically-based
failure models and criteria for laminated fibrereinforced
composites with emphasis on fibre kinking: part II: FE
implementation,” Composites Part A, vol. 37, no. 5, pp. 766–
777, 2006.
Z. P. Bazant and B. H. Oh, “Crack band theory for fracture of
concrete,” Materiaux et Constructions, vol. 16, no. 93, pp. 155–
177, 1983.
A. Puck and H. Schürmann, “Failure analysis of FRP laminates
by means of physically based phenomenological models,”
Composites Science and Technology, vol. 58, no. 7, pp. 1045–
1067, 1998.
A. Puck and H. Schürmann, “Failure analysis of FRP laminates
by means of physically based phenomenological models,”
Composites Science and Technology, vol. 62, no. 12-13, pp.
1633–1662, 2002.
M. Vural and G. Ravichandran, “Transverse failure in thick S2Glass/epoxy fiber-reinforced composites,” Journal of Composite Materials, vol. 38, no. 7, pp. 609–623, 2004.
M. V. Donadon and L. Iannucci, “Bird strike modelling using
a new woven glass failure model,” in Proceedings of the 9th
International LS-DYNA Users Conference, Dearborn, Mich,
USA, June 2006.
J. Oliver, “Consistent characteristic length for smeared cracking models,” International Journal for Numerical Methods in
Engineering, vol. 28, no. 2, pp. 461–474, 1989.
M. V. Donadon, The structural behaviour of composite structures manufactured using Resin Infusion under Flexible Tooling,
Ph.D. thesis, Department of Aeronautics, Imperial College,
London, UK, 2005.