Comments
Description
Transcript
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.