International Journal of Impact Engineering 178 (2023) 104624

Contents lists available at ScienceDirect

International Journal of Impact Engineering

journal homepage: www.elsevier.com/locate/ijimpeng

Rigid-plastic membrane response of thin plates under impulsive blast loads

using the extended Hamilton principle
Saud A.E. Alotaibi a,b ,∗,1 , Samuel E. Rigby a , Maurizio Guadagnini a , Andrew Tyas a
a
Department of Civil & Structural Engineering, University of Sheffield, Mappin Street, Sheffield, S1 3JD, UK
b
Department of Civil Engineering, Qassim University, Buraydah, Saudi Arabia

ARTICLE INFO ABSTRACT

Keywords: The response of plates subjected to blast loads is of considerable scientific interest. The loading imparted to
Near-field blast a structure following a close-in detonation of a high explosive is typically high in magnitude, near-impulsive,
Impulsive load spatially non-uniform, with localised variability and a high dependence on factors such as charge shape,
Rigid-perfectly plastic
position, and composition. The resulting structural response may induce large displacements in materials
Analytical solution
whose properties may not be fully characterised. In order to properly account for the effects of such intrinsic
Membrane behaviour
LS-DYNA
and extrinsic uncertainties, modelling approaches must balance the competing demands of accuracy and
Extended Hamilton’s principle low computational demand. This article applies the extended Hamilton’s principle to rigid-plastic thin plates
subjected to impulsive blast loads to derive the governing equation of motion without a prior assumption of
the initial specific impulse distribution. Closed-form solutions to predict the plastic response are derived for
rectangular and circular plates. The analytical models for uniform specific impulses are found to be in good
agreement with high-fidelity numerical simulations performed using LS-DYNA and experimental data available
in the literature.

1. Introduction it is more appropriate not to perform a single simulation but rather

a large number of analyses, to fully characterise the likely response
There is significant scientific interest in studying the structural and associated confidence intervals. What, then, suits the analyst in
response in the inelastic range when structures are subjected to extreme this regard is the availability of models that can be evaluated more
dynamic loadings, such as missile impacts and explosions. Extremely efficiently.
high magnitudes and short durations characterise these loads, see [1– The present work aims to develop an analytical model to predict
4], and are expected to cause mainly plastic deformations in the the response of thin ductile plates under impulsive blast loads. To this
structures. end, the current knowledge about the physical problem (see Section 2)
Tyas and Rigby and their collaborators [1–3,5] provide experimen- is utilised while attempting to minimise the complexity of the problem.
tal evidence that near-field blasts generate spatially non-uniform spon- First, the load is assumed to be perfectly impulsive and, thus,
taneous pressures that, subsequently, decay exponentially. Such intense specified in terms of an initial velocity field derived from the prescribed
loading could induce large displacement, strain rate-dependent, and blast-induced specific impulse. According to the UFC 3-340-02 design
material failure effects in the structural response, see, e.g., Jones [6]. manual [7], the impulsive regime applies when the time to maximum
For these reasons, in addition to geometric complexities and bound- response to load duration ratio is greater than three, see Figure 1–
ary effects, the blast analyst typically resorts to making use of hy- 7 of the manual. Secondly, the material is idealised as rigid-perfectly
drocodes and Finite Element (FE) programs, attempting to simulate plastic, which obeys von Mises’s yield function and its associated flow
the complete problem. However, such an approach is computationally rule. Finally, the plate, which is restrained along its outer periphery, is
expensive. The challenge becomes more critical if we recognise that the assumed to deform in membrane mode without in-plane displacements;
input data of the blast threat are never known in advance, as pointed hence, flexural effects are ignored.
out by Tyas [3]. Furthermore, the material characteristics are prone, to The first assumption, due to Rigby et al. [8], is fundamental to
some extent, to potential variabilities. Thus, to design robust structures, the present analysis (the absence of externally applied forces ensures
a monotonic deformation path in a rigid-plastic structure, as will be
discussed in Section 3). The last assumption is based on the membrane
∗ Corresponding author at: Department of Civil & Structural Engineering, thinness and loading intensity.
University of Sheffield, Mappin Street, Sheffield, S1 3JD, UK.
E-mail address: salotaibi2@sheffield.ac.uk (S.A.E. Alotaibi).
1
PhD candidate.

https://doi.org/10.1016/j.ijimpeng.2023.104624
Received 4 January 2023; Received in revised form 28 March 2023; Accepted 17 April 2023
Available online 2 May 2023
0734-743X/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

The well-known extended Hamilton’s principle, see, e.g., [9,10], 2.2. Analytical models
is applied to the above described problem to obtain the appropri-
ate equation of motion governing the transverse displacement of the A subset of earlier and recent publications that studied the dynamic
thin membrane. The equation is found as a two-dimensional linear
response of rigid-plastic structures under extreme dynamic loadings
wave equation. Therefore, it is solved by the modal decomposition
(e.g. blast and impact) is summarised in Table 1. For more discussion
technique while enforcing Drucker’s [11,12] postulate of the plastic
of the analytical method and for additional references, the reader is
work nonnegativity. Analytical solutions are provided for rectangular
membranes subjected to arbitrary distribution of specific impulse and referred to [21–25], and [26].
for axisymmetric circular membranes. Martin and Symonds [41] developed the mode approximation tech-
The closed-form solutions for rectangular and circular membranes nique to predict the response of rigid-plastic structures. Martin [40]
loaded by uniform impulses are compared to experimental data found showed that the motion of an impulsively loaded structure eventually
in Nurick, Martin and Pearce [13], Gharababaei and Darvizeh [14], converges to a mode-form response. The current practice of simplified
and Nurick, Gelman, and Marshall [15]. In addition, the analytical blast analysis, which is adopted in [7] and applied in Rigby et al. [52],
solution is compared to results from LS-DYNA [16] simulations, where is also a mode approximation technique. This single-degree-of-freedom
the loads are prescribed through the initial velocity field according SDOF technique is originally developed by Biggs [53].
to Rigby et al. [8] and using relevant material parameters adopted In an analytical and experimental work, Nurick and Martin [20,
from [8,17,18]. Finally, the present model is compared to an existing 54] derived a non-dimensional impulse (or a damage number) that
model proposed by Yuen et al. [19] as a modification to Nurick and
correlates linearly with the normalised permanent displacements of
Martin’s [20] model. The uniform and non-uniform impulse solutions
rectangular and circular thin plates obtained from a large set of exper-
are found to be reasonably accurate with very low computational
imental data. Their model was revised (25 years later) by Nurick and
expense.
his colleagues [19]. The model applies to plates loaded by impulses of
2. Literature review constant amplitudes.
Cloete and Nurick [46] hypothesised that a quadratic displacement
2.1. Background field can be assumed to analyse a uniformly loaded thin plate. Using
such field, the authors solved the problem of a circular membrane
The field of investigating the response of engineering structures loaded with uniform impulse. They found that the central residual
under intense dynamic loads has attracted the interest of many re- transverse displacement varies linearly with Nurick and Martin’s non-
searchers. By far, the idealisation of a rigid-plastic material behaviour dimensional impulse. Moreover, it was shown that such displacement
has been the basic framework based on which several approximate is independent of the ‘‘in-plane’’ components of displacement.
closed-form solutions were obtained. In earlier works, responses to
In studying the response of thin plates under non-uniform blast
uniform dynamic pressure and uniform impulse were the main topic.
loads, Tyas, Rigby, and co-authors [8,55] established two important
Both small and large displacement regimes were considered in the
response of circular and rectangular plates and beams with various relationships. Firstly, the local initial velocity, 𝑤̇ 0 , of the target (with
boundary conditions. density 𝜌, thickness ℎ, and exposed area 𝐴) is confirmed to be a linear
Again, structures were idealised as rigid-plastic, and the particular function of the imparted non-uniform specific impulse field, 𝑖,
spatial form of the dynamic load is assumed a priori. The first step 𝑖
𝑤̇ 0 = ,
is to find the ‘‘quasi-static’’ limit load by means of the upper and 𝜌ℎ
lower limit analysis theorems and to assume the incipient collapse which is based on the upper bound kinetic energy. That is, non-
mechanism. The mechanism should be consistent with the underlying uniform specific impulse generates non-uniform initial velocity, see
yield function, its associated flow rule, and boundary conditions. At
Fig. 1(c). This initial velocity field contrasts with the ‘‘lower-bound’’
the static collapse load, motion follows the ‘‘quasi-static’’ mechanism.
initial velocity field,
As the dynamic load intensity increases further, determination of the
resulting dynamic mechanism becomes highly involved, and hence it ∫𝐴 𝑖 𝑑𝐴
𝑤̇ 0,𝑙𝑏 = .
is typically assumed. Then, the exact dynamic equilibrium equation is 𝜌ℎ𝐴
formulated in which the acceleration derives from the assumed collapse which is uniform and is interpreted as the rigid-body velocity (i.e. total
mechanism. The equation is then solved while avoiding any violation
impulse divided by total mass), see Fig. 1(b).
of the yield condition and the flow rule. If yield violation cannot be
Secondly, the maximum permanent displacement is a linear func-
avoided, then the initially assumed mechanism is wrong. In this case,
another initial mechanism must be chosen, and the process is repeated. tion of the energy-equivalent uniform total impulse, 𝐼𝑘 , which is also
The solution is valid when equilibrium, yield conditions, flow rule, and based on the upper bound kinetic energy uptake, 𝐾𝑢𝑏 , developed by
boundary conditions are all satisfied. Thin plates are found to initially Tyas and Pope [5], see Fig. 1,
respond in flexure, and as displacement increases further, membrane √
1
effects evolve and become dominant. The load-free membrane response 𝐾𝑢𝑏 = 𝑖2 𝑑𝐴, and 𝐼𝑘 = 𝐴 𝑖2 𝑑𝐴.
is typically driven by some initial displacement and velocity conditions, 2𝜌ℎ ∫𝐴 ∫𝐴
which are the final states of the initial flexural phase of response. The authors showed that the blast load can be replaced by a
Typically, the eventual correct solution could be associated with prescribed initial velocity field, which remarkably simplifies the anal-
a ‘‘dynamic’’ collapse mechanism that is substantially different from ysis. This observation highlights the importance of the actual spatial
the ‘‘quasi-static’’ mechanism. Since the material is rigid-plastic, then
variation of the specific impulse. According to the authors, the imparted
motion will eventually end, which happens when the external loads are
energy from the blast onto a thin plate cannot be higher than 𝐾𝑢𝑏 or
removed and the instantaneous velocity is zero everywhere (i.e. current
lower than the kinetic energy computed with 𝑤̇ 0,𝑙𝑏 , provided the regime
kinetic energy and external work are simultaneously zero).
is impulsive.
The methodology behind the earlier works has been continued
recently where non-uniform blast loads are also studied. However, Pannell et al. [56] proposed a model that predicts the non-uniform
the so-obtained solutions apply to the initially specified spatial forms distribution of specific impulses arising from near-field blasts with
of the loading function. Response to arbitrary loading has not been scaled distances of 0.11–0.55 m/kg1/3 . Determining the specific im-
obtained. Furthermore, no solution for membranes’ responses due to pulse profile is critical in analysing structures under near-field blast
initial conditions imposed directly by the blast load is available. loading, as discussed in [8].

2
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Table 1
Analytical works to study the dynamic response of rigid-plastic structures under extreme dynamic loadings.
Circular plate Rectangular plate Beam Theory
Small displacement (flexure) [27–31] [32,33] [34–37] [38–41]
Large displacement (membrane) [5,8,42–46] [14,26,47–49] [50,51]

mentioned, some discuss the theory and provide analytical solutions,

while others give experimental observations to strengthen or refine the
theory.
Commercially available FE programs can simulate the nonlinear
structural response under specific blast settings [8,17,63,66,67]. How-
ever, such approach is computationally expensive and requires signif-
icant user expertise. Furthermore, extra computations might be neces-
sary to ensure the solution is insensitive to any artificial parameters
that are not pertinent to the physical problem (e.g. those related
to hourglassing, artificial damping, structural and material locking,
extra contact and leakage controls, etc.) Therefore, the approach is
not directly suitable for probabilistic-based blast analyses that require
many repeated calculations. The situation takes us back to the interest
in building analytical solutions to idealised mechanical problems as
they are more computationally efficient and fast-running. This is the
rationale of the present work.
As discussed in Section 1, the present work aims to develop an
Fig. 1. Non-uniform specific impulse distribution (a) and the generated profile of initial analytical solution to predict the transient and permanent displace-
velocity assuming a lower-bound (b) and an upper-bound (c) kinetic energy uptake,
ment of thin ductile plates subjected to impulsive blast loads. The
reproduced after Tyas and Pope [5].
model utilises the following simplifying idealisations: impulsive blast
load, rigid-perfectly plastic material behaviour, and membrane mode
of deformation, according to observations from the already discussed
2.3. Experimental results
literature. We apply the extended Hamilton’s principle to obtain a
general equation of motion that applies to any target’s geometry with-
Nurick and his collaborators [13,15,54,57,58] presented extensive
out a prior assumption on the initial field. The then-derived equation
experimental data for circular and rectangular plates under uniform
and solution will accommodate any distribution of the blast-induced
impulsive loads.
specific impulse.
As discussed in Section 1, the team of the Blast and Impact Engi-
neering research group at the University of Sheffield [1–3,5] confirmed
experimentally that near-field blasts produce spatially non-uniform 3. Theory
specific impulses. Furthermore, the researchers and their collaborators
presented experimental and numerical studies of the transient defor- 3.1. Problem definition
mation of thin ductile plates subjected to near-field blast [55]. In their
study, both the blast load and the full-field transient displacement of
the target were measured locally: (a) the dynamic pressures histories at Consider a thin membrane made of a rigid-perfectly plastic metal
different spatial points were measured by the Characterisation of Blast that is subjected to a prescribed initial velocity field obtained from a
Loading (CoBL) apparatus [2,59], providing the spatial variations of the specific impulse distribution according to Tyas and Pope [5] and Rigby
specific impulses; (b) full-field transient displacements of the targets et al. [8].
were measured using high-speed imaging, described further in [18]; The membrane can be of a rectangular or circular geometry, and it
and (c) the total impulse was measured using a ballistic pendulum. is supported along its outer periphery. For the rectangular geometry, 𝐿𝑥
Also, Langdon and her collaborators [18,60,61] studied the response and 𝐿𝑦 are the sides’ lengths, as indicated in Fig. 2, and 𝑅 is the radius
of steel plates under near-field blast experimentally. It was observed of the circular membrane. The specific impulse is denoted by 𝑖(𝑥, 𝑦)
that a localised blast load induces central dishing (or bulging) in the or 𝑖(𝑟), where 𝑥 and 𝑦 are the rectangular undeformed coordinates for
targets. Gharababaei and Darvizeh [14] performed blast experiments rectangular geometry, whereas 𝑟 is the radial undeformed coordinate
on steel, aluminium, and copper thin circular plates which were loaded for the circular membrane. Let 𝑡 denotes time, and 𝜌, ℎ, and 𝜎0 denote,
by detonating thin cylindrical charges of high explosives. The authors respectively, density, thickness, and characteristic yield strength of the
used rigid tubes to guide the propagation of the shock waves as they membrane.
travel towards the specimens. Large plastic deformations were observed The membrane is assumed to respond in pure membrane mode.
in the tests.
Further, it is assumed that the components of displacement along the
In addition, Aune and his team presented experimental and numer-
undeformed in-plane coordinates are negligible. That is, the only non-
ical studies on the response of thin plates under free-air blasts [62,63]
zero displacement is the one along the original out-of-plane coordinate,
and blasts produced in a shock tube facility [64]. Additional air-blast
which is denoted by 𝑤. In other words, every particle of the membrane
experiments are also reported by Spranghers et al. [65,66]. The plates
displaces vertically.
are made of ductile materials (structural steel and aluminium) and were
observed to respond impulsively and plastically. In line with the perfect plasticity assumption, the characteristic
yield strength 𝜎0 is assumed constant. Thus, for materials exhibiting
2.4. Aims of the present study substantial work-hardening [64,68], 𝜎0 could be defined such that a
rectangular stress-plastic strain curve preserves the total area under the
The above discussion serves as a short account of some notable actual uniaxial stress-plastic strain curve. This approach is adopted in
works investigating the response of simple structures to blast loads. As the UFC 3-340-02 [7] manual.

