Numerical simulation of the dynamic response in pulse-loaded fibre

Abstract

This article presents a three-dimensional constitutive model to replicate the dynamic response of blast-loaded fibre–metal laminates made of 2024-0 aluminium alloy and woven composite (glass fibre–reinforced polypropylene). Simulation of the dynamic response is challenging when extreme localised loads are of concern and requires reliable material constitutive models as well as accurate modelling techniques. It is well known that back layers in a fibre–metal laminate provide structural support for front layers; thus, proper modelling of constituent failure and degradation is essential to understanding structural damage and failure. The improved developed model to analyse damage initiation, progression and failure of the composite is implemented in finite element code ABAQUS, and a good correlation is observed with experimental results for displacements of the back and front faces as presented by other researchers. The model was also able to predict accurately the tearing impulses. Finally, the concepts of the ‘efficiency of the charge’ and ‘effectiveness of the target’ are proposed in the context of localised blast loading on a structure. Dimensionless parameters are introduced to quantify these parameters.

Keywords

Fibre–metal laminate localised blast mesh objectivity damage mechanics composite

Introduction

Fibre–metal laminates (FMLs) are a class of hybrid materials that have attracted attention and interest due to their improved impact resistance and fatigue behaviour. A number of research works have been conducted comparing FMLs to monolithic metallic plates under dynamic loads (Karagiozova et al., 2010) which signify different aspects behaviour and requirements of modelling and experimentation on this class of structures. Glare, an FML patented by a joint venture of Aluminium Company of America (ALCOA) and Akzo Nobel N.V., has been used in commercial applications in Airbus 380 (Wang, 2009). It comprises thin aluminium 2024-T3 sheets and a unidirectional or a biaxial glass fibre–reinforced epoxy composite interleaved in a periodic structure through the thickness (Langdon et al., 2007a; Vlot, 1996; Vo et al., 2012, 2013).

Fleisher (1996) studied blast resistance of light luggage containers based on Glare and reported the capacity to withstand a bomb blast greater than that of the Lockerbie air disaster. The FML container tested with a Lockerbie-type explosion shows the enormous potential of glare in the design of Blast-Resistant Aircraft Baggage Containers for in-service aircraft use (Fleisher, 1196; Vo et al., 2012, 2013).

More recently, comprehensive investigations on the blast resistance of FMLs have been carried out at the Impact Research Centre of the University of Liverpool and the Blast, Impact & Survivability Research Unit (BISRU) of the University of Cape Town. In a set of the experiments (Langdon et al., 2007a, 2007b; Lemanski et al., 2006, 2007), FMLs composed of 2024-0 aluminium alloy and glass fibre–reinforced polypropylene (GFPP) subjected to localised blast loading were investigated. The test sample consisted of an FML based on 12 different configurations, having total laminate thickness of 2–15 mm. Failure modes in this work were characterised as Mode I failure: large inelastic deformation of the back face of the panel (furthest away from the blast), Mode II failure: complete tearing of the back face and Mode II*: transition between the two failure modes – this occurs at threshold impulses. It was observed that thicker panels exhibit smaller displacements for a given impulse than thinner panels, and that the Mode I non-dimensional displacement of both the front and back faces are linearly related to the non-dimensional impulse. It was also shown that thinner panels exhibit behaviour closer to that of a monolithic plate. As the panel thickness increases, behaviour deviates from the monolithic plate response, that is, debonding and delamination become more pronounced. Further numerical studies for Mode I failure show that when expressed in terms of non-dimensional parameters, the difference between the mean displacements of the front and back faces falls within one plate thickness, following a linear trend, and the threshold impulse for the onset of tearing for Mode II failure was shown to be linearly dependent on panel thickness (Langdon et al., 2007a; Vo et al., 2012, 2013).

Karagiozova et al. (2010) studied the response of FMLs to localised blast loading in order to interpret the deformation mechanisms due to highly localised pressure pulses causing permanent deformations and damage observed. They went on to scrutinise the influence of the loading and material parameters on the final deformation characteristics examined. In addition, a comparison between the responses of relatively thin FML panels, monolithic aluminium alloy plates of equivalent mass and a foam-core panel to localised blast were also presented in their work (Karagiozova et al., 2010). Modelling parameters of GFPP and adhesives used in this work were derived from the work of Karagovizova et al. (2010a).

As observed by Vo et al. (2012, 2013), experimental modelling and subsequent optimisation are quite expensive and there is need for predictive capabilities through finite element (FE) models able to depict the behaviour of FMLs a priori rather than as post-test corroboration tools. As virtual testing proves that valuable attempts in this line have been made to model FMLs by, for example, Vo et al. (2012, 2013). However, the focus of this study is to introduce additional parameters that would increase the confidence of the analyst in modelling FMLs as well as proposing a robust constitutive model to be used in blast loading scenarios usually characterised and accompanied by occurrence of damage and fracture.

Composite laminates composed of fibre-reinforced plies are integral parts of FMLs. Thus, a thorough modelling of composites is an integral part of modelling FMLs. Failure and damage in laminated structures (composites) can be studied using micro-mechanical approach that considers failure and damage at the consistent level, mini-mechanical approach which focuses on the representative volumetric element (RVE) and statistical homogeneity and meso-mechanical approach which looks primarily into delamination or a fully macro-mechanical approach. A Continuum Damage Mechanics (CDM) approach in which material properties of the composite are homogenised and failure and damage are studied at the ply/lamina level (Batra et al., 2012) through several damage parameters is adopted in this work. It is important to note that damage studied at constituent level is both computationally expensive for practical structures and experimentally immoderate, if at all possible. Furthermore, extensive characterisation is required to determine material parameters in the damage model.

One common approach for studying composites is to use a micromechanics approach to deduce effective properties of a ply and CDM to study failure and damage at laminar level. CDM theories capture effects of microscopic damage using the thermodynamics-based theory of internal variables (Coleman and Gurtin, 1967). Ladeveze and LeDantec (1992) used this approach to degrade elastic properties of the composite due to fibre breakage and matrix cracking and plasticity theory to account for permanent deformations induced under shear loading. Hassan and Batra (2008) used three internal variables to characterise damage due to fibre breakage, matrix cracking and fibre/matrix debonding. In their work, delamination between adjacent plies was analysed using a failure surface quadratic (and thus convex) in transverse normal and the transverse shear stresses. Puck and Schürmann (1998) generalised Hashin’s (1980) stress-based failure criteria and proposed a technique to degrade elastic parameters of the lamina subsequent to the initiation of a failure mode. Xiao et al. (2007) validated experimental quasi-static punch-shear test results carried out on plain weave S-2 glass/SC-15 epoxy composite laminates. The numerical modelling was carried out using a material model called MAT162, which was incorporated into LS-DYNA. MAT162 uses damage mechanics principles for progressive damage and material degradation. In their work, the calibrated material properties have been used to simulate plates of other thicknesses, and the simulated results show good agreement with experimental results. Among the earlier works on composite degradation is the work of Matzenmiller et al. (1995). They proposed that when one of the Hashin failure criteria is satisfied at a point in the composite structure, damage ensues at that point and it is depicted by introducing damage variables for fibre breakage in tension and compression, matrix cracking in tension and compression and in-plane shear. It is shown that the evolution of these internal variables depends upon the values of stresses in Hashin’s failure criteria which are expressed in terms of stress invariants for a transversely isotropic body and strength parameters for the composite (Batra et al., 2012).

This work follows the philosophy of Matzenmiller as adopted also by Batra et al. (2012). In Batra and his team’s work, Matzenmiller damage evolution criteria are used to study three-dimensional (3D) deformations of a 16-ply laminate subjected to normal incidence impact by a slow-moving rigid sphere, and the effective elastic–plastic properties of the composite were derived using Paley and Aboudi’s (1992) method of cells but with the continuity of shear tractions across cell boundaries relaxed. A user-defined subroutine (VUMAT) was developed which updates stresses as its main duty. It has to deal with material degradation to enable accurate stress calculations. The subroutine takes as input from the commercial FE software ABAQUS values of six strain components at an integration point of a generic FE, then computes based on intact constitutive coefficients stresses in each ply using the constituent level properties, computes effective stresses at integration points, checks for Hashin’s failure criteria, computes damage if necessary, modifies stresses due to the induced damage and supplies them back to ABAQUS. As the procedure is explicit, there is no need for iterations to update the constitutive model at every increment and the base state for each increment is that of the previous increment. The constitutive tensor is only updated subsequent to calculation of stresses at the end of each increment.

