Composites Part A 115 (2018) 31–45

Contents lists available at ScienceDirect

Composites Part A
journal homepage: www.elsevier.com/locate/compositesa

Finite element modelling of Dyneema® composites: From quasi-static rates T

to ballistic impact
⁎
Mark K. Hazzarda,b, , Richard S. Traska, Ulrich Heissererb, Mirre Van Der Kampb,
Stephen R. Halletta
a
Advanced Composites Centre for Innovation and Science, University of Bristol, Queen's Building, Bristol BS8 1TR, UK
b
DSM Dyneema, P.O. Box 1163, 6160 BD Geleen, The Netherlands

A R T I C LE I N FO A B S T R A C T

Keywords: A finite element methodology to predict the behaviour of Dyneema® HB26 fibre composites at quasi-static rates
Impact behaviour of deformation, under low velocity drop weight impact, and high velocity ballistic impact has been developed. A
Finite element analysis (FEA) homogenised sub-laminate approach separated by cohesive tied contacts was employed. The modelling ap-
Dyneema® proach uses readily available material models within LS-DYNA, and is validated against experimental ob-
servations in literature. Plane-strain beam models provide accurate mechanisms of deformation, largely con-
trolled through Mode II cohesive interface properties and kink band formation. Low velocity drop weight impact
models of HB26 give force-deflection within 10% of new experimental observations, with in-plane shear strain
contour plots from models directly compared with experimental Digital Image Correlation (DIC). Ballistic impact
models utilising rate effects and damage showed similar modes of deformation and failure to that observed in
literature, and provide a good approximation for ballistic limit under 600 m/s impact speed.

1. Introduction lamina to initiate delamination. Nazarian & Zok [12] modelled the
shear response of a ± 45° tension shear coupon using shell elements
Ultra-high molecular weight polyethylene (UHMWPE) fibre com- that were tied to rebar elements, representing the matrix and fibres
posites are increasingly being used in impact scenarios due to their respectively. This captured the scissoring of fibres under large in-plane
extremely high specific tensile strength and stiffness [1]. Capturing the shear deformation, however at large strains the response was too stiff.
material behaviour through numerical modelling techniques over a Chocron et al. [13] modelled UD strips as well as laminates of Dy-
range of loading conditions provides the prospect of being able to neema® HB80 by bundling the fibres together into larger macro-scale
predict current performance for use as a design tool, and also provide strips of solid elements representing the fibre and adjoining matrix
insight into the complex deformation and failure mechanisms observed material. Capturing the meso-microscale architecture through an in-
during highspeed impact [2,3] that can be difficult to fully capture crease in scale allowed the model to be more computationally efficient,
analytically [4]. Several dimensional hierarchies exist within UHMWPE however there were limitations in the overall response of the laminate,
composites, with the scale of modelling used being dependent on a such as exhibiting a larger back-face deflection (BFD) following impact.
trade-off between the detail of phenomena required to be captured and Continuum models can provide a reasonable alternative at reduced
the computational expense. Validation of current modelling techniques computational expense, however there has been great difficulty pro-
is also just starting to come to fruition as increased test data from viding an accurate response, particularly when incorporating failure
mechanical tests [2,5,6] and ballistic impacts becomes available and damage, largely due to the poor load transfer at the fibre-fibre level
[7–10]. caused by the low mechanical properties of the adjoining matrix, which
Liu et al. [5,11] performed a static analysis of short beam shear and provides attractive ballistic performance [14]. This inherent property
long beam bending both experimentally and numerically. Numerically, contrasts with finite element theory, where perfect load transfer is as-
a finite element plane strain 2D meso-scale model captured deformation sumed within each element. Grujicic et al. [15] modelled [0°/90°]n
mechanisms under quasi-static loading conditions to investigate the impacts with homogenised solid elements using a custom continuum
flexural response, utilising cohesive zones between uni-directional (UD) unit cell response based on micro-scale modelling of a single [0°/90°]

⁎
Corresponding author at: DSM Dyneema, P.O. Box 1163, 6160 BD Geleen, The Netherlands.
E-mail address: [email protected] (M.K. Hazzard).

https://doi.org/10.1016/j.compositesa.2018.09.005
Received 3 May 2018; Received in revised form 24 July 2018; Accepted 4 September 2018
Available online 07 September 2018
1359-835X/ © 2018 Elsevier Ltd. All rights reserved.
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

yarn cross over. The difficulty in capturing delamination and BFD fol- beam response was investigated using the sub-laminate approach,
lowing impact was observed. Iannucci et al. [16] proposed a con- showing the dependence on the number of elements per sub-laminate
stitutive model for Dyneema® [0°/90°]n laminates, homogenising cross- and the cohesive interface properties (Fig. 2b). In a separate model,
ply UD layers to a shell model, providing similar deformation me- experimental and model comparisons of the strain field observed in
chanisms with clear drawing-in of the laminate due to the low shear drop weight impact simulations were made, investigating the response
properties. In a detailed study with extensive material characterisation, prior to perforation of the laminate (Fig. 2c). Finally, 3D ballistic impact
Lässig et al. [2] produced a non-linear orthotropic hydrocode model investigations were performed, making use of the damage and failure
using a solid element continuum model for impacts on Dyneema® HB26 criteria, comparing performance in terms of ballistic limit as well as
in ANSYS Autodyn. Whilst results were in good agreement for the BFD with literature experimental results (Fig. 2d).
ballistic limit, BFD was unable to be captured due to the difficulty in
capturing imperfect out-of-plane load transfer. Nguyen et al. [3,14] 2.3. Material model
recently advanced this model to better capture the ballistic limit and
BFD by using a sub-laminate approach, with sub-laminates separated by MAT162 is a commercially available composite damage and failure
tied contacts in order to represent delamination. A similar sub-laminate model produced by Materials Science Corporation [18] and accessible
approach has been used here, also in conjunction with an adjoining within LS-DYNA. It was selected to model the behaviour of the homo-
fracture mechanics based contact to represent inter-laminar matrix genised cross-ply sub-laminates. All 6 stress components contribute
adhesion in LS-DYNA. For the first time, this work validates a numerical towards damage of the elements, based on progressive Hashin failure
model for HB26 over a range of loading conditions, initially at low rates criteria [19] followed by strain softening, controlled by damage vari-
comparing beam stiffness to that observed in literature, followed by ables. A list of the implemented failure criteria is given in Table 1,
drop weight impact tests, and finally, ballistic impact. where subscripts a, b, and c are local element axes with a and b being in
the laminate plane and c being through the thickness. Subscripts T and
2. Finite element modelling philosophy C denote tension and compression, and S a delamination scaling factor,
and r7–13 are cross-ply damage thresholds [20].
2.1. Material Individual damage variables are coupled through the matrix qij:

The material selected, due to the large amount of data available was ⎡1 0 1 0 1 0 0⎤
Dyneema® HB26, a cross-ply uni-directional laminate with approxi- ⎢0 1 0 1 1 0 0⎥
0 0 0 0 1 0 1⎥
mately 80% fibre volume fraction connected by a polyurethane matrix. qij = ⎢
⎢1 1 1 1 1 1 0⎥
SK76 fibres are used with a strength of 3.6 GPa and a modulus of up to ⎢0 1 0 1 1 0 1⎥
⎢1 0 1 0 1 0 1⎥ (8)
130 GPa at ε ̇ > 1. Plies of material, which are provided in a [0°/90°/0°/ ⎣ ⎦
90°] configuration, are stacked to continue the cross-ply configuration
where i represents the 6 loading directions (damage variable) and j
to the required areal density and hot-pressed to adhere plies to one
represents the 7 damage modes (7–13 for cross-ply arrangement). A
another and produce a hard ballistic panel. Filament diameter is un-
non-linear damage progression model is used for all damage types but is
changed by the pressing process, with each uni-directional component
dependent on 4 pre-defined damage parameters, am1–4, with subscripts
of the cross ply approximately 67.5 μm in thickness.
being 1 for fibre damage in the a direction, 2 for fibre damage in the b
direction, 3 for crush damage in the c direction, and 4 for matrix shear
2.2. Model description damage. The logarithmic damage function is given by Eq. (9):