3
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Eq. (3.2) allowed a transition from the flow to the total theory
of plasticity, and it holds only when the path is monotonic. In fact,
Drucker [12] assures that when loading is monotonic, the two plasticity
theories are identical. (√ )
Using Eq. (3.2), it can be shown that 𝜆 = 3∕2 𝜀𝑝𝑖𝑗 𝜀𝑝𝑖𝑗 ∕𝜎0 , and the
total plastic work, which is

𝑊𝑝 = 𝜎𝑖𝑗 𝜀𝑝𝑖𝑗 𝑑𝑉 , (3.3)

∫𝑉
( )
can be expanded as 𝑊𝑝 = ∫𝑉 𝑠𝑖𝑗 𝜀𝑝𝑖𝑗 𝑑𝑉 = ∫𝑉 𝜆𝑠𝑖𝑗 𝑠𝑖𝑗 𝑑𝑉 = ∫𝑉 𝜆 2 2
𝜎
3 0
𝑑𝑉 ,
or
√
𝑊𝑝 = 𝜎0 2∕3 𝜀𝑝𝑖𝑗 𝜀𝑝𝑖𝑗 𝑑𝑉 . (3.4)
∫𝑉
where the plastic incompressibility, 𝜀𝑝𝑖𝑖 = 𝜆𝑠𝑖𝑖 = 0, of the von-Mises
material was utilised.
From now on, the superscript 𝑝 will be omitted as total strain and
plastic strain are identical in line with the rigid-plastic assumption. To
evaluate the quantity 𝜀𝑖𝑗 𝜀𝑖𝑗 , the four non-vanishing strain components
as functions of the transverse displacement, 𝑤(𝑥, 𝑦, 𝑡), will be utilised.
These are
( ) ( )2
1 𝜕𝑤 2 1 𝜕𝑤
𝜀𝑥 = , 𝜀𝑦 = ,
2 𝜕𝑥 2 𝜕𝑦
𝜕𝑤 𝜕𝑤 √
𝛾𝑥𝑦 = = 2 𝜀𝑥 𝜀𝑦 , (3.5)
Fig. 2. Problem definitions for the rectangular membrane: (a) plan view showing 𝜕𝑥 𝜕𝑦
undeformed geometry and (b) side view showing a typical spatial distribution of specific
and from the incompressibility condition, one has
impulse 𝑖.
𝜀𝑧 = −(𝜀𝑥 + 𝜀𝑦 ). (3.6)

Then, through some algebra simplifications, it can be shown that

3.2. Development of the equation of motion
𝜀𝑖𝑗 𝜀𝑖𝑗 = 2(𝜀𝑥 + 𝜀𝑦 )2 . (3.7)
The extended Hamilton’s principle is applied to the present problem
Therefore, the total plastic work under a monotonic deformation
to systematically derive the equation of motion. In applying the princi-
path becomes
ple, treatment of kinetic energy and external work (if any) terms are not
new and will not be shown for brevity. However, the internal energy 2
𝑊𝑝 = √ 𝜎0 (𝜀𝑥 + 𝜀𝑦 ) 𝑑𝑉
term needs special consideration. First, there is no elastic strain energy ∫𝑉 3
[ ( )2 ]
in the system, and the only allowed internal energy (arising from the ( )
2 1 𝜕𝑤 2 1 𝜕𝑤
accumulation of plastic deformation) is dissipative in nature. To obey = √ 𝜎0 + 𝑑𝑉 . (3.8)
∫𝑉 3 2 𝜕𝑥 2 𝜕𝑦
such irreversible behaviour, the rate of plastic work must always be
non-negative [11,12]. That is, when strain rate tends to change sign, It should be noted that the stress 𝜎𝑖𝑗 , appearing in Eq. (3.3), is the
stress must instantly do so in a rigid-perfectly plastic structure. second Piola–Kirchhoff stress since it is the work-conjugate to Green–
In Hamilton’s principle, the total internal energy is the one that Lagrange strain. However, for consistency of the formulation, the small
should be included, not the rate of energy. The total energy is the strain assumption implies that this mentioned stress can be replaced
time integral of the energy rate. However, the time integration can be with the true Cauchy stress. Hence, we refer to 𝜎𝑖𝑗 throughout as the
carried out beforehand (i.e. explicitly) when the total plastic strain is Cauchy stress.
monotonic since the stress would be constant in terms of sign (recall Then, the extended Hamilton’s principle is applied, which reads
that its magnitude is already constant from the perfect plasticity). Now, ( 𝑡2 )
since there are no external forces, the total plastic strain is guaranteed 𝛿 𝐻 𝑑𝑡 = 0, (3.9)
∫𝑡1
to be monotonic (i.e. the sign of strain rate is fixed) when the transverse
where 𝛿 is the variational operator, and 𝑡1 and 𝑡2 are arbitrary times.
velocity is monotonic.
In Eq. (3.9), the Hamiltonian, 𝐻, of the system (in absence of external
In the following, 𝜎𝑖𝑗 , 𝑠𝑖𝑗 , 𝜀𝑖𝑗 , and 𝜆 are the stress, deviatoric stress,
work and elastic strain energy) is given by
(Green–Lagrange) strain, and the plastic multiplier, respectively. A dot
[ ]
(. ) over a symbol denotes time differentiation, a repeated subscript 𝐻 = 𝐾 − 𝑊𝑝 . (3.10)
implies summation, and a superscript (𝑝 ) denotes the plastic part.
in which 𝐾 is the total kinetic energy of the membrane, which is
With reference to the von Mises’s yield function, 𝑓 (𝜎𝑖𝑗 ) = 12 𝑠𝑖𝑗 𝑠𝑖𝑗 ,
𝐾 = ∫𝑉 12 𝜌𝑤̇ 2 𝑑𝑉 , and 𝑊𝑝 was given in Eq. (3.8). Notice that 𝐻 becomes
yielding of the material occurs when 𝑓 (𝜎𝑖𝑗 ) = 13 𝜎02 . Associated with such
a functional of 𝑤 only.
yielding condition, the incremental flow rule reads
Finally, by (i) applying 𝛿 on 𝑤, (ii) carrying out spatial and temporal
𝜕𝑓 1 2 integration by parts, (iii) imposing constraints on 𝑤 at the membrane
𝜀̇ 𝑝𝑖𝑗 = 𝜆̇ ̇ 𝑖𝑗
= 𝜆𝑠 𝜀̇ 𝑝𝑖𝑗 ≡ 0 if 𝑓 (𝜎𝑖𝑗 ) < 𝜎 (3.1)
𝜕𝜎𝑖𝑗 3 0 edges and anywhere at times 𝑡1 and 𝑡2 where 𝛿𝑤 vanishes identically,
Now, if the material has actually yielded, i.e. 𝑠̇ 𝑖𝑗 ≡ 0, then under a and (iv) requiring 𝛿𝑤 to be otherwise arbitrary, then one obtains the
monotonic deformation regime, the flow rule can be integrated by parts (Euler–Lagrange) equation of motion governing the response of the
(while taking advantage of 𝑠̇ 𝑖𝑗 = 0 and assuming 𝜆|𝑡=0 = 0) to obtain rigid-perfectly plastic membrane as
[ 2 ]
2 𝜕 𝑤 𝜕2 𝑤
𝜀𝑝𝑖𝑗 = 𝜀̇ 𝑝𝑖𝑗 𝑑𝑡 = ̇ 𝑖𝑗 𝑑𝑡 = 𝜆𝑠𝑖𝑗 −
𝜆𝑠 𝜆𝑠̇ 𝑖𝑗 𝑑𝑡 = 𝜆𝑠𝑖𝑗 . (3.2) √ 𝜎0 + = 𝜌𝑤.
̈ (3.11)
∫ ∫ ∫ 3 𝜕𝑥2 𝜕𝑦2