Donadon et al. (2008) used a smear crack approach to develop a failure model for predicting damage in 3D deformations of a composite structure.

Iannucci and colleagues (2008, 2006) have worked extensively on the progressive failure of woven composites. Iannucci (2006) proposed a simple damage mechanics-based progressive failure model for thin woven carbon composites under impact loads (i.e. ballistic damage, bird strike and fragmentation attacked). This model can be appreciated when one considers the difficulty in characterising material constants for existing stress-based scalar functions-based techniques of modelling composite failure. The proposed approach is based on an unconventional thermodynamic maximum energy dissipation approach, which entails controlling damage evolution and hence energy dissipation per second, rather than damage. The method was implemented into explicit dynamic FE code, DYNA3D.

An energy-based damage mechanics approach was also postulated by Iannucci and Willows (2006; 2007) for woven composite materials. Five damage variables were introduced for in-plane damage per ply. The damage variables were directly related to stiffness degradation within the composite lamina and ultimately within the laminate. The evolution of damage was controlled by a series of damage–strain equations. This allows the total energy dissipated for each damage mode to be set as a material parameter. This approach was implemented into both the LLNL and LS versions of DYNA3D for plane stress (shell) elements.

The principle of smeared cracking approach highlighted by Donadon et al. (2007a, 2007b, 2008) and Bažant and Oh (1983) was used in conjunction with some assumptions to ensure its validity in this work. The fundamental principle of smeared crack approach, used in CDM, assumes that the criteria for crack propagation and the prediction of the direction of propagation are based on failure criteria expressed in terms of stresses and strains and the stress intensity factors. In this approach, the cracked material is assumed to remain a continuum, and the mechanical properties (stiffness and strength) are modified to account for the effect of cracking, according to the evolving states of stresses and strains. Thus, re-meshing is in effect not needed (Cervera and Chiumenti, 2006), and volumetric model crack density is linked to a global modal damage parameter at the element level.

Although it is possible to model composites without course for re-meshing using the extended finite element method (XFEM), this method is beyond the scope of this study. XFEM is a fracture mechanics-based approach which allows for crack propagation without re-meshing at the expense of tracking the crack front through the FE mesh and progressively enriching the displacement field with functions that represent displacement jumps across the crack and the developed singular fields at the tip of an advancing crack (Belytschko and Black, 1999; Moës and Belytschko, 2002; Sukumar et al., 2000).

Donadon et al. (2008b) and Bažant and Oh (1983) related volumetric energy, defined as area underneath the stress–strain curve to the fracture energy of the composite material (which depends upon modal stress–displacement curves) in question. This method assumes a strain-softening constitutive law for modelling the gradual stiffness reduction due to micro-cracking process within the cohesive or process zone of the material and translates this damage width to the FE domain. Thus, the fracture energy, that is, energy dissipated per unit area within a fully failed element can be written in terms of the specific energy by multiplying the specific energy by a geometric quantity defined as characteristic length; for example, for a unidirectional stress state, it should have the direction aligned with the loading direction and at most times, equals the side of the element.

This article presents an improved numerical modelling for predicting the response of blast-loaded FMLs. A simplified model is presented to analyse the damage initiation, damage progression and failure of the 3D solid woven GFPP material within the FMLs. The response of the GFPP significantly influences the overall response of the FMLs. The simplified model for the GFPP incorporates strain rate effects and an objectivity algorithm for strain softening to control energy dissipation associated with each failure mode regardless of mesh refinement and topology. This is done within the framework of CDM and using thermodynamically consistent damage variables. Cracking is treated as an irreversible thermodynamic process, and mesh objectivity and inclusion of strain rate effects in 3D are the improvements upon previous works in the literature and encompass the novel aspects of this work. To ensure the appropriateness of the model proposed in this work, the following assumptions are made: (a) the crack growth direction in the modelled composite is parallel to the edge of the FEs and (b) the composite meshes are structured. Experimental and numerical results are compared and discussed, and in all cases, good correlation is observed.

Geometry and modelling

Various 400 mm × 400 mm panels of FMLs were manufactured and tested under localised blast at the University of Cape Town from 0.025 in thick sheets of 2024-0 aluminium alloy and a woven glass–glass fibre/propylene composite. The FML panels tested had exposed areas of 300 mm × 300 mm unanimously and were labelled AXTYZ-#, where A = aluminium, X = number of aluminium layers, T = GFPP, Y = number of blocks of GFPP, Z = number of plies of GFPP per block and # indicates the panel number (Karagiozova et al., 2010; Langdon et al., 2007b; Nurick et al., 2009).

Due to the symmetrical in-plane architecture of the panel, only one-quarter of the panel was modelled in ABAQUS 6.12 with appropriate boundary conditions, that is, symmetry and fully clamped on internal and external edges, respectively. The FML is modelled as a four-part structure, that is, comprising aluminium alloy, composite, interface adhesive layers (cohesive layers representing the adhesive between composite and aluminium alloy) and laminate adhesive layers (cohesive layers representing the adhesive layers between plies of composite). In this work, we introduce a laminate adhesive layer between every two consecutive plies to reduce computational costs while retaining accuracy. Figure 1 illustrates the through-thickness structure of the FML. The aluminium part is meshed using linear brick elements, that is, C3D8R elements, which are eight-noded, linear hexahedral elements with reduced integration formulation and hourglass control. The adhesive layers in the composite and interface between composite and aluminium are modelled with 3D cohesive element with direct traction-separation formulation (COH3D8).

Figure 1.

Schematic representation of through-thickness architecture in an FML plate.

Material modelling

Composite modelling

CDM has been adopted in this work which assumes a linear elastic orthotropic response up to the point of damage initiation. Modified Hashin’s damage initiation criteria must be satisfied for damage to initiate as this model is adopted. These damage criteria are similar to the ones proposed by Xiao et al. (2007) for multi-axially loaded composites. The initial elastic constants of the undamaged material are the elastic moduli E₁, E₂ and E₃; shear moduli G₁₂, G₂₃ and G₃₁ and Poisson’s coefficients υ₁₂, υ₂₃ and υ₃₁, where 1, 2 and 3 denote local axes of the material in the direct in-plane fill (warp), transverse in-plane (weft) and out-of-plane directions, respectively.

The constitutive relation for undamaged composite ply is therefore

σ = D_{0} ε

(1)

or by defining

C_{0} = D_{0}^{- 1}

(2)

ε = C_{0} σ

(3)

where σ is the stress tensor with respect to the material principal axes, $D_{0}$ is the constitutive matrix, and the compliance matrix $C_{0} = D_{0}^{- 1}$ is expressed as follows in this system of coordinates

C 0 = [\begin{matrix} \frac{1}{E_{1}} & \frac{- υ_{21}}{E_{2}} & \frac{- υ_{31}}{E_{3}} & 0 & 0 & 0 \\ \frac{- υ_{12}}{E_{1}} & \frac{1}{E_{2}} & \frac{- υ_{32}}{E_{3}} & 0 & 0 & 0 \\ \frac{- υ_{13}}{E_{1}} & \frac{- υ_{23}}{E_{2}} & \frac{1}{E_{3}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{12}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{23}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G_{31}} \end{matrix}]

(4)

The components of the stress tensors σ and ϵ are, respectively

σ = {(σ_{11}, σ_{22}, σ_{33}, σ_{12}, σ_{23}, σ_{31})}^{T}, ε = {(ε_{11}, ε_{22}, ε_{33}, ε_{12}, ε_{23}, ε_{31})}^{T}

(5)

It is obvious that any rotation of the axes would result in a corresponding change in constitutive matrix followed by affine transformation of the coordinates, that is, a transformation which is both covariant and contravariant.

Equation (3) is further simplified in the relationship shown in equation (6)

\begin{array}{l} σ_{11} = \frac{(- 1 + υ_{23} υ_{32}) E_{1} ε_{11}}{T} + \frac{(υ_{21} + υ_{23} υ_{31}) E_{2} ε_{22}}{T} - \frac{(υ_{21} + υ_{23} υ_{31}) E_{2} ε_{22}}{T} \\ σ_{22} = \frac{(υ_{12} + υ_{13} υ_{32}) E_{2} ε_{22}}{T} + \frac{(- 1 + υ_{13} υ_{31}) E_{2} ε_{22}}{T} - \frac{(υ_{32} + υ_{12} υ_{31}) E_{2} ε_{22}}{T} \\ σ_{33} = \frac{(υ_{13} + υ_{23} υ_{12}) E_{3} ε_{33}}{T} - \frac{(υ_{23} + υ_{13} υ_{21}) E_{3} ε_{33}}{T} + \frac{(- 1 + υ_{12} υ_{21}) E_{3} ε_{33}}{T} \\ σ_{12} = 2 G_{12} ε_{12} \\ σ_{23} = 2 G_{23} ε_{23} \\ σ_{31} = 2 G_{31} ε_{31} \end{array}