All finite element analysis (FEA) was performed in the non-linear ⎧ ⎛ ⎛ ri

ami
⎞ ⎞⎫
⎜− ami ⎟
explicit solver LS-DYNA, with all modelling utilising version mpp971 Ei' = 1−⎜1−exp ⎝ ⎠ ⎟ Ei
⎨ ⎜ ⎟⎬
revision 7.1.2. All models were ran using 16 processors of 2.6 GHz
⎩ ⎝ ⎠⎭ (9)
cores, with large models requiring high memory nodes and job times
variable dependent on simulation type. Modelling HB26 at the micro- where Ei is the initial stiffness, ami is the damage factor, and ri is the
scale and meso-scale models for HB26 can be computationally ex- relevant damage threshold. Strain rate effects are also considered, and
pensive due to the low ply thicknesses combined with the low density of are captured through modification of stiffness and strength. There are 4
Dyneema®, which can lead to relatively thick laminates for ballistic rate parameters, Crate1–4, where 1 represents strain rate dependent in-
applications containing many plies. As such a cross-ply sub-laminate plane strength, 2 represents the a direction axial moduli, 3 represents
homogenisation approach, separated by cohesive tied contacts to cap- shear moduli, and 4 represents the b direction axial moduli. The loga-
ture inter-laminar matrix behaviour, was utilised (Fig. 1). rithmic law implemented is given in Eq. (10).
A hierarchical approach has been taken to investigate the model
response, first selecting parameters for a homogenised response and ε̇
SRT = S0 ⎛1 + Crate1ln ⎞
⎜ ⎟

performing single element studies (Fig. 2a). Following this the laminate ⎝ ε0̇ ⎠ (10)

Fig. 1. Schematic of the sub-laminate cross-ply homogenisation approach to capture laminate behaviour.

32
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

