Next Article in Journal
Bayesian Maximum Entropy Based Algorithm for Digital X-ray Mammogram Processing
Next Article in Special Issue
Automated Modelling of Evolving Discontinuities
Previous Article in Journal
Computer-Aided Diagnosis in Mammography Using Content-Based Image Retrieval Approaches: Current Status and Future Perspectives
Previous Article in Special Issue
Probabilistic Upscaling of Material Failure Using Random Field Models – A Preliminary Investigation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Failure Assessment of Layered Composites Subject to Impact Loadings: a Finite Element, Sigma-Point Kalman Filter Approach

Politecnico di Milano, Dipartimento di Ingegneria Strutturale, Piazza Leonardo da Vinci 32, 20133 Milano, Italy
Algorithms 2009, 2(2), 808-827; https://doi.org/10.3390/a2020808
Submission received: 1 April 2009 / Accepted: 27 May 2009 / Published: 4 June 2009
(This article belongs to the Special Issue Numerical Simulation of Discontinuities in Mechanics)

Abstract

:
We present a coupled finite element, Kalman filter approach to foresee impact-induced delamination of layered composites when mechanical properties are partially unknown. Since direct numerical simulations, which require all the constitutive parameters to be assigned, cannot be run in such cases, an inverse problem is formulated to allow for modeling as well as constitutive uncertainties. Upon space discretization through finite elements and time integration through the explicit α - method, the resulting nonlinear stochastic state model, wherein nonlinearities are due to delamination growth, is attacked with sigma-point Kalman filtering. Comparison with experimental data available in the literature and concerning inter-laminar failure of layered composites subject to low-velocity impacts, shows that the proposed procedure leads to: an accurate description of the failure mode; converged estimates of inter-laminar strength and toughness in good agreement with experimental data.

1. Introduction

Composite structures are often exposed to impact loadings during their life cycle [1]; this kind of loading conditions may become a serious issue, since composites can be very sensitive to imperfections beyond a critical threshold load. Depending on the impact energy, i.e. on the momentum of the impactor when it strikes the structural component, two different failure scenarios can be expected: (i) in case of so-called low velocity (and, therefore, low energy) impacts, stress waves propagate inside the composite leading to diffused damage processes that progressively decrease or even annihilate the residual strength and stiffness of the structure; (ii) in case of high velocity (and, therefore, high energy) impacts, local perforation of the laminate occurs, giving rise to very localized failure mechanisms. In the former case the response of a laminated plate to the impact is characterized by forced vibrations of the whole structure; in the latter case, the response of a plate and the damaging processes are very localized around the impact location, and deformation modes of the structure as a whole are only marginally excited.
Having as a target the development of a real-time health monitoring (and, hopefully, control) strategy for composite structures, in this work we focus on low velocity impact conditions. In such cases, the aforementioned diffused degradation processes consist in inter-laminar debonding (i.e. delamination) and intra-laminar damage, the latter being related to fibers failure, matrix cracking and fiber/matrix decoherence. Typically, a major role is played by delamination; therefore, throughout the whole paper we shall assume the laminae of the composite to behave elastically, while reduction of the structural load-carrying capacity may occur only because of inter-laminar damage.
To simplify modeling, the small ratio between the thickness of each inter-laminar phase and the thickness of the whole laminate is exploited by numerically treating inter-laminar phases as lumped interfaces. The relevant constitutive laws locally link the tractions acting upon the surfaces of the two adjacent layers to the displacement jump occurring across the interface, allowing for both opening and sliding (or shear) deformation modes. This lumped interface formulation has to be considered as a coarse grained description of the actual inter-phase behavior, as obtained via homogenization along the thickness direction (i.e. in the direction perpendicular to the stacking plane); the displacement jump thus accounts for both the elastic response of the inter-phase material and the possible discontinuity in the displacement field due to micro-cracking. To model micro-cracking processes and the relevant reduction of the inter-laminar strength, nonlinear constitutive laws are adopted; beyond the attainment of peak tractions, delamination inception is modeled through a softening regime. Since inter-phase materials are usually quasi-brittle, dissipative phenomena due to yielding can be disregarded. Moreover, viscous and rate-dependent properties of the composite phases are not considered here; therefore, only rate-independent constitutive laws will be adopted.
Within this framework, laminated plates are treated as a stacking sequence of laminae and interfaces. This scheme, sometimes referred to as meso-scale modeling (see, e.g. [2]), was proved to furnish accurate outcomes at the length scale of the whole laminate (also in terms of delamination prediction), provided the interface constitutive laws are appropriately calibrated. Since they are typically formulated on a pure phenomenological basis and contain several parameters to be tuned, this calibration task turns out to be difficult, also because of the lower dimensionality of interfaces (two-dimensional loci in a three-dimensional setting) and because of shadowing effects due to the global laminate response, see e.g. [3].
To accurately predict the failure mode of layered composites subject to low velocity impacts and simultaneously calibrate interface constitutive laws, a coupled finite element (FE) Kalman filter approach is here proposed. Upon space discretization (through displacement-based finite elements) and time integration (through the explicit α - method, in order to damp spurious high frequency oscillations in the fully discretized solution), Kalman filtering is adopted to fully track the structural state in a stochastic frame [4,5,6], exploiting system evolution to continuously improve estimates. Because of the nonlinearities induced by inter-laminar softening, the customarily adopted extended Kalman filter (EKF) is not able to provide robust state tracking and accurate estimations of inter-laminar properties like e.g. strength and toughness [7]. We therefore attack the problem with the newly proposed sigma-point Kalman filter (SPKF) [8,9,10], which avoids the occurrence of unstable evolution of the estimates (see also [11]). While the EKF linearizes the system evolution equations (through a Taylor series expansion of the relevant mapping, thought to be analytic everywhere, arrested at first order) and accordingly updates the statistics of the structural state, the SPKF samples the probability distribution of the state via an appropriately chosen set of sigma-points, which are then evolved according to the actual nonlinear system dynamics to obtain the updated statistics of the state itself. Among alternative approaches to system identification and, specifically, to model calibration, it is worth quoting here [12].
In this work the capability of the proposed FE-SPKF approach is assessed through comparison with experimental data available in the literature and concerning failure of laminated plates subject to plane impacts [13,14]; in such tests, the composites are mainly stressed in the region surrounding the impact location by dilatational waves propagating in the through-the-thickness direction and causing mode I delamination. Outcomes show that the failure mode is always correctly captured, and estimated values of inter-laminar strength and toughness well agree with available data.
Plane impact tests do not allow a complete calibration of the interface constitutive law, since only model parameters relevant to mode I failure can be estimated. Anyway, forthcoming results show that under very adverse conditions, as due to fast completion of the whole debonding process, the SPKF proves to be robust and accurate; on the contrary, under similar conditions the EKF was shown in [7,15] to either diverge or provide biased estimates. An extensive assessment of the SPKF performance in the presence of diffused damage/delamination patterns is beyond the scope of this paper and will be presented in future works.
As far as notation is concerned, a matrix one will be adopted throughout, with uppercase and lowercase bold symbols respectively denoting matrices and vectors. A superscript T will stand for transpose, | | will denote the norm of vector □, and a superposed dot will represent time rates.