where

T = - 1 + υ_{12} υ_{21} + υ_{13} υ_{31} + υ_{23} υ_{32} + υ_{23} υ_{12} υ_{31} + υ_{13} υ_{21} υ_{32}

(6)

Composite material properties

The mechanical properties of composite shown in Table 1 are measurable using standard test procedures (see Karagiozova et al., 2010; Nurick et al., 2009; Paley and Aboudi, 1992). These encompass simple tensile and compressive tests as well as in-plane shear tests and variety of their combinations when convex failure hypersurfaces are relevant. However, values for fracture energies of composites are not readily available in most cases. Donodan et al. (2007b) have highlighted various procedures for measuring modal fracture energies. The directional strengths of the tested composites (GFPP) were obtained from the work of Karagiozova et al. (2010). In order to determine the fracture energy associated with each fracture mode, tests can be undertaken using three different pre-cracked geometries, that is, Overhead Compact Tension (OCT) (tension/compression), double edge notch (DEN) (tension) and four-point-bending specimens. Procedures and techniques employed for the determination of these parameters were highlighted by Donadon et al. (2007b). Due to lack of availability of values for intra-laminar fracture energies in open literature, values presented in Table 2 are calibrated to validate the experimental results for the test samples.

Table 1.

Mechanical properties of GFPP (Karagiozova et al., 2010; Vo et al., 2012, 2013).

E₁ (GPa)	E₂ (GPa)	E₃ (GPa)	v ₁₂	v ₂₃	v ₃₁	G₁₂ (GPa)	G₂₃ (GPa)	G₃₁ (GPa)
13	13	4.8	0.1	0.3	0.3	1.72	1.72	1.69

GFPP: glass fibre–reinforced polypropylene.

Table 2.

Intra-laminar fracture energies.

$G_{f 11}^{t}$ (intra-laminar fracture in tension in direction 1, warp)	30 kJ/m²
$G_{f 22}^{c}$ (intra-laminar fracture in tension in weft direction 2, weft)	30 kJ/m²
$G_{f 11}^{t}$ (intra-laminar fracture in compression in direction 1, warp)	15 kJ/m²
$G_{f 22}^{c}$ (intra-laminar fracture in compression in direction 2, weft)	15 kJ/m²
$G_{f 33}^{c}$ (intra-laminar fracture through-thickness direction)	10 kJ/m²
$G_{f 12}^{s}$ (in-plane shear intra-laminar fracture)	9 kJ/m²
$G_{f 23}^{s}$ (out-plane shear intra-laminar fracture)	9 kJ/m²
$G_{f 31}^{s}$ (out-plane shear intra-laminar fracture)	9 kJ/m²

Elastic damage energy and dissipation

Using a strain-equivalent damage mechanics formulation, a damage elastic compliance matrix S is assumed with the general form