Fig. 2. Hierarchical modelling approach, (a) homogenisation, (b) beam stiffness, (c) drop weight impact, and (d) ballistic impact. (For interpretation of the references
to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1 Table 2
Failure criteria in MAT162, composite MSC [20]. MAT162 Dyneema® HB26 material card.
Failure mode Criterion Property Symbol Value Units Ref.

Fibre tension/shear (a direction) 2 2 (1) Density ρ 0.00097 g/mm 3

–
( ) + ( ) −r
σa
SaT
τca
SaFS
2
7 Young’s Modulus aa Eaa 34257 MPa [23]
Fibre tension shear (b direction) σ
2
τ
2 (2) Young’s Modulus bb Ebb 34257 MPa [23]
⎛ b ⎞ + ⎛ cb ⎞ −r82 = 0 Young’s Modulus cc Ecc 3260 MPa [2]
⎝ SbT ⎠ ⎝ SbFS ⎠
Fibre compression (a direction) 2 (3) Poisson’s Ratio ba νba 0 – –
( ) −r
σ 'a
SaC
2
9 =0 Poisson’s Ratio ca
Poisson’s Ratio cb
νca
νcb
0.013
0.013
–
–
[2]
[2]
Fibre compression (b direction) σ'
2 (4)
⎛ b ⎞ −r10
2
=0 Shear Modulus ab Gab 173.8 MPa [24]
⎝ SbC ⎠ Shear Modulus ca Gca 547.8 MPa [13]
Fibre Crush (c direction) 2 (5)
( ) −r
σc
SFC
2
11 =0 Shear modulus cb
Longitudinal tensile strength a
Gcb
SaT
547.8
1250
MPa
MPa
[13]
[25]
Matrix shear (ab direction) τ
2 (6) Longitudinal compressive strength a SaC 1250 MPa –
⎛ ab ⎞ −r12
2
=0
⎝ Sab ⎠ Transverse tensile strength b SbT 1250 MPa [25]
Matrix delamination 2 2 2 (7) Transverse compressive strength b SbC 1250 MPa –
S2 ⎧
⎨
⎩
( )
σc
ScT
τ
+ ⎛ bc ⎞ +
⎝ Sbc ⎠ ( ) ⎫⎬⎭−r
τca
Sca
2
13 =0
Through thickness tensile strength c ScT 1E20 MPa -
Crush strength c SFC 1250 MPa [26]
Fibre mode shear strength SFS 625 MPa –
Matrix mode shear strength ab Sab 1.8 MPa [24]
Matrix Mode Shear Strength bc Sbc 1.8 MPa [24]
Matrix mode shear strength ca Sca 1.8 MPa [24]
Residual compressive strength factor SFFC 0.1 – –
Failure Model (2 = fabric cross-ply) AMODEL 2 – –
Coulomb friction angle PhiC 0 degrees –
Delamination scale factor SDELM 1 – –
Limiting damage factor OMGMX 0.999 – –
Eroding axial strain E_LIMT 0.06 – –
Eroding compressive volume strain ECRSH 0.05 – –
Eroding volumetric strain EEXPN 4 – [27]
Coefficient of strain rate fibre strength Crate1 0.0287 – [23]
properties
Coefficient of strain rate for axial moduli Crate2 0.1163 – [23]
Coefficient of strain rate for shear moduli Crate3 0.225 – [23]
Coefficient of strain rate for transverse Crate4 0.1163 – [23]
moduli
Coefficient of softening for axial fibre am1 20 – –
damage
Coefficient of softening for transverse fibre am2 20 – –
Fig. 3. The effect of the in-plane damage parameter for MAT162 with HB26 damage
properties. (For interpretation of the references to colour in this figure legend, Coefficient of softening for crush damage am3 20 – –
Coefficient of softening for matrix failure am4 −0.8 – [3]
the reader is referred to the web version of this article.)

where SRT is the updated strength, S0 is the original input strength, ε ̇ is references given for the source of the data. Further description and
the current strain rate, and ε0̇ = 1 is the reference strain rate. Sub- breakdown of how particular parameters were derived are given within
stituting stiffness for strength in Eq. (10) provides rate dependent this section and investigated using single element testing. As MAT162 is
moduli. Element erosion can also occur, triggered by a limiting strain in restricted to single integration point elements, non-physical modes of
directions a and b, by reaching a compressive volume limit through deformation can occur. To counter this an hour glassing methodology of
crushing, or through reaching a volumetric limit. Element erosion type 6, Belytschko-Bindeman for solid elements [22], based on elastic
strains were selected to avoid eroding elements whilst they still con- constants with default coefficients, was used. Minimal constraints (1
tained large internal energy to limit energy loss within the system. This node fully fixed) were placed on single elements, allowing for Poisson’s
was done by using eroding strains that follow damage (see Fig. 3 ex- contraction and expansion.
ample for tension) and using high volumetric limits. The material
model has previously been fully validated against ballistic testing of
S2–glass composites [21], allowing for a direct comparison with HB26. 2.4.1. In-Plane tension & compression
In-plane UD properties were approximated from quasi-static fibre
2.4. Single element response and matrix properties due to the difficulty in measuring laminate
properties from non-uniform loading during mechanical tests [24]. The
The MAT162 material card information is given in Table 2 with rule of mixtures assuming an 83% fibre volume fraction [23] was used

33
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

Table 3 1.6
Homogenising in-plane properties of HB26 under quasi-static conditions. 1.5

Dimensionless Stiffness
Property Value Ref. 1.4
Ef (GPa) 82 [23] 1.3
Em (MPa) 60 [12] 1.2
vf (–) 0.83 [24]
E11 (GPa) 68 – 1.1
E22 (MPa) 410 –
1 Experimental
Russell DataData
et al Exp' [23][23]
G12 (MPa) 86.9 [24]
ν12 (–) 0.2832 [28] 0.9 y = 0.1163ln(x) + 1
*Ea (GPa) 34.257 –
*Eb (GPa) 34.257 –
0.8
*Gab (MPa) 173.8 [24] 1 10 100
*νab (–) 0 – Strain Rate (1/ms)
* = cross-ply
Fig. 5. The effect of Crate2 and Crate4 on the tensile longitudinal and transverse
moduli. (For interpretation of the references to colour in this figure legend, the
to determine the UD response (Table 3), where E11 and E22 principal reader is referred to the web version of this article.)
stiffnesses for a UD lamina, and ν12 the Poisson’s ratio, and v and E with
subscripts f and m represent fibre and matrix volume fraction and node numbering which allows the local elemental axis to rotate.
stiffness respectively. This was then transformed to a laminate response However, this does not capture scissoring leading to an angle reduction
using classical laminate analysis (CLA) which provided an equivalent between principal axes (Fig. 6).
homogenised [0°/90°]n response. As such, fibre scissoring must be included within the material model
In-plane tensile strength was approximated from single ply tests and itself to account for additional energy absorption. A negative damage
digital image correlation of laminate tensile tests [24,25]. The in-plane parameter in Eq. (9) has been utilised allowing elements to absorb
damage parameter (am1 = am2 = 20) was selected to minimise energy strain energy following in-plane shear damage. Shear testing of ± 45°
absorption post damage, similar to that observed in single fibre, yarn, specimens highlighted a large dependency of strength on the gauge
and ply tests [23,25], represented in Fig. 3. region geometry [3,17,24]. As this model is to be used in impact sce-
In-plane compression of laminates has been shown to initiate micro- narios with large laminate widths (complete in-plane shear failure not
buckling and kink-band formation occurring at a low stress compared observed), the data for the largest available gauge width was utilised
with tension [29]. The compressive in-plane strength was assumed [3]. MAT162 negative damage parameters allow stiffness to be altered
equal to that in tension, and the kink band formation and buckling was and approach a new asymptote, however it is limited in terms of fully
captured at the laminate level through deformation of individual sub- capturing the non-linear shear deformation observed in testing. As such
laminates. the experimental strain energy absorbed up to γab = 1.2 [3] was ap-
The effects of strain rate were estimated from yarn tests by Russell proximated by using a damage parameter of am4 = −0.8. Comparisons
et al. [23]. As the reference strain rate cannot be altered from unity, the to experimental observations are given in Fig. 7.
unit time of the simulation must be set to a suitable unit (milliseconds). A strain rate parameter for shear stiffness was estimated in the same
Crate1–4 were determined from a logarithmic fit of non-dimensionalised manner as the in-plane response, but with only a limited data set
experimental results which highlight the effect of the Crate1 on in-plane available [24]. The rate effect was larger than in-plane tension as this is
tensile strength (Fig. 4). In-plane axial moduli are controlled through a matrix dominated property, with Crate3 = 0.225 (Fig. 8).
Crate2 and Crate4 in the longitudinal and transverse direction respectively
(Fig. 5).
2.4.3. Through thickness shear
2.4.2. In-plane shear The through thickness shear properties of HB26 are difficult to
The in-plane shear behaviour during ± 45° tensile testing was measure, largely due to the inability to transmit load by the low shear
dominated by fibre re-orientation [24], which for a homogenised [0°/ strength matrix [2]. This mode of laminate deformation is a micro-re-
90°]n cross ply must be included in its response. An initial modulus sponse that may not be able to be captured in a single homogenised
based on experimental results followed by matrix shear damage was element. Here, the MAT162 stiffness has been approximated, based
implemented [24]. At large strains, non-linearities caused by fibre re- upon fibre volume fraction, assuming negligible contribution of the
alignment and scissoring become important [24], requiring invariant matrix and transverse fibres (Table 4). Transverse fibre stiffness (G32f of
90° fibres loaded in shear) was assumed to have negligible contribution
1.25 in Gca for [0°/90°]n homogenised laminates due to the difficulty in
Dimensionless Strength

1.2 loading fibres in this manner. As the element is a homogenised cross

1.15 ply, Gca = Gcb.
1.1 The through thickness matrix shear strength has been set equal to
1.05 the in-plane matrix shear strength as this is a matrix dominated prop-
1 erty. Therefore Sca = Scb = Sab = 1.8 MPa [14,24], with the response
0.95 given in Fig. 9. A coulomb friction angle can be implemented in
0.9 Experimental DataData
Russell et al Exp' [23] MAT162, however φ = 0 was selected due to the low coefficient of
0.85 y = 0.0287ln(x) + 1 friction of Dyneema® HB26. Delamination scale factor S which is used
0.8 to control delamination zone size has been left as default value of 1. A
1 10 100 1000 negative shear damage parameter am4 = −0.8 was also used for cal-
culating shear stress in the through thickness direction after damage.
Strain Rate (1/ms)
Through thickness shear stress also contributes to in-plane fibre failure
Fig. 4. Longitudinal strength determination of strain rate parameter Crate1. (For through the quadratic interaction of stresses. Maximum shear stress was
interpretation of the references to colour in this figure legend, the reader is approximated through principal stress theory based on in-plane tension,
referred to the web version of this article.) which was equal to half of in-plane tensile strength, i.e. 625 MPa.

34
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

Fig. 6. Example of the element coordinate system

(ECS) with no invariant node numbering and with
invariant node numbering [30]. (For interpretation
of the references to colour in this figure legend, the
reader is referred to the web version of this article.)

Fig. 7. The effect of a negative damage parameter on the shear response of a Fig. 9. Through thickness shear response. (For interpretation of the references
single element. (For interpretation of the references to colour in this figure to colour in this figure legend, the reader is referred to the web version of this
legend, the reader is referred to the web version of this article.) article.)

1.6
1.5
Dimensionless Stiffness

1.4
1.3
1.2
1.1
1 Measured Exp'
Experimental Stiffness
Data [24]
0.9 y = 0.225ln(x) + 1
0.8
1 10 Fig. 10. The effect of crush damage parameter. (For interpretation of the re-
Strain Rate (1/ms)
ferences to colour in this figure legend, the reader is referred to the web version
Fig. 8. The effect of strain rate on shear stiffness. (For interpretation of the of this article.)
references to colour in this figure legend, the reader is referred to the web
version of this article.)
2.5. Cohesive interface model

Table 4 Sub-laminates were joined together through cohesive tied contacts

Homogenised through-thickness properties derived from fibre shear prop- using automatic surface to surface tiebreak contact with option 9 to
erties.
represent limited load transfer capability. The sub-laminate cohesive
Property Value Ref. model utilises a mixed mode rate independent fracture mechanics ap-
proach with Mode I and Mode II interaction with a bi-linear traction
G31f (MPa) 1320 [13]
separation law [30]. Due to the low bending stiffness of Dyneema®, a
νf (–) 0.83 [24]
G31 (MPa) 1095.6 – custom experimental methodology is required to accurately capture
Gca (MPa) 547.8 – cohesive properties [31]. Cohesive interface initial stiffness was ap-
Gcb (MPa) 547.8 – proximated by assuming the cohesive interface size was equal to the
sub-laminate height (distance between cohesive interfaces). This was
equivalent to assuming poor matrix load transfer throughout the sub-
2.4.4. Through thickness compression and tensile response laminate (Eq. (11)).
Through thickness stiffness was approximated from plate flyer im-
pact tests by Lässig et al. [2], with Ec = 3620 MPa. Assuming an in- Em
KI =
direct tension mechanism drives failure [26], through-thickness com- tSL (11)
pressive strength was estimated to be equal to cross-ply tensile strength
σc = 1250 MPa. The crush damage parameter, am3 = 20, was selected where KI is Mode I traction stiffness, Em is the Young’s modulus of the
to provide an identical failure response to that of in-plane tension; once matrix material, and tSL is the thickness of the sub-laminate. KII was
fibres have failed little energy absorbing capability remains, with the approximated through using a stiffness ratio of KI/KII = 0.6 which
single element response shown in Fig. 10. provided an adequate stiffness response during short beam shear in-
vestigations based on KI input. A graphical response is given in Fig. 11

35
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

1.5 Mode I 2 Mode II

Tensile Stress (MPa)

Shear Stress (MPa)

1.5
1
GIC 1 GIIC
0.5
0.5

0 0
0 0.5 1 0 0.5 1 1.5
Displacement I (mm) Displacement II (mm)

Fig. 11. Cohesive interface model for Mode I and Mode II cohesive contact with
a sub-laminate thickness of 1 mm, at quasi-static rates. (For interpretation of the
Fig. 13. Representative short beam shear model from experimental tests [5].
references to colour in this figure legend, the reader is referred to the web
(For interpretation of the references to colour in this figure legend, the reader is
version of this article.)
referred to the web version of this article.)

Table 5
beams were fully clamped at one end during the experiment, re-
Quasi-static tied contact cohesive interface model input card, 1 mm sub-lami-
presented here by fully fixed boundary conditions.
nate thickness.
Property Symbol Value Ref.

**
3.1. Short beam shear
Mode I tensile strength (MPa) T 1.2 [2]
Mode I initial stiffness* (N/mm3) KI 60 –
Mode I Fracture Energy (N/mm) GIc 0.544 [32] Short beam shear models with a relatively low L/h and an initial
Mode I traction max’ displacement (mm) DIF 0.908 – sub-laminate thickness, tSL, of 0.25 mm (Fig. 13) were investigated first.
Initial stiffness ratio (–) KII/KI 0.6 – The number of elements through the thickness of each sub-laminate
Mode II Shear Strength** (MPa) S 1.8 [14,24] was varied as 1, 2 and 3 elements (Fig. 14). Initial stiffness was identical
Mode II fracture energy (N/mm) GIIc 1.088
Mode II traction max’ displacement (mm) DIIF 1.21 –
to that observed experimentally, however, after reaching the maximum
Exponent of Powerlaw (–) PARAM 1 [30] Mode II cohesive interface strength multiple elements per layer were
too stiff. This was due to the build-up of bending stress, in contrast a
* Initial stiffness dependency upon element size. single element layer effectively provides a direct bending stress of zero
** Quasi-static strength. due to its single integration point. Normally multiple elements through
each sub-laminate thickness would be required, however the inability
and the full cohesive interface contact card given in Table 5. of Dyneema® fibres to take compressive loads and instead micro-buckle
limits the capability for bending stress to build up in the microstructure.
Therefore, the inability of single element layers to capture bending
3. Beam deformation investigations
stress distribution provides a reasonable approximation to the short
beam post-cohesive shear damage response. Under loading this pro-
Cantilever loaded beams were modelled to investigate the effective
moted buckling of sub-laminates and limited shear stress transfer
bending stiffness of HB26 as a means of validating the deflection fol-
within the sub-laminate (Fig. 15), comparable to that observed during
lowing impact. The modelling approach was also investigated by
experimental observations [5]. One element thick sub-laminates were
varying the number of elements per sub-laminate as well as the thick-
therefore used for all modelling techniques following this investigation.
ness of each sub-laminate. Two experiments performed by Liu et al. [5]
The inability to carry further transverse load was triggered by da-
were used for comparison, with the experimental setup given in Fig. 12.
mage to the cohesive interface, controlled by the Mode II inter-laminar
The width of all beams was 20 mm, the loading cylinder diameter d was
response. Prior to reaching Mode II shear strength of the cohesive in-
6.4 mm at a distance s of 10 mm from the free end and the total lami-
terface (Fig. 16a) the response was controlled by a combination of the
nate height h was 6 mm. Loading was performed with a constant dis-
shear stiffness of the cohesive interface and the shear stiffness in the
placement rate of 2 mm/min, for beam lengths L of 10 mm and 100 mm,
sub-laminate. Nodes connected by the cohesive interface develop dis-
giving a short beam shear test and a long beam bending test respec-
placements relative to one another in the elastic region of the cohesive
tively. All modelling used one solid element with a single integration
contact rather than pure shear of elements themselves (Fig. 16b). The
point across the width, with plane strain boundary conditions. For
lack of load transfer between sub-laminates was representative of the
beam studies, rate effects were set to zero and the displacement rate
was increased to 0.05 mm/ms to reduce computational expense. All 30
1 Element per SL
3 Elements per SL
25
Load per unit width (kN/m)

2 Elements per SL

20
2 Elements per SL
15
Single Element per SL
10

3 Elements per SL 5 Experimental Data [5]

0
0 1 2 3 4 5
Deflection (mm)

Fig. 12. Experimental setup for beam studies performed by Liu et al [5]. (For Fig. 14. Effect of short beam shear response on the number of elements per
interpretation of the references to colour in this figure legend, the reader is sub–laminate (SL). (For interpretation of the references to colour in this figure
referred to the web version of this article.) legend, the reader is referred to the web version of this article.)

36
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

Fig. 15. Comparison of axial stress (MPa) and through thickness shear stress (MPa) build up in sub-laminates for 1 element per sub-laminate (SL) and 3 elements per
sub-laminate at 2 mm deflection. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

lack of load transfer between fibres, with fibres having a higher shear below 0.125 mm. Smaller sub–laminate thicknesses, approaching a
stiffness compared with the matrix. After reaching maximum Mode II single UD layer thickness, were beyond the scope of homogenised cross-
cohesive strength, gradual failure of the cohesive interface can ensue, ply investigations.
with total failure of some cohesive interfaces occurring just past 2 mm
vertical displacement, allowing sub-laminates to separate (Fig. 16c). 3.2. Long beam buckling
The initialisation of reduced load bearing capability was investigated
by varying mode II cohesive interface strength. Increased mode II shear Longer beam bending experimental results (L = 100 mm, Fig. 12)
strength postponed inter-laminar damage and increased the load re- [5] were modelled with the same parameter set, maintaining a sub-
quired to delaminate the beam, with the effect on load deflection curves laminate thickness of 0.25 mm. As in experimental observations, the
indicated in Fig. 17. load bearing capability was controlled by kink-band formation at the
Variation of sub-laminate thickness, and hence number of sub-la- root of the beam (Fig. 19). The size and shape of the kink-band was
minates in a given beam, with different cohesive interface stiffness was larger in modelling than in experimental observations, thought to be
investigated. This showed that with increasing sub-laminate thickness a due to the inability to capture micro-scale phenomena such as micro-
proportional reduction in cohesive interface stiffness was required to buckling of fibres and individual lamina deformation. This effectively
provide the same beam initial stiffness (Fig. 18). The sub-laminate led to an increase in the load required to buckle the beam at the root
thickness also influenced the post cohesive damage response, reducing (Fig. 20). Numerical noise post-buckle was induced due to the relatively
load bearing capability when the sub–laminate thickness was reduced coarse element size compared with that required to capture this micro-

Fig. 16. Finite element modelling of short beam shear deformation, (a) as displacement progresses, (b) cohesive elastic deformation, and (c) cohesive interface
failure. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

37
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

14 10
Cohesive Interface Experimental Data [5]
12 Damage Onset of Cohesive

Load per unit width (kN/m)

Model
Load per unit width (kN/m)
8
Interface Failure Model with Waviness
10
6
8
4
6

4 Experimental Data [5] 2

Mode II Strength 1.6 MPa
2 Mode II Strength 1.8 MPa 0
Mode II Strength 2.0 MPa 0 5 10 15
0 Tip Deflection (mm)
0 1 2 3 4 5
Deflection (mm) Fig. 20. Longer beam bending load-displacement response, and the effect of
introducing ply waviness. (For interpretation of the references to colour in this
Fig. 17. The effect of Mode II cohesive shear strength on short beam shear figure legend, the reader is referred to the web version of this article.)
response. (For interpretation of the references to colour in this figure legend,
the reader is referred to the web version of this article.)

10
Load per unit width (kN/m)

5
Liu et al Exp' Data [5]
Experimental
TSL
TSL ==0.25
0.25mm,
mm,KKII ==240
240
TSL
TSL ==0.5
0.5mm,
mm,KKI
I ==240240
TSL
TSL ==0.5
0.5mm,
mm,KKI
I ==120120
TSL
TSL ==0.125
0.125mm,
mm,KKI I ==240
240 Fig. 21. Sinusoidal ply waviness approximation overlaid with optical micro-
TSL ==0.125
TSL 0.125mm,
mm,KKI I ==480
480 scopy [24]. (For interpretation of the references to colour in this figure legend,
0
0 1 2 3 4 5 the reader is referred to the web version of this article.)
Deflection (mm)

Fig. 18. The effect of varying sub-laminate thickness (TSL) and cohesive stiff- laminate cross section (Fig. 21). Incorporating ply waviness reduced the
ness (KI related to KII by a factor of 0.6). (For interpretation of the references to initial buckling load by approximately half, bringing it closer to that
colour in this figure legend, the reader is referred to the web version of this observed in experiments [5] and reduced the magnitude of oscillations
article.) in the signal following buckling (Fig. 20).
Although the post buckle response was similar, the beam still
driven phenomenon. However, the post-buckle average load does pla- buckled at a higher a load compared with experimental results. The
teau to a similar value to that from the experiments. effect of micro-buckling of the fibres was difficult to capture at this
One conclusion from experimental observations was that ply wavi- scale, and needs to be investigated further to accurately capture this
ness may be a contributing factor, helping to initiate the buckling re- phenomenon. It was thought that this will have little effect on the
sponse through out of plane fibre variations observed in the micro- modelling of ballistic impact on UHMWPE fibre laminates for validation
structure. To replicate this, waviness was input into the element co- purposes as the projectile has come to rest prior to the transverse shear
ordinate system by altering the local material axes in the element wave reaching the boundary (boundary effects have limited effect on
keyword card. The degree of ply waviness was approximated from the ballistic performance for a large enough plate).
microstructural investigations [24], assuming a sinusoidal wave input
in the form: 4. Drop weight impact modelling
y (x ) = Asin(ωx ) (12)
Drop weight impact testing, without failure of the fibres, was in-
where A is the amplitude of the sin wave, and ω is the frequency of vestigated by Hazzard et al. [6]. These results were used to compare the
oscillation, equal to 0.14 mm and 0.93 rad/mm respectively. The phy- in-plane shear response to that predicted by the finite element model.
sical representation is similar to an optical microscope image of a HB26 The experiment impacted HB26 laminates 32 UD plies thick (2.2 mm

Fig. 19. Root buckle of long beam bending. (a) Modelling at deflection δ = 10 mm and (b) experimental observations from Liu et al [5]. (For interpretation of the
references to colour in this figure legend, the reader is referred to the web version of this article.)

38
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

a) Fig. 22. (a) Isometric view of drop weight impact

model, (b) cross-section highlighting clamp con-
ditions, (c) cross section highlighting mesh den-
sity through the thickness, and (c) element size
utilised within the laminate. (For interpretation
of the references to colour in this figure legend,
the reader is referred to the web version of this
article.)

b) c) d)
Elem’ L = 1 mm