2. Forward problem: layered composite dynamics

Let us consider a layered body Ω, possessing a smooth outer boundary with unit outward normal m (see Figure 1). Ω be crossed by n Γ non-intersecting surfaces, or interfaces Γ j , j = 1 , . . . , n Γ . Even though the local orientation n j of these interfaces may vary in Ω, for ease of notation we shall assume it to be constant and equal to n along every Γ j : interfaces are thus flat and parallel to each other.
Figure 1. Notation adopted for layered bodies (here n Γ = 2 ).
Figure 1. Notation adopted for layered bodies (here n Γ = 2 ).
Algorithms 02 00808 g001
The equilibrium of Ω at time t is governed by:
C T σ + b ¯ = ϱ u ¨ in Ω \ Γ
N σ = - τ j on Γ j + , N σ = τ j on Γ j - , j = 1 , . . . , n Γ
M σ = τ ¯ on Γ τ
where: Γ = j = 1 n Γ Γ j ; Γ τ is the portion of the outer boundary of Ω where tractions τ ¯ are assigned. Treating each Γ j as a locus where the displacement field u may suffer jumps (alike a crack), Γ j + and Γ j - are the two sides of Γ j ; these two sides can not be distinguished in the initial configuration (at t = 0 ) but, as soon as delamination occurs, they experience a relative movement. In the above equations, according to Voigt notation: σ is the stress vector; τ j is the traction vector acting along surface Γ j ; b ¯ is the assigned bulk load in Ω \ Γ ; ϱ is the bulk mass density; u ¨ is the acceleration field in Ω \ Γ ; C is the differential compatibility operator; M and N are the matrices collecting the components of unit vectors m and n , respectively.
Assuming linearized kinematics, compatibility reads:
ε = C u in Ω \ Γ
[ u ] j = u | Γ j + - u | Γ j - j = 1 , . . . , n Γ
Figure 2. Piecewise linear (PWL) and linear-exponential (L-E) effective traction-displacement discontinuity relationships.
Figure 2. Piecewise linear (PWL) and linear-exponential (L-E) effective traction-displacement discontinuity relationships.
Algorithms 02 00808 g002
u = u ¯ on Γ u
where: ε is the strain vector; u ¯ is the assigned displacement along portion Γ u of the outer boundary of Ω; [ u ] j is the displacement discontinuity along surface Γ j .
The body is assumed to be initially at rest, in an undeformed and unstressed state:
u 0 = 0 , u ˙ 0 = 0 in Ω
In case of low velocity impacts, dissipation due to intra-laminar damage is assumed to be negligible when compared to that due to delamination. Hence, laminae are thought to behave elastically according to:
σ = E Ω ε in Ω \ Γ
where E Ω is the elasticity matrix of the bulk.
As for inter-laminar phases, to be treated as lumped surfaces along which debonding is the result of complex micro-scale damaging processes, a phenomenological coarse-grained modeling is adopted. Strength reduction preceding delamination is governed by softening constitutive laws [3,16,17,18,19,20,21], which locally link tractions τ j (see Eq. (2) to the displacement jump [ u ] j (see Eq. (5). Since both opening/closing (along n ) and sliding (in the s 1 j - s 2 j plane) discontinuities are expected to take place under general loading conditions, interface kinematics is effectively represented by the displacement discontinuity [ u ] (in what follows, for ease of notation, we shall drop subscript j), defined as [22,23,24]:
[ u ] = [ u ] n 2 + κ 2 [ u ] s 2
where: [ u ] n = [ u ] T n and [ u ] s = [ u ] - [ u ] n n are the opening and sliding displacement discontinuities, respectively; κ is a coupling coefficient, which phenomenologically accounts for the local ongoing mixed-mode damaging in the inter-laminar phase. The effective traction τ, work-conjugate to [ u ] , reads (see [22,23]):
τ = τ n 2 + τ s 2 κ 2
where τ n and τ s are the normal and shear traction components.
Since quasi-brittle materials (like inter-phase resins) show a high ratio between compressive and tensile strengths, inter-laminar damage in compression due to low velocity impacts can be disregarded. On the other hand, under tensile stress pulses the response of inter-phase materials is characterized by softening beyond the attainment of a peak traction. Two simple effective descriptions of this behavior are here envisaged: a piecewise linear (PWL) law
τ = K [ u ] if [ u ] [ u ] e τ = τ M + Q [ u ] - [ u ] e if [ u ] e < [ u ] [ u ] U τ = 0 if [ u ] > [ u ] U
and a linear-exponential (L-E) law
τ = K [ u ] if [ u ] [ u ] e τ = τ M exp - ς ( [ u ] - [ u ] e ) if [ u ] > [ u ] e
In these cohesive-like envelopes: [ u ] e is the effective displacement discontinuity at the attainment of the peak traction τ M = K [ u ] e (elastic limit); K is the stiffness in the elastic regime; Q is the (negative) slope of the softening branch in the PWL case; [ u ] U = 1 - K Q [ u ] e is the effective displacement discontinuity beyond which the interaction of opened crack faces is annihilated in the PWL case; ς defines the slope of the softening branch in the L-E case. The two effective models (11) and (12) are compared, for [ u ] > 0 , in Figure 2 at assigned effective fracture toughness, or work of separation G = 0 τ d [ u ] . It is worth noting that, since the L-E law assumes the interaction between crack faces to continue up to [ u ] , delamination growth is modeled via a breakdown threshold; according to [25,26], we assume this interaction to suddenly vanish as soon as the current traction τ becomes smaller than a pre-defined fraction (say 5-10 %) of the peak value τ M .
To allow for irreversibility of the micro-mechanical damage processes, unloading from the tensile envelope follows a radial path to the origin of the τ - [ u ] plane. In rate form, constitutive modeling for interface Γ j can thus be written:
τ ˙ j = E Γ [ u ˙ ] j
where E Γ is the interface tangent stiffness matrix. For additional details, readers are referred to [23,24].
Upon space discretization via displacement-based finite elements, the weak form of the equilibrium equations for the layered body reads:
M u ¨ h + K Ω u h + j = 1 n Γ R j = F
where u h is the vector of nodal displacements. In (14) the mass matrix M , the bulk stiffness matrix K Ω , the internal force vectors R j , j = 1 , . . . , n Γ , and the external load vector F are respectively given by:
M = Ω \ Γ ϱ Φ T Φ d Ω K Ω = Ω \ Γ B Ω T E Ω B Ω d Ω R j = Γ j B Γ j T τ j d Γ j F = Ω \ Γ Φ T b ¯ d Ω + Γ τ Φ T τ ¯ d Γ τ
Here: Φ is the matrix of nodal shape functions; B Ω and B Γ j are, respectively, the compatibility matrices for the bulk Ω \ Γ and surface Γ j [7,27].
In the space-discretized formulation here above we have implicitly assumed that the displacement field u h may be discontinuous only along inter-element boundaries. Smarter formulations [23,28,29,30,31] allow also for intra-element discontinuities by exploiting the partition of unity property of standard nodal shape functions [32]; since delamination in layered continua occurs along the a-priori known inter-laminar surfaces, the description of the discontinuity topology is simple and the aforementioned feature of the smart (usually termed extended) finite element formulation is of no help.
We now need to adopt a time stepping scheme for (14), able to reduce to a minimum or even avoid spurious high frequency oscillations of the velocity field in the fully-discretized solution. In [33] it was shown that local corrective flux terms can be adopted to accurately track discontinuities in the velocity field due to the propagation of stress wave fronts. Since this time stepping scheme is computationally expensive while our target is a low-cost scheme, a less expensive (and, indeed, less accurate) explicit α - method is adopted [34,35]. Having partitioned the time interval of interest according to [ t 0 t N ] = i = 0 N - 1 [ t i t i + 1 ] , the solution of (14) at time t i + 1 in terms of nodal displacements, velocities and accelerations is obtained according to:
  • prediction:
    u ˜ i + 1 h = u i h + Δ t u ˙ i h + Δ t 2 ( 1 2 - β ) u ¨ i h
    u ˙ ˜ i + 1 h = u ˙ i h + Δ t ( 1 - γ ) u ¨ i h
    where Δ t = t i + 1 - t i is the time step size;
  • integration:
    u ¨ i + 1 h = M - 1 F i + 1 + α - ( 1 + α ) K u ˜ i + 1 h + j R ˜ i + 1 j + α K u i + j R i j
    where: F i + 1 + α = F ( t i + ( 1 + α ) Δ t ) ; R i j and R ˜ i + 1 j are respectively computed (see Equation 15 3 ) adopting τ j , i = τ j ( [ u i h ] j ) and τ ˜ j , i + 1 = τ j ( [ u ˜ i + 1 h ] j ) ;
  • correction:
    u i + 1 h = u ˜ i + 1 h + Δ t 2 β u ¨ i + 1 h
    u ˙ i + 1 h = u ˙ ˜ i + 1 h + Δ t γ u ¨ i + 1 h
To get oscillations of small amplitude in the velocity field, the following values of the algorithmic parameters α, β and γ have been adopted: α = - 0 . 3 ; β and γ finely tuned around values allowing second-order accuracy in linear elastodynamics [36].
The time step size Δ t has been set so as to fulfill the Courant condition in the bulk. In case of growing delamination, interface elements in the process zone (where damage is locally reducing the inter-laminar strength and stiffness) need to be sized so as to accurately resolve the traction field; this requirement is satisfied if the characteristic size of interface elements is a fraction (say 1 3 - 1 5 ) of the process zone length, estimated e.g. through Irwin’s law [37]. Since the surrounding FEs in the bulk are associated to interface elements via matching grids, it turns out that inter-laminar strength and toughness, which affect the process zone length, actually affect the time step size Δ t too; for details, readers are referred to [38,39].

3. Inverse problem: constrained sigma-point Kalman filter

In our study the inverse problem has to provide real-time forecasts of the delamination processes, along with estimates of uncertain model parameters. Allowing for the above described explicit time integration scheme, a state vector x gathering all the information on the evolving system is obtained by joining the nodal displacement, velocity and acceleration vectors to a vector ϑ gathering model parameters in need of calibration:
x = u h u ˙ h u ¨ h ϑ
Since laminate failure due to delamination is represented by displacement discontinuities arising along the inter-laminar surfaces, vector x in (21) also provides the requested details of the failure mode. State vectors shaped in this way, with model parameters explicitly accounted for, are typically adopted in extended Kalman filtering [5,40]. While the current structural state (here u h , u ˙ h and u ¨ h ) is always at least partially observed, model parameters to be calibrated can not be directly measured; by joining structural and parameter vectors in x , state tracking can consistently enhance model calibration.
Allowing for modeling and measurement errors in a stochastic framework, the discretized state-space model of the system within the time interval [ t i t i + 1 ] reads:
x i + 1 = f i x i + v i y i = H x i + w i
where: y i is the observation vector, which collects measurables; v i and w i are the process and measurement noises. These noises are assumed to be additive, uncorrelated white and Gaussian processes, with zero mean and time-invariant covariance matrices V and W [6,40]. Because of inter-laminar strength degradation, mapping f i turns out to be highly nonlinear; on the other hand, since either surface displacements or velocities are measured during testing, the observation equation (22 2 ) shows up as a linear relation between y i and x i (see Eq. (21).
Table 1. Sigma-point Kalman filter.
Table 1. Sigma-point Kalman filter.
  • Initialization at t 0 :
    x ^ 0 = E [ x 0 ] P 0 = E [ x 0 - x ^ 0 x 0 - x ^ 0 T ]
  • At t i , for i = 0 , . . . , N
    • Prediction:
      χ ^ i , j = x ^ i + Δ χ i , j j = 0 , . . . , N χ χ ^ i + 1 , j - = f i ( χ ^ i , j ) x ^ i + 1 - = j = 0 N χ ω j χ ^ i + 1 , j - P i + 1 - = R i + 1 - + V
      where
      R i + 1 - = j = 0 N χ ω j χ ^ i + 1 , j - - x ^ i + 1 - χ ^ i + 1 , j - - x ^ i + 1 - T
    • Correction:
      x ^ i = x ^ i - + G i U y i - H x ^ i - P i = P i - - G i U H R i -
      where
      G i U = R i - H T H R i - H T + W - 1
The mentioned nonlinearities of mapping f i cause estimates obtained through a standard EKF to be affected by biases [41,42]; to get accurate results, the recently proposed SPKF [8,9,43,44] is instead here adopted. While the EKF linearizes the mapping f i , the SPKF deterministically samples the probability distribution of x at time t i (i.e. at the beginning of the time step) through a set of properly chosen sigma-points χ ^ i , j , j = 0 , . . . , N χ . These sigma-points are allowed to evolve according to the nonlinear mapping f i ; the statistics of x at time t i + 1 (i.e. at the end of the time step) are then obtained through a weighted averaging procedure performed on the evolved sigma-points [8]. The filtering procedure is detailed in Table 1, where E [ ] represents the expected value of □. In this scheme, the number N χ of sigma-points, terms Δ χ i , j , and weights ω j and ω j need to be set so as to achieve maximum accuracy; such issue is treated in Appendix A, where fundamental results are summarized.
In Kalman filtering it is usually assumed that all the state vector components are Gaussian variables; hence, the relevant probability distributions cannot be compliant with physical side constraints (like, e.g. positiveness of inter-laminar strength and toughness). In this work we do not aim to modify this filtering frame with the adoption of alternative probability distributions; instead, we define terms Δ χ i , j so that all the sigma-points be compliant with the aforementioned side constraints (see Appendix A). Through this procedure, which was proposed in [27], we can assure that all the sigma-points do not follow unphysical paths in the state vector space. It is worth noting that this procedure plays a role only when the sigma-points are well scattered around the mean x ^ i (see the prediction phase in Table 1), i.e. mainly at the beginning of filtering.

4. Failure assessment of laminates exposed to impacts

Figure 3. Impact on a 7-layer composite [13]: space-time diagram (the vertical dashed line here represents the possibly debonding surface of a brittle, homogeneous material).
Figure 3. Impact on a 7-layer composite [13]: space-time diagram (the vertical dashed line here represents the possibly debonding surface of a brittle, homogeneous material).
Algorithms 02 00808 g003
To assess the capability of the proposed FE-SPKF approach to provide accurate failure forecasts for layered composites subject to low velocity impacts and to simultaneously calibrate the inter-laminar constitutive laws, two experiments reported in the literature have been considered: an impact test on a 7-layer composite plate (experiment FY06001 in [13]), and an impact test on a 11-layer composite plate backed by another 11-layer composite plate (experiment 1 in [14]). In both cases, plane impacts caused the propagation of compressive waves behind the impact location; upon reflection at the rear surface of the plates, a spall-like failure governed by tensile tractions perpendicular to the inter-laminar surfaces took place. Because of the experimental set-up, shear tractions do not play any role in the failure events; in what follows, for ease of notation we therefore avoid using subscript n to denote (mode I) opening/closing displacement discontinuities along interfaces.
Figure 4. Impact on a 7-layer composite [13]: (a) comparison between experimental and tracked free surface velocity records; (b-c) estimated values of inter-laminar openings [ u ] j , j = 1 , . . . , 6 (b: PWL law; c: L-E law).
Figure 4. Impact on a 7-layer composite [13]: (a) comparison between experimental and tracked free surface velocity records; (b-c) estimated values of inter-laminar openings [ u ] j , j = 1 , . . . , 6 (b: PWL law; c: L-E law).
Algorithms 02 00808 g004
Figure 5. Impact on a 7-layer composite [13]: evolution in time of estimated inter-laminar (a) strength τ M and (b) toughness G.
Figure 5. Impact on a 7-layer composite [13]: evolution in time of estimated inter-laminar (a) strength τ M and (b) toughness G.
Algorithms 02 00808 g005
As far as space discretization is concerned, owing to the major role played by dilatational waves the full three-dimensional problem can be approached as a one-dimensional one featuring constrained deformation in the plane perpendicular to the wave propagation direction. Accordingly, one-dimensional meshes of constant strain elements have been adopted for the bulk; the characteristic size (length) of FEs in the through-the-thickness direction was determined by a trade-off between accuracy requirements and computational economy (leading to around 10 elements to discretize each lamina).
In the first experiment, each lamina was 1.37 mm in thickness, and was made of a balanced 5-harness satin weave E-glass and LY564 epoxy. The overall wave speed in the through-the-thickness direction was measured as 3.34 km/s, while the mass density was ϱ = 1885 kg/m 3 ; additional details and a thorough discussion of the results from an experimental perspective can be found in [13]. The specimen was stricken by a 12.5 mm thick aluminum impactor, flying at velocity v = 71 m/s, which led to dynamic delamination. The free surface velocity profile u ˙ r was measured during the test with a velocity interferometer for any reflector (VISAR).
Figure 6. Impact on a 7-layer composite [13]: evolution in time of estimated inter-laminar strength τ M and toughness G, and of the relevant covariances (a: PWL law; b: L-E law). In the plots the squares denote the common initial guess ϑ ^ 0 , whereas the circles denote the final estimates at time t N = 10 μs.
Figure 6. Impact on a 7-layer composite [13]: evolution in time of estimated inter-laminar strength τ M and toughness G, and of the relevant covariances (a: PWL law; b: L-E law). In the plots the squares denote the common initial guess ϑ ^ 0 , whereas the circles denote the final estimates at time t N = 10 μs.
Algorithms 02 00808 g006
The space-time diagram of this test is shown in Figure 3: in this diagram, wave reflection and dispersion due to the inter-laminar phases [45] have been neglected. Moreover, since failure is actually caused by delamination, the possibly debonding surface highlighted in the diagram is here to be considered as a reference one: even in the presence of wave dispersion, the actual failure of the laminate is expected to occur close to this reference surface.
Figure 4 shows a comparison between the estimated and the experimentally observed free surface velocity u ˙ r , and the time evolution of the tracked inter-laminar openings [ u ] j , j = 1 , . . . , 6 ; filter results are shown for both linear and exponential softening models. Before t 7 μs information on the mechanical properties of the failing inter-laminar surface are not filtered, since they are not yet available through measurables in y (see Figure 3), and some discrepancies between experimental and numerical graphs of u ˙ r are shown (see Figure 4(a)). Delamination is foreseen to occur along the third inter-laminar surface away from the impact location; in fact, only [ u ] 3 is diverging in Figure 4(b) and Figure 4(c), while all the other interface openings remain bounded. This latter outcome, which is independent of the shape of the softening envelope, well agrees with the state-space diagram of Figure 3 once wave dispersion caused by inter-laminar phases and other inner composite inhomogeneities is allowed for.
Previous results are accompanied by model calibration; in Figure 5 the time evolution of estimated inter-laminar strength τ M and toughness G are reported at varying initialization of ϑ within the domain:
C ϑ = 50 τ M , 0 200 (MPa) 0 . 5 G 0 2 . 0 (N/mm)
These estimates are not modified by the filter (apart from small fluctuations due to the stochastic environment) till the so-called pull back signal is processed at t 7 μs. Since [ u ] 3 is already diverging at that time (see Figure 4), the SPKF cannot obtain much information on the inter-laminar properties; hence, depending on the initialization in x ^ 0 , it sometimes happens that parameter estimates are eventually affected by biases. Anyway, it is worth noting that, under the same conditions, the standard EKF is not able to provide results featuring the same degree of accuracy [7]. In fact, the average converged estimates of τ M are in good agreement with the tensile strength of 119.5 MPa reported in [13], while final estimates of G well represent the behavior of this kind of composites. On the basis of the similar accuracy in the calibration of the softening envelopes here considered, it can be said that the shape of the softening regime does not play a major role in this kind of laminate failure.
To provide further insights into the performance of the FE-SPKF approach, Figure 6 shows an example of the time evolution of estimated model parameters along with the relevant statistical dispersion (here represented through error bars) around the expected values. Because of the aforementioned fast opening of the only decohering inter-laminar surface, entries in the covariance matrix P are not enhanced much by the SPKF.
Figure 7. Impact on a 11+11-layer composite [14]: space-time diagram (the vertical dashed line here represents the possibly debonding surface of a brittle, homogeneous material).
Figure 7. Impact on a 11+11-layer composite [14]: space-time diagram (the vertical dashed line here represents the possibly debonding surface of a brittle, homogeneous material).
Algorithms 02 00808 g007
In the second experiment a glass fiber reinforced plastic (GRP) specimen, 7 . 02 mm thick, was backed by a 6.91 mm thick GRP plate. The overall wave speed in the through-the-thickness direction of the plates amounted to 3.19 km/s, and the mass density was ϱ = 1867 kg/m 3 [14]. The specimen was stricken by a 5-layer GRP flyer, 2 . 96 mm in thickness and flying at velocity v = 85 m/s, which caused dynamic delamination. During the test the free surface velocity profile was again measured via a VISAR. The reference space-time diagram is depicted in Figure 7; it is here shown that, because of the test set-up, release waves interact causing delamination inside the back plate.
Results of the FE-SPKF approach are reported in Figure 8 in terms of comparison between experimental and tracked free surface velocity u ˙ r , and in terms of estimated displacement jumps along all the inter-laminar surfaces of the back plate (sequence starting at the specimen/back plate contact surface).
Figure 8. Impact on a 11+11-layer composite [14]: (a) comparison between experimental and tracked free surface velocity records; (b) estimated values of inter-laminar openings [ u ] j , j = 1 , . . . , 10 .
Figure 8. Impact on a 11+11-layer composite [14]: (a) comparison between experimental and tracked free surface velocity records; (b) estimated values of inter-laminar openings [ u ] j , j = 1 , . . . , 10 .
Algorithms 02 00808 g008
Tracking capabilities turn out again to be independent of the interface law: in Figure 8 discrepancies due to the shape of the softening regime cannot be detected. Such results have been obtained for any initialization values of τ M and G inside the domain:
C ϑ = 25 τ M , 0 100 (MPa) 0 . 1 G 0 0 . 6 (N/mm)
Delamination is foreseen to take place along the 7-th inter-laminar surface, in good agreement with the space-time diagram of Figure 7 and with the results of [7,14]. As for model calibration, outcomes are qualitatively similar to those reported for the previous test in terms of accuracy and convergence of the estimates.

5. Concluding remarks

In this paper we have proposed a finite element, sigma-point Kalman filter approach to the failure assessment of layered composites subject to low velocity impacts. Since under such loading conditions the damage and energy dissipation are mainly linked to the propagation of inter-laminar cracks (delamination), softening interface constitutive laws have been adopted to locally link the tractions transmitted between adjacent laminae to the opening/sliding displacement discontinuities occurring across the debonding surfaces.
To track in real time the whole state of the laminate and predict failure, a Kalman filter has been adopted. The filter has also allowed to calibrate the interface constitutive laws, whose parameters can be difficult to estimate due to the lower dimensionality of the inter-laminar surfaces. To properly deal with the nonlinearities induced by strength degradation phenomena, the sigma-point version of the Kalman filter has been implemented: this filter was already proved to be robust in highly nonlinear regimes, and to furnish more accurate outcomes than the customarily used extended Kalman filter (see [27,35]).
Through analysis of failures induced by plane waves propagating in the through-the-thickness direction of stricken laminates, it has been shown that debonding events are accurately tracked by the offered procedure. Furthermore, inter-laminar strength and toughness, which have been assumed unknown in the simulations, turned out to be accurately estimated.
In future works, to track and monitor diffused damage/delamination events, even featuring more complicated patterns than those here investigated, order reduction techniques will be coupled to the present methodology; this coupling will allow to get closer to the design of health monitoring systems for complex heterogeneous structures, able to detect in real time possible strength degradation phenomena induced by loading and environment.

A. Sigma-point transformation

In this Appendix we summarize the fundamentals of sigma-point transformation, with focus on its performance in representing the probability distribution of a random variable vector undergoing a nonlinear transformation [8,9].
Let x be a Gaussian random vector, with mean x ^ i and covariance P i at time t i . If x undergoes an analytic, nonlinear transformation f i within the time interval [ t i t i + 1 ] , its expected value at time t i + 1 reads:
x ^ i + 1 - = E [ x i + 1 - ] = f i ( x ^ i ) + n = 1 1 n ! E D Δ x n , i
where:
D Δ x n , i = 1 N x f i x | x = x ^ i ( x - x ^ i , ) n
N x being the dimension of x . The relevant error covariance matrix is:
P i + 1 - = E x i + 1 - x ^ i + 1 - x i + 1 - x ^ i + 1 - T = n = 1 m = 1 1 n ! 1 m ! E D Δ x n , i D Δ x m , i T - E D Δ x n , i E D Δ x m , i T
Through the sigma-point transformation, the probability distribution of x at time t i is sampled via a set of sigma-points χ ^ i , j scattered around x ^ i according to (see Table 1):
Δ χ i , 0 = 0 Δ χ i , k = + ψ P i 1 k k = 1 , . . . , N χ 2 Δ χ i , N χ 2 + k = - ψ P i 1 k
where: ψ is a scaling factor to be determined; P i is the square root of the (positive definite) covariance matrix P i ; 1 k is a unit vector pointing towards the k - th component of the state vector in the relevant space, whose role in (28) is to sample the k - th column of matrix P i ; N χ = 2 N x .
Information conveyed by the evolved sigma-points are merged at time t i + 1 via a weighted averaging scheme, according to:
x ^ i + 1 , SPT - = j = 0 N χ ω j χ ^ i + 1 , j - = j = 0 N χ ω j f i ( x ^ i ) + n = 1 1 n ! j = 0 N χ ω j D Δ χ i , j n , i
and:
P i + 1 , SPT - = j = 0 N χ ω j χ ^ i + 1 , j - - x ^ i + 1 , SPT - χ ^ i + 1 , j - - x ^ i + 1 , SPT - T = j = 0 N χ ω j n = 1 m = 1 1 n ! 1 m ! D Δ χ i , j n , i - r = 0 N χ ω r D Δ χ i , r n , i D Δ χ i , j m , i - s = 0 N χ ω s D Δ χ i , s m , i T
where:
D Δ χ i , j n , i = 1 N x f i x | x = x ^ i Δ χ i , j n
Weights ω j in (29) can be determined by matching the series expansions (25) and (29) up to third order; by assuming ω j = ω for j = 1 , . . . , N χ this requirement results in:
ω 0 + N χ ω = 1 2 ψ 2 ω = 1
To also set the value of the scaling factor ψ, a further condition needs to be established. Usually, it is suggested to match the diagonal entries of the fourth order terms in (25) and (29): this leads to ψ = 3 . In [27], partially exploiting the features of the so-called scaled sigma-point transformation [46], we proposed to adaptively set ψ so as to always fulfill the physical constraints:
ϑ m ϑ ϑ M
where ϑ m and ϑ M are, respectively, the minimum and maximum (if applicable) allowed values of model parameters. These side constraints need to be fulfilled by all the sigma-points χ ^ i , j : for j = 0 , conditions (33) are automatically satisfied; for j = 1 , . . . , N χ , conditions (33) are satisfied if:
ψ min ϑ ^ - ϑ m a k ϑ M - ϑ ^ a k , k = 1 , . . . , N χ 2 ; = 1 , . . . , N ϑ
where: a k = B ϑ P i 1 k ; B ϑ is a Boolean matrix, defined as ϑ = B ϑ x ; N ϑ is the dimension of vector ϑ.
Weights ω j in (30) are set by first letting ω j = ω j = ω for j = 1 , . . . , N χ . Value of ω 0 is then obtained by matching the fourth order terms involving D 2 , i D 2 , i T in (27) and (30), thus getting ω 0 = 4 - N χ ψ 2 - ψ 2 [27].

References and Notes

  1. Abrate, S. Impact on composite structures; Cambridge University Press: Cambridge, UK, 1998. [Google Scholar]
  2. Ladevéze, P.; Lubineau, G.; Marsal, D. Towards a bridge between the micro- and mesomechanics of delamination for laminated composites. Compos. Sci. Technol. 2006, 66, 698–712. [Google Scholar] [CrossRef]
  3. Corigliano, A.; Mariani, S. Simulation of damage in composites by means of interface models: parameter identification. Compos. Sci. Technol. 2001, 61, 2299–2315. [Google Scholar] [CrossRef]
  4. Kalman, R.E. A new approach to linear filtering and prediction problems. Trans. ASME J. Basic Eng. 1960, 82D, 35–45. [Google Scholar] [CrossRef]
  5. Bittanti, S.; Maier, G.; Nappi, A. Inverse problems in structural elastoplasticity: a Kalman filter approach. In Plasticity today. Modelling, Methods and Applications; Sawczuk, A., Bianchi, G., Eds.; Elsevier Applied Science Publishers: London, UK - New York, NY, USA, 1984; pp. 311–329. [Google Scholar]
  6. Gelb, A. Applied optimal estimation; MIT Press: Cambridge, USA, 1974. [Google Scholar]
  7. Mariani, S.; Corigliano, A. Impact induced composite delamination: state and parameter identification via joint and dual extended Kalman filters. Comput. Meth. Appl. Mech. Eng. 2005, 194, 5242–5272. [Google Scholar] [CrossRef]
  8. Julier, S.; Uhlmann, J.; Durrant-Whyte, H. A new method for the nonlinear transformation of means and covariances in filters and estimators. IEEE Trans. Automat. Control 2000, 45, 477–482. [Google Scholar] [CrossRef]
  9. Wan, E.; van der Merwe, R. The unscented Kalman filter. In Kalman filtering and neural networks; Haykin, S., Ed.; John Wiley & Sons, Inc.: New York, NY, USA, 2001; pp. 221–280. [Google Scholar]
  10. Lefebvre, T.; Bruyninckx, H.; Schutter, J.D. Kalman filters for nonlinear systems: a comparison of performance. Int. J. Contr. 2004, 77, 639–653. [Google Scholar] [CrossRef]
  11. Khalil, M.; Sarkar, A.; Adhikari, S. Nonlinear filters for chaotic oscillatory systems. Nonlinear Dynamics 2009, 55, 113–137. [Google Scholar] [CrossRef]
  12. Nguyen, H.-M.; Allix, O.; Feissel, P. A robust identification strategy for rate-dependent models in dynamics. Inverse Problems 2008, 24, 065006. [Google Scholar] [CrossRef]
  13. Yuan, F.; Tsai, L.; Prakash, V.; Rajendran, A.; Dandekar, D. Spall strength of glass fiber reinforced polymer composites. Int. J. Solids Struct. 2007, 44, 7731–7747. [Google Scholar] [CrossRef]
  14. Espinosa, H.; Dwivedi, S.; Lu, H.-C. Modeling impact induced delamination of woven fiber reinforced composites with contact/cohesive laws. Comput. Met. Appl. Mech. Eng. 2000, 183, 259–290. [Google Scholar] [CrossRef]
  15. Corigliano, A.; Mariani, S. Parameter identification in explicit structural dynamics: performance of the extended Kalman filter. Comput. Met. Appl. Mech. Eng. 2004, 193, 3807–3835. [Google Scholar] [CrossRef]
  16. Allix, O.; Ladevéze, P. Interlaminar interface modelling for the prediction of delamination. Compos. Struct. 1992, 22, 235–242. [Google Scholar] [CrossRef]
  17. Corigliano, A. Formulation, identification and use of interface elements in the numerical analysis of composite delamination. Int. J. Solids Struct. 1993, 30, 2779–2811. [Google Scholar] [CrossRef]
  18. Corigliano, A.; Allix, O. Some aspects of interlaminar degradation in composites. Comput. Met. Appl. Mech. Eng. 2000, 185, 203–224. [Google Scholar] [CrossRef]
  19. Alfano, G.; Crisfield, M. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. Int. J. Num. Met. Eng. 2001, 50, 1701–1736. [Google Scholar] [CrossRef]
  20. Corigliano, A.; Mariani, S.; Pandolfi, A. Numerical modeling of rate-dependent debonding processes in composites. Compos. Struct. 2003, 61, 39–50. [Google Scholar] [CrossRef]
  21. Corigliano, A.; Mariani, S.; Pandolfi, A. Numerical analysis of rate-dependent dynamic composite delamination. Compos. Sci. Technol. 2006, 66, 766–775. [Google Scholar] [CrossRef]
  22. Camacho, G.; Ortiz, M. Computational modelling of impact damage in brittle materials. Int. J. Solids Struct. 1996, 33, 2899–2938. [Google Scholar] [CrossRef]
  23. Mariani, S.; Perego, U. Extended finite element method for quasi-brittle fracture. Int. J. Num. Met. Eng. 2003, 58, 103–126. [Google Scholar] [CrossRef]
  24. Comi, C.; Mariani, S. Extended finite element simulation of quasi-brittle fracture in functionally graded materials. Comput. Met. Appl. Mech. Eng. 2007, 196, 4013–4026. [Google Scholar] [CrossRef]
  25. Needleman, A.; Rosakis, A. The effect of bond strength and loading rate on the conditions governing the attainment of intersonic crack growth along interfaces. J. Mech. Phys. Solids 1999, 47, 2411–2449. [Google Scholar] [CrossRef]
  26. Landis, C.; Pardoen, T.; Hutchinson, J. Crack velocity dependent toughness in rate dependent materials. Mech. Mater. 2000, 32, 663–678. [Google Scholar] [CrossRef]
  27. Mariani, S. Failure of layered composites subject to impacts: constitutive modeling and parameter identification issues. In Strength of Materials; Mendes, G., Lago, B., Eds.; Nova Publishers. (in press)
  28. Moës, N.; Dolbow, J.; Belytschko, T. A finite element method for crack growth without remeshing. Int. J. Num. Met. Eng. 1999, 46, 131–150. [Google Scholar]
  29. Wells, G.; Sluys, L. A new method for modelling cohesive cracks using finite elements. Int. J. Num. Met. Eng. 2001, 50, 2667–2682. [Google Scholar] [CrossRef]
  30. Comi, C.; Mariani, S.; Perego, U. An extended FE strategy for transition from continuum damage to mode I cohesive crack propagation. Int. J. Num. Anal. Met. Geomech. 2007, 31, 213–238. [Google Scholar] [CrossRef]
  31. de Borst, R.; Remmers, J.; Needleman, A. Mesh-independent discrete numerical representations of cohesive-zone models. Eng. Fract. Mech. 2006, 73, 160–177. [Google Scholar] [CrossRef]
  32. Melenk, J.; Babusˇka, I. The partition of unity finite element method: basic theory and applications. Comput. Met. Appl. Mech. Eng. 1996, 139, 289–314. [Google Scholar] [CrossRef]
  33. Mariani, S.; Martini, R.; Ghisi, A. A finite element flux-corrected transport method for wave propagation in heterogeneous solids. Algorithms 2009, 2, 1–18. [Google Scholar] [CrossRef] [Green Version]
  34. Hilber, H.; Hughes, T.; Taylor, R. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Eng. Struct. Dyn. 1977, 5, 283–292. [Google Scholar] [CrossRef]
  35. Mariani, S.; Ghisi, A. Unscented Kalman filtering for nonlinear structural dynamics. Nonlinear Dyn. 2007, 49, 131–150. [Google Scholar] [CrossRef]
  36. Hughes, T. The finite element method. Linear static and dynamic finite element analysis; Dover Publications: Mineola, NY, USA, 2000. [Google Scholar]
  37. Irwin, G.R. Structural aspects of brittle fracture. Appl. Mat. Res. 1964, 3, 65–81. [Google Scholar]
  38. Espinosa, H.; Zavattieri, P. A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. Part I: theory and numerical implementation. Mech. Mat. 2003, 35, 333–364. [Google Scholar] [CrossRef]
  39. Mariani, S.; Pandolfi, A.; Pavani, R. Coupled space-time multiscale simulations of dynamic delamination tests. Mat. Sci. 2005, 23, 509–519. [Google Scholar]
  40. Wan, E.; Nelson, A. Dual extended Kalman filter methods. In Kalman filtering and neural networks; Haykin, S., Ed.; John Wiley & Sons, Inc.: New York, NY, USA, 2001; pp. 123–173. [Google Scholar]
  41. Ljung, L. Asymptotic behavior of the extended Kalman filter as a parameter estimator for linear systems. IEEE Trans. Autom. Control 1979, AC-24, 36–50. [Google Scholar] [CrossRef]
  42. Reif, K.; Gu¨nther, S.; Yaz, E.; Unbehauen, R. Stochastic stability of the discrete-time extended Kalman filter. IEEE Trans. Autom. Control 1999, 44, 714–728. [Google Scholar] [CrossRef]
  43. Julier, S.; Uhlmann, J. A general method for approximating nonlinear transformations of probability distributions. Technical report, RRG. Department of Engineering Science, University of Oxford, 1996. [Google Scholar]
  44. Julier, S.; Uhlmann, J.; Durrant-Whyte, H. A new approach for filtering nonlinear systems. In Proceedings of the American Control Conference, Seattle, USA, June 1995; pp. 1628–1632.
  45. Zhuang, S.; Ravichandran, G.; Grady, D. An experimental investigation of shock wave propagation in periodically layered composites. J. Mech. Phys. Solids 2003, 51, 245–265. [Google Scholar] [CrossRef]
  46. Julier, S. The scaled unscented transformation. In Proceedings of the American Control Conference, Anchorage, USA, May 2002; Vol. 6, pp. 4555–4559.

Share and Cite

MDPI and ACS Style

Mariani, S. Failure Assessment of Layered Composites Subject to Impact Loadings: a Finite Element, Sigma-Point Kalman Filter Approach. Algorithms 2009, 2, 808-827. https://doi.org/10.3390/a2020808

AMA Style

Mariani S. Failure Assessment of Layered Composites Subject to Impact Loadings: a Finite Element, Sigma-Point Kalman Filter Approach. Algorithms. 2009; 2(2):808-827. https://doi.org/10.3390/a2020808

Chicago/Turabian Style

Mariani, Stefano. 2009. "Failure Assessment of Layered Composites Subject to Impact Loadings: a Finite Element, Sigma-Point Kalman Filter Approach" Algorithms 2, no. 2: 808-827. https://doi.org/10.3390/a2020808

Article Metrics

Back to TopTop