Mariani-2009.pdf

FAILURE OF LAYERED COMPOSITES SUBJECT TO
IMPACTS: CONSTITUTIVE MODELING AND
PARAMETER IDENTIFICATION ISSUES
Stefano Mariani
Politecnico di Milano, Dipartimento di Ingegneria Strutturale
Piazza L. da Vinci 32, 20133 - Milano (Italy)
E-mail: [email protected]
Layered composites subject to impacts can fail by delamination, i.e. by debonding between laminae, if the stress waves cause damaging phenomena to take place mainly within
the resin-enriched interlaminar phases. To simulate delamination at the structural level, processes dissipating energy are lumped onto fictitious zero-thickness interlaminar surfaces,
and softening interface constitutive laws are adopted to describe the progressive failure of
the interlaminar phases.
Since delamination occurs inside very narrow regions, results of experimental testing
on whole composites need to be accurately and reliably filtered to calibrate the interface
constitutive laws. To this aim, here we propose a sigma-point Kalman filter approach. The
performances of the proposed methodology, in terms of constitutive parameter estimations
and dynamic delamination tracking, are assessed through pseudo-experimental testings on
a two-layer composite, and real testings on multi-layer glass fiber reinforced plastic composites.
Keywords: composites, delamination, impact loading, interface constitutive modeling, parameter identification, sigma-point Kalman filter.
1
2
Stefano Mariani
1 Introduction
Foreign objects striking the outer surface of composite structures may cause permanent
damage, or even sudden failure [1]. On the basis of the velocity of the striker, impacts
can be roughly distinguished into two main classes: low-velocity impacts, characterized
by small or negligible effects of inertial forces on the damage/failure mode; high-velocity
impacts, whereby inertial forces strongly affects the damage/failure event. In case of a highvelocity impact, the eventual failure of layered composites can be due to the propagation
of interlaminar cracks only (delamination), or to a local perforation accompanied also by
intralaminar damage [2].
In this work we focus on impact-induced delamination and provide a review of a classical approach to debonding in layered bodies [3–8]. Within this approach it is assumed
that laminae always behave elastically (i.e. intralaminar damaging phenomena are disregarded, or thought to be negligible), and interlaminar resin-enriched phases are lumped
onto zero-thickness surfaces, along which debonding can occur because of the impinging
shock waves. To model delamination, softening interface constitutive laws are adopted
to link the tractions acting upon each interlaminar surface with the displacement jumps
occurring across it. If these laws are able to phenomenologically describe the micromechanical processes leading to debonding, the above approach furnishes accurate results at
the structural level, see e.g. [4]. However, calibration of the interface laws is still an open
issue: since damaging and cracking phenomena linked to delamination take place inside
very narrow regions, direct testing on the interlaminar phases can not be devised. Instead,
inverse analysis procedures can be adopted to efficiently manage experimental data in order
to estimate uncertain constitutive parameters [6, 9–11]; needles to say, the aforementioned
experimental data have to be informative, i.e. they have to carry information on the current
response of interphases to the impact loading.
As for model calibration purposes, standard filtering procedures have been proven accurate enough in case of static loadings [12, 13]; in case of impacts the inverse problem becomes stiff since composite failure, once incepted, usually occurs almost instantaneously [2, 14]. To deal with this issue, we recently adopted the extended Kalman filter
(EKF) [11, 15]; when compared to alternative approaches (like, e.g., neural networks and
least squares), the EKF has the great advantages of being able to work in real-time and of
being explicitly linked to the physics of the ongoing delamination processes. Moreover, the
EKF exploits the evolution in time of the measured fields in order to continuously enhance
model calibration.
As pointed out by several authors (see, e.g., [16, 17]), the EKF looses stability when
nonlinearities become dominant in the equations governing the inverse problem. This is due
to the fact that the EKF replaces nonlinear functions with their relevant tangent surfaces,
leading to possibly biased or even divergent parameter estimates. An alternative approach
to treat nonlinearities in a stochastic framework recently led to the proposal of the so-called
sigma-point Kalman filter (SPKF), also termed unscented Kalman filter [16]. Instead of
linearizing the governing equations, thereby introducing approximations, the SPKF samples
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
3
the statistics of the current state of the system and of model parameters to draw a set of
sigma-points. These sigma-points are then let to evolve according to the actual, nonlinear
dynamics of the problem. The filter estimates are eventually computed by averaging the
information conveyed by all the evolved sigma-points through an ad-hoc numerical scheme
[16, 18]. When compared to the EKF, the SPKF can achieve a much higher accuracy in
the model calibration task [19], often preventing the aforementioned bias and divergence
occurrences.
If used to deal with the nonlinearities accompanying delamination in layered composites, the SPKF has to furnish accurate estimates of uncertain interface model parameters
while tracking the state of the laminate. In this Chapter we analyze this framework in
details, trying to point out the main strengths of the SPKF. To this purpose, in Section
2 we present the equations governing the dynamics of a layered, possibly delaminating
body. Constitutive models for laminae and resin-enriched interphases are then discussed
so as to recognize the basic parameters, in need of an accurate estimation. An explicit
time integration scheme for the equations of motion of the composite, and an explanation
on how to handle model calibration via the SPKF are eventually offered. Section 3 deals
with sigma-point Kalman filtering: physical arguments are used to propose a slightly new
scheme for the drawing of the sigma-points. In Section 4 the performances of the SPKF
are assessed: pseudo-experimental, i.e. fictitious impact testings on a two-layer composite
are first considered; hence, real experiments on multi-layer glass fiber reinforced plastic
(GRP) composites are used to testify the robustness of the SPKF in promptly detecting
delamination.
Throughout the whole Chapter, a matrix notation is adopted, with uppercase and lowercase bold symbols respectively denoting matrices and vectors; a superscript T stands for
transpose, while a superposed dot represents time rates.
2 Dynamics of layered composites
2.1 Governing relations
Let Ω be a three-dimensional, layered body; its smooth outer boundary, with unit outward
normal m, be constituted by the two disjoint parts Γu and Γτ , where displacement and
traction fields are respectively assigned. Ω is assumed to be crossed by nΓ non-intersecting
surfaces, or interfaces Γj , j = 1, ..., nΓ , see Figure 1. Each resulting portion Ωj , j =
1, ..., nΓ + 1, of the bulk will be here referred to as layer, or lamina. Interfaces in laminates
are usually flat and parallel to each other, with a common orientation defined by the unit
vector n; opening (mode I) and sliding (mode II and mode III) displacement discontinuities
along each Γj take place along the n direction and in the sj1 − sj2 plane, respectively.
The equilibrium of Ω at time t is governed by the following equations:
C T σ + b̄ = ̺ü
Mσ = τ̄
in Ω\Γ
on Γτ
(1)
(2)
4
Stefano Mariani
m
1
W
n
W
s11
s21
2
G1
x3
G2
W3
x2
x1
Figure 1: geometry of a three-dimensional layered continuum (here nΓ = 2), and notation.
(
on Γ+
j
on Γ−
j
N σ = −τ
Nσ = τ
j = 1, ..., nΓ
(3)
nΓ +1 j
−
Γ
where: Ω = ∪j=1
Ω and Γ = ∪nj=1
Γj ; Γ+
j and Γj are the sides of Γj respectively
belonging to layers Ωj and Ωj+1 , according to the notation of Figure 1. The two sides
can not be distinguished in the initial configuration (at t = 0) but, as soon as delamination
is incepted, they can experience a relative movement. In the equations above, adopting a
standard Voigt notation: σ is the stress vector, which gathers the independent components of
the stress tensor; b̄ and τ̄ are the assigned loads in the bulk Ω\Γ and along Γτ , respectively;
̺ is the bulk mass density; ü is the acceleration field in the bulk Ω\Γ; C is the differential
compatibility operator; M and N are the matrices respectively collecting the components
of the unit vectors m and n. To locally ensure equilibrium along each surface, sides Γ+
j
and Γ−
are
respectively
acted
upon
by
the
traction
vectors
−τ
and
τ
.
j
Since we aim at modeling phenomena occurring within very short time intervals after
the impact, linearized kinematics proves sufficient. Compatibility thus reads:
ε = Cu
in Ω\Γ
(4)
u = ū
on Γu
(5)
where: ε is the strain vector, which gathers the independent components of the strain tensor;
u is the displacement field in the bulk Ω\Γ; ū is the assigned displacement along Γu .
Across each interface, the displacement discontinuity [u] is defined as:
[u] = u − u
j = 1, ..., nΓ
(6)
Γ+
j
Γ−
j
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
5
The body is assumed to be initially at rest, in an undeformed and unstressed state, such
that:
(
u0 = 0
in Ω
(7)
u̇0 = 0
2.2 Constitutive modeling
Energy dissipation in layered composites subject to impacts arises from: the spreading of
damage and the subsequent propagation of cracks within the interlaminar, resin-enriched
phases; the spreading of damage inside laminae. This latter dissipation mechanism can be
disregarded, as done here, if the impact energy is not high enough to cause a penetration
of the impactor inside the laminate. In such a case, laminae can be assumed to behave
elastically, according to:
σ = E Ωε
in Ω\Γ
(8)
where E Ω is the elasticity matrix of the bulk. Each lamina is usually an orthotropic body;
though not explicitly shown in (8), E Ω can change from lamina to lamina because of a
different orientation of the axes of elastic symmetry.
To simulate strength reduction in the interlaminar phases, eventually leading to
debonding, softening interface constitutive laws prove efficient [3–8, 20, 21]. A detailed,
micromechanics-based representation of the damaging processes caused by the constrained
deformation field inside each interphase is not looked for; instead, a phenomenological re−
lationship is adopted to link the tractions τ acting upon the interface sides Γ+
j and Γj to
the displacement jump [u] occurring across Γj . This constitutive law might be conceived
as the macroscopic (homogenized through the interphase thickness) behavior of the interphase. Since the ratio between the thickness of each interlaminar phase and the thickness of
the whole laminate is usually quite small, this approach can furnish accurate results at the
laminate length scale; on the contrary, an accurate representation of the micromechanical
phenomena preceding delamination is in need of a multi-scale approach (see, e.g., [22–24]).
Along each Γj , both opening/closing (along n) and sliding (in the sj1 − sj2 plane) discontinuities take place under general loading conditions. By way of a simplified scheme
adopted in [25–28], the local kinematics of an interface is governed by the effective displacement discontinuity [u], defined as:
p
(9)
[u] = [u]2n + κ2 [u]2s
where: [u]n = [u]T n and [u]s = |[u] − [u]n n| are, respectively, the opening and sliding
displacement discontinuities; κ is a coupling coefficient, which accounts for the interaction
between stretching and shearing deformation modes inside the interphase. Through an
incremental work equivalence the effective traction τ , work-conjugate to [u], turns out to
be (see [25, 27]):
r
τ2
(10)
τ = τn2 + s2
κ
6
Stefano Mariani
(a)
(b)
(c)
Figure 2: effective traction-displacement discontinuity relationships under tensile loading.
(a) piecewise linear law (12); (b) linear-exponential law (13); (c) exponential law (14).
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
7
where τn and τs are defined like [u]n and [u]s . The mechanical behavior of the interface
can now be furnished in terms of an effective τ − [u] relationship.
For quasi-brittle materials, like interphase resins, the tensile strength is typically much
smaller than the compressive one. Therefore, when dynamic delamination occurs without
intralaminar damage, the response of the interface under compressive stress states along n
can be assumed linear elastic, according to:
τ = K[u]
if [u] < 0
(11)
where K is the stiffness of the interface. On the other hand, the response of the interface
under tensile stress states along n is characterized by strength reduction, leading to softening, beyond the attainment of a peak traction. A simple way to model the transition from the
initial elastic regime to the subsequent softening one, is through a piecewise linear (PWL)
law (see Figure 2(a)):


if [u] ≤ [u]e
τ = K[u]
(12)
τ = τM + Q ([u] − [u]e ) if [u]e < [u] ≤ [u]U


τ =0
if [u] > [u]U
where: [u]e is the effective displacement discontinuity corresponding to the
peak traction
τM = K[u]e ; Q is the (negative) slope of the softening branch; [u]U = 1 − K
Q [u]e is
the effective displacement discontinuity at which the interaction between the crack faces
ceases.
The linear softening in (12) can be a too crude approximation of the post-peak regime
for some materials, whose response features an initial steep descent followed by a long,
much less steep tail (see, e.g., [29]). In such a case, the softening regime can be modeled
through an exponential law (hereafter referred to as linear-exponential, L-E law because of
the pre-peak linear elastic phase) according to (see Figure 2(b)):
(
τ = K[u]
if [u] ≤ [u]e
(13)
τ = τM exp (−ς([u] − [u]e )) if [u] > [u]e
where ς is a model parameter that allows to match the slope of the softening branch just
beyond the attainment of the peak traction τM .
Sometimes, a smooth transition from the elastic regime to the softening one turns out to
be more representative of the actual interphase response. The nonlinear binding model,
originally proposed in [30, 31] for metals and bimetallic compounds and later adopted
also in nonlinear fracture mechanics [20, 21, 32], allows to describe such smooth transition through the following exponential (EXP) law (see Figure 2(c)):
[u]
(14)
τ = K[u] exp −
[u]e
8
Stefano Mariani
Besides the effective stiffness K and strength τM , a full characterization of the nonlinear behavior of the interface has to match the fracture energy, or work of separation G. In
terms of effective quantities, G is defined as the amount of energy required to annihilate the
interaction between the opening/sliding crack faces, i.e.:
G=
Z∞
τ d[u]
(15)
0
From a model calibration perspective, parameters Q and ς in (12) and (13) can be tuned
to accurately match the actual G, since they do affect only the softening branch of the
interface law. On the other hand, after having assigned K in (14), only [u]e can be adjusted:
therefore, both τM and G can not be accurately matched. In the exponential law, in fact, the
following constraint holds:
s
KG
(16)
2 =e
τM
e being the Nepero number. To avoid problems related to this fictitious constraint, a modified exponential law is here formulated as follows (see Figure 3):
[u] q
τ = K[u] exp −
(17)
[u]e
where q shows up as an additional constitutive parameter. In law (17): K still represents
the initial elastic stiffness; [u]e is a reference effective displacement discontinuity, while
peak traction τM is attained at [ū] ≡ (1/q)1/q [u]e . The effective peak traction and fracture
energy are affected by q, according to:
1/q
1
1
exp −
τM =K[u]e
q
q
2
K
G = [u]2e γf
q
q
(18)
γf being the gamma function. The dependence of τM and G on the parameter q is depicted
in Figure 4: it can be seen that G is a monotonically decreasing function of q, whereas τM is
lower-bounded by the value corresponding to q = 1. Having tuned K, this law thus allows
the calibration of both τM and G.
All the above laws but the piecewise linear one, assume that the interaction between
the opening interface sides continues up to [u] → ∞, which seems not physical at the
macroscale. To simulate delamination growth a breakdown threshold therefore needs to be
introduced [33, 34]: as soon as the current traction τ reduces to a small fraction (say 5-10
%) of the peak value τM , the interaction is suddenly assumed to vanish.
When unloading from the tensile envelope occurs, i.e. when [u̇] < 0, the above interface
models can be viewed as either reversible, if τ always belongs to the envelope (leading to
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
9
Figure 3: effect of q on the modified exponential traction-displacement discontinuity law
(17).
(a)
(b)
Figure 4: modified exponential law (17). Effects of q on (a) the effective peak traction τM
and on (b) the effective fracture energy G.
interface healing if softening has already started), or irreversible, if τ decreases following a
radial path to the origin of the τ − [u] plane. These two alternative constitutive assumptions
lead to different entries in the interface tangent stiffness matrix E Γ , linking rates of τ and
[u] in the local n − sj1 − sj2 reference frame according to:
τ̇ = E Γ [u̇]
For additional details, readers are referred to [27, 28].
(19)
10
Stefano Mariani
2.3 Finite element formulation
The weak form of the equilibrium equations (1)-(3) reads:
Z
Ω\Γ
T
εv σ dΩ =
Z
T
v (b̄ − ̺ü) dΩ +
Z
T
v τ̄ dΓτ −
Γτ
Ω\Γ
nΓ Z
X
[v]T τ dΓj
j=1 Γ
∀v ∈ U0 (20)
j
where: v is the test function; εv = Cv; U is the trial solution space, collecting displacement fields u continuous in Ω\Γ, possibly discontinuous along each Γj and fulfilling the
boundary condition (5) on Γu ; U0 is the relevant variation space, with zero prescribed
displacements on Γu . In (20), in view of the assumed linearized kinematics, the relation
−
Γj ≡ Γ+
j ≡ Γj for each interface has been exploited.
Allowing for the elastic bulk constitutive law (8), the following variational statement is
arrived at:
find u ∈ U :
Z
T
̺ v ü dΩ +
Z
T
εv E Ω ε dΩ +
=
Z
[v]T τ dΓj
j=1 Γ
Ω\Γ
Ω\Γ
nΓ Z
X
j
T
v b̄ dΩ +
Z
T
v τ̄ dΓτ
(21)
∀v ∈ U0
Γτ
Ω\Γ
Now, let the finite element approximation of the displacement and deformation fields in
Ω\Γ be (see [35] for the notation):
u∼
= Φuh
(22)
ε∼
= CΦuh = B Ω uh
(23)
where matrix Φ gathers the nodal shape functions, and vector uh collects the nodal displacements.
If delamination is allowed to occur only along element boundaries, the discrete displacement jump field can be written:
[u]
Γj
∼
= B Γj uh
j = 1, ..., nΓ
(24)
Owing to the discrete interpolation fields defined above, the semi-discretized equations
of motion of the composite turn out to be:
h
h
M ü + K Ω u +
nΓ
X
j=1
Rj = F
(25)
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
11
where the mass matrix M , the bulk stiffness matrix K Ω , the internal force vectors Rj and
the external load vector F are, respectively:
Z
̺ ΦT Φ dΩ
M=
Ω\Γ
KΩ =
Z
B TΩ E Ω B Ω dΩ
Ω\Γ
j
R =
Z
(26)
B TΓj τ dΓj
Γj
F =
Z
Ω\Γ
T
Φ b̄ dΩ +
Z
ΦT τ̄ dΓτ
Γτ
Smarter finite element formulations, like the extended or generalized ones [36,37], have
been recently formulated to simulate mixed-mode crack growth in homogeneous solids,
see e.g. [27, 38, 39]. These methodologies allow cracks to propagate not only along interelement edges, but also inside the elements; possible constraints imposed by the mesh layout on crack trajectories, evidenced e.g. in [40], can be therefore alleviated. When dealing with delamination in layered continua, where debonding occurs only along the a-priori
known interlaminar surfaces, crack description looks simple and the aforementioned feature
of the extended finite element method looses much of its advantages.
2.4 Time integration
In our previous works [11, 15] it was shown that the time integration scheme can strongly
affect the stability of the filtering procedure and, therefore, the accuracy of model parameter
estimates.
In case of impact loadings, which cause the propagation of shock waves inside the composite, the time marching algorithm has to dump the spurious high frequency oscillations
linked to space discretization. Otherwise, the numerically computed displacement and velocity fields do not prove reliable enough to be compared to the experimental data.
We adopt here the explicit α−method [41, 42] to advance in time the solution of the
equations of motion (25), see also [19]. After having partitioned the time interval of interest
Nt −1
according to [t0 tN ] = ∪i=0
[ti ti+1 ], at the end of the generic time step [ti ti+1 ]
the solution to (25) is obtained according to the following predictor-integrator-corrector
splitting:
• predictor:
1
ũi+1 =ui + ∆t u̇i + ∆t2 ( − β)üi
2
˜ i+1 =u̇i + ∆t(1 − γ)üi
u̇
(27)
(28)
12
Stefano Mariani
where ∆t = ti+1 − ti ;
• explicit integrator:


üi+1 = M −1 F i+1+α − (1 + α) K ũi+1 +
where: F i+1+α = F (ti + (1 + α)∆t);
R T
B Γj τ̃ i+1 dΓj ;
X
Rji
j
j


R̃i+1  + α Kui +
=
R
Γj
B TΓj τ i dΓj
and
X
j

Rji 
(29)
=
j
R̃i+1
Γj
• corrector:
ui+1 =ũi+1 + ∆t2 β üi+1
˜ i+1 + ∆tγ üi+1
u̇i+1 =u̇
(30)
(31)
In the above equations α, β and γ are algorithmic coefficients. To get a non-oscillatory
velocity field, α = −0.3 has been adopted in all the forthcoming simulations; furthermore,
β and γ have been finely tuned around the values allowing second-order accuracy in linear
elasto-dynamics.
To ensure accuracy of the filtering procedure, the time step size ∆t has been always set
so as to fulfill the Courant condition in the bulk of the composite. Moreover, to speed up
the explicit integrator phase (29), the mass matrix M has been diagonalized by means of a
standard row-sum lumping procedure [42].
Account taken of the explicit format of the integrator stage, the space-time discretized
equations of motion of the laminate (state equations) can be formally written:


ui+1 
z i+1 = u̇i+1 = f zi (z i )
(32)


üi+1
where z is the structural state vector, and mapping f z turns out to be nonlinear because of
the softening interface behavior.
3 Constrained sigma-point Kalman filtering
3.1 Parameter identification via joint Kalman filtering
According to a standard methodology [43], the calibration of constitutive laws can be pursued by Kalman filtering if a state vector x is obtained by joining the structural state vector
z (see Eq. 32) with a vector ϑ gathering all the model parameters to be tuned. At time ti
this can be written:
zi
(33)
xi =
ϑi
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
13
While the current structural state z is always at least partially observed, model parameters to be identified can not be directly measured; by joining z and ϑ, state tracking can
consistently improve model calibration.
In case of irreversible constitutive laws, internal state variables must be gathered by x
too, see e.g. [12, 19].
Allowing for model and measurement errors, the state-space model describing the evolution within the time interval [ti ti+1 ] of the joint state vector and its link with observations turns out to be:
(
xi+1 = f i (xi ) + v i
(34)
yi
= Hxi + wi
where: y is the observation vector, which collects the measured components of the state
vector; v is the process noises; w is the measurement noise. v, w are assumed to be
additive, uncorrelated white and Gaussian processes, with zero mean and covariances V
and W [44, 45]. Since z is defined according to Eq. (32), the observation equation in (34)
shows up as a linear relation between y and x. On the contrary, the interface behavior
renders the evolution equation f nonlinear.
By way of the EKF [12, 46], within the time step the nonlinear mapping f is expanded
in Taylor series, up to the first order, around the current estimates of the state vector and
of model parameters. Bounds on the required accuracy of the initialization of x, and on
the statistics of noises v and w to assure filter stability were provided for linear systems
in [47] and, more recently, for nonlinear systems in [48]. Even in the absence of filter
instabilities, the softening response of the interlaminar surfaces does not always guarantee
the achievement of an accurate model calibration, see [15, 19].
To improve the results when nonlinearities become dominant, the SPKF has been recently proposed [16, 49–51]. At the beginning of the time step, the probability distribution
of x is deterministically sampled through a set of sigma-points χ̂j , j = 0, ..., Nχ . These
sigma-points are then allowed to evolve according to the nonlinear mapping f . The statistics of x at the end of the time step are finally obtained through a proper weighted averaging
scheme [18]. This filtering procedure is detailed in Table 1, where E[2] represents the expected value of 2.
The number of sigma-points and their location in the state vector space are accurately
chosen, so as to achieve high accuracy in the estimated probability distribution of x at the
end of each time step; when compared to the EKF, a better performance of the SPKF, also
in terms of model calibration, is therefore expected [16]. The enhanced accuracy of the
SPKF is discussed next; even though these results have been already presented elsewhere,
they are here collected to show how possible constraints on parameter estimates, not dealt
with by the standard SPKF, can be managed.
3.2 Accuracy of a constrained sigma-point transformation
In this Section we focus on the time interval [ti
i + 1 to simplify the notation.
ti+1 ], but we avoid using indexes i and
14
Stefano Mariani
• Initialization at t0 :
x̂0 =E[x0 ]
P 0 =E[(x0 − x̂0 ) (x0 − x̂0 )T ]
• At ti , for i = 0, ..., N
1. Predictor phase:
χ̂i,j
−
χ̂i+1,j
x̂−
i+1
P−
i+1
=x̂i + ∆χi,j
j = 0, ..., Nχ
=f i (χ̂i,j )
=
Nχ
X
ωj χ̂−
i+1,j
j=0
=R−
i+1
+V
where
R−
i+1 =
Nχ
X
j=0
T
−
−
−
−
x̂
χ̂
−
x̂
ωj χ̂−
i+1
i+1
i+1,j
i+1,j
2. Corrector phase:
U
−
x̂i =x̂−
i + Gi y i − H x̂i
U
−
P i =P −
i − Gi HRi
where
−
T
T
GU
HR−
i =Ri H
i H +W
Table 1: sigma-point Kalman filter.
−1
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
15
Let x be a random vector, featuring at the beginning of the time step mean x̂ and
covariance P . We conceive x as the sum of the mean x̂ and a zero-mean disturbance
∆x = x − x̂. If x undergoes a nonlinear transformation, governed by a mapping f
analytic everywhere so that it can be expanded in Taylor series about x̂, at the end of the
time step we get:
∞
X
1 n
x = f (x) = f (x̂ + ∆x) =f (x̂) +
D
n! ∆x
n=1
−
(35)
n in the series expansion is:
where, with a slight abuse in notation, the n−th order term D∆x
n
D∆x
≡
!n
Nx
X
∂f ∆xℓ
∂xℓ x=x̂
(36)
ℓ=1
Nx being the number of components of the state vector x. Since the derivatives of f in (36)
are evaluated at x = x̂, they are not random variables. The expected value of x− therefore
reads:
#
"∞
X 1
Dn
x̂− = E[x− ] =E [f (x̂)] + E
n! ∆x
n=1
(37)
∞
X
1 n =f (x̂) +
E D∆x
n!
n=1
The relevant error covariance matrix is:
h
T i
P − =E x− − x̂− x− − x̂−

!
∞
X
1
n
n
− E D∆x
D∆x
=E 
n!
n=1
=
∞
X
m 1
m
− E D∆x
D∆x
m!
m=1
∞ X
∞
i
X
n h m T i
1 1 h n
m T
− E D∆x
E D∆xD∆x
E D∆x
n!
m!
n=1 m=1
!T 

(38)
Now, let us suppose to sample the probability distribution of x through a set of sigmapoints χ̂j , j = 0, ..., Nχ , chosen around the current mean x̂ according to:
χ̂j = x̂ + ∆χj ,
j = 0, ..., Nχ
(39)
where the terms ∆χj need to be determined.
Similarly to x, within the time step each sigma-point undergoes the transformation:
χ̂−
j = f (χ̂j ) = f (x̂) +
∞
X
1 n
D
n! ∆χj
n=1
(40)
16
Stefano Mariani
where:
Dn
∆χj
≡
!n
Nx
X
∂f ∆χ
∂xℓ x=x̂ j ℓ
(41)
ℓ=1
At the end of the time step, the information in the evolved sigma-points are collected
via a weighted averaging procedure to obtain:
x̂−
SPT
=
Nχ
X
ωj χ̂−
j
j=0




Nχ
Nχ
∞
X
X
X
1 

=
ωj D n
ωj  f (x̂) +
∆χj
n!
n=1
j=0
(42)
j=0
where ωj are the weights of the sigma-point transformation relevant to the mean of x. The
corresponding error covariance matrix is given by:
P−
SPT
=
Nχ
X
j=0
=
Nχ
X
j=0
T
−
−
−
−
x̂
χ̂
−
x̂
ωj⋆ χ̂−
SPT
SPT
j
j


T
Nχ
Nχ
∞ X
∞
X
X
X
1 1  n
 D m

ωr D n
ωs D m
D ∆χj −
ωj⋆
∆χr
∆χs
∆χj −
n!
m!
r=0
s=0
n=1 m=1
(43)
where ωj⋆ are the weights of the sigma-point transformation relevant to the covariance of x.
If x is a Gaussian random vector, its probability distribution is symmetric with respect
n , n = 1, 3, 5, ..., are zero. To be
to the mean x̂; therefore, all the odd central moments D∆x
compliant with this condition, couples of sigma-points are symmetrically placed around x̂,
according to [16]:


=0

∆χ0
√
N
∆χk
= +ψ P 1k
k = 1, ..., 2χ
(44)
√


∆χ Nχ
= −ψ P 1k
+k
√
2
Here: P represents the square root of matrix P , computed e.g. through a Cholesky
factorization; ψ is a scaling parameter; 1k is a unit vector aligned with component k in the
state vector space. The series expansions (37) and (42) agree up to third order if:
P
Nχ

ωj = 1


 j=0

P

χ

1
 N
j=0 ωj D ∆χj = 0
(45)
2 PNχ
2


ω
D
=
E
D
j

∆χ
∆x
j=0
j




PNχ ωj D 3 = 0
∆χj
j=0
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
17
To simplify the matter, let us assume ωj = ω for j = 1, ..., Nχ ; relations involving the first
and third order terms in (45) are then automatically satisfied. Relations involving the zeroth
and second order terms in (45) then furnish:
(
ω 0 + Nχ ω = 1
(46)
2ψ 2 ω = 1
A further condition to set ω0 , ω and ψ can be furnished by matching the diagonal entries of
the fourth order terms (kurtoses) in (37) and (42). This leads to [19]:
ω0 = 1 −
Nχ
,
6
1
ω= ,
6
ψ=
√
3
(47)
Here we propose an alternative condition to determine ω0 , ω and ψ, partially exploiting
the features of the so-called scaled unscented transformation [52]. Let us assume that model
parameters have to satisfy the constraints:
ϑm ≤ ϑ ≤ ϑM
(48)
where ϑm and ϑM respectively gather the minimum and maximum (if any) allowed values
of model parameters. This requirement must be fulfilled by each sigma-point χ̂−
j , j =
0, ..., Nχ . For j = 0 the conditions (48) are automatically satisfied, since x̂ (and, therefore,
ϑ̂) is computed at the end of the previous time step by averaging sigma-points all fulfilling
the constraints. Further, if ϑ = B ϑ x, B ϑ being a Boolean matrix, conditions (48) are
satisfied by all the sigma-points if:
)
(
ϑM
ϑ̂ℓ − ϑm
Nχ
ℓ − ϑ̂ℓ
ℓ
,
k = 1, ...,
; ℓ = 1, ..., Nϑ
(49)
,
ψ ≤ min
ℓ
ℓ
2
ak
ak
√
where ak = B ϑ P 1k , and Nϑ is the√number of model parameters in ϑ. In the forthcoming
examples, we initially assume ψ = 3 (according to what reported in 47) and reduce its
value if necessary, according to relation (49).
⋆
As for the error covariance matrix P −
SPT , by letting ωj = ωj = ω for j = 1, ..., Nχ , we
get:


Nχ
Nχ
Nχ
∞ X
∞
X
X
X
X
1
1
T
⋆
m T

ωj D n
ωk D m
ωj D n
P−
∆χj
∆χk
∆χj D ∆χj + (ω0 + Nχ ω − 2)
SPT =
n!
m!
n=1 m=1
j=0
j=0
k=0
(50)
In case of Gaussian random variables, independently of the value of ω0⋆ , (38) and (50) agree
up to the third order. Weight ω0⋆ can be set by matching part of the fourth order terms
T
(specifically those involving D 2∆χj D 2∆χk in 50), thereby obtaining (see also [52]):
ω0⋆ = 4 −
Nχ
− ψ2
ψ2
(51)
18
Stefano Mariani
Figure 5: impact on a two-layer composite. Space-time diagram (the vertical dashed line
here represents the possibly debonding surface when a brittle, homogeneous material is
subject to the same impact).
In what precedes we have assumed the mapping f to be analytic everywhere in the x
space. If f is not differentiable, low order terms of the Taylor series expansions of mean
and covariance get affected. It is difficult to quantify the discrepancies with respect to
the analytic case, because they depend on whether the sigma-points sample the loci of nondifferentiability. However, it can be generally said that the order of accuracy is detrimentally
affected.
4 Results
To assess the performances of the proposed filtering approach in calibrating the interlaminar constitutive law while detecting impact-induced delamination, we first study a simple
problem consisting of a two-layer composite stricken by a homogeneous impactor. Hence,
two different impact tests on GRP composites [14, 53] are considered to mainly show the
accuracy in detecting delamination in real-time.
In all the cases, it is assumed that the contact between specimen and impactor is perfect (i.e. distributed all over their approaching surfaces) and that failure of the laminate
occurs because of the propagation of dilatational plane waves in the through-the-thickness
direction: the interlaminar surfaces are therefore subject to pure mode I loading.
4.1 Pseudo-experimental testings
As a starting benchmark, a pseudo-experimental testing condition is conceived. The
pseudo-experimental response of the laminate to the impact loading has been computed
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
19
by adding a white noise (of assigned variance) to sampled outcomes of finite element analyses. Even though this approach is sometimes criticized, being not clear whether one is
testing with it the filtering approach or its implementation, it helps in getting insights into
the performance of the filter in terms of stability and convergence rate. Indeed, calibration
of interlaminar constitutive models may become difficult if delamination occurs almost instantaneously: the filter has to be highly sensitive to model parameters to promptly react
to the information conveyed by measurements. This requires a careful setting of filter parameters, like P 0 and V . It is worth mentioning that, to further complicate the problem,
the effects of the shape of the tensile envelope of interlaminar laws (at assigned strength
and toughness) on the overall response of a laminate subject to impacts have not been thoroughly understood yet: depending on the loading and boundary conditions, on the composite geometry and on the stacking sequence, in some circumstances the shape affects the
response, whereas in others it does not [54].
The capabilities of the SPKF are therefore first assessed through a simple test: a twolayer composite is stricken by a homogeneous impactor. The laminae and the impactor are
assumed to be isotropic and elastic, featuring Young’s modulus E = 10 GPa, Poisson’s ratio
ν = 0.35 and mass density ̺ = 1500 kg/m3 (see also [11]). Each lamina and the impactor
are 0.75 mm in thickness. Target mechanical properties of the interlaminar surfaces are
assumed:
K =277.09
(N/mm3 )
τM =75
G =0.15
(MPa)
(52)
(N/mm)
According to the space-time diagram of Figure 5, failure can occur only along the interlaminar surface because of the interaction of the two release waves propagating inwards from
the free surfaces of impactor and specimen.
Two different values of the velocity v of the impactor are considered. In a first case
v = 10.19 m/s leads to the propagation in the through-the-thickness direction of a compressive/tensile wave of amplitude τ̄ = 50 MPa, which does not cause interface failure. In
a second case v = 20.38 m/s causes laminate failure, i.e. whole delamination, because of
the propagation of a compressive/tensile wave of amplitude τ̄ = 100 MPa. Outcomes of
the two tests are respectively reported in Figures 6 and 7, in terms of time evolution of the
free surface velocity velocity u̇r at the rear laminate surface, of the opening displacement
discontinuity [u] and of the normal traction τ (here and in what follows the subscript n has
been dropped to simplify the notation). Results are shown for all the constitutive models
described in Section 2.2, having assumed q = 1 for the modified exponential interface law
(17). For comparison purposes, the response of an interface-free specimen is reported too;
in such a case, the purely elastic behavior of the material leads to the propagation of sharp
fronts of a shock wave.
If delamination is not incepted (τ̄ = 50 MPa), the response is almost independent of
the shape of the interface law. Only in the presence of an interface that behaves according
to the exponential model, the pre-peak nonlinearity of the constitutive law (see Figure 2(a))
leads to a larger opening [u] of the interlaminar surface when subject to a tensile stress. For
20
Stefano Mariani
(a)
(b)
(c)
Figure 6: impact on a two-layer composite, τ̄ = 50 MPa. Effects of the interlaminar laws on
the time evolution of (a) the free surface velocity u̇r , of (b) the displacement discontinuity
[u] and (c) relevant traction τ along the interface.
any constitutive model, the signature of the interface is shown in Figure 6(a) by the delay in
the sudden changes of u̇r with respect to the reference, interface-free solution. This delay,
which grows in time, is caused by the compliance of the interface, that is additional to the
bulk one.
In case of failure (τ̄ = 100 MPa), the free surface velocity u̇r is affected by the interface
model only when the waves, traveling across the interlaminar surface while softening is
taking place, reach the rear surface; this occurs in the present case around 1µs after the
impact, see Figure 7(a). After failure, the rear lamina detaches from the front one and
freely flies off, as testified by the diverging [u] history in Figure 7(b). Part of the shock
waves then get confined inside the back lamina: this explains the subsequent doubling of
the drops in the u̇r evolution. It is worth noting that the time elapsed between softening
inception (t ∼
= 0.72 µs) and whole failure of the interface (t ∼
= 0.9 µs), see Figure 7(c), is
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
21
(a)
(b)
(c)
Figure 7: impact on a two-layer composite, τ̄ = 100 MPa. Effects of the interlaminar laws
on the time evolution of (a) the free surface velocity u̇r , of (b) the displacement discontinuity [u] and (c) relevant traction τ along the interface.
very short; only the information on this failure event, constituting the so-called pull-back
signal (PBS), conveyed by the shock waves to the free laminate surface, can be used by the
filter to calibrate the interface law.
The effects on the PBS of the shape of the tensile envelope and of the strength τM and
toughness G values need to be assessed. In the absence of any dissipation mechanisms, in
this test the free surface velocity u̇r would drop to zero at t = 0.917 µs (see Figure 6(a));
in case of delamination, the minimum attained velocity after the arrival of the unloading
tensile wave and the shape of the PBS do furnish information on the interface response. To
understand the roles of interface law, strength and toughness, the results of a parametric
analysis are shown in Figure 8: for any interface model, τM and G are varied by 20%
at most with respect to the target values (52). The piecewise linear and linear-exponential
laws, having a common initial elastic phase in tension, lead to a common descending branch
in the PBS. At variance, the local tangent stiffness of the exponential law in the hardening
22
Stefano Mariani
(a)
(b)
(c)
(d)
(e)
(f)
Figure 8: impact on a two-layer composite, τ̄ = 100 MPa. Effects of the interface strength
τM (left column) and toughness G (right column) on the pull-back signal. (a-b) piecewise
linear law (12); (c-d) linear-exponential law (13); (e-f) exponential law (14).
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
23
phase is affected by τM and G: the slope of the descending branch of the PBS is therefore
affected by τM and G too. Independently of the interface law, τM turns out to affect the
starting stage of debonding, whereas G affects its tail and the time needed to complete it;
this is clearly evidenced by the PBS, since τM modifies the minimum attained velocity,
whereas G influences only the ascending branch with no effects on the pull-back velocity.
From a model calibration perspective, it is clear that the SPKF has chances to improve
parameter estimates only in the time interval 0.9 ≤ t ≤ 1.2 µs. It is therefore hard to
figure out from the previous plots whether the effects of tensile envelope, τM and G can
be actually interpreted by the filter to improve model calibration. The performances of
the SPKF are hence tested here not only looking at parameter estimates, but also checking
its capability to detect whether a laminate is failing and, in case of delamination, where it
actually takes place.
Typical results of the filtering procedure are depicted in Figures 9 and 10; in this case
the pseudo-experimental data, which consists in the free surface velocity alone, have been
supposed very accurate, featuring a standard deviation δ = 0.33 m/s (the measurement
error covariance matrix W becomes scalar-valued, with entry W = δ2 ). As far as the
process covariance matrix V is concerned, in case of pseudo-experimental testing it can
be assumed to be vanishing, since the filter employs the same structural model adopted to
get the pseudo-experimental data. Components of P 0 instead need to be finely tuned to
enhance filter convergence [11, 15].
Figure 9 shows the obtained estimates of τM and G as a function of their initial guess in
x̂0 (here respectively denoted by τM,0 and G0 ), for all the interface models. These estimates
evolve from τM,0 and G0 once the PBS is processed by the filter; after that, they become
stationary. The tracked state of the laminate is shown in Figure 10 in terms of predicted
interface opening [u] and free surface velocity u̇r . Knowing the target response of the
composite to the impact, here denoted by the dashed lines, allows to certify stability and
convergence of the SPKF, no matter if displacement is diverging in a part of the system and
what kind of interface constitutive law has been adopted. It can be seen that estimates get
enhanced as soon as the filter senses the PBS: in fact, the sudden changes in the estimate of
[u] show up only while processing the PBS, starting from t ∼
= 0.9 µs.
These outcomes testify that the SPKF is very efficient in tracking the state of the laminate. i.e. in understanding whether the structure is failing or not. Model calibration is
instead less accurately accomplished: independently of the interface law, τM is quite precisely estimated, provided the initial guess τM,0 is not too far from the target value, whereas
G can be hardly inferred. These conclusions are in agreement with the results of the parametric analysis: while τM affects the whole PBS, G affects only its ascending branch.
Therefore, filtering out from u̇r the effects of G alone turns out to be extremely difficult.
In case of a much higher scattering of pseudo-experimental data (δ = 3.3 m/s), see
Figure 11, results loose accuracy as for the calibration task. Contrariwise, state tracking
maintain accuracy: even though measurements contain poor information, the SPKF is again
able to provide the evolution of the free surface velocity in the PBS.
24
Stefano Mariani
(a)
(b)
(c)
(d)
(e)
(f)
Figure 9: impact on a two-layer composite, τ̄ = 100 MPa (W=10−1 m2 /s2 ). Effects of the
initialization values τM,0 and G0 on the converged estimates of τM (left column) and G
(right column). (a-b) piecewise linear law; (c-d) linear-exponential law; (e-f) exponential
law.
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
25
(a)
(b)
(c)
(d)
(e)
(f)
Figure 10: impact on a two-layer composite, τ̄ = 100 MPa (W=10−1 m2 /s2 ). Evolution in
time of interface opening [u] (left column) and free surface velocity u̇r (right column), and
comparison among tracked state (orange squares), actual state (dashed lines) and pseudoexperimental data (blue circles). (a-b) piecewise linear law; (c-d) linear-exponential law;
(e-f) exponential law.
26
Stefano Mariani
(a)
(b)
(c)
(d)
(e)
(f)
Figure 11: impact on a two-layer composite, τ̄ = 100 MPa (W=10 m2 /s2 ). Evolution in
time of interface opening [u] (left column) and free surface velocity u̇r (right column), and
comparison among tracked state (orange squares), actual state (dashed lines) and pseudoexperimental data (blue circles). (a-b) piecewise linear law; (c-d) linear-exponential law;
(e-f) exponential law.
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
27
(a)
(b)
Figure 12: impact on a 7-layer composite [53]. (a) space-time diagram (the vertical dashed
line here represents the possibly debonding surface when a brittle, homogeneous material
is subject to the same impact), and (b) experimentally measured free surface velocity.
4.2 Actual experimental testings
To finally check the performances of the SPKF when dealing with multi-layered composites, we consider two of the experiments reported in [53] and [14].
In the first experiment (experiment FY06001 in [53]), the specimen is a 7-layer composite plate; each lamina is 1.37 mm in thickness, and is made of a balanced 5-harness satin
weave E-glass and LY564 epoxy. The wave speed in the through-the-thickness direction is
3.34 km/s, while the mass density is ̺ = 1885 kg/m3 [53]. This laminate was subject to a
plane impact, stricken by an aluminum impactor (12.5 mm thick) flying at velocity v = 71
m/s. The relevant space-time diagram is shown in Figure 12, along with the free surface
velocity profile measured via a velocity interferometer for any reflector (VISAR).
Account taken of the high accuracy of the experimentally measured u̇r , the identifica-
28
Stefano Mariani
(a)
(b)
(c)
(d)
(e)
(f)
Figure 13: impact on a 7-layer composite [53]. Evolution in time of the estimated values of
τM (left column) and G (right column). (a-b) piecewise linear law; (c-d) linear-exponential
law; (e-f) exponential law.
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
29
(a)
(b)
(c)
(d)
(e)
(f)
(g)
Figure 14: impact on a 7-layer composite [53]. Evolution in time of (a) free surface velocity
u̇r and (b-g) estimated interface openings [u]1 − [u]6 .
30
Stefano Mariani
tion procedure have furnished the results reported in Figure 13 in terms of time evolution
of estimates of τM and G as a function of their initialization values within the domain:
Cϑ = {50 ≤ τM,0 ≤ 250 (MPa), 0.5 ≤ G0 ≤ 2.5 (N/mm)}
(53)
Converged estimates of τM are in good agreement with the spall strength of 119.5 MPa
reported in [53]; on the other hand, final estimates of G are well representative for this kind
of composites. Figure 14 reports the estimated state of the specimen: the capability to track
the measured free surface velocity and to foresee delamination along the third interlaminar
surface away from the impact plane, is evidenced. This latter result, allowing also for wave
dispersion caused by interlaminar surfaces and inner inhomogeneities of the composite,
well agrees with the state-space diagram of Figure 12(a).
In the second experiment (experiment 1 in [14]), a GRP specimen, 7.02 mm thick, is
backed by another GRP plate, 6.91 mm thick; both laminates are made of 11 plies. The
wave speed in the through-the-thickness direction now amounts to 3.19 km/s, and the mass
density to ̺ = 1867 kg/m3 [14]. The specimen is stricken by a 5-layer GRP flyer, 2.96 mm
in thickness, flying at velocity v = 85 m/s. The corresponding space-time diagram, and the
free surface velocity profile measured via a VISAR are reported in Figure 15. Because of
the test set-up, the release waves interact causing delamination inside the back plate.
Results of the filtering process are reported in Figure 16 is terms of tracked free surface
velocity and estimated displacement jumps along all the interfaces inside the back plate
(sequence starts at the specimen-back plate contact surface). These estimations turn out
once again to be independent of the interface law and of the initialization values of τM and
G inside the domain:
Cϑ = {25 ≤ τM,0 ≤ 100 (MPa), 0.1 ≤ G0 ≤ 0.6 (N/mm)}
(54)
While the free surface velocity is accurately tracked, delamination is foreseen to take place
along the 7-th interlaminar surfaces, in agreement with the results of [11, 14]. As far as
model calibration is concerned, outcomes turn out to be qualitatively in agreement with
those already reported for the previous tests.
5 Conclusion
In this Chapter we have addressed some issues related to constitutive modeling and parameter identification in finite element simulations of layered composites subject to impacts.
Assuming the impact energy to be high enough to cause damage spreading inside the interlaminar resin-enriched phases, but not high enough to result in penetration of the impactor
accompanied by intralaminar damage, a numerical scheme for structural-level analyses has
been revised. Within this scheme the laminae are assumed to behave elastically, whereas
dissipation mechanisms are lumped onto zero-thickness interlaminar surfaces. Along these
interlaminar surfaces strength reduction, eventually leading to delamination, is governed by
softening interface constitutive laws linking tractions to displacement jumps.
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
31
(a)
(b)
Figure 15: impact on a 11+11-layer composite [14]. (a) space-time diagram (the vertical
dashed line here represents the possibly debonding surface when a brittle, homogeneous
material is subject to the same impact), and (b) experimentally measured free surface velocity.
32
Stefano Mariani
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
Figure 16: impact on a 11+11-layer composite [14]. Evolution in time of (a) free surface
velocity u̇r and (b-k) estimated interface openings [u]1 − [u]10 in the back plate.
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
33
Interface laws are known to be difficult to calibrate, since direct testing on a single
interlaminar phase can not be devised. Here we have offered a sigma-point Kalman filtering
approach to estimate uncertain model parameters of the aforementioned interface laws. This
technique outperforms most of the customarily adopted ones, since it efficiently deals with
nonlinearities, which are a result of interlaminar strength degradation in the case under
study.
The performances of the filtering procedure have been assessed through pseudoexperimental testings on a two-layer composite, and through real testings on multi-layer
glass fiber reinforced plastic composites. It has been shown that the state of the composite,
including delamination inception and growth, is always tracked with a noteworthy level of
accuracy. Due to the fast failure events, model calibration is instead less accurately performed and sometimes requires initialization values of uncertain model parameters to be
appropriately chosen to avoid biased estimates.
References
[1] S. Abrate. Impact on composite structures. Cambridge University Press, 1998.
[2] H.D. Espinosa, S. Dwivedi, P.D. Zavattieri, and G. Yuan. A numerical investigation of
penetration in multilayered material/structure systems. International Journal of Solids
and Structures, 35:2975–3001, 1998.
[3] O. Allix and P. Ladevéze. Interlaminar interface modelling for the prediction of delamination. Composite Structures, 22:235–242, 1992.
[4] A. Corigliano. Formulation, identification and use of interface elements in the numerical analysis of composite delamination. International Journal of Solids and Structures, 30:2779–2811, 1993.
[5] A. Corigliano and O. Allix. Some aspects of interlaminar degradation in composites.
Computer Methods in Applied Mechanics and Engineering, 185:203–224, 2000.
[6] A. Corigliano and S. Mariani. Simulation of damage in composites by means of interface models: parameter identification. Composites Science and Technology, 61:2299–
2315, 2001.
[7] G. Alfano and M.A. Crisfield. Finite element interface models for the delamination
analysis of laminated composites: mechanical and computational issues. International
Journal for Numerical Methods in Engineering, 50:1701–1736, 2001.
[8] A. Corigliano. Damage and fracture mechanics techniques for composite structures. In
I. Milne, R.O. Ritchie, and B. Karihaloo, editors, Comprehensive Structural Integrity,
volume 3, chapter 9, pages 459–539. Elsevier Science, 2003.
34
Stefano Mariani
[9] A. Corigliano and S. Mariani. Parameter identification of a time-dependent elasticdamage interface model for the simulation of debonding in composites. Composites
Science and Technology, 61:191–203, 2001.
[10] A. Corigliano and S. Mariani. Identification of a constitutive model for the simulation of time-dependent interlaminar debonding processes in composites. Computer
Methods in Applied Mechanics and Engineering, 191:1861–1894, 2002.
[11] S. Mariani and A. Corigliano. Impact induced composite delamination: state and parameter identification via joint and dual extended Kalman filters. Computer Methods
in Applied Mechanics and Engineering, 194:5242–5272, 2005.
[12] S. Bittanti, G. Maier, and A. Nappi. Inverse problems in structural elastoplasticity:
a Kalman filter approach. In A. Sawczuk and G. Bianchi, editors, Plasticity today.
Modelling, Methods and Applications, pages 311–329. Elsevier Applied Science Publishers, 1984.
[13] A. Corigliano, S. Mariani, and B. Orsatti. Identification of Gurson-Tvergaard material
model parameters via Kalman filtering technique - I. Theory. International Journal of
Fracture, 104:349–373, 2000.
[14] H.D. Espinosa, S. Dwivedi, and H.-C. Lu. Modeling impact induced delamination of
woven fiber reinforced composites with contact/cohesive laws. Computer Methods in
Applied Mechanics and Engineering, 183:259–290, 2000.
[15] A. Corigliano and S. Mariani. Parameter identification in explicit structural dynamics:
performance of the extended Kalman filter. Computer Methods in Applied Mechanics
and Engineering, 193:3807–3835, 2004.
[16] S.J. Julier, J. Uhlmann, and H.F. Durrant-Whyte. A new method for the nonlinear
transformation of means and covariances in filters and estimators. IEEE Transactions
on Automatic Control, 45:477–482, 2000.
[17] T. Lefebvre, H. Bruyninckx, and J. De Schutter. Kalman filters for nonlinear systems:
a comparison of performance. International Journal of Control, 77:639–653, 2004.
[18] K. Ito and K. Xiong. Gaussian filters for nonlinear filtering problems. IEEE Transactions on Automatic Control, 45:910–927, 2000.
[19] S. Mariani and A. Ghisi. Unscented Kalman filtering for nonlinear structural dynamics. Nonlinear Dynamics, 49:131–150, 2007.
[20] A. Corigliano, S. Mariani, and A. Pandolfi. Numerical modeling of rate-dependent
debonding processes in composites. Composite Structures, 61:39–50, 2003.
[21] A. Corigliano, S. Mariani, and A. Pandolfi. Numerical analysis of rate-dependent
dynamic composite delamination. Composites Science and Technology, 66:766–775,
2006.
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
35
[22] R. Lund and J. Fish. Space-time multiscale laminated theory. International Journal
of Multiscale Computational Engineering, 2:26–47, 2004.
[23] S. Mariani, A. Pandolfi, and R. Pavani. Coupled space-time multiscale simulations of
dynamic delamination tests. Materials Science, 23:509–519, 2005.
[24] R.R. Talreja Y.W. Kwon, D.H. Allen, editor. Multiscale Modeling and Simulation of
Composite Materials and Structures. Springer, 2007.
[25] G.T. Camacho and M. Ortiz. Computational modelling of impact damage in brittle
materials. International Journal of Solids and Structures, 33:2899–2938, 1996.
[26] Y. Arun Roy and R.H. Dodds. Simulation of ductile crack growth in thin alluminum panels using 3-D surface cohesive elements. International Journal of Fracture,
110:21–45, 2001.
[27] S. Mariani and U. Perego. Extended finite element method for quasi-brittle fracture.
International Journal for Numerical Methods in Engineering, 58:103–126, 2003.
[28] C. Comi and S. Mariani. Extended finite element simulation of quasi-brittle fracture
in functionally graded materials. Computer Methods in Applied Mechanics and Engineering, 196:4013–4026, 2007.
[29] S. Li, M.D. Thouless, A.M. Waas, J.A. Schroeder, and P.D. Zavattieri. Mixed-mode
cohesive-zone models for fracture of an adhesively bonded polymer matrix composite.
Engineering Fracture Mechanics, 73:64–78, 2006.
[30] J.H. Rose, J. Ferrante, and J.R. Smith. Universal binding energy curves for metals and
bimetallic interfaces. Physical Review Letters, 47:675–678, 1981.
[31] J.H. Rose, J.R. Smith, and J. Ferrante. Universal features of bonding in metals. Physical Review B, 28:1835–1845, 1983.
[32] X.-P. Xu and A. Needleman. Numerical simulation of fast crack growth in brittle
solids. Journal of the Mechanics and Physics of Solids, 42:1397–1434, 1994.
[33] A. Needleman and A.J. Rosakis. The effect of bond strength and loading rate on
the conditions governing the attainment of intersonic crack growth along interfaces.
Journal of the Mechanics and Physics of Solids, 47:2411–2449, 1999.
[34] C.M. Landis, T. Pardoen, and J.W. Hutchinson. Crack velocity dependent toughness
in rate dependent materials. Mechanics of Materials, 32:663–678, 2000.
[35] O. C. Zienkiewicz and R. L . Taylor. The finite element method, volume II. McGrawHill, London, 1989.
36
Stefano Mariani
[36] N. Moës, J. Dolbow, and T. Belytschko. A finite element method for crack growth
without remeshing. International Journal for Numerical Methods in Engineering,
46:131–150, 1999.
[37] T. Strouboulis, K. Copps, and I. Babuška. The generalized finite element method.
Computer Methods in Applied Mechanics and Engineering, 190:4081–4193, 2001.
[38] G.N. Wells and L.J. Sluys. A new method for modelling cohesive cracks using finite
elements. International Journal for Numerical Methods in Engineering, 50:2667–
2682, 2001.
[39] C. Comi, S. Mariani, and U. Perego. An extended FE strategy for transition from
continuum damage to mode I cohesive crack propagation. International Journal for
Numerical and Analytical Methods in Geomechanics, 31:213–238, 2007.
[40] M. Jirasek and T. Zimmerman. Embedded crack model: II. Combination with smeared
cracks. International Journal for Numerical Methods in Engineering, 50:1291–1305,
2001.
[41] H.M. Hilber, T.J.R. Hughes, and R.L. Taylor. Improved numerical dissipation for time
integration algorithms in structural dynamics. Earthquake Engineering and Structural
Dynamics, 5:283–292, 1977.
[42] T.J.R. Hughes. The finite element method. Linear static and dynamic finite element
analysis. Dover Publications, 2000.
[43] L. Ljung. System identification. Theory for the user. Prentice Hall, second edition,
1999.
[44] A. Gelb. Applied optimal estimation. The M.I.T. Press, 1974.
[45] E.A. Wan and A.T. Nelson. Dual extended Kalman filter methods. In S. Haykin,
editor, Kalman filtering and neural networks, pages 123–173. John Wiley & Sons,
Inc., 2001.
[46] R. E. Kalman and R. S. Bucy. New results in linear filtering and prediction theory.
ASME Journal of Basic Engineering, 83D:95–108, 1961.
[47] L. Ljung. Asymptotic behavior of the extended Kalman filter as a parameter estimator
for linear systems. IEEE Transactions on Automatic Control, AC-24:36–50, 1979.
[48] K. Reif, S. Günther, E. Yaz, and R. Unbehauen. Stochastic stability of the discretetime extended Kalman filter. IEEE Transactions on Automatic Control, 44:714–728,
1999.
[49] S.J. Julier and J.K. Uhlmann. A general method for approximating nonlinear transformations of probability distributions. Technical report, RRG, Department of Engineering Science, University of Oxford, 1996.
Failure of layered composites subject to impacts: constitutive modeling and parameter
identification issues
37
[50] S.J. Julier, J. Uhlmann, and H.F. Durrant-Whyte. A new approach for filtering nonlinear systems. In Proceedings of the American Control Conference, pages 1628–1632,
Seattle, June 1995.
[51] E.A. Wan and R. van der Merwe. The unscented Kalman filter. In S. Haykin, editor,
Kalman filtering and neural networks, pages 221–280. John Wiley & Sons, Inc., 2001.
[52] S.J. Julier. The scaled unscented transformation. In Proceedings of the American
Control Conference, volume 6, pages 4555–4559, Anchorage, May 2002.
[53] F. Yuan, L. Tsai, V. Prakash, A.M. Rajendran, and D.P. Dandekar. Spall strength
of glass fiber reinforced polymer composites. International Journal of Solids and
Structures, 44:7731–7747, 2007.
[54] G. Alfano. On the influence of the shape of the interface law on the application of
cohesive-zone models. Composites Science and Technology, 66:723–730, 2006.