TSL = 0.1375 mm

12
Experimental Data [6]
HB26 Experiment
10
Model
Contact Force (kN)

0
0 10 20 30 40 50
Displacement (mm)
Fig. 25. Projectile and target mesh. (For interpretation of the references to
Fig. 23. Force-deflection curve comparison from HB26 drop weight impact test colour in this figure legend, the reader is referred to the web version of this
and FEA model. (For interpretation of the references to colour in this figure article.)
legend, the reader is referred to the web version of this article.)

of the impactor material selected to simulate a mass of 26.3 kg, totalling

thickness) with a 150 J blunt hemispherical impactor. Purely frictional 150 J impact energy. Clamps were modelled with a single element
clamping of the 200 mm2 laminates around a 125 mm diameter aper- through the thickness (1 mm thick), with both the impactor and clamps
ture constrained the laminate, replicated in LS-DYNA by modelling the using a rigid material, as deformation of the impactor and clamp did not
entire laminate, clamp, and blunt hemispherical impactor (Fig. 22a). An occur. The lower clamp was fully constrained whilst the top clamp had
initial velocity was given to the impactor of −3.377 m/s, with a density an equivalent pressure load on its top surface to provide 2 MPa pressure

Fig. 24. HB26 drop weight impact in-plane shear strain contour map comparison of experimental results recorded via DIC compared with model results with varying
time after impact. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