4
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624
√
𝜋𝑐
Eq. (3.11) is a two-dimensional scalar (plastic) wave equation in 𝜔𝑚𝑛 = (𝐿𝑦 𝑚)2 + (𝐿𝑥 𝑛)2 (4.3)
𝐿𝑥 𝐿𝑦
𝑤(𝑥, 𝑦, 𝑡) in a rectangular coordinate system and with a wave speed ( ) ( )
[ √ ]1∕2 𝐿𝑥 𝐿𝑦
2𝜎0 ∕( 3 𝜌) . Note that Eq. (3.11) is a field equation so that it applies 𝑚𝜋𝑥 𝑛𝜋𝑦
𝐼𝑚𝑛 = 𝑖(𝑥, 𝑦) sin sin 𝑑𝑦 𝑑𝑥. (4.4)
∫0 ∫0 𝐿𝑥 𝐿𝑦
to any membrane geometry with restrained edges.
Later, the equation is solved for a rectangular membrane, as defined The pair (𝑚, 𝑛) defines a particular mode with mode shape 𝜙𝑚𝑛 (𝑥, 𝑦).
in Fig. 2, under the following kinematic conditions It is known that the modes are orthogonal over the membrane domain,
and hence they are independent. Thus, to fulfill the requirement of
𝑤(0, 𝑦, 𝑡) = 𝑤(𝐿𝑥 , 𝑦, 𝑡) = 𝑤(𝑥, 0, 𝑡) = 𝑤(𝑥, 𝐿𝑦 , 𝑡) = 0, a monotonic deformation path (and hence obey the plastic work’s
𝑤(𝑥, 𝑦, 0) = 0, (3.12) non-negativity), a strategy is adopted to enforce the termination of a
particular mode when the associated (modal) velocity reaches zero for
and the dynamic condition the first time. This occurs when
𝑖(𝑥, 𝑦) 𝜋
̇
𝑤(𝑥, 𝑦, 0) = . (3.13) 𝑡 = 𝑡𝑚,𝑛 = . (4.5)
𝜌ℎ 2𝜔𝑚𝑛
Eq. (3.13) is experimentally shown to hold for thin plates under non- From 𝜔𝑚𝑛 , see Eq. (4.3), it is clear that the sequence of turning off
uniform specific impulse [5,8,55], which was derived from the balance the modal contributions is ordered from highest to lowest modes in
of linear momentum of shear non-rigid thin plates. terms of frequency. The notion of sequential terminations of the modes
For a circular membrane and under axisymmetric conditions (which was previously used in [43,45]. The last contributing mode is the first,
will be assumed throughout), it can be shown, through the standard i.e. with (𝑚 = 𝑛 = 1). Thus, the whole membrane ceases motion at no
transformation from rectangular to polar coordinates, that the equation later than 𝑡 = 𝑡1,1 , which is given by
of motion is given as 𝐿∗
[ 2 ] 𝑡1,1 = (4.6)
2 𝜕 𝑤 1 𝜕𝑤 2𝑐
√ 0 𝜎 + = 𝜌𝑤.
̈ (3.14)
3 𝜕𝑟2 𝑟 𝜕𝑟 where 𝐿∗ , defined for convenience, is the ratio of the membrane area
to the length of its diagonal and given by
where a term within the brackets on the left-hand side (1∕𝑟2 ) 𝜕 2 𝑤∕𝜕𝜃 2
has been omitted due to the axisymmetric assumption. Eq. (3.14) will 𝐿𝑥 𝐿𝑦
𝐿∗ ≡ √ . (4.7)
be solved under the following conditions 𝐿2𝑥 + 𝐿2𝑦
𝑖(𝑟)
𝑤(𝑅, 𝑡) = 𝑤(𝑟, 0) = 0, ̇ 0) =
𝑤(𝑟, . (3.15) Note that if 𝑡1,1 is less than three times the duration of blast load, the
𝜌ℎ
global response is less likely to be impulsive, based on [7], since 𝑡1,1 is
Although the actual problem involves plastic deformations, the an upper bound on the response time. The actual time of maximum
obtained equations of motion are linear. Hence, they can be solved by response is the modal time 𝑡𝑚,𝑛 of the dominant mode (whose total
Fourier’s decomposition, which is based on the principle of superposi- modal impulse is 𝐼𝑚𝑛 ) of the membrane under a particular distribution
tion. of specific impulse 𝑖. If there are several dominant modes, then the
It should be recalled that the above equations of motion are valid maximum time is the largest 𝑡𝑚,𝑛 among these modes.
as long as the deformation path remains monotonic so that Eq. (3.2) is The permanent shape, 𝑤𝑝 (𝑥, 𝑦), of the membrane is given by 𝑤(𝑥, 𝑦, 𝑡)
not violated. Hence, the solutions (to be presented later) of Eqs. (3.11) when 𝑡 ≥ 𝑡1,1 , or
and (3.14) are valid up to the instant of time when the transverse
4 ∑
∞
velocity, 𝑤,
̇ tends to change sign, or simply when velocity reaches zero. 𝐼𝑚𝑛
In other words, a component of the solution must terminate whenever
𝑤𝑝 (𝑥, 𝑦) = √ 𝜙𝑚𝑛 (𝑥, 𝑦). (4.8)
𝜋𝜌𝑐ℎ 𝑚,𝑛=1
(𝐿𝑥 𝑛)2 + (𝐿𝑦 𝑚)2
the velocity associated with it reaches zero for the first time. Otherwise,
plastic work would decrease and thereby violating its irreversibility or For all cases in which the specific impulse distribution is symmetric
dissipating nature, and thus the solution becomes non-physical. about the membrane’s centre, the peak displacement is located at the
Many textbooks, e.g. references that treat the elastic free vibration centre. The central permanent displacement, 𝑤𝑐 ≡ 𝑤𝑝 (𝐿𝑥 ∕2, 𝐿𝑦 ∕2), is
of pre-tensioned membranes under small displacements, such as [69,
( ) ( )
4 ∑
∞
Ch. IX] or [70, Sec. 69], show how equations similar to Eqs. (3.11) 𝐼𝑚𝑛 𝑚𝜋 𝑛𝜋
𝑤𝑐 = √ sin sin . (4.9)
and (3.14) can be solved. Excellent stepwise derivations are presented 𝜋𝜌𝑐ℎ 𝑚,𝑛=1 2 2 2 2
(𝐿𝑥 𝑛) + (𝐿𝑦 𝑚)
in [71]. Hence, the derivation steps of the solution will be omitted for
brevity. Instead, the solutions themselves are given. According to Pannell et al. [56], the specific impulse distribution
from a near-field spherical charge blast is of a Gaussian form as a
4. Rectangular membrane function of the angle of incidence. However, it was not possible to
evaluate 𝐼𝑚𝑛 symbolically for a specific impulse distribution, 𝑖(𝑥, 𝑦), as
4.1. Response of rectangular membrane predicted by Pannell et al. Hence, numerical integration is needed.
A practical Matlab code for calculating 𝐼𝑚𝑛 using the Fast Fourier
The rectangular membrane equation of motion, Eq. (3.11), was Transform (FFT) is given in Appendix B.
solved by the modal decomposition technique under the prescribed Fig. 3 shows the normalised permanent displacement, 𝑤𝑝 , profiles
geometric conditions, Eqs. (3.12) and (3.13). Its solution is along 𝑦 = 𝐿𝑦 ∕2 due to three impulse distributions with constant ampli-
∑∞
𝐼𝑚𝑛 ( ) tudes applied over varying central areas of a rectangular membrane. In
4
𝑤(𝑥, 𝑦, 𝑡) = 𝜙 (𝑥, 𝑦) sin 𝜔𝑚𝑛 𝑡 (4.1) the figure, the legends indicate the ratios of the loaded to total areas. It
𝜌ℎ𝐿𝑥 𝐿𝑦 𝑚,𝑛=1 𝜔𝑚𝑛 𝑚𝑛
is, thus, evident that a localised impulse induces localised displacement
with the wave speed 𝑐, mode shape 𝜙𝑚𝑛 (𝑥, 𝑦), modal angular frequency shape, i.e. with central dishing, while the case of uniform impulse
𝜔𝑚𝑛 , and total modal impulse 𝐼𝑚𝑛 given by applied over the whole area of the membrane results in global uniform
√ dishing.
2 𝜎
𝑐= √ 0 (4.2) In practice, a finite number of modes is used to numerically evaluate
3 𝜌 the permanent, 𝑤𝑝 (𝑥, 𝑦), and permanent central, 𝑤𝑐 , displacements.
( ) ( )
𝑚𝜋𝑥 𝑛𝜋𝑦 To maintain sufficient accuracy while truncating the infinite series in
𝜙𝑚𝑛 (𝑥, 𝑦) = sin sin Eqs. (4.8) and (4.9), the error estimate given in Appendix A can be
𝐿𝑥 𝐿𝑦

5
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Fig. 3. Permanent displacement profile for rectangular membrane under three sets of
impulses of constant amplitudes applied over central rectangular regions with loaded-
to-total area ratios of (1) 0.062, (2) 0.25, and (3) 1.0. Curve (1) is associated with
localised impulse, and (3) results from a uniform impulse over the entire membrane.

Fig. 4. The converged value of truncated sum, 𝑆0 , for the rectangular membrane
used. The error relates the sum of modal impulses (included in the associated with uniform impulse case as a function of the membrane aspect ratio,
𝐿𝑦 ∕𝐿𝑥 , where 𝐿𝑦 ≤ 𝐿𝑥 .
approximation) to the energy-equivalent total impulse, 𝐼𝑘 which is
√
𝐼𝑘 = 𝐴 𝑖2 𝑑𝐴
∫𝐴 5. Circular membrane
in which 𝐴 is the area of the membrane. The absolute importance
of a particular mode is shown to be indicated by the measure 0 ≤ 5.1. Response of circular membrane
(2𝐼𝑚𝑛 ∕𝐼𝑘 )2 ≤ 1, which can also be used to identify the dominant
mode(s); dominant modes have values closer to unity. As stated in Section 3, axisymmetric conditions are assumed for
the circular membrane problem. Eq. (3.14) was solved by the modal
4.2. Uniform specific impulse case — rectangular membrane decomposition technique. The solution is presented below, Eq. (5.1),
which gives the displacement of a circular rigid-perfectly plastic mem-
For the case where the specific impulse has a constant distribution, brane of radius 𝑅, mass density 𝜌, characteristic yield strength 𝜎0 , and
i.e. 𝑖(𝑥, 𝑦) = 𝑖0 , the modal impulse, 𝐼𝑚𝑛 , simplifies to thickness ℎ, due to a specific impulse (impulse per unit area), 𝑖(𝑟).
{ 4𝑖 𝐿 𝐿
2 ∑
∞
0 𝑥 𝑦
when (𝑚, 𝑛) are odd 𝐼𝑚 ( )
𝐼𝑚𝑛 = 𝜋 2 𝑚𝑛 𝑤(𝑟, 𝑡) = 𝜙 (𝑟) sin 𝜔𝑚 𝑡 , (5.1)
0 otherwise 𝜌𝑐𝑅ℎ 𝑚=1 𝑗0,𝑚 𝐽1 (𝑗0,𝑚 )2 𝑚

which suggests that the first mode is the most dominant. where
( )
The peak displacement is located at the plate’s centre and is given 𝑗0,𝑚
𝜙𝑚 (𝑟) = 𝐽0 𝑟 , (5.2)
by 𝑅
( ) ( ) 𝑐𝑗0,𝑚
16𝑖0 𝐿𝑦 ∑ ∞ sin 𝑚𝜋
2
sin 𝑛𝜋2
𝜔𝑚 = , (5.3)
𝑅
𝑤𝑐 = √ . (4.10) 𝑅 ( )
𝜋 3 𝜌𝑐ℎ 𝑚,𝑛=1,3,5 𝑚𝑛 𝑛2 + 𝑚2 (𝐿 ∕𝐿 )2 𝑗0,𝑚
𝑦 𝑥 𝐼𝑚 = 𝑖(𝑟) 𝐽0 𝑟 𝑟𝑑𝑟 (5.4)
∫0 𝑅
In Eq. (4.10), the summand depends only on the membrane’s aspect
𝐼𝑚 is the total modal impulse per unit radian, 𝑐 is the wave speed
ratio 𝐿𝑦 ∕𝐿𝑥 . Denoting the sum by 𝑆0 , it was observed to converge. The
given in Eq. (4.3), 𝜙𝑚 (𝑟) is the 𝑚th mode shape, and 𝜔𝑚 is the corre-
value of 𝑆0 for any aspect ratio in the range [0.1 − 1.0] can be read from
sponding frequency. In above, 𝐽0 (𝑥) and 𝐽1 (𝑥), respectively, are Bessel
Fig. 4, in which 𝐿𝑦 ≤ 𝐿𝑥 .
functions of the first kind of order zero and one, while the scalar value
With 𝑆0 known, the central displacement of a membrane due to a
𝑗0,𝑚 is the 𝑚th root of 𝐽0 (𝑥) = 0, i.e. 𝐽0 (𝑗0,𝑚 ) ≡ 0. The solution does
uniform impulse (with intensity 𝑖0 ) is given by
not involve Bessel functions of the second kind to avoid infinite (non-
16𝑖0 𝐿𝑦 physical) response at the origin (plate’s centre). Furthermore, the modal
𝑤𝑐 = 𝑆0 . (4.11)
𝜋 3 𝜌𝑐ℎ solution depends only on the zeroth order Bessel function due to the
Or, using the total impulse 𝐼0 = ∫𝐴 𝑖(𝑥, 𝑦) 𝑑𝐴 = 𝐴 𝑖0 , the last expression axisymmetry of the problem. Again, the 𝑚th mode shape is 𝐽0 (𝑗0,𝑚 𝑟∕𝑅),
[ ]2
becomes and the square of its norm (per unit radian) is (1∕2)𝑅2 𝐽1 (𝑗0,𝑚 ) .
16𝑆0 The modes, 𝜙𝑚 (𝑟), are orthogonal to each other, and thus they are
𝑤𝑐 = 𝐼 0 ≡ 𝑘0 𝐼 0 (4.12) independent.
𝜋 3 𝜌𝑐ℎ𝐿𝑥
Similar to the rectangular case, the contribution from a given mode
where a structural parameter 𝑘0 was introduced and defined as
in the solution is valid until the corresponding modal velocity reaches
16𝑆0 zero when 𝑡 ≥ 𝑡𝑚 , where 𝑡𝑚 = 𝜋𝑅∕(2𝑐𝑗0,𝑚 ). Thus, the modes will be
𝑘0 = . (4.13)
𝜋 3 𝜌𝑐ℎ𝐿𝑥 switched off sequentially in descending order with respect to frequency.
It should be re-emphasised that 𝐿𝑦 is the shorter side’s length, i.e. 𝐿𝑦 ∕ Thus, the whole membrane motion terminates at or before 𝑡 = 𝑡1 =
𝐿𝑥 ≤ 1, when 𝑆0 is estimated from Fig. 4. 𝜋𝑅∕(2𝑐𝑗0,1 ), where 𝑗0,1 = 2.405.

6
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Table 2
The first seven roots of the zeroth-order Bessel function 𝐽0 (𝑥) and their related quantities, computed using [72].
𝑚 1 2 3 4 5 6 7
𝑗0,𝑚 2.4048 5.5201 8.6537 11.7915 14.9309 18.0711 21.2116
𝐽1 (𝑗0,𝑚 )2 0.2695 0.1158 0.0737 0.0540 0.0427 0.0352 0.0300
𝐽1 (𝑗0,𝑚 )∕𝑗0,𝑚 0.2159 −0.0616 0.0314 −0.0197 0.0138 −0.0104 0.0082

in which the numerical value of the infinite sum is denoted with 𝑆0

and evaluates to 0.2674. Hence, the central permanent displacement
of a circular membrane due to a uniform specific impulse (of intensity
𝑖0 ) is
2𝑅𝑖0
𝑤𝑐 = × 0.2674. (5.8)
𝜌𝑐ℎ
Alternatively, in terms of total impulse 𝐼0 = 𝜋𝑅2 𝑖0 , the central
displacement becomes
2 × 0.2674
𝑤𝑐 = 𝐼0 ≡ 𝑘0 × 𝐼0 (5.9)
𝜋𝑅𝜌𝑐ℎ
with the circular structural parameter 𝑘0 defined as
2 × 0.2674
𝑘0 = . (5.10)
𝜋𝑅𝜌𝑐ℎ
The model predicts the normalised permanent membrane profile
Fig. 5. Permanent displacement profile of an axisymmetric circular membrane sub-
for the uniform specific impulse case as shown by curve (2) in Fig. 5,
jected to a uniform impulse which is applied over a localised region (1) with a
loaded-to-total area ratio of 0.062 and over the whole area (2) of the membrane as presented earlier.
predicted by the present solution.
6. Model verification

The permanent membrane profile, 𝑤𝑝 (𝑟), is obtained from Eq. (5.1) 6.1. Rectangular model verification
for 𝑡 ≥ 𝑡1 and is given by
6.1.1. Uniform impulse
2 ∑
∞
𝐼𝑚
𝑤𝑝 (𝑟) = 𝜙 (𝑟). (5.5) The solution 𝑤𝑐 for the uniform impulse case, as given by Eq. (4.12),
𝜌𝑐𝑅ℎ 𝑚=1 𝑗0,𝑚 𝐽1 (𝑗0,𝑚 )2 𝑚
is validated using experimental data obtained from the literature as
Since the problem is axisymmetric, the central displacement is the detailed in this section. Nurick, Martin, and Pearce [13] carried out
peak. Using the fact that 𝐽0 (0) = 1, the central permanent displacement, 82 experiments using rectangular and square thin plates made of steel,
𝑤𝑐 ≡ 𝑤𝑝 (𝑟 = 0), is and the specimens were loaded impulsively by distributed sheets of ex-
plosives. The authors suggest that the impulse distributions are uniform
2 ∑
∞
𝐼𝑚
𝑤𝑐 = . (5.6) over the specimens’ exposed surfaces. The sides measured 113 mm and
𝜌𝑐𝑅ℎ 𝑚=1 𝑗0,𝑚 𝐽1 (𝑗0,𝑚 )2
70 mm for the rectangular plates, and a side length of 89 mm was given
Some available tools, e.g. Matlab native function besselj() and the for the square plates. The thickness and static yield strength were given
user-built Matlab function in [72], can be utilised to evaluate Bessel as 1.6 mm and 296 MPa, respectively; the mass density is assumed to
quantities appearing in the above expressions. Table 2 is provided for be 7830 kg/m3 . In the tests, the amounts of explosives were varied,
quick estimation purposes. resulting in different values of total impulses measured using a ballistic
Similar to the rectangular case, it was not possible to evaluate pendulum. The permanent central displacements were measured and
the integral 𝐼𝑚 for a specific impulse distribution, 𝑖(𝑟), of the type given in the paper, and further details can be found there.
predicted by Pannell et al.’s [56] model. Thus, 𝐼𝑚 needs to be computed Nurick et al. [13] data is used to validate the solution, Eq. (4.12).
numerically. The static yield strength reported in the experiments is taken as 𝜎0 in
Experiments indicate that when a circular membrane is subjected to the model. The results of the comparisons are shown in Figs. 6 and 7.
localised impulse, say, over the membrane’s central region, then central Nurick et al. experiments were simulated numerically using LS-
dishing results, see Curry and Langdon [18]. Thus, to qualitatively test DYNA [16], in which the plates were subjected to uniform initial
the developed solution, such a loading case was simulated that results velocity fields. The steel material was modelled using the *Mat_Simplif-
in the permanent shape depicted in Fig. 5 by curve (1), in which a ied_Johnson_Cook model, available in LS-DYNA [73], which accounts for
central bulging can be seen. strain-hardening and strain-rate effect on the current yield stress. Ther-
mal softening and strain-based failure were neglected in the analyses.
5.2. Uniform specific impulse case — circular membrane The material parameters, except the static yield strength, were taken
from [8,17]. The rectangular plates were modelled as fully integrated
When the specific impulse (impulse per unit area) is spatially uni- shell elements using *Element_Shell and *Section_Shell keywords with
form with intensity 𝑖0 , the total modal impulse 𝐼𝑚 simplifies to 𝐼𝑚 = Elform = 16 (free of hourglass modes) and three through-thickness in-
[𝑖0 𝑅2 𝐽1 (𝑗0,𝑚 )]∕𝑗0,𝑚 . tegration points, NIP = 3, (to incorporate flexural effects). The uniform
In this case, the permanent central displacement of a circular mem- initial velocities were prescribed using the *Initial_Velocity_Node key-
brane is word. Nodes on the plates’ peripheries were restrained in all (including
2𝑅𝑖0 ∑ [ 2 ]−1
∞
the rotational) degrees of freedom.
𝑤𝑐 = 𝑗 𝐽 (𝑗 )
𝜌𝑐ℎ 𝑚=1 0,𝑚 1 0,𝑚 The peak displacement at the plate’s centre was used to determine
2𝑅𝑖0 an appropriate element mesh density, which was then held fixed in
= 𝑆 (5.7) subsequent analyses. The permanent displacement was determined by
𝜌𝑐ℎ 0

7
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Furthermore, the slight discrepancy can also be linked to the perfect

plasticity assumption on the material behaviour. The model excludes
strain-hardening and strain-rate effects, and for simplicity the static
yield strength is taken as the perfectly plastic limit. However, many
practical ductile materials exhibit various levels of yield enhancement
due to these plastic effects. Heat-treated metals, e.g., typically possess
substantial work-hardening where the current yield strength progresses
gradually from the (less distinct) initial yield stress up to the ultimate
strength that could be as twice as the initial yield strength [68]. A cer-
tain class of steel, which has been used in blast environments, exhibits
different hardening-ductility characteristics depending on the material
composition, heat-treatment, and manufacturing processes [64]. Ideal-
ising plates made of hardening materials as perfectly plastic would then
result in overestimations of the responses under large dynamic loads.
Therefore, it is important that such plastic characteristics are incorpo-
rated in computer numerical analyses since mathematical complexity
is not a barrier.
The perfect plasticity assumption is adopted in the present work
to obtain a ‘‘first-order approximation’’ model. Addition of the two
mentioned effects into the model was found, by the authors, to lead
to a nonlinear equation of motion. Therefore, the solution given by the
Fig. 6. Comparison of model predictions to experimental data [13] and LS-DYNA
numerical results in terms of central residual displacement 𝑤𝑐 of rectangular and square present perfectly plastic model should be regarded as an upper-bound
membranes under uniform impulse of total magnitude of 𝐼0 ; 𝑘0 is the membrane’s solution for structures made of materials that (in practice) deviate from
parameter defined in Eq. (4.13). the perfect plasticity behaviour.
If desirable, we propose that the yield strength to be used in the
model might be adjusted (e.g. amplified) to compensate for large strain-
averaging (through time integration) the displacement time history hardening strengthening. For example, the actual area under the plastic
beyond the first peak over a small number of vibration cycles. part of the engineering stress–strain curve could be converted to a rect-
LS-DYNA results are compared to the experiments and the analytical angular area, and by maintaining the ultimate strain, the characteristic
predictions for the rectangular and square tests (separately) as shown (or effective) yield strength can then be determined. This procedure
in Fig. 7. should give better predictions while maintain the upper-bound sense.
Aune and his team at SIMlab presented experimental and numerical However, in the foregoing validation work, the model is evaluated
studies on the response of thin plates under free-air blasts [62,63] and using the ‘‘static’’ yield strength because it was found to agree with
blasts produced in a shock tube facility [64]. Additional air-blast exper- the data.
iments are also reported by Spranghers et al. [65,66]. The rectangular A simplified study to assess the effect of work-hardening on the
plates are made of ductile materials (structural steel and aluminium). response of a single-degree-of-freedom (SDOF) due to blast-type loading
is presented in Appendix C. Few representative numerical cases are
Some specimens are seen to respond impulsively and plastically. From
given in which the ratios of the constant hardening moduli to elastic
the span-to-thickness ratios of some experiments and the intensity
stiffness are in a range of practical values. Furthermore, a numerical
of the generated blast loads, the targets are thought to respond in
parametric study using LS-DYNA of the effect of strain-hardening on
membrane mode. Therefore, their experiments can, in general, be used
the residual response of thin plates is also given in Appendix C; in
to further assess the accuracy of the present model. However, few tests
addition, the corresponding predictions by the present analytical model
should be excluded from the validation analysis in which the plates:
are compared to LS-DYNA results.
(1) experienced complete fractures, (2) had pre-formed cracks or holes,
The experiments in [62–66], discussed earlier, can be used to
or (3) experienced a (counterintuitive) bifurcation or rebound buckling
quantify the effect of neglecting the work-hardening in the model
due to their ultra thinness when combined with low blast intensity. The
as some uniaxial tensile specimens show a significant presence of
validation of the model against the above experiments is not presented
strain-hardening.
herein because the original works, cited above, did not provide quanti-
Finally, it should be re-emphasised that values of the material
tative information about the distribution of the blast-generated specific
parameters, except the static yield strength, used in the LS-DYNA sim-
impulses, which the present model requires. ulations are assumed. The source papers reported neither the Johnson–
Cook (JC) parameters nor re-usable stress–strain data to enable iden-
6.1.2. Discussion and model limitations tifications of the constitutive plastic parameters. For example, Nurick
Figs. 6 and 7 exhibit a reasonable accuracy of the analytical model, et al. [13] present the engineering stress–strain data at different strain
Eq. (4.12), although the model accounts only for membrane behaviour, rates. Using their data, it was attempted to obtain digitally converted
rigid-perfectly plastic material model, and idealised impulsive load. true stress–strain curves and determine the material parameters by
However, as can be seen in the figures, the analytical model ap- curve-fitting the data to JC model. However, this was not possible
pears to be vertically offset when compared to the data. The offset due to precision issues associated with the strain axis resolution. Thus,
is attributed to the fact that the model equation passes through the practical JC parameters (except the static yield strength) were adopted
origin of 𝑤𝑐 versus 𝐼0 , which is a shortcoming of the model since one from [8,17], as mentioned already. Curry, in [17], pointed out that
expects a negative intercept on the 𝑤𝑐 -axis (note that 𝑤𝑐 is a plastic when quasi-static stress–strain data are fitted to JC model, the resulting
displacement). static yield stress (𝐴𝐽 𝐶 ) underestimates the value observed experimen-
From the other hand, it should be noted that Nurick et al. [13] tally, due to the hardening power-law in JC model. Although, the
measured (indirectly) the total impulse, and the actual specific impulse used material parameters were found to give reasonable predictions
distribution (which is assumed uniform herein) was not reported. The when compared to the experiments. Example techniques to determine
slight discrepancy between the present model prediction and the ex- such parameters are explained in [17,74], and validated JC material
perimental data could then be attributed to the implied lack of precise parameters for some commonly used ductile materials are given in [8,
data. 26,61,63,66].

8
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Fig. 7. Comparison of model predictions to experimental data [13] and LS-DYNA numerical results in terms of central residual displacement 𝑤𝑐 of (a) rectangular and (b) square
membranes under uniform impulse of total magnitude of 𝐼0 . LS-DYNA data are for permanent displacements. The data was separated for visual convenience.

using beam elements along the radial axis of the plates using *Elem-
ent_Beam and *Section_Beam keywords with Elform = 8 and three through-
thickness integration points, IR/QR = 3. With this set-up, it is unneces-
sary to specify the conditions at the axis of symmetry. Because the ax-
isymmetric solver was utilised, the *Mat_Modified_Johnson_Cook model
was used; it should be noted that the ‘‘modified’’ and ‘‘simplified’’
Johnson–Cook models differ in describing the strain-rate sensitivity.
Again, applicable material parameters (except the static yield strength)
were adopted from [8,17]. Motion is induced by a uniform initial
velocity field calculated from the uniform specific impulse, and the
periphery node was fully restrained. Finally, the central displacement
was chosen to carry out mesh convergence study to determine the
appropriate mesh density, which was then maintained throughout.
The finite element results for the permanent displacement are com-
pared to the present analytical model and the experimental data of
Nurick et al. [13], as shown in Fig. 8.
Gharababaei and Darvizeh [14] report results from 86 experiments
on circular plates. The authors measured the central permanent dis-
placements of steel, copper, and aluminium thin plates and the total
impulses they were subjected to using a ballistic pendulum. All spec-
Fig. 8. Comparison of model predictions to experimental data [13] and LS-DYNA
imens had circular exposed areas with a fixed diameter of 100 mm.
results in terms of central residual displacement 𝑤𝑐 of circular membranes under Among the experiments, 42 tests are assumed to generate spatially
uniform impulse of total magnitude of 𝐼0 ; 𝑘0 is the circular structural parameter defined uniform specific impulse based on the following. The blast loads were
in Eq. (5.10). generated by detonating thin (cylindrical) disks of C4 explosives lo-
cated at a stand-off distance of 300 mm from the plates’ centres for
the 42 tests. The smallest scaled distance was Z = 1.12 m/kg1/3 , and a
6.2. Circular model verification 9.5◦ angle of incidence at the plates’ periphery was held constant for
the tests. Furthermore, the authors used a rigid circular tube of equal
diameter as that of the specimens to guide the propagation of the shock
6.2.1. Uniform impulse waves along its axis. Further details can be found in the original paper.
In this section, we compare the predictions of the circular model for The predictions from the present model, Eq. (5.9), and LS-DYNA
uniform specific impulse, as given by Eq. (5.9), to experimental data simulations are compared to the experimental data of Gharababaei and
and numerical LS-DYNA predictions. Nurick, Martin, and Pearce [13], Darvizeh. Again, the static yield strengths reported in experiments are
discussed in Section 6.1.1, also report experimental data for the per- taken as the characteristic yield strengths in the model calculations.
manent displacements of circular membranes of fixed diameters of The data for the aluminium plates were excluded due to numerical
100 mm when subjected to uniform impulses of varying amplitudes. difficulties in simulating their behaviour as the material is not strain-
The membrane material properties and thickness are as described for rate sensitive, and there is no available material data given in [14]
regarding its strain-hardening parameters. The results are given in
the rectangular membrane, see Section 6.1.1. The data are compared
Fig. 9.
to predictions from the model for the uniformly loaded circular mem-
As shown in Fig. 9, the model does not accurately predict the
branes, Eq. (5.9). The characteristic yield strength 𝜎0 is taken as the
experimental outcomes of Gharababaei and Darvizeh [14], in particular
static yield strength in [13]. The results are shown in Fig. 8. for 𝑘0 × 𝐼0 larger than 0.03 m. The model overpredicts the permanent
In addition, LS-DYNA was used to replicate the tests of Nurick displacement by at most a factor of two.
et al. [13] for the circular membrane case. Axisymmetric conditions Fig. 10 compares the model to data from LS-DYNA alone, which
were assumed in the simulations, and hence the problems were solved combines data already shown in Figs. 8 and 9. In the validation data

9
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Fig. 9. Comparison of model predictions to experimental data [14] and LS-DYNA Fig. 11. Comparison of model predictions to experimental data [15] in terms of
results in terms of central residual displacement 𝑤𝑐 of circular membranes under central residual displacement 𝑤𝑐 of circular membranes under uniform impulse of total
uniform impulse. In [14], the total impulse 𝐼0 was measured by ballistic pendulum; magnitude of 𝐼0 . 𝑘0 is the circular structural parameter defined in Eq. (5.10).
however, it was generated from thin disks of explosives with charge-to-target radius
ratios of 0.2 and 0.30, but with a stand-off distance to plate’s diameter ratio of 3.0.
Experimental data for specimens with diameter-to-thickness ratio, (𝐷∕ℎ), less than 50
are highlighted with filled red markers. 𝑘0 is the circular structural parameter defined MPa, respectively. These specimens were clamped with sharp edges.
in Eq. (5.10).
Two additional sets of specimens with filleted clamping supports had
diameters of 100 mm and yield strength of 251 MPa each. All data,
including the total impulses and central residual displacements, are
tabulated in Nurick et al. [15], and further details can be found
there. Their data, excluding five tests for which displacements are not
reported, are used to verify the circular membrane solution, Eq. (5.9),
as shown in Fig. 11.
So far, the experimental tests used for the validation involve limited
ranges of plates’ thicknesses. To further assess the performance of the
model under combined variations of total impulse and plate thickness
in broader ranges, the input data presented in Rigby et al.’s [8] was
used in additional LS-DYNA simulations. The set-ups are similar to
those described earlier. From a sample simulation run, the central
displacement time history is shown in Fig. 12. The results of the latter
investigation are depicted in Fig. 13, which demonstrates the accuracy
of the present model under varying plate thicknesses.

6.2.2. Discussion and limitations

The discussion and limitations presented in Section 6.1.2 for the
rectangular case are also applicable to the validation of the uniform
circular model against the experiments of Nurick et al. [13,15]. In
general, the model slightly deviates from the experimentally observed
Fig. 10. Comparison of model predictions to results from LS-DYNA simulations measurements.
corresponding to the experimental set-ups in [13,14] in terms of central residual However, there is a pronounced discrepancy when comparing the
displacement 𝑤𝑐 of circular membranes under uniform impulse of total magnitude uniform model to Gharababaei and Darvizeh’s [14] data. This required
of 𝐼0 . Numerical material parameters (except the static yield strength) were taken
a reappraisal of the model, the input data, and the testing set-up to
from [8] for the steel plates. Notice that LS-DYNA simulations account for elasticity,
strain-hardening, strain-rate sensitivity, bending and shear effects, which the present trace the source of discrepancy. As previously discussed, the authors
analytical model completely ignores. 𝑘0 is the circular structural parameter defined in used a rigid tube as a blast wave guide that (when combined with the
Eq. (5.10). moderate stand-off distance) could produce wave reflection effects due
to the interaction near the tube’s wall. However, the assumption of
uniformity of the impulse is ruled out. If the resulting specific impulse
in Fig. 10, yield strength, plate thickness and density are variables in was indeed non-uniform, then according to [5,8], the displacement
addition to the total impulse. As the model properly captures the trend would then be even larger than that induced by a uniform impulse.
of LS-DYNA predictions, it is concluded that the functional form of 𝑤𝑐 Although the model is built to predict membrane behaviour of thin
in terms of these structural parameters is accurate. plates, it performs better for the relatively thicker plates in [14] than
Nurick, Gelman and Marshall [15] present additional 162 exper- it does for specimens with diameter-to-thickness ratios greater than 50.
iments on circular steel plates, of varying diameters, loaded with This can be seen in Fig. 9, in which the markers for thicker specimens
uniform impulses. All had 1.6 mm thickness and assumed density of are filled in red to highlight such observation. Furthermore, the trend
7850 kg/m3 . For plates with diameters 60, 80, 100, 120 mm, the of the experimental data deviates from the expected behaviour that
corresponding static yield strengths were given as 251, 220, 270, 220 moderately thicker plates are stiffer. This is because: (1) extra resisting

10
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Table 3
Comparison of the present model against the modified Nurick and Martin’s model
proposed by Yuen et al. [19].
Geometry: 𝛼 𝛽
Yuen et al. Present Yuen et al. Present
Circular 0.241 0.281 0.298 0.0
Rectangular 0.253 0.480 × 𝑆0 (𝐿𝑦 ∕𝐿𝑥 )1∕2 −0.158 0.0
Square 0.253 0.270 −0.158 0.0

Furthermore, similar to the discussion in Section 6.1.2, the authors

in [14] measured the total impulse indirectly using a ballistic pendu-
lum. Rigby et al. [55] pointed out that targets typically experience
∼ 67% of the total impulses calculated from the ballistic pendulum
measurements. Their conclusion can explain the observed discrepancy
between the model predictions and Gharababaei and Darvizeh’s [14]
data.

6.3. Comparison of uniform circular and rectangular models to the modified

Fig. 12. Time history of the central displacement 𝑤𝑐 from LS-DYNA simulation for a
Nurick and Martin model
sample problem using input parameters given in [8]. Negative values on the ordinate
𝑦-axis correspond to downward displacements. It is seen that the resulting displacement
is mainly plastic. Yuen et al. [19] proposed modifications to Nurick and Martin’s [54]
model. In their model, the normalised permanent displacement (i.e.
displacement-thickness ratio), 𝑤𝑐 ∕ℎ, of circular or rectangular thin
plates with exposed areas, 𝐴, is linearly related to the non-dimensional
̂
impulse, 𝜙,
𝐼0
𝜙̂ = √ ,
ℎ2 𝜌𝜎0 𝐴
by the following empirical relation, see [19],
𝑤𝑐
= 𝛼 𝜙̂ + 𝛽
ℎ
where 𝛼 and 𝛽 are correlation coefficients given in Yuen et al. [19] and
presented in Table 3 for convenience. Note that we divided 𝛼 (reported
√
in Yuen et al. [19]) by 𝜋 and 2, respectively, for the circular and
rectangular models to unify the form of 𝜙̂ for the two geometries, as
given above.
Their empirical models were shown to reasonably predict the re-
sponse of the thin plates compared to 699 (circular) and 356 (rect-
angular) experiments. Thus, it is of interest to compare the present
model, developed herein, to their findings. For the present model,
the coefficients 𝛼 and 𝛽 are obtained from the previously developed
̂
solutions by rewriting them in terms of 𝜙.
Fig. 13. Comparison of model predictions to LS-DYNA results in terms of central
The comparisons are summarised in Table 3. The values 𝑆0 = 0.5627
residual displacement 𝑤𝑐 of circular membranes under uniform impulse of total √
magnitude of 𝐼0 . In the numerical simulations, input parameters were obtained from and 𝐿𝑦 ∕𝐿𝑥 = 1 were used to calculate 𝛼 for the square geometry using
the parametric study in [8] to investigate the model performance for a wide range of the expression of 𝛼 for the rectangular model, given in Table 3. Note
plate’s thickness. Peak (p) and residual (r) displacements from LS-DYNA analyses are
that 𝛽 is zero for all cases of the present model. While Yuen et al.’s 𝛼
shown.
and 𝛽 are identical for the rectangular and square plates, the present
model’s 𝛼 depends on the aspect ratio, 𝐿𝑦 ∕𝐿𝑥 , (the parameter 𝑆0 also
depends on the aspect ratio). Since the value of 𝛽 from Yuen et al. is
modes are involved, i.e. bending and shear, that also contribute to
absorb (or dissipate) the initial kinetic energy; and (2) due to in- small and will be multiplied by the plate’s thickness to give a fraction
creased mass per unit area, thicker plates attenuate the initial velocity of 𝑤𝑐 , 𝛽 can be neglected for our purpose of comparison. Overall, the
generated by the (externally) imparted impulse according to Rigby present model overestimates the displacement by around 7∽17% (based
et al.’s [8]. on 𝛼) compared to the predictions of Yuen et al. [19].
However, since an increase in the impulse leads to increase in From the comparisons of the present model against data from
initial velocity (for a given plate), the dynamic yield stress increases experiments and LS-DYNA simulations, presented in Sections 6.1.1 and
(due to strain-rate effects), which in turn would reduce the permanent 6.2.1, it is concluded in overall that the model is reasonably accurate
displacement. Thus, the discrepancy between the model predictions and in predicting the permanent displacements of rectangular and circular
the experimental displacements for the range of larger impulses can membranes under uniformly distributed specific impulses. The compar-
be partly attributed to neglecting the strain-rate sensitivity of the yield ison of the uniform model to the already validated predictions of Nurick
strength and possibly, as discussed in Section 6.1.2, the work-hardening and Martin [54] and Yuen et al. [19], presented in this section, also
of the specimens (in particular made of copper). supports this conclusion.

11
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

7. Limitations 8. Summary and conclusions

Analytical solutions were developed to predict the profile and peak

Despite the relative accuracy of the present model, it has several permanent displacement of thin plates under impulsive blast loads.
shortcomings. The plates’ materials were assumed rigid-perfectly plastic, obeying
First, during the development of the equation of motions, the mem- von Mises’s yield criteria, and their motions were initiated from the
brane is assumed to be already yielding as long as there is motion. That blast-induced specific impulse. This problem set-up led to a monotonic
is, there is no rigid body motion. The whole membrane experiences deformation path, which was exploited in deriving the governing equa-
(non-uniformly) distributed plasticity. tion of motion systematically through the application of the extended
Hamilton’s principle. The equations of motion apply to thin plates
Second, ‘‘plastic’’ motion is assumed to be initiated by initial ve- deforming mainly in membrane modes and plastically. Although the
locity, which is directly given by the imposed specific impulse. Hence, obtained equations are general, solutions were given for two membrane
the model cannot determine the amount of impulse to be elastically geometries, rectangular and circular, by the modal decomposition tech-
absorbed. Therefore, any amount of given impulse leads immediately nique. The modal decomposition was supplemented with sequential
to an onset of ‘‘plastic’’ deformation. The model indicates that kinetic mode terminations, a strategy that is justified herein by Drucker’s
energy is at maximum at the initial time state, which gradually and requirement of plastic work non-negativity.
monotonically decays as plasticity evolves. The model assumes com- The rectangular solution applies to any spatial distribution of spe-
cific impulses. A practical MATLAB code is proposed to efficiently
plete removal of the external load (that induced the specific impulse)
compute the total modal impulses involved in the general solutions for
before motion, and hence the decreasing rate of kinetic energy is
non-uniform impulses. Moreover, a procedure to estimate the errors
identical to the increasing rate of plastic work; the time history of from truncating the infinite series in the solutions was discussed in
plastic work is of an inverted form of the kinetic energy history. connection to the concept of the upper bound kinetic energy of Rigby
Third, the yield strength is independent of the evolution of the mem- et al. [8]. The circular solution is restricted to axisymmetric specific
brane’s motion. That is, enhancement (or hardening) of the yield stress impulses.
due to strain-hardening and strain-rate sensitivity are not (correctly) However, it was possible to provide simple closed-form solutions
incorporated. This is due to the perfect plasticity assumption, and a for the case of uniform specific impulse. The closed-form solutions of
rectangular and circular thin plates were verified against experimental
detailed discussion was given earlier in Section 6.1.2.
data in the literature and results from LS-DYNA simulations performed
Fourth, the displacement components along the in-plane coordi- by the authors.
nates of the undeformed membrane are assumed negligible as com- The present models for the rectangular and circular membranes
pared to the out-of-plane component. In fact, the in-plane components were shown to be reasonably accurate in comparison to the experi-
are set identically to zero throughout. Hence, the longitudinal Green– mental and numerical results. It should be emphasised that the models
Lagrange strains are merely due to the quadratic terms of the gradient account only for simplified idealisations: an impulsive blast load, a
of the out-of-plane displacement. This, in turn, leads to zero membrane rigid-perfectly plastic material behaviour, and a membrane mode of de-
strains at points where zero gradients occur, such as at the centre formation. As a result, the obtained models are believed to be (already
validated) simple and fast-running tools. Thus, they can be used by
of a symmetrically loaded membrane. Experiments often show that
structural blast engineers for probabilistic-based analyses. Although the
this is not true; in fact, the peak membrane strains in uniformly and
analytical models were compared against cases with uniform imparted
impulsively loaded membranes are located near their central regions. impulsive loads, the general relations derived herein can readily model
Despite this shortcoming, the model is able to predict thinning of non-uniform impulsive loads.
the membrane near its restraints (i.e. where it is supported). This is
directly due to the plastic incompressibility, which necessarily gives CRediT authorship contribution statement
the transverse longitudinal strain as the negative of the sum of the
Saud A.E. Alotaibi: Formal analysis, Validation, Data curation,
other two longitudinal strains (which are always non-negative); thus,
Writing – original draft, Visualization. Samuel E. Rigby: Conceptu-
thickness shortens in strained regions.
alization, Methodology, Supervision, Software, Writing – review &
Finally, the model focuses on thin plates subjected to extreme editing, Validation. Maurizio Guadagnini: Supervision, Methodology,
out-of-plane blast loads. The model excludes flexural and direct-shear Writing – review & editing. Andrew Tyas: Conceptualization, Method-
mechanisms. The model should, then, be limited to thin plates with ology, Supervision, Resources.
high span-to-thickness ratios and subjected to intense lateral forces.
Declaration of competing interest
Therefore, the model can only predict failure in the form of thinning-
induced splitting. Thus, the model does not apply to cases where
Saud A. E. Alotaibi reports financial support was provided by Saudi
direct shear failures near the supports occur. Note that as the load Arabian Cultural Bureau - London.
intensity becomes very low or the load duration increases significantly,
the model is not suitable since elastic and flexural effects become Data availability
important.
Data will be made available on request.
Attempts by the authors to incorporate as many sophistications as
needed to replicate the actual behaviour of the studied membrane Acknowledgements
response would eventually lead to an incremental analysis of a non-
linear dynamic problem. We did not attempt to re-discover an already The authors are grateful to the contributions of others (cited
existing problem which can be tackled by existing techniques (e.g. FE throughout the text) that gave technical enlightenment to accomplish
analysis using LS-DYNA). Instead, we aimed to simplify the problem the work presented herein. The authors thank the Department of
Civil and Structural Engineering and the Blast and Impact Engineering
using overall observations that capture the ‘‘dominant’’ behaviour to
Research Group at the University of Sheffield. The first author also
ultimately obtain a model that rationally balances accuracy, simplicity,
thanks Qassim University, Saudi Arabia, and the Saudi Arabian Cultural
and efficiency. Bureau, London, UK, for their financial support. Finally, the authors are
Further justifications of the model’s assumptions and the aim of the grateful to the dedicated effort and technical enrichment provided by
work were given earlier in Sections 1 and 2. the reviewers to enhance the article.

12
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Appendix A. Assessment of error due to series truncation then, one reaches the following general condition
∑
∞
Unless the specific impulse distribution identically matches the cos2 𝜃𝑚𝑛 = 1. (A.8)
shape of a particular mode, the exact solution, as given by Eqs. (4.8) 𝑚,𝑛=1
and (4.9), requires taking an infinite number of terms in the sum.
It is vital to recognise that there is no single component in the
However, in practice, a finite number of terms is used to approximate
infinite series above, Eq. (A.8), with a magnitude larger than unity since
the solution within reasonable accuracy. An appropriate measure to
each term is always positive while the right-hand side is one. Hence, if
evaluate the sufficiency of the approximation is the total kinetic energy.
one term is identically one, then all other terms must vanish, which is
Suppose that the kinetic energy computed by including a finite number
of modes is close to the ‘‘exact’’ kinetic energy, which is a strict upper the situation when the specific impulse field matches the shape of the
bound due to the non-negativity of kinetic energies. Then, the discarded surviving (or resonating) mode.
modes will be insignificant as their total contribution is bounded from Finally, if only a finite number of terms is used to calculate 𝑤𝑝 or 𝑤𝑐 ,
above by the implied error (or difference). It is the initial kinetic energy then it is sufficient to verify that the sum of cos2 𝜃𝑚𝑛 , for those modes
that is referred to, which is included, is close to unity (from below). In other words, the error, 𝜖,
𝐿𝑥 𝐿𝑦 due to truncation at 𝑚 = 𝑀 and 𝑛 = 𝑁, can be estimated as
1
𝐸𝑘 = 𝜌ℎ𝑤̇ 20 𝑑𝑦 𝑑𝑥 ∑
𝑀,𝑁 ∑
∞
∫0 ∫0 2
𝐿𝑥 𝐿𝑦
𝜖𝑀𝑁 = 1 − cos2 𝜃𝑚𝑛 ≡ cos2 𝜃𝑚𝑛 . (A.9)
1 2 𝑚,𝑛=1 𝑀,𝑁
= 𝑖(𝑥, 𝑦) 𝑑𝑦 𝑑𝑥. (A.1)
2𝜌ℎ ∫0 ∫0
The quantity cos 𝜃𝑚𝑛 , on its own, gives the absolute physical im-
Again, this is taken as the ‘‘exact’’ kinetic energy at time 𝑡 = 0, and it is portance of the (𝑚, 𝑛)th mode in relation to all other modes since it
termed the ‘‘upper bound kinetic energy’’ uptake in Tyas and Pope [5] compares the modal energy (through 𝐼𝑚𝑛 ) to the strict upper bound
and Rigby et al. [8]. kinetic energy, 𝐸𝑘 , (through 𝐼𝑘 ).
Due to the assumption of deformation monotonicity, the total plastic
work is associated with the total strain at the final time (as the problem Appendix B. A practical method to compute 𝑰𝒎𝒏
is path-independent). That is, after 𝑡 ≥ 𝑡1,1 , the plastic work 𝑊𝑝∗ is
evaluated using the final displacement field, 𝑤𝑝 (𝑥, 𝑦), and hence is given
For smoothly varying distribution of specific impulse, 𝐼𝑚𝑛 decreases
by
[ ( as 𝑚 and 𝑛 increase due to cancellations associated with high spatial
𝐿𝑥 𝐿𝑦 )2 ( )2 ]
2 1 𝜕𝑤𝑝 1 𝜕𝑤𝑝 oscillations. Thus, practically a finite number of modes suffices to
𝑊𝑝∗ = √ 𝜎0 ℎ + 𝑑𝑦 𝑑𝑥. (A.2)
∫0 ∫0 2 𝜕𝑥 2 𝜕𝑦 estimate the displacement accurately.
3
In the expression for 𝑤𝑐 , instead of carrying out the numerical
By exploiting the modal orthogonality property and after some lengthy integrations directly by the trapezoidal rule, it is observed that the two-
algebraic simplifications, the plastic work evaluates to dimensional Fast Fourier Transform (FFT) could be utilised to reduce
1 4 ∑∞ the computational time. However, according to the form of 𝐼𝑚𝑛 , the
𝑊𝑝∗ = 𝐼2 . (A.3) specific impulse distribution should be slightly manipulated first. The
2𝜌ℎ 𝐿𝑥 𝐿𝑦 𝑚,𝑛=1 𝑚𝑛
procedure is quite simple to derive, and it is briefly described below
The above, Eq. (A.3), is precisely the expression for the (initial) and followed by a practical Matlab code for the implementation.
kinetic energy if we would evaluate it by time differentiating the From the actual specific impulse, 𝑖(𝑥, 𝑦), defined on the actual
general solution 𝑤(𝑥, 𝑦, 𝑡), Eq. (4.1), then set 𝑡 = 0, which confirms that membrane that spans the domain [0, 𝐿𝑥 ] × [0, 𝐿𝑦 ], construct a fictitious
the initial kinetic energy is converted into (plastic) internal energy. specific impulse 𝑖∗ (𝑥, 𝑦) that covers an extended rectangular region
Now, when the expressions of the exact initial kinetic energy, 𝐸𝑘 , [−𝐿𝑥 , 𝐿𝑥 ] × [−𝐿𝑦 , 𝐿𝑦 ], which is defined as
from Eq. (A.1), and the final plastic work, 𝑊𝑝∗ , are set equal, the
following condition, known as Parseval’s formula [75], is obtained ⎧𝑖(𝑥, 𝑦), (𝑥, 𝑦) ∈ [0, 𝐿𝑥 ] × [0, 𝐿𝑦 ],
⎪
∑
∞
𝐿𝑥 𝐿𝑦 𝐿𝑥 𝐿𝑦 ⎪−𝑖(−𝑥, 𝑦), (𝑥, 𝑦) ∈ [−𝐿𝑥 , 0] × [0, 𝐿𝑦 ],
2
𝐼𝑚𝑛 = 𝑖(𝑥, 𝑦)2 𝑑𝑦 𝑑𝑥 = ‖ ‖2
‖𝜙𝑚𝑛 ‖ ⋅ ‖𝑖‖
2
(A.4) 𝑖∗ (𝑥, 𝑦) = ⎨ (B.1)
4 ∫0 ∫0 ⎪−𝑖(𝑥, −𝑦), (𝑥, 𝑦) ∈ [0, 𝐿𝑥 ] × [−𝐿𝑦 , 0],
𝑚,𝑛=1
√ ⎪𝑖(−𝑥, −𝑦), (𝑥, 𝑦) ∈ [−𝐿𝑥 , 0] × [−𝐿𝑦 , 0].
⎩
where ‖𝑓 (𝑥, 𝑦)‖ = ∫𝐴 𝑓 (𝑥, 𝑦)2 𝑑𝐴, is the norm of a function 𝑓 , and
√ Now, the real components of the two-dimensional discrete Fourier
‖𝜙𝑚𝑛 ‖ is (𝐿𝑥 𝐿𝑦 ∕4).
‖ ‖ transform of 𝑖∗ (𝑥, 𝑦) is denoted by 𝑏𝑚𝑛 and given by
Therefore,
( ) ( )
( )2 1
𝐿𝑥 𝐿𝑦
2 𝑚𝜋𝑥 2𝑛𝜋𝑦
∑∞
𝐼𝑘 𝑏𝑚𝑛 = 𝑖∗ (𝑥, 𝑦) sin sin 𝑑𝑦 𝑑𝑥. (B.2)
2
𝐼𝑚𝑛 = (A.5) 𝐿𝑥 𝐿𝑦 ∫−𝐿𝑥 ∫−𝐿𝑦 2𝐿𝑥 2𝐿𝑦
𝑚,𝑛=1
2
It is important to note that the 𝑥-interval is 2𝐿𝑥 and that along 𝑦 is 2𝐿𝑦 .
in which 𝐼𝑘 is the energy-equivalent impulse due to Rigby et al. [8],
Next, 𝑏𝑚𝑛 = 𝑏(𝑚, 𝑛) is related to the complex Fourier coefficients by
which, for a rectangular target with loaded area 𝐴 = 𝐿𝑥 𝐿𝑦 , is given by
√ 𝑏𝑚𝑛 = −𝑐(𝑚, 𝑛) + 𝑐(−𝑚, 𝑛) + 𝑐(𝑚, −𝑛) − 𝑐(−𝑚, −𝑛) (B.3)
𝐼𝑘 = 𝐴 𝑖(𝑥, 𝑦)2 𝑑𝐴 (A.6) where 𝑚 and 𝑛 are indices corresponding to positive integers, and −𝑚
∫𝐴
√ and −𝑛 are indices corresponding to negative integers. Then, using the
= 𝐴 × ‖𝑖‖ . piecewise definition of 𝑖∗ (𝑥, 𝑦), it can be shown that 𝑏𝑚𝑛 reads
It should be noted that in [8], 𝐼𝑘 was derived directly from the 𝐿𝑥 𝐿𝑦 ( ) ( )
4 𝑚𝜋𝑥 𝑛𝜋𝑦
physical problem using equivalency between kinetic energies due to 𝑏𝑚𝑛 = 𝑖(𝑥, 𝑦) sin sin 𝑑𝑦 𝑑𝑥
𝐿𝑥 𝐿𝑦 ∫0 ∫0 𝐿𝑥 𝐿𝑦
non-uniform and uniform specific impulses.
4𝐼𝑚𝑛
Now, defining an angular parameter, which measures how much = . (B.4)
the actual specific impulse field is along the direction of one particular 𝐿𝑥 𝐿𝑦
mode (in the inner product sense), by Finally,
2𝐼 𝐿𝑥 𝐿𝑦
cos 𝜃𝑚𝑛 = 𝑚𝑛 (A.7) 𝐼𝑚𝑛 = [−𝑐(𝑚, 𝑛) + 𝑐(−𝑚, 𝑛) + 𝑐(𝑚, −𝑛) − 𝑐(−𝑚, −𝑛)] . (B.5)
𝐼𝑘 4

13
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

The values of 𝑐(𝑚, 𝑛), 𝑐(−𝑚, 𝑛), 𝑐(𝑚, −𝑛), and 𝑐(−𝑚, −𝑛) are the stan-
dard outputs of an FFT, to within a constant multiplier. FFT gives the
amplitudes of the modes, e.g. 𝑐(0, 0) reflects the amplitude of the mode
associated with 𝑚 = 0 and 𝑛 = 0. In Matlab, this is achieved using
the built-in function fft2(). Note that in Matlab, the output sorts the
number of modes in each direction in a special order: first, the zeroth
mode (which we do not need), followed by positive modes in ascending
order, and lastly negative modes in descending order.
The specific procedure to be implemented in Matlab is what follows.
Let the actual specific impulse, 𝑖(𝑥, 𝑦), be stored in 2D array I. Further,
let the plate lengths along 𝑥 and 𝑦 be Lx and Ly, respectively.
Then, denote 𝑖∗ (𝑥, 𝑦) by Istar, which can easily be formed in Matlab
using the built-in function flip(). Finally, say the total modal impulse,
𝐼𝑚𝑛 , will be stored in the array Imn and be evaluated using FFT. The
Matlab procedure is given in Script 1.

%% start of script Fig. C.1. Normalised residual displacement 𝑥𝑟 ∕𝑥𝑦,0 of elasto-plastic SDOF as function
% I =specific impulse matrix (2D array)...
of the ratios of load duration to elastic period 𝑡𝑑 ∕𝑇𝑒 and hardening modulus to elastic
% rows of I --> variation along x
stiffness 𝐻𝑘 . The grey curves are for intermediate ratios of 𝐻∕𝑘, which are bounded
% cols of I --> variation along y
by the values of the blue and red curves.
Istar=[flip(flip(I,1),2),-flip(I,1);-flip(I,2),I];

C0=fft2(Istar);
C0=real(C0); where 𝑘 and 𝐻 are the elastic stiffness and the hardening modulus;
C0=1/(numel(Istar))*C0; %undo multiplier 𝑅𝑦,0 , 𝑥𝑦,0 , and 𝑡𝑦,0 are the initial yield force, the corresponding yield
C1=C0(2:end,2:end); %skip zeroth mode
displacement, and the time at which first yielding occurs. The maxi-
B=zeros(ceil((size(C1)-1)/2));
mum displacement is 𝑥𝑚 . The last expression for 𝑅 describes the elastic
for i=1:size(B,1) unloading for sufficiently small-time interval after 𝑡 = 𝑡𝑚 .
for j=1:size(B,2) The applied force 𝐹 is assumed to be a rectangular pulse with
B(i,j) = -C1(i,j)+C1(end-i+1,j)+C1(i,end-j+1)-C1(end-i+1,end-j+1);
end
amplitude 𝐹0 and duration 𝑡𝑑 .
end The elasto-plastic response of the SDOF was solved numerically
using an explicit time integration scheme. The solution is terminated
Imn=Lx*Ly*B/4;
after completing a half cycle from the first occurrence of maximum
% Imn stores total modal impulses, Imn. displacement (i.e. after sufficient time during the elastic unloading/re-
% Use this directly in expressions for displacement w_p or w_c. bound), in order to obtain sufficient response to compute the residual
%% end of script.
displacement, 𝑥𝑟 .
Script 1: Matlab script to compute 𝐼𝑚𝑛 efficiently using FFT √We considered the following input data. The elastic period 𝑇𝑒 =
2𝜋 𝑚∕𝑘 is 0.0811 s. The ratio of the external force to the initial yield
limit 𝐹0 ∕𝑅𝑦,0 is 5.0. The ratio of load duration to the elastic period 𝑡𝑑 ∕𝑇𝑒
By comparing the results (not shown herein) from computing 𝐼𝑚𝑛 and the ratio of hardening modulus to elastic stiffness 𝐻∕𝑘 are varied
via the trapezoidal rule for a large number of modes in each direction independently; some practical ratios of 𝐻∕𝑘 are considered. The results
to those using the FFT, the observations are FFT is superiorly efficient for the normalised residual displacement 𝑥𝑟 ∕𝑥𝑦,0 are given in Fig. C.1.
and very reasonably accurate. Hence, the use of FFT to evaluate 𝐼𝑚𝑛 The Rigid-plastic response of the SDOF is obtained by substituting
is recommended. This computational strategy is rarely pointed out in 𝑅 = 𝑅𝑦,0 in the equation of motion, which is valid since 𝐹0 > 𝑅𝑦,0 , and
the literature as an efficient method to compute the modal amplitudes hence it is expected that 𝑥̇ ≥ 0. Therefore, the maximum displacement
appearing in the displacement response of plates. is
𝑅𝑦,0
Appendix C. Effect of work-hardening on the response of an SDOF 𝑥𝑚,𝑟𝑝 = 𝑥𝑡𝑑 + 𝑥̇ 𝑡𝑑 (𝑡𝑚 − 𝑡𝑑 ) − (𝑡 − 𝑡𝑑 )2
and plate 2𝑚 𝑚
where,
( )
𝐹0 − 𝑅𝑦,0 2
C.1. Forced response - SDOF 𝑥 𝑡𝑑 = 𝑡𝑑 , (C.1)
2𝑚
𝐹 𝑡 𝑅𝑦,0 𝑡𝑑
In this section, an assessment of the effect of work-hardening is 𝑥̇ 𝑡𝑑 = 0 𝑑 − , (C.2)
presented based on the response of a single-degree-of-freedom (SDOF). 𝑚 𝑚
𝑥̇ 𝑡
The system consists of a mass 𝑚, a massless spring with resistance 𝑅, 𝑡𝑚 = ( 𝑑 ) + 𝑡𝑑 (C.3)
𝑅
and an applied dynamic force, and its general equation of motion is 𝑦,0
𝑚
𝑚𝑥̈ + 𝑅 = 𝐹 The maximum rigid-perfectly plastic displacement 𝑥𝑚,𝑟𝑝 is plotted
For simplicity, we are concerned with motion up to a half cycle in Fig. C.1 as a function of the ratio of load duration to elastic period
beyond the first maximum response. The SDOF is assumed as initially 𝑡𝑑 ∕𝑇𝑒 and compared to the residual displacement of the elasto-plastic
at rest. problem.
The resistance 𝑅 is of bilinear form during loading and has an elastic It can be seen, from the figure, that the rigid-perfectly plastic solu-
unloading tion gives reliable predictions as compared to the elasto-plastic solution
with various work-hardening, in particular as the load duration to
⎧𝑘𝑥 𝑡 ≤ 𝑡𝑦,0 ,
⎪ elastic period ratio 𝑡𝑑 ∕𝑇𝑒 becomes very small. In Fig. C.1, the difference
𝑅 = ⎨𝑅𝑦,0 + 𝐻(𝑥 − 𝑥𝑦,0 ) 𝑡𝑦,0 ≤ 𝑡 ≤ 𝑡𝑚 , between the responses for the rigid-perfectly plastic (black) and elasto-
⎪ plastic with 𝐻∕𝑘 = 0.0001 (blue) is attributed to elastic deformations
⎩𝑅𝑦,0 + 𝐻(𝑥𝑚 − 𝑥𝑦,0 ) + 𝑘(𝑥 − 𝑥𝑚 ) 𝑡𝑚 ≤ 𝑡

14
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

of the elasto-plastic SDOF. Such difference should decrease as the the residual displacement is
external force to initial yield limit (𝐹0 ∕𝑅𝑦,0 ) increases while small 𝑡𝑑 𝑇
𝑡 + 2𝑒
𝑚
is maintained. 2
𝑥𝑟 = 𝑥(𝑡) 𝑑𝑡
𝑇𝑒 ∫𝑡𝑚
which simplifies to
C.2. Impulsive response - SDOF [ ]
𝑅𝑦,0 + 𝐻(𝑥𝑚 − 𝑥𝑦,0 ) − 𝑘𝑥𝑚
𝑥𝑟 = −
𝑘
In this section, we present additional assessment of the effect of
work-hardening on the response of the same SDOF as in the previous The initial kinetic energy, 𝐸𝑘,0 , of the SDOF is
section. However, we consider the system where the external force is 1
𝐸𝑘,0 = 𝑚 𝑥̇ 20
absent and assume the response to be driven by initial velocity. The 2
initial displacement is assumed zero. and, the maximum (initial) elastic energy, 𝐸𝑒,𝑚 , is
We consider the initial velocity, 𝑥̇ 0 , to be large enough to cause 1
𝐸𝑒,𝑚 = 𝑘 𝑥2𝑦,0
initial yielding. 2
Denoting the time of initial yielding by 𝑡𝑦,0 , the response for (0 ≤ 𝑡 ≤ For a given ratio of 𝐸𝑘,0 ∕𝐸𝑒,𝑚 and a ratio of 𝐻∕𝑘, one can study the
)
𝑡𝑦,0 is governed by 𝑚𝑥̈ + 𝑘𝑥 = 0 with initial conditions 𝑥0 = 0 and response in terms of ratio of maximum (or residual) displacement 𝑥𝑚
𝑥̇ 0 > 0. The response is (or 𝑥𝑟 ) to the initial yield displacement 𝑥𝑦,0 . This gives direct assessment
√ of the influence of work-hardening.
𝑘
𝜔𝑒 = , The rigid-perfectly plastic solution is characterised by the maximum
𝑚
( ) response time 𝑡𝑚,𝑟𝑝 (when velocity reaches, and subsequently held
𝑥̇
𝑥(𝑡) = 0 sin 𝜔𝑒 𝑡 constant at, zero) and the corresponding maximum displacement 𝑥𝑚,𝑟𝑝 .
𝜔𝑒
These are defined by
From which, the state at 𝑡 = 𝑡𝑦,0 is 𝑥̇
𝑡𝑚,𝑟𝑝 = ( 0 ) ,
𝑅𝑦,0 𝑅 𝑦,0
𝑥𝑦,0 = , 𝑚
𝑘 (𝑥 )
𝑦,0 𝜔𝑒 1 𝑅𝑦,0 2
sin−1 𝑥𝑚,𝑟𝑝 = 𝑥̇ 0 𝑡𝑚,𝑟𝑝 − 𝑡
𝑡𝑦,0 =
𝑥̇ 0
, 2 𝑚 𝑚,𝑟𝑝
𝜔𝑒 The rigid-perfectly plastic response 𝑥𝑚,𝑟𝑝 can be compared to the
( )
𝑥̇ 𝑦,0 = 𝑥̇ 0 cos 𝜔𝑒 𝑡𝑦,0 response of the elasto-plastic with hardening solutions, the maximum
𝑥𝑚 or residual 𝑥𝑟 displacements, to assess the effects of both work-
Denoting the time at maximum response (i.e. just prior to elastic
( ) hardening and elasticity.
unloading) by 𝑡𝑚 , the response for 𝑡𝑦,0 ≤ 𝑡 ≤ 𝑡𝑚 is governed by 𝑚𝑥̈ + We consider a particular case for which the ratio of initial kinetic
𝑅𝑦,0 + 𝐻(𝑥 − 𝑥𝑦,0 ) = 0 with initial conditions 𝑥(𝑡𝑦,0 ) = 𝑥𝑦,0 and 𝑥(𝑡
̇ 𝑦,0 ) = 𝐸
energy to maximum elastic energy is large, 𝐸 𝑘,0 = 9.0. The elastic
𝑥̇ 𝑦,0 . The response is 𝑒,𝑚

√ period is 𝑇𝑒 = 0.0811 s. Also, the ratio of the hardening modulus to

𝐻 the elastic stiffness 𝐻∕𝑘 is taken as (0.01). The comparison of the time
𝜔𝐻 = ,
𝑚 history responses of the elasto-plastic with hardening SDOF and the
[ ( )] 𝑥̇ 𝑦,0 [ ( )] corresponding rigid-perfectly plastic SDOF is shown for the present case
𝑥(𝑡) =𝑥𝑦,0 cos 𝜔𝐻 𝑡 − 𝑡𝑦,0 + sin 𝜔𝐻 𝑡 − 𝑡𝑦,0 +
𝜔𝐻 in Fig. C.2. The displacement 𝑥 is normalised by 𝑥𝑦,0 , and the time axis
( )
𝑅𝑦,0 − 𝐻𝑥𝑦,0 { [ ( )]} is normalised by 𝑡𝑦 . The figure indicates that the rigid-perfectly plastic
(−1) 1 − cos 𝜔𝐻 𝑡 − 𝑡𝑦,0 .
𝐻 solution is suitable and highly efficient for response of the SDOF for the
considered ratio of 𝐸𝑘,0 ∕𝐸𝑒,𝑚 = 9.0.
By definition 𝑡𝑚 is the time when velocity becomes zero for the first
As a second example, we consider a problem with 𝐸𝑘,0 ∕𝐸𝑒,𝑚 = 25.0,
time, i.e. 𝑥(𝑡
̇ 𝑚 ) = 0. From which, the state at 𝑡 = 𝑡𝑚 is
𝑇𝑒 = 0.1147 s, and 𝐻∕𝑘 = 0.002. The responses from the elasto-
⎡ ⎤ plastic with hardening model and the rigid-perfectly plastic model are
⎢ 𝑥̇ 𝑦,0 ⎥
𝛽 = tan−1 ⎢ ( 𝑅 +𝐻𝑥 ) ⎥, compared in Fig. C.3.
⎢ 𝜔𝐻 𝑥𝑦,0 + 𝑦,0 𝑦,0
⎥
⎣ 𝑚𝜔𝐻 ⎦
𝛽 C.3. Plate response - LS-DYNA
𝑡𝑚 = + 𝑡𝑦,0 ,
𝜔𝐻
( ) In this section, results from LS-DYNA simulations are presented to
𝑥̇ 𝑦,0 𝑅𝑦,0 − 𝐻𝑥𝑦,0
𝑥𝑚 = 𝑥𝑦,0 cos (𝛽) + sin (𝛽) − [1 − cos (𝛽)] . assess the influence of work-hardening on the response of a ductile thin
𝜔𝐻 𝐻
plate loaded by a uniform impulse.
( )
For 𝑡 ≥ 𝑡𝑚 , the response is an elastic rebound, which is governed Flexural, shear, and membrane effects are all considered. Further-
by 𝑚𝑥̈ + 𝑅𝑦,0 + 𝐻(𝑥𝑚 − 𝑥𝑦,0 ) + 𝑘(𝑥 − 𝑥𝑚 ) = 0, with initial conditions more, a general elasto-plastic material behaviour is adopted. Plasticity
𝑥(𝑡𝑚 ) = 𝑥𝑚 and 𝑥(𝑡 ̇ 𝑚 ) = 0. The rebound response, for 𝑡 ≥ 𝑡𝑚 , is follows a von-Mises yield function in which the yielding stress is given
[ ] by Johnson–Cook (JC) model. In this study, the hardening saturation
𝑥(𝑡) =𝑥𝑚 cos 𝜔𝑒 (𝑡 − 𝑡𝑚 ) + stress 𝐵 is varied while all other material (and geometry and loading)
( )
𝑅𝑦,0 + 𝐻(𝑥𝑚 − 𝑥𝑦,0 ) − 𝑘𝑥𝑚 { [ ]} parameters fixed.
(−1) 1 − cos 𝜔𝑒 (𝑡 − 𝑡𝑚 ) .
𝑘 The plate is impulsively loaded using a prescribed uniform initial
The residual (plastic) displacement, 𝑥𝑟 , is obtained from the above transverse velocity 𝑤̇ 0
rebound solution when the vibration terms (i.e. the cosine terms) are 𝑖0
𝑤̇ 0 =
eliminated. With the elastic period 𝜌ℎ
√ where the uniform specific impulse, density, and uniform thickness are
𝑚
𝑇𝑒 = 2𝜋 denoted by 𝑖0 , 𝜌, and ℎ.
𝑘

15
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Table C.4
JC material parameters.
Parameter Value Unit
𝐸 206 × 109 Pa
𝜌 7830 kg/m3
𝜈 0.29 –
𝐴 296 × 106 Pa
𝑛 0.5597 –
𝑚 0 –
𝑐 0.032 –
𝜀̇ 0 1.4 × 10−6 1/s

Fig. C.2. Normalised displacement 𝑥(𝑡)∕𝑥𝑦,0 of elasto-plastic SDOF (with work-

hardening) as function of normalised time 𝑡∕𝑡𝑦,0 ; hardening modulus to elastic stiffness
ratio is 𝐻∕𝑘 = 0.01. The response is for an impulsively loaded SDOF, i.e. with nonzero
( )
initial velocity and zero external force. The initial kinetic energy 𝐸𝑘,0 = 1∕2 𝑚 𝑥̇ 20 is
( )
nine times larger than the maximum elastic energy 𝐸𝑒,𝑚 = 1∕2 𝑘 𝑥2𝑦 . Horizontal (blue)
dashed lines are drawn at 𝑥𝑦,0 , 𝑥𝑟 , and 𝑥𝑚 for the elasto-plastic response, and the
additional (red) dashed line corresponds to 𝑥𝑚,𝑟𝑝 for the rigid-plastic response.

Fig. C.4. Variation of JC current yield stress 𝜎𝑦 versus effective (von-Mises) plastic
strain 𝜀p,eff at the quasi-static plastic strain rate 𝜀̇ 0 . Curves correspond to different
amount of strain-hardening in terms of 𝐵∕𝐴 ratio.

rates. On the other hand, two softening mechanisms (thermal softening

when 𝑚 is non-zero and/or damage-induced softening) can lower 𝜎𝑦 .
We assume that softening is absent for simplicity (i.e. 𝑚 = 0 and damage
threshold = ∞).
In the study, we considered the material input values as shown
in Table C.4. Three values of 𝐵 as multiples of 𝐴 are considered. In
particular, we consider ratios of 𝐵∕𝐴 of 1.7, 0.89, and 0.05. These ratios
Fig. C.3. Normalised displacement 𝑥(𝑡)∕𝑥𝑦,0 of elasto-plastic SDOF (with work- correspond to the (true) yield stress-effective plastic strain curves in
hardening) as function of normalised time 𝑡∕𝑡𝑦,0 ; hardening modulus to elastic stiffness Fig. C.4.
ratio is 𝐻∕𝑘 = 0.002. The response is for an impulsively loaded SDOF, i.e. with non-zero The other input parameters are the following. The plate is square
( )
initial velocity and zero external force. The initial kinetic energy 𝐸𝑘,0 = 1∕2 𝑚 𝑥̇ 20 is
( ) with a side length of 0.089 m and thickness of 1.6 mm. Additionally,
25 times larger than the maximum elastic energy 𝐸𝑒,𝑚 = 1∕2 𝑘 𝑥2𝑦 . Horizontal (blue)
a total impulse of −16.1 N m. is applied uniformly to the plate to give
dashed lines are drawn at 𝑥𝑦,0 , 𝑥𝑟 , and 𝑥𝑚 for the elasto-plastic response, and the
additional (red) dashed line corresponds to 𝑥𝑚,𝑟𝑝 for the rigid-plastic response.
a uniform initial transverse velocity of −162.24 m∕s. The input values
correspond to the square plate tests of Nurick et al. cited and discussed
in Section 4.2. All boundary nodes are restrained in all translational
and rotational degrees of freedom. A fine element mesh is selected. The
According to JC model without Voce’s hardening, the current yield plate is modelled with two-dimensional shell elements using ELFORM
stress 𝜎𝑦 is influenced by current effective plastic strain 𝜀eff , the instan- = 7 (selectively-reduced integrated Hughes–Liu co-rotational formula-
taneous effective plastic strain rate 𝜀̇ eff , and the absolute temperature tion) and 10 Gauss integration points to better represent the extent of
𝑇 using plasticity through the plate’s thickness.
[ ( )]
( ) 𝜀̇ eff 𝑐 Nodal, element, and energy databases are requested with fine sam-
𝜎𝑦 = 𝐴 + 𝐵 𝜀𝑛eff 1 + (1 − 𝑇 ∗𝑚 )
𝜀̇ 0 pling rates. The transverse (z-) displacement and energy histories are
post-processed in MATLAB. The residual displacement is computed by
where 𝑇 ∗ is
time integration averaging beyond the first peak response.
𝑇 − 𝑇𝑟
𝑇∗ = The history results obtained from LS-DYNA for the three ratios of
𝑇𝑚 − 𝑇𝑟 𝐵∕𝐴 are compared to the predictions of the present model. Namely,
in which 𝑇𝑚 and 𝑇𝑟 are the (absolute) melting and room temperatures. the displacements at the plate’s centre, the time to maximum response,
𝐵 and 𝑛 are power-law hardening parameters, 𝑐 is a strain-rate sensi- and the evolution of different energies are analysed. In the model
tivity parameter. 𝐴 is the initial quasi-static yield stress at the threshold predictions, the characteristic yield strength 𝜎0 is taken as the static
strain rate 𝜀̇ 0 . Note that 𝑛 is typically in the range 0 ≤ 𝑛 ≤ 1, and then 𝐵 yield strength 𝐴.
could be regarded as an increase in 𝜎𝑦 (from 𝐴) due to strain-hardening. The results are shown in Figs. C.5–C.8. It can be concluded that the
A further strengthening in the strain-hardened 𝜎𝑦 occurs at high strain higher the degree of strain-hardening the smaller the peak and residual

16
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

Fig. C.5. Central transverse displacement 𝑤 time history from LS-DYNA simulation (for Fig. C.7. Effect of the degree of strain-hardening (in terms of 𝐵∕𝐴 ratio) on the
𝐵∕𝐴 = 1.7) and the corresponding model prediction. transient and residual central displacements from LS-DYNA. The model prediction
assumes 𝜎0 = 𝐴.

Fig. C.6. Histories of global energies from LS-DYNA (for 𝐵∕𝐴 = 1.7), and the
corresponding kinetic (KE) and plastic internal (IE) energies predicted by the model. Fig. C.8. LS-DYNA maximum (peak) and residual (res.) displacements at the plate’s
From LS-DYNA, the internal (IE), kinetic (KE), total (TE), and hourglass energies and centre as functions of the degree of hardening (in terms of 𝐵∕𝐴).
external work (EW) are shown. The maximum response time 𝑡𝑚𝑎𝑥 predicted by the
model is shown as a vertical dashed line.

//dx.doi.org/10.3390/ICEM18-05364, URL https://www.mdpi.com/2504-3900/

2/8/471.
displacements are, as clearly depicted in Figs. C.7 and C.8. In addition, [5] Tyas A, Pope DJ. The energy take-up of panels subjected to near-field blast
Figs. C.5–C.7 indicate that the model predictions capture both the time loading. Int Conf Response Struct Extrem. Load. 2003;1(C).
evolution of residual displacement and internal and kineitic energies, [6] Jones N. Inelastic response of structures due to large impact and blast load-
ings. J Strain Anal Eng Des 2010;45(6):451–64. http://dx.doi.org/10.1243/
in particular as compared to the case with the largest 𝐵∕𝐴 ratio.
03093247JSA611, URL https://journals-sagepub-com.sheffield.idm.oclc.org/doi/
abs/10.1243/03093247JSA611.
[7] U.S. Army CoE. Structures to resist the effects of accidental explosions UFC
References 3-340-02. In: Unified facil. criteria 3-340-02. Tech. rep., U.S. Army Corps of
Engineers; 2008, p. 1796, URL http://dod.wbdg.org/.
[1] Rigby SE, Knighton R, Clarke SD, Tyas A. Reflected near-field blast pressure [8] Rigby SE, Akintaro OI, Fuller BJ, Tyas A, Curry RJ, Langdon GS, Pope DJ.
measurements using high speed video. Exp Mech 2020;60(7):875–88. http://dx. Predicting the response of plates subjected to near-field explosions using
doi.org/10.1007/s11340-020-00615-3. an energy equivalent impulse. Int J Impact Eng 2019;128(January):24–
[2] Rigby SE, Tyas A, Clarke SD, Fay SD, Reay JJ, Warren JA, Gant M, 36. http://dx.doi.org/10.1016/j.ijimpeng.2019.01.014, URL https://doi.org/10.
Elgy I. Observations from preliminary experiments on spatial and temporal 1016/j.ijimpeng.2019.01.014.
pressure measurements from near-field free air explosions. Int J Prot Struct [9] Reddy JN. Energy principles and variational methods in applied mechanics.
2015;6(2):175–90. http://dx.doi.org/10.1260/2041-4196.6.2.175. John Wiley & Sons, Inc.; 2002, p. 592. http://dx.doi.org/10.1016/0307-904x(85)
[3] Tyas A. Blast loading from high explosive detonation: What we know and don’t 90131-3.
know. In: 13th int. conf. shock impact loads struct. SILOS 2019, 14-15 December [10] Washizu K. Variational methods in elasticity and plasticity. Int. ser. monogr.
2019, Guangzhou, China, no. December. Guangzhou, China; 2019, p. 65–76. aeronaut. astronaut. div. I solid struct. mech., 2nd ed.. Vol. 9, (2):Exeter,
[4] Tyas A. Experimental measurement of pressure loading from near-field blast Great Britain: Pergamon Press; 1975, p. 70–1. http://dx.doi.org/10.1002/zamm.
events: Techniques, findings and future challenges. Proceedings 2018;2(8). http: 19840640121.

17
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

[11] Drucker DC. A definition of stable inelastic material. Tech. rep., Providence, [40] Martin JB. Convergence to mode form solutions in impulsively loaded piecewise
Rhode Island: Division of Engineering, Brown University; 1957, p. 23, URL linear rigid-plastic structures. Tech. rep., 2, Great Britain: Pergamon Press Ltd.;
https://apps.dtic.mil/sti/pdfs/AD0143756.pdf. 1983, p. 125–41.
[12] Drucker DC. Variational principles in the mathematical theory of plasticity. [41] Martin JB, Symonds PS. Mode approximations for impulsively loaded rigid plastic
Tech. rep., Providence, R. I.: DIVISION OF APPLIED MATHEMATICS, BROWN structures. 1965, p. 61.
UNIVERSITY; 1956, p. 24. [42] Hudson GE. A theory of the dynamic plastic deformation of a thin diaphragm.
[13] Nurick GN, Pearce HT, Martin JB. The deformation of thin plates subjected J Appl Phys 1951;22(1):1–11.
to impulsive loading. In: Bevilacqua L, Feijno R, Valid R, editors. Inelast. [43] Jones N. Impulsive loading of a simply supported circular plate. Tech. rep.,
behav. plates shells iutam symp. Rio Janeiro. Rio de Janeiro; 1985, p. 597–616. February, Division of Engineering Brown University Providence, Rhode Island;
http://dx.doi.org/10.1007/978-3-642-82776-1_29. 1967.
[14] Gharababaei H, Darvizeh A. Experimental and analytical investigation of large [44] Mihailescu-Suliciu M, Wierzbicki T. Wave solution for an impulsively loaded
deformation of thin circular plates subjected to localized and uniform impulsive rigid-plastic circular membrane. Arch Mech 2002;54(5–6):553–75.
loading. Mech Based Des Struct Mach 2010;38(2):171–89. http://dx.doi.org/10. [45] Wierzbicki T, Nurick GN. Large deformation of thin plates under localised
1080/15397730903554633. impulsive loading. Int J Impact Eng 1996;18(7–8):899–918. http://dx.doi.org/
[15] Nurick GN, Gelman ME, Marshall NS. Tearing of blast loaded plates with clamped 10.1016/s0734-743x(96)00027-9.
boundary conditions. Int J Impact Eng 1996;18(7–8):803–27. http://dx.doi.org/ [46] Cloete TJ, Nurick GN. On the influence of radial displacements and bending
10.1016/s0734-743x(96)00026-7. strains on the large deflections of impulsively loaded circular plates. Int J Mech
[16] LST. LS-DYNA | Livermore Software Technology Corp., URL https://www.lstc. Sci 2014;82:140–8. http://dx.doi.org/10.1016/j.ijmecsci.2014.02.026.
com/products/ls-dyna. [47] Mehreganian N, S. Fallah A, Louca LA. Inelastic dynamic response of square
[17] Curry R. Response of plates subjected to air-blast and buried explosions (Ph.D. membranes subjected to localised blast loading. Int J Mech Sci 2018;148:578–95.
thesis), University of Cape Town; 2017, p. 318. http://dx.doi.org/10.1016/j.ijmecsci.2018.09.017.
[18] Curry RJ, Langdon GS. Transient response of steel plates subjected to close [48] Mehreganian N, Fallah AS, Louca LA. Nonlinear dynamics of locally pulse loaded
proximity explosive detonations in air. Int J Impact Eng 2017;102:102–16. square Föppl–von Kármán thin plates. Int J Mech Sci 2019;163(September).
http://dx.doi.org/10.1016/j.ijimpeng.2016.12.004. http://dx.doi.org/10.1016/j.ijmecsci.2019.105157.
[19] Chung Kim Yuen S, Nurick GN, Langdon GS, Iyer Y, Yuen SCK, Nurick GN,
[49] Jones N. Dynamic inelastic response of strain rate sensitive ductile plates due
Langdon GS, Iyer Y. Deformation of thin plates subjected to impulsive load:
to large impact, dynamic pressure and explosive loadings. Int J Impact Eng
Part III – an update 25 years on. Int J Impact Eng 2016;107:108–17. http:
2014;74:3–15. http://dx.doi.org/10.1016/J.IJIMPENG.2013.05.003.
//dx.doi.org/10.1016/j.ijimpeng.2016.06.010.
[50] Ploch J, Wierzbicki T. Bounds for large plastic deformations of dynamically
[20] Nurick GN, Martin JB. Deformation of thin plates subjected to impulsive loading
loaded continua and structures. Int J Solids Struct 1981;17(2):183–95. http:
- A review Part I: Theoretical considerations. Int J Impact Eng 1989;8(2):159–70.
//dx.doi.org/10.1016/0020-7683(81)90074-3.
[21] Jones N. Structural impact. 2nd ed.. New York: Cambridge University Press;
[51] Ronter ARS, Martin JB. Bounds for impulsively loaded plastic structures. J Eng
2012, p. 584. http://dx.doi.org/10.1155/1999/124792.
Mech Div 1972;98(EM1):107–19.
[22] Chakrabarty J. Applied plasticity. 2nd ed.. l: Springer Science+Business Media;
[52] Rigby SE, Tyas A, Bennett T. Elastic-plastic response of plates subjected to
2010, p. 544. http://dx.doi.org/10.1007/b22134, arXiv:arXiv:1011.1669v3.
cleared blast loads. Int J Impact Eng 2014;66:37–47. http://dx.doi.org/10.1016/
[23] Jones N. Recent studies on the dynamic plastic behavior of structures - An
J.IJIMPENG.2013.12.006.
update. Appl Mech Rev 1996;49(10):S112—-S117. http://dx.doi.org/10.1115/1.
[53] Biggs JM. Introduction to structural dynamics. New York, NY, USA: McGraw-Hill
3101962.
Book Company; 1964, p. 355.
[24] Jones N. Recent studies on the dynamic plastic behavior of structures. Appl Mech
[54] Nurick GN, Martin JB. Deformation of thin plates subjected impulsive loading–A
Rev 1989;42(4):95–115.
review PART II: Experimental studies. Int J Impact Eng 1989;8(2):171–86.
[25] Jones N. Recent progress in the dynamic plastic behavior of structures. 1978, p.
[55] Rigby SE, Tyas A, Curry RJ, Langdon GS. Experimental measurement of specific
44.
impulse distribution and transient deformation of plates subjected to near-
[26] Lomazzi L, Giglio M, Manes A. Analytical and empirical methods for the
field explosive blasts. Exp Mech 2019;59(2):163–78. http://dx.doi.org/10.1007/
characterisation of the permanent transverse displacement of quadrangular metal
s11340-018-00438-3.
plates subjected to blast load: Comparison of existing methods and development
of a novel methodological approach. Int J Impact Eng 2021;154(March):103890. [56] Pannell JJ, Panoutsos G, Cooke SB, Pope DJ, Rigby SE. Predicting specific
http://dx.doi.org/10.1016/j.ijimpeng.2021.103890. impulse distributions for spherical explosives in the extreme near-field using
[27] Hopkins HG, Prager W. On the dynamics of plastic circular plates. Z Angew a Gaussian function. Int J Prot Struct 2021;1–23. http://dx.doi.org/10.1177/
Math Phys Zamp 1954;5(4):317–30. http://dx.doi.org/10.1007/BF01587827. 2041419621993492.
[28] Wang AJ, Hopkins HG. On the plastic deformation of built-in circular plates [57] Nurick GN, Shave GC. The deformation and tearing of thin square plates
under impulsive load. J Mech Phys Solids 1954;3(1):22–37. http://dx.doi.org/ subjected to impulsive loads - An experimental study. Int J Impact Eng
10.1016/0022-5096(54)90036-8. 1996;18(1):99–116. http://dx.doi.org/10.1016/0734-743X(95)00018-2.
[29] Wang AJ. The permanent deflection of a plastic plate under blast loading. J Appl [58] Teeling-Smith RG, Nurick GN. The deformation and tearing of thin circular
Mech 1955;22(3):375–6. http://dx.doi.org/10.1115/1.4011092. plates subjected to impulsive loads. Int J Impact Eng 1991;11(1):77–91. http:
[30] Micallef K, Fallah AS, Pope DJ, Louca LA. The dynamic performance of simply- //dx.doi.org/10.1016/0734-743X(91)90032-B.
supported rigid-plastic circular steel plates subjected to localised blast loading. [59] Clarke SD, Fay SD, Warren JA, Tyas A, Rigby SE, Elgy I. A large scale
Int J Mech Sci 2012;65(1):177–91. http://dx.doi.org/10.1016/j.ijmecsci.2012.10. experimental approach to the measurement of spatially and temporally localised
001. loading from the detonation of shallow-buried explosives. Meas Sci Technol
[31] Micallef K, Fallah AS, Pope DJ, Louca LA. Dynamic performance of simply 2015;26(1):1–10. http://dx.doi.org/10.1088/0957-0233/26/1/015001.
supported rigid plastic circular thick steel plates subjected to localized blast [60] McDonald B, Bornstein H, Langdon GS, Curry R, Daliri A, Orifici AC. Experimen-
loading. J Eng Mech 2014;140(1):159–71. http://dx.doi.org/10.1061/(asce)em. tal response of high strength steels to localised blast loading. Int J Impact Eng
1943-7889.0000645. 2018;115(January):106–19. http://dx.doi.org/10.1016/j.ijimpeng.2018.01.012.
[32] Cox AD, Morland LW. Dynamic plastic deformations of simply-supported square [61] Mehreganian N, Louca LA, Langdon GS, Curry RJ, Abdul-Karim N. The response
plates. J Mech Phys Solids 1959;7(4):229–41. http://dx.doi.org/10.1016/0022- of mild steel and armour steel plates to localised air-blast loading-comparison of
5096(59)90022-5. numerical modelling techniques. Int J Impact Eng 2018;115(May 2017):81–93.
[33] Mehreganian N, Fallah AS, Louca LA. Plastic dynamic response of simply http://dx.doi.org/10.1016/j.ijimpeng.2018.01.010.
supported thick square plates subject to localised blast loading. Int J Impact [62] Aune V, Fagerholt E, Hauge KO, Langseth M, Bø rvik T. Experimental study
Eng 2019;126:85–100. http://dx.doi.org/10.1016/j.ijimpeng.2018.12.010. on the response of thin aluminium and steel plates subjected to airblast loading.
[34] Parkes EW. The permanent deformation of a cantilever struck transversely at its Int J Impact Eng 2016;90:106–21. http://dx.doi.org/10.1016/j.ijimpeng.2015.11.
tip. Proc R Soc Lond A 1955;228(1175):462–76. 017.
[35] Parkes EW. The permanent deformation of an encastre beam struck transversely [63] Aune V, Valsamos G, Casadei F, Larcher M, Langseth M, Børvik T. Numerical
at any point in its span. Proc Inst Civ Eng 1958;10(3):277–304. study on the structural response of blast-loaded thin aluminium and steel plates.
[36] Wang XD, Yu TX. Parkes revisited: Effect of elastic deformation at the root of a Int J Impact Eng 2017;99:131–44. http://dx.doi.org/10.1016/j.ijimpeng.2016.08.
cantilever beam. Int J Impact Eng 1991;11(2):197–209. 010.
[37] Symonds PS, Fleming WT. Parkes revisited: On rigid-plastic and elastic-plastic [64] Elveli BS, Iddberg MB, Bø rvik T, Aune V. On the strength–ductility trade-
dynamic structural analysis. Int J Impact Eng 1984;2(1):1–36. http://dx.doi.org/ off in thin blast-loaded steel plates with and without initial defects — An
10.1016/0734-743X(84)90013-7. experimental study. Thin-Walled Struct 2022;171(November 2021):108787. http:
[38] Hopkins HG. On the plastic theory of plates. Proc R Soc A 1957;241(1225):153– //dx.doi.org/10.1016/j.tws.2021.108787.
79, URL http://www.jstor.org/stable/2414705. [65] Spranghers K, Vasilakos I, Lecompte D, Sol H, Vantomme J. Full-field deforma-
[39] Martin JB. Impulsive loading theorems for rigid-plastic continua. J Eng Mech tion measurements of aluminum plates under free air blast loading. Exp Mech
Div 1964;90(5):27–42. 2012;52(9):1371–84. http://dx.doi.org/10.1007/s11340-012-9593-5.

18
S.A.E. Alotaibi et al. International Journal of Impact Engineering 178 (2023) 104624

[66] Spranghers K, Vasilakos I, Lecompte D, Sol H, Vantomme J. Numerical simulation [72] Nicholson J. Bessel zero solver - file exchange - MATLAB central. 2022,
and experimental validation of the dynamic response of aluminum plates under URL https://uk.mathworks.com/matlabcentral/fileexchange/48403-bessel-zero-
free air explosions. Int J Impact Eng 2013;54:83–95. http://dx.doi.org/10.1016/ solver.
j.ijimpeng.2012.10.014. [73] Livermore Software Technology (LST). Keyword user’s manual volume ii: mate-
[67] Rigby S, Fuller B, Tyas A. Validation of near-field blast loading in LS-DYNA, 1, rial models R12, Vol. II. r12 ed.. Livermore, CA: Livermore Software Technology
(August) (2018) 1–8. (LST), An ANSYS Company; 2020, p. 3099.
[68] Granum H, Aune V, Bø rvik T, Hopperstad OS. Effect of heat-treatment on [74] Spranghers K, Lecompte D, Sol H, Vantomme J. Material characterization of
the structural response of blast-loaded aluminium plates with pre-cut slits. Int blast loaded plates. Meca.Rma.Ac.Be 2012;(May):9–10, URL http://meca.rma.ac.
J Impact Eng 2019;132(February):103306. http://dx.doi.org/10.1016/j.ijimpeng. be/nctam/ExperimentalTechniques/4{_}Spranghers.pdf.
2019.05.020. [75] Weisstein EW. Parseval’s Theorem - from Wolfram MathWorld, URL https://
[69] Rayleigh J. The theory of sound. The theory of sound, Vol.. 1, Macmillan; 1894. mathworld.wolfram.com/ParsevalsTheorem.html.
[70] Timoshenko S, Young DH. Vibration problems in engineering. Van Nostrand;
1955.
[71] Rao SS. Vibration of membranes. In: Vib. contin. syst.. John Wiley & Sons, Ltd;
2019, p. 427–63. http://dx.doi.org/10.1002/9781119424284.ch13, URL https:
//onlinelibrary.wiley.com/doi/abs/10.1002/9781119424284.ch13, (Chapter 13).

19