S = [\begin{matrix} \frac{1}{(1 - d_{1}) E_{1}} & - v_{21} / E_{2} & - v_{31} / E_{3} & 0 & 0 & 0 \\ - v_{12} / E_{1} & \frac{1}{(1 - d_{2}) E_{2}} & - v_{32} / E_{3} & 0 & 0 & 0 \\ - v_{13} / E_{1} & - v_{23} / E_{2} & \frac{1}{(1 - d_{3}) E_{3}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{(1 - d_{4}) G_{12}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{(1 - d_{5}) G_{23}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{(1 - d_{6}) G_{31}} \end{matrix}]

(7)

The above compliance matrix has six scalar damage parameters d₁, d₂, d₃, d₄, d₅ and d₆. Intrinsic to this formulation is the assumption that symmetry of the compliance matrix is preserved throughout the analysis and is never violated. These six damage parameters have values between 0 and 1 corresponding to intact (virgin) and fully damaged (fractured) material, respectively. Values of d₁, d₂ and d₃ are associated with damage in the warp, weft and through-thickness directions, respectively, while d₄, d₅ and d₆ correspond to the in-plane and out-of-plane shear failures. Parameters d₁–d₆ can be viewed as invariant functions that represent physical parameters. Note there are no additional damage parameters to model independent degradation in Poisson’s ratios. The general damage mechanics formulation is based on an internal energy function ϕ for an orthotropic solid. The function ϕ is a scalar field function shown in the relationship in equation (8)

φ = \frac{1}{2} σ^{T} S σ = \frac{1}{2} S_{i j k l} σ_{i j} σ_{k l}

(8)

We introduce thermodynamic forces ( Y = Y₁, Y₂, Y₃, Y₄, Y₅ and Y₆) which act as driving forces for damage development. We show that the strain tensor ϵ and the thermodynamic forces Y can be derived from the internal energy function as

ε = \frac{\partial φ}{\partial σ} = S σ (ε_{i j} = \frac{\partial φ}{\partial σ_{i j}} = S_{i j k l} σ_{k l}) and Y = \frac{\partial φ}{\partial d} (Y_{i} = \frac{\partial φ}{\partial d_{i}})

(9)

which implies that for elastic material

\begin{array}{l} Y_{1} = \frac{σ_{11}^{2}}{2 {(1 - d_{1})}^{2} E_{1}} Y_{2} = \frac{σ_{22}^{2}}{2 {(1 - d_{2})}^{2} E_{2}} \\ Y_{3} = \frac{σ_{33}^{2}}{2 {(1 - d_{3})}^{2} E_{3}} Y_{4} = \frac{σ_{12}^{2}}{2 {(1 - d_{4})}^{2} G_{12}} \\ Y_{5} = \frac{σ_{23}^{2}}{2 {(1 - d_{5})}^{2} G_{23}} Y_{6} = \frac{σ_{31}^{2}}{2 {(1 - d_{6})}^{2} G_{31}} \end{array}

(10)

The variable Y_i is referred to as damage energy rates in particular failure modes. It can be inferred, as we will see later that

d_{i} = f_{i} (Y_{1}, Y_{2}, Y_{3}, Y_{4}, Y_{5} and Y_{6})

(11a)

where f_i (f₁, f₂, f₃, f₄, f₅ and f₆) is to be determined from a multi-axial failure or interaction between damage states. Further model assumptions are as follows. (a) The six damage modes in all stress directions are decoupled and determined by Y₁, Y₂, Y₃, Y₄, Y₅ and Y₆. This assumption was made by A.F. Johnson (2001) in modelling two-dimensional (2D) plane elements for reinforced composites. (b) Ply material is non-healing, thus on unloading, subsequent to being damaged, the damage parameter remains constant until a larger damaging load is applied; in other words, each damage parameter is an absolutely ascending function of time irrespective of loading history and multi-axial pattern of damage. Thus, the evolution function depends on the maximum values of Y_i attained. (c) The inelastic behaviour of the matrix is ignored and was found in this particular case not to affect the correlation between the experiments and the numerical models.

In order to implement assumption (b), we introduce a parameter ${\bar{Y}}_{i}$ related to maximum value of damage forces reached during the previous loading history, that is

{\bar{Y}}_{i} (t) = \max (Y_{i} (τ)) where τ ⩽ t

(11b)

Thus, imagining an elastic domain without damage at the outset of loading and subsequent evolution to a lower than unity value of the modal damage initiation parameter and full separation or cut-off at an upper damage threshold, we deduce that

\begin{matrix} d_{1} = 0 for {\bar{Y}}_{1} < Y_{10}; d_{1} = f ({\bar{Y}}_{1} - Y_{10}) for {\bar{Y}}_{10} < Y_{1} < Y_{1 f} or d_{1} = 1 \\ d_{2} = 0 for {\bar{Y}}_{2} < Y_{20}; d_{2} = f ({\bar{Y}}_{2} - Y_{20}) for {\bar{Y}}_{2} < Y_{20} < Y_{2 f} or d_{2} = 1 \\ d_{3} = 0 for {\bar{Y}}_{3} < Y_{30}; d_{3} = f ({\bar{Y}}_{3} - Y_{30}) for {\bar{Y}}_{3} < Y_{30} < Y_{3 f} or d_{3} = 1 \\ d_{4} = 0 for {\bar{Y}}_{4} < Y_{40}; d_{4} = f ({\bar{Y}}_{4} - Y_{40}) for {\bar{Y}}_{4} < Y_{40} < Y_{4 f} or d_{4} = 1 \\ d_{5} = 0 for {\bar{Y}}_{5} < Y_{50}; d_{5} = f ({\bar{Y}}_{5} - Y_{50}) for {\bar{Y}}_{5} < Y_{50} < Y_{5 f} or d_{5} = 1 \\ d_{6} = 0 for {\bar{Y}}_{6} < Y_{60}; d_{5} = f ({\bar{Y}}_{6} - Y_{60}) for {\bar{Y}}_{6} < Y_{60} < Y_{6 f} or d_{6} = 1 \end{matrix}

(12)

In this model, we assume nonlinear (exponential) forms for function f, that is, for d₁, d₂, d₃, d₄, d₅ and d₆. Thus, the evolution equations of a 3D composite element requires the relation between d_i and $({\bar{Y}}_{i} - Y_{i 0}),$ where i = 1, 2, 3, 4, 5 and 6 as shown in the equation above. In the above equation, the threshold parameters, defined as Y₁₀, Y₂₀, Y₃₀, Y₄₀, Y₅₀, Y₆₀, Y_1f, Y_2f, Y_3f, Y_4f, Y_5f and Y_6f, determine the bounds for the damage parameters.

Damage initiation criteria

As mentioned earlier, the composite is assumed to behave elastically in an orthotropic manner until the onset of damage. The modified Hashin criteria take into consideration more modes of failure than the traditional one and are applicable to woven as well as unidirectional fibrous composites. This is shown in equation (13a) to (13h) which determine the point of damage initiation. Damage initiates in any of the failure modes when r_m = 1.

Fill and warp fibre tensile/shear failure mode

The failure surface that characterises this mode is given by the quadratic interaction between the associated axial and through-thickness shear strains; that is, for fill direction (direction 1)

{(f_{1})}^{2} = {[\frac{E_{1} 〈 ε_{1} 〉}{S_{1 T}}]}^{2} + {[\frac{G_{31} ε_{31}}{S_{1 F S}}]}^{2}

(13a)

For warp direction (direction 2)

{(f_{2})}^{2} = {[\frac{E_{2} 〈 ε_{2} 〉}{S_{2 T}}]}^{2} + {[\frac{G_{23} ε_{23}}{S_{2 F S}}]}^{2}

(13b)

where the sign <x> is equal to x for positive values of x and is 0 otherwise. For the woven composite plies, 1, 2 and 3 represent in-plane fill (direction 1), in-plane warp (direction 2), and out-of-plane (direction 3) directions, respectively. E and G are tensile and shear moduli, respectively. S_1T and S_2T are tensile strengths in the fill and warp directions. S_1FS and S_2FS are fibre shear strengths in 1–3 and 2–3 directions, ϵ₁ and ϵ₂ are failure tensile strains in a and b directions, ϵ₃₁ and ϵ₂₃ are shear strains in 1–3 and 2–3 planes. We assume that S_1FS = S_FS and S_2FS = S_FS*S_2T/S_1T, where S_FS is the fibre shear strength in Table 3.

Table 3.

Strength of composite (GFPP) (Vo et al., 2013).

S_1T (MPa)	S_2T (MPa)	S_1C (MPa)	S_2C (MPa)	S_3T (MPa)	S_FC (MPa)	S_FS (MPa)
300	300	200	200	300	200	140
S₁₂ (MPa)	S₂₃ (MPa)	S₃₁ (MPa)
140	140	140

Table 4.

Mechanical properties of aluminium (Al 2024-0) (Vo et al., 2013).

ρ (kg/m³)	E (GPa)	A (MPa)	B (MPa)	n	C	ϵ₀ (s⁻¹)
2700	73.4	85	325	0.4	0.001	0.0083
	UTS (MPa)	d ₁	d ₂	d ₃	d ₄	d ₅
	186	0.13	0.13	1.5	0.011	0

Fibre compressive failure mode

The in-plane compressive damage in directions 1 and 2 are given by the following failure criterion

For fill direction (direction 1)

{(f_{3})}^{2} = {[\frac{E_{1} {ε^{'}}_{1}}{S_{1 C}}]}^{2} where {ε^{'}}_{1} = - ε_{1} - 〈 ε_{3} 〉 \frac{E_{3}}{E_{1}}

(13c)

For warp direction (direction 2)

{(f_{4})}^{2} = {[\frac{E_{2} {ε^{'}}_{2}}{S_{2 C}}]}^{2} where {ε^{'}}_{2} = - ε_{2} - 〈 ε_{3} 〉 \frac{E_{3}}{E_{2}}

(13d)

where S_1C and S_2C are in-plane compressive strengths in directions 1 and 2, respectively.

Fibre crush failure mode

The crush damage because of high through-thickness compressive pressure from blast waves can be modelled with the hypersurface of equation (13e)

{(f_{5})}^{2} = {[\frac{E_{3} ε_{3}}{S_{F C}}]}^{2}

(13e)

where S_FC is the fibre crush strength.

Fibre in-plane shear failure mode (planes 1 and 2)

A woven layer can damage under in-plane shear stress without occurrence of fibre breakage. The in-plane matrix damage mode is given by

{(f_{6})}^{2} = {[\frac{G_{12} ε_{12}}{S_{12}}]}^{2}

(13f)

where S₁₂ is the layer shear strength due to matrix shear failure.

Fibre out-of-plane shear failure mode (planes 2 and 3)

The in-plane matrix damage mode is given by

{(f_{7})}^{2} = {[\frac{G_{b c} ε_{b c}}{S_{23}}]}^{2}

(13g)

where S₂₃ is the corresponding layer shear strength due to matrix shear failure.

Fibre in-plane shear failure mode (planes 1–3)

The in-plane matrix damage mode is given by

{(f_{8})}^{2} = {[\frac{G_{31} ε_{31}}{S_{31}}]}^{2}

(13h)

where S₃₁ is the corresponding layer shear strength due to matrix shear failure.

Strain rate effect

The effect of strain rate on the strength values of composite failure modes is modelled by multiplying the associated strength values {S₀} by a scale factor

{S_{R T}} = {S_{0}} [1 + C_{r a t e} \ln (\frac{{\dot{\bar{ε}}}}{{\dot{ε}}_{0}})]

(14)

where C_rate is the dimensionless strain rate constant, and {S₀} is the available strength value of {S_RT} at the reference strain rate ${\dot{ε}}_{0} = 1 s^{- 1}$ . The values of the strain rate used for the different strength values of the GFPP in this work were obtained from Vo et al. (2012, 2013; Al-Ani and Habibi, 2013). In this work C_rate = 3.5.

It has been observed that the through-thickness elastic modulus of GFPP used in this work is strain rate dependent. A through-thickness elastic modulus of 4.8 GPa used was extracted from an experiment under a strain rate of 10³s⁻¹, while the static through-thickness elastic modulus was found to be half of the measured value under such strain rate. However, E₁ and E₂ for GFPP were observed to be strain rate independent (Karagiozova et al., 2010).

Damage progression criteria

The evolution function for the damage proposed in the relationship of equation (12) is empirically determined following a procedure similar to the one adopted by Johnson (2001) and Matzenmiller et al. (1995) (Batra et al., 2012). The damage evolution at a material point is defined in terms of an internal variable Q_m (m = 1, 2, 3, 4, 5, 6, 7, 8) associated with failure index f_m by the empirical relationship in equation (13). The threshold values of the failure indexes, f_m, (represented as f₀) in this context are analogous to the terms Y₁₀, Y₂₀, Y₃₀, Y₄₀, Y₅₀, Y₆₀. The eight failure modes showed by the modified Hashin criteria in equation (13) determine the value of the failure index f₀.The value H_m measures the rate of damage evolution

Q_{m} = 1 - \frac{f_{o}}{f_{m}} \exp {- 2 H_{m} (\frac{f_{m} - f_{0}}{f_{o}})} (m = 1, 2, 3, 4, 5, 6) r_{m} ⩾ 1

(15)

Equation (16) shows how the various failure modes (internal variable, Q_m) affects the properties of the composite in all directions (i.e. the values of d_i)

\begin{array}{l} d_{1} = \max {Q_{1}, Q_{3}, Q_{5}} \\ d_{2} = \max {Q_{2}, Q_{4}, Q_{5}} \\ d_{3} = Q_{5} \\ d_{4} = \max {Q_{1}, Q_{2}, Q_{3}, Q_{4}, Q_{5}, Q_{6}} \\ d_{5} = \max {Q_{3}, Q_{4}, Q_{5}, Q_{7}} \\ d_{6} = \max {Q_{1}, Q_{2}, Q_{5}, Q_{8}} \end{array}

(16)

In the relationship above, we see that more than one failure index, Q_m, affects the various individual values of failure parameter, d_i. Fibre damage in either the fill or warp directions results in the reduction in stiffness in the loading direction and in the related shear direction. For failure in shear direction, Q₅ influences all six components of the damage vector, d_i. Notice that Q₅ is related to fibre crushing mode, that is, as a resultant of high through-thickness compressive stress as shown in equation (13e).

In this model, we degrade elastic properties as damage evolves, and the damaged and undamaged constitutive equations are checked and updated in the model per increment. The irreversibility of the internal variables is accounted for by requiring that dQ_m = 0 whenever dr_m ⩽ 0, where dr_m represents an increment in r_m for an increment in the applied load.

Fracture energy (energy release rate) and mesh objectivity

Equation (16) proposed above defines the softening evolution in composites after the damage is initiated. The finite area under the, that is, individual stress–displacement curve (i.e. each representing an individual damage mode) represents the energy dissipated for each mode. For validity of this proposed model, we assume that the crack growth direction is parallel to one edge of the FE, and the composite meshes are structured.

The specific energy dissipated per unit volume is

\int_{t = 0}^{t = \infty} \dot{D} d t = \frac{G_{m}}{l^{*}} = \int_{t = 0}^{t = \infty} Y \dot{Q} d t = \frac{1}{2 E_{m}} \int_{r = r_{0}}^{r = \infty} S_{m}^{2} r^{2} Q^{'} d r where \dot{Q} = Q^{'} \dot{r}

(17)

It is convenient to express the effective stress damage evolution as

σ (r) = S_{m} (1 - d (r)) r r_{o} ⩽ r

(18)

In equation (17), E_m corresponds to corresponding shear or elastic modulus, and S_m is the effective strength of the composite in the particular direction in question factoring in rate dependency.

Thus, the rate of mechanical dissipation is valid if the damage indexes increase monotonically, thus

\dot{D} = Y \dot{d} ⩾ 0

(19)

We propose from the above analysis a value for H_m as

H_{m} = \frac{{\bar{H}}_{s} l_{e l m}}{1 - {\bar{H}}_{s} l_{e l m}} ⩾ 0 {\bar{H}}_{s} = \frac{S_{m}^{2}}{2 E G_{m}}

(20)

where G_m is the specific fracture energy corresponding to each damage mode, g_m is the energy dissipated per unit volume, l_elm is the characteristic length of the FE, S_m is the effective strength of the composite in the particular direction in question factoring in rate dependency and d(r) depicts that the damage parameter as a function of r, where r₀ corresponds to the point of damage initiation.

The value H_m measures the brittleness of the FE. As observed by Cervera and Chiumenti (2006), in an FE analysis, the state variables of the local model are computed at the integration points in terms of the local stress or strain history. Thus, the characteristic length l_elm is related to the volume or area of the FE. This implies that for a simple beam element, for instance, the characteristic length can be taken as the size of the element. If we restrict the configuration of our elements in this model as equilateral, we will have for equilateral plane triangular element with area A_e, and we will have the characteristic length as equation (21a). In addition, for equilateral tetrahedral and cubic elements of volume V_e, we have the characteristic length as shown, respectively, by equations (21b) and (21c) as follows

l_{e l m}^{* 2} = (\frac{4}{\sqrt{3}}) A_{e}

(21a)

l_{e l m}^{* 3} = (\frac{12}{\sqrt{2}}) V_{e}

(21b)

l_{e l m}^{* 3} = V_{e}

(21c)

Aluminium constitutive model

The aluminium component of the FML was modelled using the famous Johnson–Cook plasticity failure model, which is a special type of Mises plasticity model with analytical forms of the hardening law and rate dependency. Its suitability for high-strain rate deformation of many materials (including most metals) makes it an attractive option in the development of this model. In using this model, adiabatic transient response is usually assumed. Johnson–Cook’s hardening model (a type of isotropic hardening) is explained below.

The strain rate assumption of the Johnson–Cook material model postulates that

\bar{σ} = σ^{0} ({\bar{ε}}^{p l}, θ) R ({\dot{\bar{ε}}}^{p l}) where {\dot{\bar{ε}}}^{p l} = {\dot{ε}}_{0} [\frac{1}{C} (R - 1)]

(22)

where $\bar{σ}$ is the dynamic yield stress, that is, yield stress at non-zero strain rate; ${\dot{\bar{ε}}}^{p l}$ is the equivalent plastic strain rate; ${\dot{ε}}_{0}$ and C are material parameters at or below the transition temperature, θ_transition; $σ^{0} ({\bar{ε}}^{p l}, θ)$ is the static yield stress; and $R ({\dot{\bar{ε}}}^{p l})$ is the ratio of the yield stress at non-zero strain rate to the static yield stress (such that $R ({\dot{ε}}_{0}) = 1.0) .$ Thus, the yield stress is

\bar{σ} = [A + B {({\bar{ε}}^{p l})}^{n}] [1 + C \ln (\frac{{\dot{\bar{ε}}}^{p l}}{{\dot{ε}}_{0}})] [1 - {\hat{θ}}^{m}]

(23)

The effective plastic strain ${\bar{ε}}^{p l}$ (PEEQ) is defined as

{\bar{ε}}^{p l} = {{\bar{ε}}^{p l} |}_{0} + \int_{0}^{t} \sqrt{\frac{2}{3} {\dot{ε}}^{p l} : {\dot{ε}}^{p l}} d t

(24)

where ${{\bar{ε}}^{p l} |}_{0}$ is the initial plastic strain and is usually taken as 0. The effective stress $\bar{σ}$ is defined based on the J₂-plasticty model as follows

\bar{σ} = \sqrt{\frac{3}{2} σ_{i j} σ_{i j}}

(25)

The Johnson–Cook failure model is used in conjunction with its plasticity model in the numerical modelling of failure of the aluminium metal sheets. The failure mechanism is based on the value of the effective (or equivalent) plastic strain at element integration points. Fracture occurs in the aluminium sheet when the damage parameter ω exceeds 1.0. The evolution of ω is given by the accumulated incremental effective plastic strains divided by the current stain at fracture

ω = \sum \frac{Δ {\bar{ε}}^{p l}}{{\bar{ε}}_{f}^{p l}}

(26)

where $Δ {\bar{ε}}^{p l}$ is an increment of the equivalent plastic strain; ${\bar{ε}}_{f}^{p l}$ is the strain at failure, ${\bar{ε}}_{f}^{p l}$ is assumed to be dependent on a non-dimensional plastic strain rate; ${\dot{\bar{ε}}}^{p l} / {\dot{ε}}_{0}$ is a dimensionless pressure–deviatoric stress ratio, p/q (where p is the pressure stress (hydrostatic component of stress tensor corresponding to the spherical, purely dilatational or irrotational part of the deformation) and q is the Mises stress (related to the second invariant of the stress deviatoric tensor J₂ representing equivoluminal distortion)); and the non-dimensional temperature, $\hat{θ}$ , defined earlier in Johnson–Cook’s hardening model. The dependencies are assumed separable and are of the form

{\bar{ε}}_{f}^{p l} = [d_{1} + d_{2} \exp (d_{3} \frac{p}{q})] [1 + d_{4} \ln (\frac{{\dot{\bar{ε}}}^{p l}}{{\dot{ε}}_{0}})] [1 + d_{5} \hat{θ}]

(27)

where d₁–d₅ are non-dimensional failure parameters measured at or below the transition temperature. It is important to note that the failure parameter d₃ is reported as being negative in some literature. However, ABAQUS’ general expression for strain at fracture expects this parameter to be positive. It is worthy of mentioning that failure to properly account for the sign of d₃ will result in an inaccurate response. Choosing a value for the fracture energy, which is used as a data parameter for the damage evolution law, completes the setting of the fracture model. Elements are deleted by default upon reaching maximum degradation according to the usual rules of ABAQUS progressive damage framework. Table 2 shows the values of the plasticity and failure parameters used for Al 2024-0 in this work.

Cohesive elements and interface simulation

The constitutive behaviour of the adhesive layer in this work is described in terms of direct traction versus separation model implemented in ABAQUS (usually preferable for bonded interface where the thickness of the adhesive is negligible and the path of fracture lies essentially within the adhesive) (Luo et al., 2012). This means that the variations of stress and strain through the thickness of the adhesive are disregarded. This usually assumes an initial linear elastic model followed by the initiation and evolution of the damage. The nominal stress vector is related to the nominal strain vector across the interface with an elastic constitutive matrix. The nominal stresses are the force components divided by the original area (similar to the second Piola–Kirchhoff stress) at each integration point, while the nominal strains are the separations divided by the original thickness at each integration point (similar to the Green–Lagrange strain). However, the default choice of the constitutive thickness in terms of traction separation is 1.0 (irrespective of the thickness of the adhesive layer). Thus, the diagonal terms of elastic matrix shown in equation 30 are

K_{n n} = \frac{E_{n}}{t_{c}}, K_{s s} = \frac{E_{s}}{t_{c}}, K_{t t} = \frac{E_{t}}{t_{c}} and ρ = ρ_{c} t_{c}

(28a)

{\begin{matrix} t_{n} \\ t_{s} \\ t_{t} \end{matrix}} = [\begin{matrix} K_{n n} & K_{n s} & K_{n t} \\ K_{s s} & K_{s t} & K_{s n} \\ K_{t t} & K_{t s} & K_{t n} \end{matrix}] {\begin{matrix} ε_{n} \\ ε_{s} \\ ε_{t} \end{matrix}}

(28b)

where t_n represents the normal traction, and t_s and t_t are the two shear tractions (i.e. in the first and second shear directions). The term ρ_c is the actual density of the cohesive element, and ρ is the inputted density. The parameter t_c in equation (28a) denotes the constitutive thickness and should not be confused with a traction. The elasticity matrix in equation (28b) provides a fully coupled behaviour between all components of the traction vector and the separation vector (strain vector). Note the off-diagonal terms in the matrix can be set to 0 if the uncoupled behaviour between stress and strain is desired. This is the case for the model used in this work.

As shown by the graph of typical traction–separation response in Figure 2, the damage initiation refers to the beginning of degradation in material constitutive response at a material point. The damage initiation in damage models begins when the stress and/or strains meet some damage initiation criteria, in this case, a quadratic nominal stress criterion as shown in equation (29)

{\frac{t_{n}}{t_{n}^{0}}}^{2} + {\frac{t_{s}}{t_{s}^{0}}}^{2} + {\frac{t_{t}}{t_{t}^{0}}}^{2} = 1

(29)

where $t_{n}^{0}, t_{s}^{0} and t_{t}^{0}$ represent the peak values of the nominal stress when the deformation is either purely normal to the interface or purely in the first or the second shear direction, respectively. Equation (30) shows the dependency of the fracture energy on the mode mix, where α = 1. This law states that failure under mixed-mode conditions is governed by a power interaction of the energies required to cause failure in the individual (normal and two shears) modes.

{\frac{G_{n}}{G_{n}^{c}}}^{α} + {\frac{G_{s}}{G_{s}^{c}}}^{α} + {\frac{G_{t}}{G_{t}^{c}}}^{α} = 1

(30)

where the mixed-mode fracture energy is G_c = G_n + G_s + G_t when the above equation is satisfied. In the expression, the quantities G_n, G_s and G_t refer to the work done by the traction and its conjugate relative displacement in the normal, the first and the second shear directions, respectively. The specified quantities $G_{n}^{c}, G_{s}^{c} and G_{t}^{c}$ refer to the critical fracture energies required to cause failure in the normal, the first and the second shear direction, respectively. In this work, the thickness of the interface adhesive between aluminium and GFPP is taken as 0.001 mm, and the thickness of adhesive between the GFPP plies is taken as 0.0005 mm. The interface and laminate adhesives are assumed to have the same mechanical properties in this model. The properties of the cohesive layer are shown in Table 5.

Figure 2.

Traction–separation response of cohesive element (the parameter δ is the modal separation (difference in top and bottom displacements)).

Table 5.

Properties of cohesive layer.

Elastic and inertial properties				Damage initiation			Damage evolution
ρ_c (kg/m³)	E_n (GPa)	E_s (GPa)	E_t (GPa)	σ_n (MPa)	σ_s (MPa)	σ_t (MPa)	G_n (J/m²)	σ_s (MPa)	σ_t (MPa)
920	2.05	0.72	0.72	140	300	300	2000	3000	3000

Localised blast

The simulation of the dynamic response in localised blast-loaded FMLs is of importance owing to the differences between the response modes to this type of loading and global blasts. The apparently peculiar nature of this kind of loads renders the dynamic solution, damage pattern, thresholds and critical conditions intrinsically different. When localised blast loads are applied to plates, it is normal practice to idealise the localised blast as pressure loading relating to impulse obtained from actual experiment. For localised blast at the centre of a plate, the pressure loading is a function of time and distance from the centre of the plate (Karagiozova et al., 2010).

However, for FMLs, like for other plated systems, the response to localised blast includes two regimes: one associated with the initial through-thickness compression phase, the subsequent reflected tensile from the back face, and the other related to the overall response. These characteristics make the time history of localised blast load very important.

This work has been validated by the experimental results generated by other researchers at the University of Cape Town, the loading comprises a disk of 8-g explosive (cylindrical charge) with a diameter of d₀ = 30 mm and a leader of 1 g mass which was placed at the centre of the panel. Results from AUTODYN simulation show that the pulse is well described by an exponentially decaying function exp (−2t/t₀) in time with t₀ = 0.008 ms (Karagiozova et al., 2010). The spatial distribution of the loading shows that the leader influences the shape of the pressure pulse by producing a variable pressure within the area bounded by the disk of explosive, r ⩽ d₀/2 (Karagiozova et al., 2010), where d₀ denotes the diameter of the charge. Figure 3 shows the temporal and spatial distributions of the localised blast load.

Figure 3.

(a) Temporal pressure–time history at a generic fixed position x₀ and (b) spatial pressure–time history at a generic time t₀.

The pressure function, P (r, t) used in this work is

P (r, t) = p_{1} (r) p_{2} (t)

(31)

where

p_{1} (r) = {\begin{matrix} P_{0} \exp [- k (r - r_{0})] r_{0} < r ⩽ r_{b} \\ P_{0} \exp [- k (r - r_{0})] r_{0} < r ⩽ r_{b} \\ 0 r > r_{b} \end{matrix}

(32a)

p_{2} (t) = \exp (\frac{2 t}{t_{0}})

(32b)

In the above equation, r₀ = 15 mm which is the radius of the explosive disc used in the experiments, (k) = (m⁻¹) is an exponential decay constant, which models the pressure distribution over the exposed area of the plate, r_b < L/2, where L is the length of the panel, and t₀ is the characteristic decay time for the pulse (Karagiozova et al., 2010). The decay constant, k, depends on the ratio d₀/L through the relationship

k = 130 - 261 (\frac{d_{0}}{L}) + 948 {(\frac{d_{0}}{L})}^{2} 0.15 < \frac{d_{0}}{L} < 0.6

(33)

The range d₀/L was restricted by the experimental configurations used; however, the above expression is still valid for current configuration, d₀/L = 0.1. Thus, k ≈ 114 m⁻¹. The total impulse imparted to the plate is given by the expression

I = 2 π \int_{0}^{\infty} \int_{0}^{r_{b}} P (r, t) r d r d t

(34)

Figure 4 shows the spatial distribution of local blast on the centre of the FML plate. Table 6 shows the values of calculated peak pulse, P₀. It must be noted that the spatial distribution of the load is an axisymmetric function.

Figure 4.

Panel shows distribution of applied localised blast load.

Table 6.

Computed peak pulse pressure.

Lay-ups	No. of layers	Thickness (mm)	Impulse (N s)	Pressure, P₀ (MPa)
A2T14-2	6	3.36	5.89	741.5
A2T18-3	10	5.35	6.17	775.1
A2T18-4		5.6	7.94	1000.1
A3T22-2	7	4.09	7.57	948.2
A3T22-4		4.13	7.7	975.1
A3T24-5	11	6.08	10.58	1326.3
A3T24-7		6.27	3.76	471
A3T24-8		6.06	7.85	984
A3T26-1	15	8.49	7.8	980
A3T26-3		8.1	9.54	1196
A3T26-4		8.41	11.29	1415.5
A3T28-4	19	9.84	12.43	1560.1
A3T28-5		9.82	10.34	1296.2

Results and discussion

A subroutine (VDLOAD) was developed to depict the localised blast load illustrated in equations (32) to (34). The loads are applied to the developed FML model, and the results are compared with the experimental results presented by Langdon et al. (2007a). The significant back face deflection was observed as the first response mode of the structures in the numerical models. This can be attributed to the through-thickness-reflected tensile wave propagation in the FML as the early time response, scabbing or spalling (regime I of response). Figure 5(b) shows large localised back face deflection with no visible damage in the composites. This can be compared to the diamond-shaped deflected back face of panel A2T14-2 shown in Figure 5(a). Figure 5(c) shows the displacement time history at the midpoint of the plate. This is obviously the point of maximum deflection.

Figure 5.

(a) Photograph showing the back faces of the FML panels form experiment for A2T18-3 (Langdon et al., 2007a), (b) numerical simulation for FML panel A2T18-3 for 1 mm × 1 mm mesh at 0.15 ms and (c) displacement time history in the middle (front face) for A2T18-4.

At increased blast loads, the numerical model in Figure 6(a) and (b) shows large plastic deformation and back face delamination. The displacements of the back face and middle start to diverge and more pronounced debonding is observed. Figure 6(c) shows the displacement time history in the middle of the plate.

Figure 6.

(a) Photograph showing the back faces of the FML panels form experiment for panel A3T28-4 (Langdon et al., 2007a), (b) numerical simulation for FML panel A2T28-4 for 1.5 mm × 1.5 mm mesh at 0.15 ms and (c) displacement time history in the middle (front face).

At increased blast loads, the numerical model in Figure 7(b) and (c) captures the tearing and eventual perforation of panel A2T14-1 as shown in Figure 7(a). Figure 7(b) shows the numerical simulation of the FML with a mesh size of 1 mm × 1 mm and corresponding softening modulus, while Figure 7(c) shows a mesh size of 1.5 mm × 1.5 mm. Both figures show an average of 75% damage in warp direction in the composite patches close the centre of the panel. Figure 7(a) and (b) shows the failure of the front and back faces of panel A3T24-3 which compares favourably the experimental results in Figure 7(a). Both the 1 mm × 1 mm mesh and 1.5 mm × 1.5 mm mesh were able to capture the failure in the composites with adequate softening modulus. Figure 8(a) and (b) shows similar results as the case depicted in Figure 7.

Figure 7.

(a) Photograph showing the back faces of the FML panels form experiment for A2T14-1 (Langdon et al., 2007a), (b) numerical simulation for FML panel A2T14-1 for 1 mm × 1 mm mesh at 0.01 ms and (c) numerical simulation for FML panel A2T14-1 for 1.5 mm × 1.5 mm mesh at 0.01 ms.

Figure 8.

(a) Photograph showing the back faces of the FML panels form experiment for A3T24-3 (Langdon et al., 2007a) and (b) numerical simulation for FML panel A2T24-3 for 1.5 mm × 1.5 mm mesh at 0.01 ms.

A good correlation is observed for the maximum displacements of the back and front face panels between the numerical simulation and experimental models of the FML panels. The displacements of the FML panels correlated well as shown in Table 7.

Table 7.

Comparison of displacement of FML panels.

Panel	Mass PE 4 (g)	Mean thickness in (mm)	Impulse (N s)	Experimental back face deflection	Experimental front face deflection	Numerical back face deflection	Numerical front face deflection
AT14-1	2.5	1.01	3.38	N/A	N/A	N/A	N/A
A2T18-4	2.8	5.60	7.94	16.3	9.1	16.0	8.1
A3T28-4	5.0	9.84	12.43	26.7	4.1	26	4.9
A2T24-3	5.3	1.74	6.06	N/A	N/A	N/A	N/A

FML: fibre–metal laminate; PE: plastic explosive.

Efficiency of the charge and effectiveness of the target

The efficiency of localised blast on the target can be represented in a non-dimensional form for correlation studies and investigation into the effect of such charges on the target. The equations of spatial and temporal distributions of the above localised blast load on the FML can be rewritten as follows

\begin{array}{l} p_{s} (r) = {\begin{matrix} P_{0} & r ⩽ R_{0} \\ P_{0} \exp [- k (R_{0} - r)] & R_{0} < r ⩽ + \infty \end{matrix} \\ p_{t} (t) = {\begin{matrix} e^{- ω t} & 0 ⩽ t < t_{d} \\ 0 & t > t_{d} \end{matrix} \end{array}

(35)

Thus, for $r < R_{0},$ it can be shown that the total impulse $I_{r < R_{0}} (r, t_{d})$ is

\lim_{t_{d} \to 0^{+}} I_{r < R_{0}} (r, t_{d}) = π P_{0} r^{2} \lim_{t_{d} \to 0^{+}} (\frac{1 - e^{- ω t_{d}}}{ω}) = π r^{2} P_{0} t_{d}

(36)

In addition, for $R_{0} < r < \infty,$ it can be shown that the total impulse, $I_{R_{0} < r < \infty} (r, t_{d})$ is

\begin{matrix} \lim_{t_{d} \to 0^{+}} I_{R_{0} < r < \infty} (r, t_{d}) = (\frac{2 π P_{0}}{b^{2}} + \frac{2 π P_{0} R_{0}}{b} - \frac{2 π P_{0} e^{b R_{0}}}{b^{2} e^{b r}} - \frac{2 π P_{0}}{b e^{b r}}) \lim_{t_{d} \to 0^{+}} (\frac{1 - e^{- ω t_{d}}}{ω}) \\ = (\frac{2 π}{b^{2}} + \frac{2 π R_{0}}{b} - \frac{2 π e^{b R_{0}}}{b^{2} e^{b r}} - \frac{2 π}{b e^{b r}}) P_{0} t_{d} \end{matrix}

(37)

Introducing non-dimensional parameters, it follows that

λ = \frac{r}{R_{0}} and I_{r} (λ) = \frac{I (r, t_{d})}{π R_{0}^{2} P_{0} t_{d}}, η = \frac{I_{r} (λ)}{I_{r} (\infty)}

(38)

where η is referred to as the efficiency of the charge, and I_r is the non-dimensional total impulse on the FML plate. After simplification, η can be written as

η = \frac{I_{r} (λ)}{I_{r} (\infty)} = \frac{1 + \frac{2}{R_{0}^{2} b^{2}} + \frac{2}{R_{0} b} - \frac{2 e^{b R_{0}}}{R_{0}^{2} b^{2} e^{b λ R_{0}}} - \frac{2 e^{b R_{0}} λ}{R_{0} b e^{b λ R_{0}}}}{1 + \frac{2}{R_{0}^{2} b^{2}} + \frac{2}{R_{0} b}}

(39)

It can be shown that for the configuration studied, the efficiency of the charge is 99.5%. This indicates that in these scenarios, the full impulse of the charge is taken by the system. Effectiveness of the target, on the other hand, in a particular localised blast scenario can be defined as the ratio of the impulse corresponding to a particular damage level (most frequently, the tearing threshold damage) to the total impulse delivered to it in that blast scenario, that is

ξ = \frac{I_{c}}{I_{r} (λ)}

(40)

where the effectiveness of the target is denoted by $ξ,$ I_c is the impulse corresponding to a threshold damage level and the denominator is already defined. If the threshold corresponds to perforation, then this relation would determine how efficient the target is, that is, if $ξ > 1,$ then the target is fully functional and efficient. On the other hand, if $ξ < 1,$ the target will perforate and cannot withstand the impulse exerted to it. The lower the efficiency parameter, the higher is the degree of damage and extent of perforation. One final remark, however, must be made with regard to the definitions above. These definitions only hold if the load is fully characterised by the total impulse which is equivalent to assuming impulsive response regime for the structure. In any such sort of scenario, the specific pulse loading time history can be disregarded and total impulse rendered only the relevant loading parameter.

Conclusion

Several aspects of the study presented in this article are validated against the experimental data for localised blast-loaded FMLs. Thus, using numerical simulations, one is provided by a cost-effective, predictive and efficient way of virtual testing and preliminary study for the response of FMLs without having to go through the rigours of blast testing physical panel prototypes. The response of FMLs has been found to be highly dependent on the thickness of the loaded plate and the characteristics of underlying layers. Thus, it is of crucial significance that in the behaviour of FMLs, highly nonlinear transient dynamic phenomena (i.e. the large plastic deformation, debonding, delamination, tearing and failure of the constituent components of the plate) be captured in order to be able to describe the full response of the FML in such a loading environment.

Layers at the back in a typical FML support layers at the front providing the reason why the back face easily debonds from the rest of the plate when loaded. Thus, the failure of the intermediate composites in the FMLs needs to be adequately captured in order to be able to predict the health of the top aluminium plates. If the composite layers are modelled improperly and do not fail when they ought to fail as observed in the experiments, the maximum deflection of the system would be predicted incorrectly and thus, the overall response of the system rendered under-predicted. As a result of this, subsequent tearing of the aluminium plate might not be captured in the model. The developed model captures the state/health of the aluminium and composite in the FML and the maximum displacements of the back and front faces.

The concepts of charge efficiency and target effectiveness are defined and relevant equations are derived. These are important as far as defensive shields are concerned.

One final remark regarding the constitutive models used in this work is in order. While several nonlinear phenomena have been considered and the corresponding material models are incorporated, the hydrodynamic equations of state (EoSs) have not been used for the target. The loading is obtained as a result of a full fluid–structure interaction (FSI) simulation in AUTODYN; however, material EoSs such as those of Rankine–Hugoniot or Mie–Grüneisen are not included (Cheng et al., 2007; Hanstrom and Lazor, 2000). Such models, while necessary to simulate high-rate phenomena as ballistic perforation are of little use for charge sizes and stand-offs of this study. A more comprehensive study would be needed to show with certainty the little relevance of these EoSs in a blast loading scenario replicated here as well as to give bounds of relevance on charge sizes, stand-offs or related dimensionless parameters.

Since the constitutive models incorporate strain rate effects in metal and composite parts and an objectivity algorithm for strain softening to control energy dissipation associated with each failure mode regardless of mesh refinement and topology, the confidence of the analyst in determining the response of the intermediate composite layers as well as the global response is improved.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship and/or publication of this article.

References

Al-Ani

Habibi

(2013) A new energy optimization strategy for pumping operation in water distribution systems. In: ASME 2012 International Mechanical Engineering Congress and Exposition (vol. 6, Parts A and B), Houston, TX, 9–15 November, pp. 1613–1621. New York: ASME.

Batra

Gopinath

Zheng

(2012) Damage and failure in low energy impact of fiber-reinforced polymeric composite laminates. Composite Structures 94: 540–547.

Bažant

(1983) Crack band theory for fracture of concrete. Matériaux et Construction 16: 155–177.

Belytschko

Black

(1999) Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering 45: 601–620.

Cervera

Chiumenti

(2006) Mesh objective tensile cracking via a local continuum damage model and a crack tracking technique. Computer Methods in Applied Mechanics and Engineering 196: 304–320.

Cheng

Itoh

Jen

. (2007) A new analytical model for high-velocity impact of thick composites. International Journal of Crashworthiness 12: 57–65.

Coleman

Gurtin

(1967) Thermodynamics with internal variables. Journal of Chemical Physics 47: 597–613.

Donadon

Falzon

Iannucci

. (2007a) A 3-D micromechanical model for predicting the elastic behaviour of woven laminates. Composites Science and Technology 67: 2467–2477.

Donadon

Falzon

Iannucci

. (2007b) Intralaminar toughness characterisation of unbalanced hybrid plain weave laminates. Composites Part A: Applied Science and Manufacturing 38: 1597–1611.

10.

Donadon

Iannucci

Falzon

. (2008) A progressive failure model for composite laminates subjected to low velocity impact damage. Computers & Structures 86: 1232–1252.

11.

Fleisher

(1996) Design and explosive testing of a blast resistant luggage container. In: Proceedings of the structures under shock and impact conference IV, Udine, July 1996, pp. 51–59. New Forest, UK: WIT press.

12.

Fleisher

(1996) Design and explosive testing of a blast resistant luggage container. In: Jones

Watson

Brebbia

(eds) Structures under Shock and Impact IV. Southampton: WIT Press, pp. 51–59.

13.

Hanstrom

Lazor

(2000) High pressure melting and equation of state of aluminium. Journal of Alloys and Compounds 305: 209–215.

14.

Hashin

(1980) Failure criteria for unidirectional composites. Journal of Applied Mechanics 47: 329–334.

15.

Hassan

Batra

(2008) Modeling damage in polymeric composites. Composites Part B: Engineering 39: 66–82.

16.

Iannucci

(2006) Progressive failure modelling of woven carbon composite under impact. International Journal of Impact Engineering 32: 1013–1043.

17.

Iannucci

Willows

(2006) An energy based damage mechanics approach to modelling impact onto woven composite materials – part I: numerical models. Composites Part A: Applied Science and Manufacturing 37: 2041–2056.

18.

Iannucci

Willows

(2007) An energy based damage mechanics approach to modelling impact onto woven composite materials: part II. Experimental and numerical results. Composites Part A: Applied Science and Manufacturing 38: 540–554.

19.

Johnson

(2001) Modelling fabric reinforced composites under impact loads. Composites Part A: Applied Science and Manufacturing 32: 1197–1206.

20.

Karagiozova

Langdon

Nurick

. (2010) Simulation of the response of fibre–metal laminates to localised blast loading. International Journal of Impact Engineering 37: 766–782.

21.

Ladeveze

LeDantec

(1992) Damage modelling of the elementary ply for laminated composites. Composites Science and Technology 43: 257–267.

22.

Langdon

Lemanski

Nurick

. (2007a) Behaviour of fibre–metal laminates subjected to localised blast loading: part I – experimental observations. International Journal of Impact Engineering 34: 1202–1222.

23.

Langdon

Nurick

Lemanski

. (2007b) Failure characterisation of blast-loaded fibre–metal laminate panels based on aluminium and glass–fibre reinforced polypropylene. Composites Science and Technology 67: 1385–1405.

24.

Lemanski

Nurick

Langdon

. (2006) Understanding the behaviour of fibre metal laminates subjected to localised blast loading. Composite Structures 76: 82–87.

25.

Lemanski

Nurick

Langdon

. (2007) Behaviour of fibre metal laminates subjected to localised blast loading – part II: quantitative analysis. International Journal of Impact Engineering 34: 1223–1245.

26.

Luo

Liang

Chen

. (2012) Structural optimization of variable section slender manipulator based on sensitivity analysis method. Mechatronics and Information Technology 2–3: 291–295.

27.

Matzenmiller

Lubliner

Taylor

(1995) A constitutive model for anisotropic damage in fiber-composites. Mechanics of Materials 20: 125–152.

28.

Moës

Belytschko

(2002) Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics 69: 813–833.

29.

Nurick

Langdon

Chi

. (2009) Behaviour of sandwich panels subjected to intense air blast – part 1: experiments. Composite Structures 91: 433–441.

30.

Paley

Aboudi

(1992) Micromechanical analysis of composites by the generalized cells model. Mechanics of Materials 14: 127–139.

31.

Puck

Schürmann

(1998) Failure analysis of FRP laminates by means of physically based phenomenological models. Composites Science and Technology 58: 1045–1067.

32.

Sukumar

Moës

Moran

. (2000) Extended finite element method for three-dimensional crack modelling. International Journal for Numerical Methods in Engineering 48: 1549–1570.

33.

Vlot

(1996) Impact loading on fibre metal laminates. International Journal of Impact Engineering 18: 291–307.

34.

Guan

Cantwell

. (2012) Low-impulse blast behaviour of fibre-metal laminates. Composite Structures 94: 954–965.

35.

Guan

Cantwell

. (2013) Modelling of the low-impulse blast behaviour of fibre–metal laminates based on different aluminium alloys. Composites Part B: Engineering 44: 141–151.

36.

Wang

(2009) Calculation of stress intensity factors for interfacial cracks in fiber metal laminates. Available at: http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20090019750.pdf

37.

Xiao

Gama

Gillespie

(2007) Progressive damage and delamination in plain weave S-2 glass/SC-15 composites under quasi-static punch-shear loading. Composite Structures 78: 182–196.

Numerical simulation of the dynamic response in pulse-loaded fibre–metal-laminated plates

Abstract

Keywords

Introduction

Geometry and modelling

Material modelling

Composite modelling

Composite material properties

Elastic damage energy and dissipation

Damage initiation criteria

Fill and warp fibre tensile/shear failure mode

Fibre compressive failure mode

Fibre crush failure mode

Fibre in-plane shear failure mode (planes 1 and 2)

Fibre out-of-plane shear failure mode (planes 2 and 3)

Fibre in-plane shear failure mode (planes 1–3)

Strain rate effect

Damage progression criteria

Fracture energy (energy release rate) and mesh objectivity

Aluminium constitutive model

Cohesive elements and interface simulation

Localised blast

Results and discussion

Efficiency of the charge and effectiveness of the target

Conclusion

Footnotes

Declaration of conflicting interests

Funding

References