39
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

Table 6 and matched by that of the impactor, and an eroding surface contact
MAT 98 material model for 20 mm diameter FSP [3]. was used between in the impactor and laminate in case of element
Property Symbol Value Units failure (none expected from experiments). Identical material properties
within each sub-laminate were used as in the beam test modelling, with
Density ρ 0.0078 g/mm3 rate effects now turned on (Table 2).
Young’s modulus E 207000 MPa
Poisson’s ratio ν 0.33 –
Yield stress A 1030 MPa
4.1. Drop-weight model results
Hardening constant B 477 MPa
Hardening exponent N 0.18 – The predicted force-deflection results from the finite element ana-
Strain rate constant C 0.012 – lysis were similar to experimental results (Fig. 23). Maximum BFD for
Effective plastic strain at failure ƐPF 1E17 MPa
the [0°/90°]16 laminate was almost identical; however the maximum
Maximum stress before rate effects σx 1E20 MPa
Saturation Stress σSAT 1E20 MPa contact force was approximately 10% larger in the model. Energy ab-
Reference strain rate ε 0̇ 1 s−1 sorbed was also within 10%, however a larger rebound is observed
within the model as gravity has not been applied. Initially the contact
force was lower, thought to be caused through excess in-plane shear
Table 7 strain of the model compared with that of the experiment, visualised by
DOP mesh study input parameters at 500 m/s impact velocity of a 50 mm la- a full field strain comparison with time in Fig. 24. Although strains were
minate target. of the same magnitude, locally higher in-plane shear strains were ob-
Impact zone No. Sub- KI* (N/ Outer DOP % Diff’ from served around the boundary and close to the point of impact. It is
element size laminates mm3) element (mm) next largest thought that the inability to fully capture the in-plane shear curve can
(mm) aspect ratio mesh
lead to locally higher shear strains, whilst the relatively simple
5 10 12 1 5 – boundary model does not capture the potential variation in coefficient
3.33 15 18 1.5 10 50 of friction in static and dynamic loading with the dependence on
2 25 30 2.5 14 28.6 boundary pressure. Modelling was also able to show that the largest in-
1 50 60 5 16 12.5 plane shear deformations tended to occur under the boundary itself,
something that couldn’t be captured by digital image correlation during
* From Eq. (11).
experiments. Larger scale wrinkling was also visible in the model
-500 (Fig. 24), thought to be due to the relatively large in-plane element
1 mm
dimensions compared with the length scale of the wrinkles observed
Projectile Velocity (m/s)

-400 during testing. The out-of-plane curvatures in experiments with a scale

2 mm
3.33 mm of less than 1 mm in length could therefore not be captured. Instead the
-300
5 mm wrinkles form at a larger scale due to the larger elements and a higher
-200 compressive strength/stiffness within the sub-laminates. The formation
of the wrinkles was thought to limit shear deformation in certain areas,
-100
effectively increasing laminate stiffness at greater deflection values and
0 increasing contact force (Fig. 23).

100 5. Ballistic impact modelling

0.00 0.10 0.20 0.30 0.40 0.50
Time (ms)
Following on from low rate modelling, ballistic impact modelling of
Fig. 26. Projectile velocity time trace for varying mesh densities during DOP Dyneema® was compared with open literature experimental data.
investigation. (For interpretation of the references to colour in this figure le- Nguyen et al. [7] impacted target laminates with in-plane dimensions of
gend, the reader is referred to the web version of this article.) 300 × 300 mm with varying thickness. The projectile used was a
20 mm diameter fragment simulating projectile (FSP) at a range of
Table 8 impact velocities, providing a ballistic limit approximation through a
Mesh sensitivity of VBL for a 10 mm laminate target. Lambert-Jonas approximation [33]. Laminates were modelled with an
identical approach as previously described, with no boundary con-
Impact zone No. Sub- KI (N/ Outer Approximate % Diff’
element size laminates mm3) element VBL (m/s) from next straints surrounding the laminate. Although the boundary was clamped
(mm) aspect ratio largest during the experiment, the effect of the clamping condition was thought
mesh to be minimal as the plate was large enough for the projectile to come
to rest prior to the transverse shear wave reaching the boundary com-
2 5 30 2.5 475 –
1.66 6 36 3 420 13.1
bined with the low coefficient of friction of Dyneema®. A biased mesh
1.25 8 48 4 409 2.7 was incorporated, increasing mesh density towards the impact zone
1 10 60 5 399 2.5 (central 20 × 20 mm), with the elements in this zone being cubic to
0.66 15 90 5 361 10.5 provide the optimum aspect ratio for failure analysis (Fig. 25). In all
ballistic analyses an identical material and contact card to previous
models was used, as in Tables 2 and 5, with one exception; shear
over the laminate contact area. This was equivalent to controlling the
strength in the material model as well as the cohesive interface strength
clamping pressure experimentally through torque tightening of sur-
was increased from 1.6 to 2.6 MPa, equivalent to that found at higher
rounding bolts during the experiment. Contact was applied between the
rate experiments, similar to those observed during ballistic impact [14].
clamps and the laminate with a static and dynamic coefficient of fric-
During analyses, it was found that there was a contact issue re-
tion of 0.1 (Fig. 22b). The laminate itself was discretised into sub-la-
garding interpenetration of sub-laminates resulting in large negative
minates with thickness of 0.1375 mm, equivalent to a single [0°/90°]
contact energies, thought to be caused by a rate effect with the surface
cross-ply of Dyneema® (Fig. 22c, d). Cohesive interfaces, with stiffness
contact algorithm related to low stiffness properties of the cohesive
adjusted according to Eq. (11), were inserted between each layer of sub-
formulation. Previous studies have found a similar issue, and modelled
laminates. The in–plane laminate element dimensions were 1×1 mm
laminates with gaps in-between sub-laminates to prevent inter-part

40
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

550
500
450
2 mm Model Data
400 2 mm (0.88, 6.12, 475.3)
350 1.66 mm Model Data
VR (m/s) 1.66 mm (0.91, 3.65, 419.8)
300 1.25 mm Model Data
250 1.25 mm (0.92, 3.98, 409.3)
1 mm Model Data
200
1 mm (0.91, 4.52, 398.7)
150 0.66 mm Model Data
100 0.66 mm (0.92, 4.66, 361.4)

50
0
350 400 450 500 550 600
VI (m/s)

Fig. 27. Numerical residual velocity prediction with Lambert-Jonas VBL fit with parameters given in the legend (a, p, VBL). (For interpretation of the references to
colour in this figure legend, the reader is referred to the web version of this article.)

Table 9 MAT98 in LS-DYNA, with identical parameters to that used by Nguyen

Mesh sensitivity of BFD for a 10 mm laminate target impacted at 300 m/s. et al. [3], given in Table 6.
Impact No. Sub- KI (N/ Outer Maximum % Diff’ DOP
zone laminates mm3) element BFD (mm) from next (mm)
element aspect largest 5.2. Mesh refinement study
size (mm) ratio mesh
A mesh refinement study was performed to determine the de-
2 5 30 2.5 42.87 – 0
1.66 6 36 3 42.31 0.1 0 pendency of the element geometry on the ballistic performance of the
1.25 8 48 4 46.13 8.3 2.5 laminate. Here, depth of penetration (DOP) has been investigated on
1 10 60 5 50.12 8.0 3 50 mm thick target laminates, followed by the ballistic limit (VBL) and
0.66 15 90 5 52.91 5.3 4 BFD of thinner 10 mm target laminates. Meshes produced in all studies
maintained an element aspect ratio of unity in the impact zone. During
beam studies, it was highlighted that single element sub-laminates more
penetration [3]. Here an additional automatic surface to surface contact
accurately represent the near zero bending stiffness of Dyneema®. This
with default options and an updated contact algorithm was applied to
meant that increasing the mesh density also increased the number of
maintain a gapless approach between sub-laminates.
sub-laminates within the model. The cohesive interface has therefore
had its initial stiffness KI updated in accordance with Eq. (11). The
5.1. Projectile model biased mesh with increasing element size approaching the edge of the
laminate had a maximum in-plane dimension of 5 mm and a maximum
The FSP had a mesh size equal to that of the target laminate impact aspect ratio of 5.
zone to avoid any contact mismatch through varying mesh density. The Firstly, DOP studies were performed on 50 mm thick laminates at an
FSP was represented using the simplified Johnson-Cook material model impact speed of 500 m/s. Cubic sized central elements varying from

1200
1100 t = 5 mm Model Data
1000 t = 5 mm (0.95, 6.83, 338.9)
t = 10 mm Model Data
900
t = 10 mm (0.91, 4.52, 398.7)
800 t = 15 mm Model Data
700 t = 15 mm (0.92, 3.54, 499.1)
VR (m/s)

t = 20 mm Model Data
600
t = 20 mm (0.90, 3.80, 634.7)
500 t = 30 mm Model Data
400 t = 30 mm (0.80, 6.03, 787.9)
300 t = 40 mm Model Data
t = 40 mm (0.78, 5.12, 889.0)
200
t = 50 mm Model Data
100 t = 50 mm (0.75, 5.40, 990.3)
0
300 400 500 600 700 800 900 1000 1100 1200 1300
VI (m/s)

Fig. 28. Lambert-Jonas ballistic approximation from model impacts with varying thickness t, with curve fit parameters given in the legend (a, p, VBL). (For inter-
pretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

41
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

1600

1400 ADt .Ap /mp = 0.08

1200 Lee Exp Spectra 1000 5.58 mm FSP 1994 [10]

Lee Exp Fit 1994 [10]
1000 Heisserer HB26 5.46 mm FSP 2012 [9]

V50 (m/s)
Heisserer HB26 20 mm FSP 2012 [9]
800 Nguyen HB26 20 mm FSP 2015 [7]
Nguyen HB26 12.7 mm FSP 2015 [7]
600 Nguyen Analytical Model 2015 [7]
Phoenix & Porwal Membrane Model 2003 [4]
400 Current Numerical Model

200

0
0 0.1 0.2 0.3 0.4 0.5
ADt.Ap/mp (-)

Fig. 29. Experimental, analytical, and numerical ballistic limit compared in terms of non-dimensioned areal density. (For interpretation of the references to colour in
this figure legend, the reader is referred to the web version of this article.)

Fig. 30. (a) Schematic of apex deflection and hinge position, progression with time for (b, c) 10 mm target at 365 m/s impact velocity and (d, e) 20 mm target at
615 m/s impact velocity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

1 mm to 5 mm resulted in the sub-laminate thickness being the same as where VI is the initial velocity, and a and p are selected through a least
the in-plane element dimensions to maintain an aspect ratio of unity, square curve fit with VBL estimated at VR = 0. Mesh parameters for the
see Table 7. A relatively thick laminate was chosen to investigate the thinner target tests are provided in Table 8. The investigation showed
progressive failure regime prior to large membrane deflection, with that reducing mesh size lowered VBL, with the smallest mesh sizes of
DOP defined as the thickness of the number of layers that failed. DOP 1 mm and 0.66 mm within ± 10% of one another. This was not a full
was found to increase with increasing mesh density, partly due to large convergence, which highlights that bulge failure at the rear of the la-
elements giving limited resolution in terms of failure of the sub-lami- minate was more susceptible to mesh variation than progressive failure
nate thickness. Convergence of the mesh was found at approximately in DOP tests. Due to the limitations in runtime, computational cost, and
1 mm in size with DOP being within 2 mm of the next largest element the physical mesh size approaching that of a UD model, smaller mesh
size (Table 7). The change in velocity of the projectile with time also sizes were not investigated. The Lambert-Jonas least square fit of results
highlights the mesh convergence at higher mesh densities (Fig. 26). The is also given in Fig. 27, highlighting the variation induced in VBL by
largest mesh produced a stiffer response, decelerating the projectile at a changes in mesh geometry.
faster rate with a lower DOP. 10 mm thick laminates were also impacted with a 20 mm FSP at
Mesh dependency of VBL and BFD of thinner laminates was also 300 m/s to investigate the unperforated response of thinner laminates.
investigated with the 20 mm FSP. All laminates were 10 mm thick, It was found that reducing the mesh size induced earlier failure of the
whilst impact velocity was varied to provide a minimum of 3 residual target strike face (Table 9), which is consistent with previous DOP
velocity (VR) values following perforation to estimate VBL using the studies in the progressive failure regime (Table 7). This increased the
Lambert-Jonas approximation in Eq. (13). BFD of the remaining laminate thickness deforming in a bulging
manner, and may play a significant role in the variation of VBL of
VR = a (VIp−VBL
p 1/ p
) (13) thinner laminates as additional failed layers account for ≈10% of the

42
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

Bottom View
t = 10 mm, V = 365 m/s t = 50 mm, V = 900 m/s

Side View

Cut Isometric View

Projectile Effective Plastic Strain

Fig. 31. Model imagery of partially perforated 10 and 50 mm thick laminates at 0.5 ms after impact at 365 and 900 m/s respectively. (For interpretation of the
references to colour in this figure legend, the reader is referred to the web version of this article.)

laminate thickness, roughly equal to the reduction in ballistic limit performed and compared with previous experimental results and ana-
(Table 8). lytical methods, with VBL approximations from the model provided in
Following the mesh studies, a 1 mm central mesh with a bias to Fig. 28. Initial results gave good correlation between experimental re-
5 mm at the outer edge was selected for use in all further analyses. This sults and the model, with increasing ballistic limit with increasing
appeared to show mesh convergence during DOP studies, whilst VBL target thickness. The different mechanisms of energy absorption were
and BFD are within 10% of surrounding mesh sizes, whilst achieving also captured, bridging the gap between analytical models at non-di-
reasonable computational run times. Typical runtimes for a 10 mm mensional areal densities ≈0.1. The model does however produce
target with a 1 mm central mesh size were approximately 4–6 h de- lower than expected VBL at higher non-dimensional areal densities
pendent on impact velocity. (ADt.Ap/mp) when compared with experimental results from literature,
shown in Fig. 29, where ADt is the areal density of the target, Ap is the
projected area of the projectile, and mp is the mass of the projectile.
5.3. Model-experimental comparison
Lassig et al. [17] incorporated the Mie-Gruneisen equation of state into
a ballistic impact model for Dyneema® and its effect on VBL at higher
Investigation at a range of non-dimensional areal densities was

43
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

Fig. 32. Projectile contact force against a 50 mm target laminate at 900 m/s impact speed with time, highlighting the arrival time of the reflected pressure relief
wave. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

impact speeds following its inclusion. At higher non-dimensional areal relief wave is only one factor in controlling the transition from pro-
densities, the impact velocity required to perforate the target rises and gressive local failure to bulge deformation. In the 50 mm target model,
shock effects become increasingly important. the projectile force at which the mode switch was observed is ap-
proximately 90 kN (≈286 MPa pressure assuming the FSP projected
circumferential area) and at a later time than shown in Fig. 32. Fur-
5.4. Deformation and failure mechanisms of laminates thermore, some sub-laminates ahead of the one in direct contact with
the projectile were also observed to be damaged by the high-pressure
BFD and shear hinge progression was compared with experimental front prior to tensile relief wave interaction (Fig. 32). This indicates
results from Nguyen et al. [3], showing good correlation to experi- that some damage may be imparted to the laminate by the compressive
mental results for a 10 mm target thickness for both BFD and shear shockwave, prior to projectile interaction in that part of the laminate.
hinge expansion (Fig. 30b, c). At a greater target thickness of 20 mm
and higher impact velocity, a larger BFD was predicted by the model
(Fig. 30d), supporting previous observations of over perforation at 6. Conclusions
larger target thicknesses and higher impact velocities. A third result at
36 mm thickness and 888 m/s impact velocity showed complete per- A finite element modelling approach within LS-DYNA has been
foration, whilst experimental results showed an unperforated response proposed to capture quasi-static stiffness behaviour, drop-weight im-
[3]. Thicker laminates (≥20 mm) impacted at higher velocities pact behaviour, and finally ballistic behaviour of Dyneema® compo-
(≥600 m/s) have a larger portion of the laminate failing in a locally sites. Beam studies highlighted the importance of the cohesive interface
progressive manner around the projectile, with thinner laminates properties in limiting load transfer between sub-laminates and de-
having a higher proportion of the laminate in bulging membrane action termining the stiffness response. Longer beam performance was de-
(Fig. 31). This suggests a lower energy absorption during localised termined by buckling at the root, with initiation load found to be de-
failure compared with experimental results. Large amounts of in-plane pendent on fibre waviness. Drop weight impact modelling provided
shear pull in was visible during bulging of the rear surface on all models BFD and impactor contact force within 10% of that observed experi-
during bulge deformation. This helped form the pyramid BFD shape mentally, with discrepancies thought to be caused by the difficulty in
observed in [0°/90°]n cross-ply ballistic impacts due to large anisotropy accurately capturing boundary conditions and the inability to capture
in the cross-ply [7]. In the transition zone from progressive local failure microscale wrinkling phenomena in a larger continuum model. The
on the front face to bulging membrane action, large amounts of dela- magnitude and distribution of in–plane shear strain contour maps was
mination was also visible, similar to that observed during experiments similar between experimental digital image correlation and model
[14]. Minimal deformation of thicker laminates on the front surface was output. Small differences arise as the complex in-plane shear response is
observed surrounding the impact zone, however thinner laminates not fully captured by the material model, and may be an area for future
tended to delaminate on the top surface following the direction of fibres improvement. Ballistic performance in terms of VBL and BFD was shown
that were in contact with the projectile. At the largest target thickness to be accurate at areal densities less than 0.15 ADtAp/mp and impact
investigated (50 mm shown in Fig. 31), increased projectile deforma- speeds less than 600 m/s. At lower target thickness and lower impact
tion was also visible, with blunting of sharp edges of the 20 mm FSP speeds, membrane and bulge deformation dominated the response. At
observed at an impact speed of 900 m/s, similar to that observed during higher areal densities, the ballistic limit from the model was lower than
experiments [3]. experimental results, thought to be due to the increasing importance of
An investigation of the mechanism triggering change in failure an equation of state, currently absent within this model. The mechan-
mode from local progressive failure where effective fibre failure occurs isms of failure were similar to that observed in ballistic impact, with
to bulge and membrane action was performed using a 900 m/s FSP large scale delamination, shear pull-in, and fibre failure all visible
impact on a 50 mm thick laminate. It was observed that the reflective within the model. Mode switch between local progressive failure and
tensile relief wave [34] arrived at the progressively penetrating pro- bulge deformation was thought to be caused by reaching a critical
jectile front in under 40 µs, substantially reducing contact force on the contact force between the projectile and the laminate, below which the
projectile from the laminate (Fig. 32). However, the projectile con- laminate can behave in a bulge deformation.
tinues to fail plies in a progressive manner, indicating that the tensile

44
M.K. Hazzard et al. Composites Part A 115 (2018) 31–45

Acknowledgements shear strength on the ballistic response of laminated composite plates. Eur J Mech A
Solids 2013;42:35–53. https://doi.org/10.1016/j.euromechsol.2013.04.002.
[15] Grujicic M, Glomski PS, He T, Arakere G, Bell WC, Cheeseman BA. Material mod-
The authors acknowledge the support of the EPSRC through the eling and ballistic-resistance analysis of armor-grade composites reinforced with
ACCIS Centre for Doctoral Training grant, EP/G036772/1, and the fi- high-performance fibers. J Mater Eng Perform 2009;18(9):1169–82.
nancial contribution of the Defence Science and Technology Laboratory [16] Iannucci L, Pope D, Dalzell M, Dstl P, Salisbury W. A constitutive model for
Dyneema UD composites. In: iccm-central.org; 2009.
(dstl) for funding for this research. We would also like to thank, [17] Nguyen L, Lässig T, Ryan S. Numerical modelling of ultra-high molecular weight
Professor Lorenzo Iannucci and Professor Paul Curtis for helpful back- polyethylene composite under impact loading. In: the 13th Hypervelocity Impact
ground discussions, Assistant Professor Bazle Gama for discussions Symposium; 2015, p. 436–43.
[18] Materials Sciences Corporation, 135 Rock Rd, Horsham, Pennsylvania, 19044.
surrounding MAT162, and Dr Galal Mohamed for general advice on [19] Hashin Z. Failure criteria for unidirectional fiber composites. J Appl Mech
numerical modelling in LS-Dyna. Dyneema® is a trademark of DSM. 1980;47(2):329–34.
[20] B.A. Gama. User manual: a progressive composite damage model for unidirectional
and woven fabric, no. 215. Materials Sciences Corporation (MSC) & University of
References
Delaware Center for Composite Materials (UD-CCM); 2015.
[21] Xiao J, Gama B, Gillespie J. Progressive damage and delamination in plain weave S-
[1] Cunniff PM. Dimensionless parameters for optimization of textile-based body armor 2 glass/SC-15 composites under quasi-static punch-shear loading. Compos Struct
systems. In: Proceedings of the 18th International Symposium on Ballistics; 1999, p. 2007.
1303–10. [22] Belytschko T, Bindeman L. Assumed strain stabilization of the eight node hexahe-
[2] Lässig T, Nguyen L, May M, Riedel W, Heisserer U, van der Werff H, et al. A non- dral element. Comput Methods Appl Mech Eng 1993;105:225–60.
linear orthotropic hydrocode model for ultra-high molecular weight polyethylene in [23] Russell BP, Karthikeyan K, Deshpande VS, Fleck NA. The high strain rate response
impact simulations. Int J Impact Eng 2015;75:110–22. of Ultra High Molecular-weight Polyethylene: from fibre to laminate. Int J Impact
[3] Nguyen L, Lässig T, Ryan S, Riedel W. A methodology for hydrocode analysis of Eng 2013;60:1–9. https://doi.org/10.1016/j.ijimpeng.2013.03.010.
ultra-high molecular weight polyethylene composite under ballistic impact. Compos [24] Hazzard MK, Curtis PT, Iannucci L, Hallett S, Trask R. An Investigation of the in-
A Appl Sci Manuf 2016;84:224–35. plane performance of ultra-high molecular weight polyethylene fibre composites.
[4] Phoenix SL, Porwal P. A new membrane model for the ballistic impact response and In: International Conference of Composite Materials, Copenhagen; 2015, p. 19–24.
V50 performance of multi-ply fibrous systems. Int J Solids Struct [25] Heisserer U, van der Werff H. Strength matters: which strength of Dyneema fiber
2003;40(24):6723–65. composites to use in hydrocode models? A Discussion. In: Internation Ballistics
[5] Liu G, Thouless M. Collapse of a composite beam made from ultra high molecular- Symposium, Edinburgh; 2016, no. Iso 527, p. 1–5.
weight polyethylene fibres. J Mech Phys Solids 2014;63:320–35. [26] Attwood JP, Khaderi SN, Karthikeyan K, Fleck NA, Omasta MR, Wadley HNG, et al.
[6] Hazzard MK, Hallett S, Curtis PT, Iannucci L, Trask RS. Effect of fibre orientation on The out-of-plane compressive response of Dyneema ® composites. J Mech Phys
the low velocity impact response of thin Dyneema® composite laminates. Int J Solids 2014;70:200–26. https://doi.org/10.1016/j.jmps.2014.05.017.
Impact Eng 2016;100:35–45. https://doi.org/10.1016/j.ijimpeng.2016.10.007. [27] Gama B, Gillespie JW. Punch shear based penetration model of ballistic impact of
[7] Nguyen L, Ryan S, Cimpoeru S. The effect of target thickness on the ballistic per- thick-section composites. Compos Struct 2008.
formance of ultra high molecular weight polyethylene composite. Int J Impact Eng [28] Levi-Sasson A, Meshi I, Mustacchi S, Amarilio I, Benes D, Favorsky V, et al.
2015;75:174–83. Experimental determination of linear and nonlinear mechanical properties of la-
[8] Ćwik T, Iannucci L, Curtis P, Pope D. Investigation of the ballistic performance of minated soft composite material system. Compos B Eng 2014;57:96–104. https://
ultra high molecular weight polyethylene composite panels. Compos Struct doi.org/10.1016/j.compositesb.2013.09.043.
2016;149:197–212. [29] Attwood J, Fleck N. The compressive response of ultra-high molecular weight
[9] Heisserer U, van der Werff H. The relation between Dyneema® fiber properties and polyethylene fibres and composites. Int J Solids Struct 2015;71:141–55.
ballistic protection performance of its fiber composites. In: 15th International [30] Hallquist J. LS-DYNA keyword user’s manual. Livermore Software Technology
Conference on Deformation, Yield and Fracture of Polymers; 2012, vol. 3, p. Corporation; 2007.
242–46. [31] Lässig T, Nolte F, Riedel W, May M. An assessment of experimental techniques for
[10] Lee BL, Song JW, Ward JE. Failure of spectra® polyethylene fiber-reinforced com- measuring the mode I fracture toughness of UHMW-PE composites, In: Proceedings
posites under ballistic impact loading. J Compos Mater 1994;28(13):1202–26. of the 17th European Conference on Composite Materials ECCM17; 2016.
[11] Liu G. Modelling microbuckling failure of a composite cantilever beam made from [32] Grujicic M, Arakere G, He T. Multi-scale ballistic material modeling of cross-plied
ultra high molecular-weight polyethylene fibres. Acta Mech 2015;226(4):1255–66. compliant composites. Compos B Eng 2009;40(6):468–82.
[12] Nazarian O, Zok F. Constitutive model for the shear response of Dyneema® fiber [33] Lambert J, Jonas G. Towards standardization in terminal ballistics testing: velocity
composites. Compos A Appl Sci Manuf 2014;66:73–81. representation, no. Technical report AD-A021 389. Ballistic Research Laboratories
[13] Chocron S, Nicholls AE, Brill A, Malka A, Namir T, Havazelet D, et al. Modeling Aberdeen Proving Ground, MD; 1976.
unidirectional composites by bundling fibers into strips with experimental de- [34] Nguyen L, Ryan S, Orifici AC. A numerical investigation on the response of thick
termination of shear and compression properties at high pressures. Compos Sci ultra-high molecular weight polyethylene composite to ballistic impact. In:
Technol 2014;101:32–40. https://doi.org/10.1016/j.compscitech.2014.06.016. Proceedings - 29th International Symposium on Ballistics, BALLISTICS 2016, vol. 2;
[14] Karthikeyan K, Russell BP, Fleck NA, Wadley HNG, Deshpande VS. The effect of 2016.

45