Aluminum Structures Exposed To Blast Loading
Aluminum Structures Exposed To Blast Loading
Aluminum Structures Exposed To Blast Loading
loading
TITLE:
BY:
SUMMARY:
Design of blast resistance in structures is an important aspect in modern society. Plated structures are used in a lot of
constructions that can be especially vulnerable to explosions. This includes e.g. protective, offshore or automotive
structures. The main objective in this study was to investigate the response of thin plates made of the aluminum alloy
1050A-H14, as well as to evaluate the available computational methods.
Material tensile tests were performed in order to determine a material model. The material model was created through
inverse modeling by the help of the computer code EUROPLEXUS.
Airblast experiments were performed by using composition C4 in experiments performed at Østøya in Horten. The
experiments were performed in collaboration with the Norwegian Defence Estates Agency (NDEA). Experiments were
performed on stiff calibration plates in order to measure the pressure from the explosions, as well as on aluminum plates in
order to determine the structural response.
Analytical and numerical computations were performed in order to evaluate them compared to the experimental results.
Numerical simulations were performed with both a Lagrangian and an embedded FSI approach.
The analytical approach underestimated the displacements. However, the results indicated that the analytical calculations
could have given better results at lower deformations.
The Lagrangian approach overestimated the displacements. Additional simulations indicate that this could be due to
overestimation of the load by the AIRB approach in EUROPLEXUS.
The embedded FSI approach underestimated the displacements. The simulations performed indicate that the description
of the load is fairly accurate, but that the displacements were underestimated due to the use of elements that were too
large.
MASTEROPPGAVE 2014
TITTEL:
UTFØRT AV:
SAMMENDRAG:
Dimensjonering mot eksplosjonslast er en viktig faktor i det moderne samfunn. Plater blir benyttet i mange konstruksjoner
som kan være spesielt utsatt for eksplosjoner. Dette inkluderer f.eks. beskyttende konstruksjoner, offshore konstruksjoner
eller i bilindustrien. Hovedmålet med dette studiet var å undersøke responsen fra tynne plater, laget av
aluminiumslegeringen 1050A-H14. I tillegg skulle tilgjengelige beregningsmetoder vurderes.
Strekktester ble utført med materialet for å etablere en materialmodell. Denne modellen ble etablert ved hjelp av
inversmodellering i elementprogrammet EUROPLEXUS.
Eksplosjonsforsøkene ble utført ved å detonere C4-ladninger i eksperimenter utført på Østøya i Horten. Forsøkene ble
utført i samarbeid med Forsvarbygg. Forsøkene ble utført på stive kalibrasjonsplater for å måle trykket fra eksplosjonene. I
tillegg ble det utført forsøk på aluminimiumsplater for å undersøke deformasjonen.
Analytiske og numeriske beregninger ble utført for å sammenligne med de eksperimentelle resultatene. Numeriske
simuleringer ble utført gjennom en Lagrange-analyse, og gjennom simuleringer som tar hensyn til interaksjon mellom fluid
og konstruksjon.
Den analytiske tilnærmingen underestimerte deformasjonene. Imidlertid indikerte resultatene at de analytiske metodene
kunne gitt bedre resultater ved mindre deformasjoner.
Analysene som tok hensyn til interakson mellom fluid og konstruksjon, underestimerte forskyvningen. Simuleringene
indikerer at lasten er beskrevet forholdsvis nøyaktig, men at forskyvningene muligens ble underestimert grunnet bruk av for
store elementer.
i
ii
Abstract
Design of blast resistance in structures is an important aspect in modern society.
Plated structures are used in a lot of constructions that can be especially vulnerable
to explosions. This includes e.g. protective, offshore or automotive structures. The
main objective in this study was to investigate the response of thin plates made of
the aluminum alloy 1050A-H14, as well as to evaluate the available computational
methods.
Material tensile tests were performed in order to determine a material model. The
material model was created through inverse modeling by the help of the computer
code EUROPLEXUS.
Airblast experiments were performed by using composition C4 in experiments per-
formed at Østøya in Horten. The experiments were performed in collaboration with
the Norwegian Defence Estates Agency (NDEA). Experiments were performed on
stiff calibration plates in order to measure the pressure from the explosions, as well
as on aluminum plates in order to determine the structural response.
Analytical and numerical computations were performed in order to evaluate them
compared to the experimental results. Numerical simulations were performed with
both a Lagrangian and an embedded FSI approach.
The analytical approach underestimated the displacements. However, the results
indicated that the analytical calculations could have given better results at lower
deformations.
The Lagrangian approach overestimated the displacements. Additional simulations
indicate that this could be due to overestimation of the load by the AIRB approach
in EUROPLEXUS.
The embedded FSI approach underestimated the displacements. The simulations
performed indicate that the description of the load is fairly accurate, but that the
displacements were underestimated due to the use of elements that were too large.
iii
iv
Contents
Acknowledgements i
Abstract iii
Nomenclature ix
1 Introduction 1
2 Theory 3
2.1 State of the Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Blast Load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.1 Groups of Explosions . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.2 Explosives Classification . . . . . . . . . . . . . . . . . . . . . 7
2.2.3 Explosion Processes . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.4 Ideal Blast Wave . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.5 Air Burst . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.6 Internal Explosions . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.7 Blast Load Design . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3 Equation of State (EoS) . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4 Failure Modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.5 Material Characterization . . . . . . . . . . . . . . . . . . . . . . . . 21
2.5.1 Strains and Stresses . . . . . . . . . . . . . . . . . . . . . . . . 21
2.5.2 Necking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.5.3 Fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.6 Constitutive Equations . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.7 Material Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.7.1 Johnson-Cook Material Model . . . . . . . . . . . . . . . . . . 27
2.7.2 Outline of the Johnson-Cook Material Model . . . . . . . . . . 30
2.8 Finite Element Formulations . . . . . . . . . . . . . . . . . . . . . . . 33
2.8.1 Lagrangian Formulation . . . . . . . . . . . . . . . . . . . . . 37
2.8.2 Eulerian Formulation . . . . . . . . . . . . . . . . . . . . . . . 40
2.8.3 Arbitrary Lagrangian-Eulerian (ALE) Formulation . . . . . . 42
2.9 Discrete Particle Method . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.9.1 Theory behind DPM . . . . . . . . . . . . . . . . . . . . . . . 45
v
CONTENTS
3 Preliminary Studies 65
3.1 Analytical Consideration . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.2 Mesh-Sensitivity Study . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.2.1 Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.3 Evaluation of Stand-off Distances . . . . . . . . . . . . . . . . . . . . 73
3.3.1 Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4 Experimental Work 77
4.1 Aluminum alloy EN AW-1050A-H14 . . . . . . . . . . . . . . . . . . 77
4.2 Material Tensile Tests . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.2.1 Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.3 Airblast Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3.1 Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3.2 Results - Calibration Tests . . . . . . . . . . . . . . . . . . . . 92
4.3.3 Results - Component Tests . . . . . . . . . . . . . . . . . . . . 102
vi
CONTENTS
7 Discussion 163
10 Bibliography 173
A Calibration Tests A1
A.1 Comparison of pressure in transducers . . . . . . . . . . . . . . . . . A1
A.2 Pressure curves used for calculating positive impulse . . . . . . . . . . A22
A.3 Results from calibration tests . . . . . . . . . . . . . . . . . . . . . . A43
B Component Tests B1
B.1 Displacement Measurements . . . . . . . . . . . . . . . . . . . . . . . B1
B.2 Deformed plates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B2
B.3 Comparison of pressure in transducers . . . . . . . . . . . . . . . . . B4
B.4 Deformation profiles at different timesteps . . . . . . . . . . . . . . . B9
C MATLAB-files C1
C.1 Script used to calculate results from calibration test R11 . . . . . . . C1
C.2 Script used to calculate results from component test A11 . . . . . . . C6
C.3 Script used to generate 3D-plot from component tests . . . . . . . . . C11
D Cast3M-files D1
D.1 File used to convert mesh from SALOME . . . . . . . . . . . . . . . . D1
D.2 File used to generate mesh for FSI-model . . . . . . . . . . . . . . . . D2
E EUROPLEXUS-files E1
E.1 Material tensile test . . . . . . . . . . . . . . . . . . . . . . . . . . . . E1
E.2 Explosion - Lagrangian approach . . . . . . . . . . . . . . . . . . . . E3
E.3 Explosion - Embedded FSI approach . . . . . . . . . . . . . . . . . . E4
E.4 Explosion - ALE approach . . . . . . . . . . . . . . . . . . . . . . . . E6
vii
CONTENTS
viii
Nomenclature
Acronyms
ALE Arbitrary Lagrangian-Eulerian
CCFV Cell-Centered Finite Volumes
CPU Central Processing Unit
DIC Digital Image Correlation
EoS Equation of State
FE Finite Element
FEA Finite Element Analysis
FSI Fluid-structure interactions
FV Finite Volume
HE High Explosives
JWL Jones-Wilkins-Lee
LE Low Explosives
MDOF Multi-Degree-of-Freedom
MPP Massively Parallel Processing
NCFV Node-Centered Finite Volumes
PDEs Partial Differential Equations
SDOF Single-Degree-of-Freedom
SIMLab Structural Impact Laboratory
SPH Smoothed Particle Hydrodynamics
Greek Letters
α Johnson’s dimensionless damage number
α The linear thermal expansion coefficient
ix
NOMENCLATURE
x
NOMENCLATURE
xi
NOMENCLATURE
xii
NOMENCLATURE
xiii
NOMENCLATURE
xiv
NOMENCLATURE
xv
NOMENCLATURE
dσ Stress increment
dε Strain increment
dσi Stress increment in the ith direction
dτij Shear stress increment in the i-j plane
xvi
Chapter 1
Introduction
1
CHAPTER 1. INTRODUCTION
2
Chapter 2
Theory
Preliminary discussion
What are the objectives of today’s research within structures subjected to blast
loading? What are the researchers looking for? And what type of approaches do
they use? The overall objective in today’s research within structures subjected to
blast loading, is to obtain realistic descriptions of material responses and of the
blast loading itself. Since full-scale testing is both expensive and hard to validate,
other alternatives are necessary. This is where the Finite Element Method (FEM)
plays an important role. By comparing experimental results with corresponding
results obtained by use of finite element codes, researchers are attempting to obtain
reliable solution methods, within reasonable costs. These methods can then be used
to predict the complex structural behaviour due to blast loading. However, full-
scale experiments will never be redundant, as this is the only way of obtaining fully
realistic results.
Several approaches exist within the field of numerical simulations of blast load
problems, both simple and quite complex ones. Recent advances in commercially
3
CHAPTER 2. THEORY
Previous observations
Menkes and Opat [2] studied the response pattern of fully clamped aluminum beams,
subjected to impulsive loading, and divided it into the following three cases:
• Mode I: large inelastic deformation of the entire beam
• Mode II: large deformations and tensile-tearing over the support
• Mode III: transverse shear failure at the support
Subsequently, these modes of failure were identified for other type of structures.
Teeling-Smith and Nurick [3] identified these modes for circular plates, while Nurick
and Shave [4] observed this for square plates. Nurick and Martin [5, 6] have reported
the associated analysis. This way of categorizing the deformation pattern are com-
monly referred to in the literature, and is therefore considered as fundamental in
this context.
Some important results obtained in the recent years will now be highlighted. Several
papers can be looked into for a thorough review of results and observations obtained
in earlier years. Nurick and Martin [5] presented an overview of theoretical and
experimental results up to 1989, primarily for uniformly loaded plates. In the paper
by Nurick and Shave [4], a study of the deformation and tearing of thin square plates
subjected to impulsive loading is presented. They have also given a review of earlier
work within this subject. Nurick et al. [7] highlighted the significant effects of the
4
CHAPTER 2. THEORY
boundary conditions for the purpose of predicting tearing, and similar discussion
were taken by Wierzbicki and Nurick [8].
Balden and Nurick [9] performed numerical simulations of uniformly loaded circular
plates. In fact they described the numerical results from two experimental studies,
one of them published by Teeling-Smith and Nurick [3], in which the deformation
and post-failure response of a plate subjected to uniform blast loading were investi-
gated. Among others, they emphasized the close correlation to what was observed
by Teeling-Smith and Nurick, for both mode I and mode II responses.
Jacob et al. [10] studied fully clamped, circular steel plates subjected to blast
loading. They investigated the effect of stand-off distance and charge mass on the
response of plates with radius of 53 mm. The blast loads were travelling through
tubular structures, giving a focused or localized loading, at least for the nearest
stand-off distances. For stand-off distances greater than 100 mm, a uniformly dis-
tributed loading over the structure was observed. The response of the plates varied
from large inelastic deformation to complete tearing at the boundary.
Neuberger et al. [11] focused on the difference between transient dynamic and
residual deflections on blast loaded, fully-clamped, circular steel plates. In this
paper they introduced the term springback, which represents the difference between
the peak transient deflection and the residual plastic deflection. In addition to
experimental results, a numerical investigation of the differences is reported. In
addition they have presented a quantitative relation between the magnitude of the
springback effect and the stand-off distance.
Børvik et al. [12] investigated edge-clamped steel plates subjected to sand-buried
blast loading. In a fully coupled approach, they used a discrete particle method to
determine the load due to the high explosive detonation products, the air shock and
the sand. The method was validated against experimental tests on square, edge-
clamped steel plates, where spherical 150 g C4 charges were detonated at various
stand-off distances. They observed good quantitative agreement between experimen-
tal data and the numerical simulations. Further, they concluded that the discrete
particle approach is able to describe the physical interactions in the problem.
Sprangers et al. [13] focused on numerical simulations of blast loaded, thin, square
aluminum plates. Numerical results were validated against small-scale blast load
experiments. To capture the dynamic response of the plate, 3D high-speed Digital
Image Correlation (DIC) was used. In order to study how different parameters af-
fect the accuracy of the solution, different numerical configurations were employed.
This included variation in element type and element size, among others. Due to
sufficiently small shear strains, employment of shell elements were found to give
acceptable results when it comes to modeling the deformation. In addition they
concluded that a reduction in computational time, without significant loss of ac-
curacy, can be observed when using an explicit integration scheme instead of an
implicit procedure.
5
CHAPTER 2. THEORY
Alia and Souli [14] described an air blast simulation by use of Eulerian multi-material
formulation. Validation of their numerical approach was done by comparing nu-
merical results with corresponding results obtained from experiments. In addition,
simulated pressure-time histories and pressure impulses were in accordance with ex-
perimental results. According to Alia and Souli [14], the physical phenomena in a
blast load problem can very well be described using a method based on Eulerian
multi-material formulations.
Olovsson et al. [15] compared the coupled Eulerian-Lagrangian formulation and
the corpuscular approach. They showed that the corpuscular method represents a
useful way of simulating close-range blast effect on structures. Numerical robustness,
fast and easy to use, were some of the characteristics given by Olovsson et al.
[15]. A good agreement was further observed between results obtained from the
corpuscular approach, corresponding ALE simulations, and available experimental
data. However, they emphasized that several aspects regarding the implementation
in the finite element code LS-DYNA should be improved, in order to make the
method even more applicable.
Counter-intuitive behaviour of plates have been reported by many authors, e.g.
Flores-Johnson and Li [16], Galiev [17] and Li et al. [18]. Flores-Johnson and Li [16]
investigated numerically the counter-intuitive response of square plates subjected to
blast loading. This effect can in some cases be seen for elastic-plastic materials
subjected to blast loading, and for plate materials this means that the deformation
end up in the opposite direction of the loading. Flores-Johnson and Li [16] defined
a counter-intuitive region for both simply and fully supported plates, in terms of
non-dimensional numbers. It was further observed that a reduction in boundary
constraints lowered the chance of the rebounding instability to occur. Asymmetrical
response of the square plates in the counter-intuitive regions was also reported.
Nowadays, a drastic increase in use of numerical simulation methods can be seen.
However, analytical approaches are still employed, and serve as a useful supple-
ment to numerical methods. Nurick and Martin [5, 6] and Jacob et al. [19] focused
on predictive, empirical methods, based on experimental observations for different
materials and experimental setups. Based on results from previous performed ex-
periments, Nurick and Martin [6] pursued a relation for prediction of mid-point
deflection for impulsively loaded plates. By putting both geometrical and material
properties, in addition to loading conditions, into a dimensionless number, an ex-
pression for the deflection-thickness ratio was obtained. Relations for both circular
and quadrangular plates were established.
6
CHAPTER 2. THEORY
Explosive materials, or explosives, are reactive materials that contain a great amount
of potential energy. Rapid release of this energy is known as an explosion, which is
characterized by production of light, heat, sound and pressure [20]. This definition
is somewhat incomplete, and more thoroughly definitions exist. However, due to the
fact that it is easily understandable and presents the basic concepts quite intuitively,
this definition was chosen.
Explosions can be caused by a countless number of different activities, but generally
they can be divided into three groups [21]:
• Natural explosions
• Intentional explosions
• Accidental explosions
Volcano eruption and lightning are examples of natural explosions. Intentional ex-
plosions are in some way arranged and to some extent controllable. Nuclear weapon
explosions and firing of guns are examples of such type of explosions. Accidental
explosions are characterized by lack of control regarding the release of the energy.
Combustion explosions in enclosures, for instance in gas containers, are most often
unintended and are therefore regarded as an accidental explosion.
Explosives can either detonate or deflagrate. Explosives that detonate are called
High Explosives (HE), while explosives that deflagrate are referred to as Low Ex-
plosives (LE). HE are characterized by supersonic combustion propagating through
high intensity shock waves, while subsonic combustion through heat transfer is char-
acteristic for LE [20]. LE are included here only for completeness, and will not be
discussed extensively in this text.
Comparison of different types of HE can be done in multiple ways. One way, which
is widely accepted as a basis for comparison, is to convert the mass of the relevant
explosive into TNT equivalent mass. This can be achieved by multiplying the mass
of the explosive by a conversion factor, based on its specific energy relative to that
of TNT [22]. Thus, a relation between a given explosive’s demolition power to that
of TNT can be established. In Table 2.1, the relative effectiveness factor for some
types of HE are presented [23].
7
CHAPTER 2. THEORY
The values for the relative effectiveness factors in Table 2.1 are not constant. Other
values may very well be used, but Table 2.1 gives at least a roughly indication of
expected values for different HE. In the case of C4 for instance, use of a relative
effectiveness factor ranging from 1.19 to 1.37 have been reported [13]. The variation
is related to the scope of the study, i.e. what parameter is to be measured.
Different methods can be used to initiate explosives, e.g. percussion or use of elec-
trical energy. The choice of method will depend on both the relevant explosive and
area of application. However, the amount of energy needed to initiate the explo-
sives is called the activation energy, EA . A diagrammatic representation is given in
Figure 2.1. The explosive material contains a given amount of energy prior to the
ignition. After the reaction process starts, the reaction energy decreases to a level
beneath the initial level. For this schematic description of an exothermic reaction,
Q is representing the release of heat energy [22].
8
CHAPTER 2. THEORY
Reaction
energy
EA
Start
Finish
Extent of reaction
In Chapter 2.2.3, a blast wave was referred to as layers of compressed air. More
exact, a blast wave includes both sonic compression waves, shock waves and rar-
efaction waves. The characteristics of a blast wave can vary, depending on how and
9
CHAPTER 2. THEORY
when the energy is released, as well as the distance from the point of detonation
[24].
Before going into details regarding the pressure-time history of an ideal blast wave,
an introduction of different measures of pressure will be given.
For ideal blast waves, it is assumed that propagation occur freely through the air.
However, modification of these waves will occur if interactions with solid or more
dense objects are encountered [25]. In order to describe the blast wave induced
changes in pressure, different measures must be introduced. Important terms in this
discussion are:
• Ambient (static) pressure, P0
• Dynamic pressure, Pdyn
• Side-on (incident) pressure, Ps
• Reflected pressure, Pr
The static pressure is the pressure level prior to the forming of blast waves, and
often this corresponds to the atmospheric pressure. The dynamic pressure is directly
related to motion of the blast waves, and can be expressed as a function of the flow
density, ρf , and velocity, u [24]:
1
Pdyn = ρf u2 (2.1)
2
The stagnation pressure, Pstag , can be given as [24]:
The stagnation pressure can be defined as the static pressure at a stagnation point
in a fluid flow [26].
Two more important definitions are the side-on pressure, Ps , and the reflected pres-
sure, Pr . The side-on pressure, or incident pressure, is measured perpendicular to
the propagation direction of the blast wave, while the reflected pressure is measured
parallel to propagation of the blast wave. Graphical descriptions of how these pres-
sures are measured are given by Bjerketvedt et al [24], and these are reproduced in
Figure 2.2.
10
CHAPTER 2. THEORY
Pr
Wall Wall
Ps
Figure 2.2: Side-on and reflected pressure [24].
By including terms like ideal or idealized, it is assumed that the explosion occurs in
a still, homogeneous atmosphere. In addition, the source is assumed to be spherical
symmetric, such that the characteristics of the blast are functions only of the distance
from the center of the source, R, and the time, t [25]. The main characteristics of an
ideal blast wave are stated below [27]. For a graphical description, see Figure 2.3.
• The arrival time of the shock wave to the point in consideration is represented
by ta . For the idealized case, a discontinuous rise in the pressure level can be
observed at this point in time.
• The peak overpressure is reached for P (t) = P (ta ) = P0 + Ps+ . The pressure
level starts decaying, and will at time t = ta +t+ intersect the ambient pressure
level, P0 .
• The positive phase have a duration of t+ . The area which the positive phase
represents is denoted the positive (incident overpressure) impulse, Is+ , given
in Equation (2.3).
• The positive phase is followed by a negative phase of duration t− . The ex-
pression for the negative impulse, Is− , is given in Equation (2.4). Maximum
negative pressure is reached for P (t) = P0 − Ps− .
• Eventually, the pressure level returns to the ambient pressure, P0 , at time
t = ta + t+ + t− .
11
CHAPTER 2. THEORY
P0 +P+
s
Positive phase
P (t)
Negative phase
P0
P0 -P−
s
t
+ + −
ta ta +t ta +t +t
Z ta +t+
Is+ = [P (t) − P0 ] dt (2.3)
ta
Z ta +t+ +t−
Is− = [P (t) − P0 ] dt (2.4)
ta +t+
Several descriptions of the pressure-time history of the ideal blast wave exist, both
simple and quite complex formulations. Expressions for this type of formulations
are given by Baker [25], where emphasis has been given on fitting the positive phase.
The simplest one only includes two parameters, where a linear decay of pressure is
assumed. This formulation is given by:
1
P (t) = P0 + Ps+ (1 − ), for 0 < t ≤ t+ (2.5)
t+
12
CHAPTER 2. THEORY
In Equation (2.5), t is the time after arrival of the shock wave. The true value for Ps+
is usually preserved when fitting this expression to data. Also, the positive phase
duration, t+ , is adjusted such that the true positive impulse, I + , is maintained [25].
A more complex formulation, called the modified Friedlander equation, is also men-
tioned by Baker [25]. This equation is defined by:
t −b +t
P (t) = P0 + Ps+ (1 − )e t (2.6)
t+
which again stems from the original Friedlander equation. The only difference is
the additional parameter, b, which describes the decay of the curve. The modi-
fied Friedlander equation is stated to be the best compromise, when both desired
accuracy and applicability of the formulation are taken into account. By allowing
adjustment to conform to the most important blast wave properties, without get-
ting too complex, the modified Friedlander equation has proven to be an adequate
approach [25].
Expressions for wave front velocity, uw , air density behind the wave front, ρs , and
the maximum dynamic pressure, Pdyn,max , are given as [22]:
6Ps + 7P0
s
uw = a0 (2.7)
7P0
6Ps + 7P0
ρs = ρ0 (2.8)
Ps + 7P0
5Ps2
Pdyn,max = (2.9)
2(Ps + 7P0 )
where ρ0 is the density of air at ambient pressure ahead of the blast wave, and a0
is the speed of sound in air at ambient pressure. By assuming that air behaves as a
perfect gas, it can be shown that the following relation for reflected peak pressure,
Pr , may be obtained [22]:
7P0 + 4Ps
Pr = 2Ps (2.10)
7P0 + Ps
The lower limit for the reflected pressure, Pr = 2Ps , represents very weak shocks,
Ps P0 , for which acoustic approximations are valid. For stronger shocks, the value
13
CHAPTER 2. THEORY
of the reflected pressure increases dramatically. The upper limit is Pr = 8Ps . This
can be obtained from the Rankine-Hugoniot relations. Note that these limits are
constructed under the assumption that air behaves as a perfect gas. Under strong
shock conditions however, where high pressures and temperatures are present, this
assumption may underestimate the values of the reflected pressure severely [25].
For explosions above ground, the distribution of blast waves will be affected by the
reflecting effects from the ground. The first wave that hits the ground is called the
incident or initial wave. The reflection of this wave will interact and merge with
the initial blast wave, resulting in a so-called Mach front. This front is a result of
the reinforcement of the incident wave by the ground surface. A simple, graphical
description is given in Figure 2.4.
Initial wave
Reflected wave
Charge
Mach front
Ground
The characteristics of the Mach front are basically the same as for the initial wave.
However, the magnitude of the blast parameters are somewhat larger [28].
Explosions that occur inside a confined space are defined as internal explosions, see
Figure 2.5. Complicated blast wave phenomena can be observed in this type of
blast environment. However, initially the characteristics of an ideal blast wave, as
described in Chapter 2.2.4, can be recognized. That is a sharp spike in the time-
pressure curve, and this phase may be referred to as the shock pressure phase. The
duration of this phase depends on the shock front velocities and the distance from
the charge to the surrounding surfaces. The duration of this phase is in general very
short [29].
14
CHAPTER 2. THEORY
Axis of symmetry
Pressure distribution
Wall
Charge
Incident shock
Reflected shock
Floor
Figure 2.5: Shock reflections from walls during internal detonation [29].
Shortly after the shock pressure phase, the gas pressure phase, characterized by a
complicated blast environment, comes into existence. Here, reflected shock waves
propagate and interact, leading to difficulties with defining the blast environment.
This process goes on for quite a long time, compared to the shock pressure phase.
The gas pressure phase, sometimes referred to as a quasi-static or a pseudo-static
phase, terminate when the pressure reaches the ambient level [29].
15
CHAPTER 2. THEORY
3-340-02, step-by-step analyses and design procedures can be found. This include
information about the blast loading, principles of non-linear dynamic analyses, and
considerations about design of reinforced concrete and steel structures. The TM
5-855-1, on the other hand, contains design and analysis procedures for structures
subjected to blast effects from conventional weapons. Information about blast effects
and blast effects on structures can be found in this technical manual. In addition,
close-form equations for generating pressure-time histories are given. However, no
blast wave-structure phenomena is taken into account, which means that this is a
highly simplified approach. Also, for close-range blast loading, these assumptions
result in poor approximations.
A commonly used basis within blast load design, is the scaled distance, Z. Given
the standoff distance, R, and the TNT equivalent mass, WE , a cube root scaling
method, given in design standards, can be employed to find the scaled blast distance
[28]. The scaled distance is given as:
R
Z= 1 (2.11)
WE3
d
Hexp
WE = WEXP (2.12)
HTd N T
where
d
Hexp heat of detonation of explosive in question
HTd N T heat of detonation of TNT
WEXP weight of the explosive in question
Based on the scaled distance, Z, diagrams given in the manual can be used to esti-
mate different blast wave parameters for three different burst environments, namely
free air burst, air burst and surface burst [32]. For free air burst, the ground does
not disturb the blast wave propagation. For the air burst group, the reflections from
the ground interact with the incident wave, according to Figure 2.4. Surface bursts
are valid for detonations close to or on the ground. In this case, the incident wave
and the reflected wave due to the ground will merge immediately after detonation.
This leads to a hemispherical shaped, reinforced wave propagating away from the
point of detonation.
CONWEP is a computer program based on the empirical equations made by Kingery
and Bulmash [33]. They developed equations for predicting air blast parameters for
both spherical air burst and hemispherical surface bursts. Graphical representations
16
CHAPTER 2. THEORY
17
CHAPTER 2. THEORY
where ω = γ − 1 and γ is the heat capacity ratio. A, B, R1 and R2 are material con-
stants that can be found experimentally. The specific internal energy at atmospheric
pressure is given by eint . The relative density, ρ̄, is defined as:
ρsol V
ρ̄ = = (2.14)
ρ Vsol
where ρsol and Vsol are the density and volume of the explosive material in solid
state, respectively. ρ is the current density and V is the current volume. In Table
2.2, JWL parameters for both TNT and C4 are given. The detonation velocity, vdet ,
which comes from an ignition law related to the release of pressure, is also included
in the table.
Unit TNT C4
A Pa 3.738 · 1011 5.98155 · 1011
B Pa 3.747 · 109 0.13750 · 1011
eint,0 J/kg 3.68 · 106 5.4341 · 106
R1 – 4.15 4.5
R2 – 0.90 1.5
ω – 0.35 0.32
ρsol kg/m3 1630 1601
vdet m/s 6930 8500
18
CHAPTER 2. THEORY
Mode I
Mode II
Mode III
Figure 2.6: Failure modes for beams defined by Menkes and Opat [2].
Nurick and Shave [4] observed several phases in the mode II failure region. According
to their work, failure mode II can be divided into three groups of tensile tearing, as
listed below.
• Mode II*: Only partial tearing
• Mode II a: Plate is totally torn. Increasing mid-point deflection is observed
for increasing impulses.
19
CHAPTER 2. THEORY
20
CHAPTER 2. THEORY
Different approaches may be used for defining strains and stresses, depending on
whether small or large deformations are likely to occur. Small deformations, which
implies infinitesimal strains, are often assumed in engineering applications. In this
type of problems, where the undeformed and deformed configurations of a body is
assumed identical, the engineering stress, s, and engineering strain, e, are defined
by the following equations [37]:
F
s= (2.15)
A0
∆L
e= (2.16)
L0
In Equation (2.15), F represents the axial force, while A0 represents the initial
cross-sectional area. In Equation (2.16), L0 represents the initial length, while ∆L
represents the incremental change in length. Engineering stress and engineering
strain are only applicable for deformations within the elastic domain of the material,
in which a linear relationship between stress and strain can be observed. This linear
relationship is given by Hooke’s law:
s = Ee (2.17)
In Equation (2.17), E represents the elastic modulus, also called the Young’s mod-
ulus. As mentioned, this relation is only valid for small deformations. For large
deformations, the geometrical changes in the material must be taken into account.
As the applied force increases, the cross-section of the specimen will decrease. To
include this effect, the Cauchy stress, or true stress, σ, and the true strain, ε, are
needed. These can be computed by use of the following equations [37]:
σ = s(1 + e) (2.18)
21
CHAPTER 2. THEORY
ε = ln(1 + e) (2.19)
It should be noted that within the elastic region, the true stress and strain are
approximately the same as the engineering stress and strain.
The true strain in Equation (2.19), consists of an elastic and a plastic part. The elas-
tic part can be developed directly from Hooke’s law. Consequently, the expression
for the plastic strain εp , can be given by:
σ
εp = ε − (2.20)
E
2.5.2 Necking
The plastic strains start to develop after reaching the yield stress limit for the
material, which marks the border between the elastic and the plastic domains. From
this point on, permanent deformations can be observed. After yielding, further load
can increase the stress. This is due to an effect called strain hardening. However,
at the ultimate stress, that is at the maximum value for an engineering stress-
strain curve, diffuse necking can be observed. Diffuse necking means that the cross-
sectional area will start to decrease in a localized region of the specimen, instead of
over its entire length [38]. This can be seen as a contraction somewhere along the
length of the specimen . For illustration purposes only, the point of diffuse necking
is indicated on an engineering stress-strain curve in Figure 2.7.
s0
Diffuse necking
Yielding
Fracture
e
Figure 2.7: Engineering stress-strain curve in tension [37].
22
CHAPTER 2. THEORY
Necking is basically an instability that can occur for tensile deformations. This
happens when the cross-sectional area starts to decrease at a higher rate than the
material strain hardens [39]. Diffuse necking occurs when the following criteria is
fulfilled [37]:
dσ
=σ (2.21)
dε
σ vs ε
dσ
dε =σ
dσ
dε vs ε
εnecking ε
Figure 2.8: Graphical illustration of necking criterion.
It is important to identify the point of necking since it marks the end of the uniform
deformations. The relations defined in Chapter 2.5.1 are based on an assumption of
uniform deformation, and diffuse necking can therefore be used to appoint the end
of the useful part of a material tensile test [37].
After further elongation, localized necking will occur. Normally, this takes place in
the region in which diffuse necking initiated.
2.5.3 Fracture
When a material loses its load-carrying capacity, material failure occurs. For ductile
materials, the term failure, is often related to initiation of yielding, and not the actual
fracture, as indicated in Figure 2.7 [38]. For brittle materials, material failure and
fracture define the same phenomenon.
23
CHAPTER 2. THEORY
In Figure 2.9 the necking phenomenon can be seen clearly as a reduction in the
cross-sectional area. Due to localized necking in the region where diffuse necking
initiated, fracture occurs in the same area.
24
CHAPTER 2. THEORY
where superscripts e and p refer to elastic and plastic, respectively. When it comes
to stress increments, only the elastic part is included [43]:
{dσ} = [E]{dε } or {dσ} = [E] {dε} − {dε }
e p
(2.23)
where
{dσ} stress increment vector
[E] elastic material property matrix
As stated by Cook et al. [43], {dσ} contains in general increments of all six stress
components. This means that the stress vector can be written as:
where
dσi stress increment in the ith direction
dτij shear stress increment in the i-j plane
In order to make the description of the elastic-plastic material behaviour possible,
the following must be obtained [43]:
• a yield criterion
• a flow rule
• a hardening rule
25
CHAPTER 2. THEORY
An yield criterion relates the state of stress to the onset of yielding. When an
increment of plastic flow occurs, a flow rule is required to give a relation between
the state of stress, given by {σ}, and the corresponding six increments of plastic
strain, {dεp }. In order to take the material hardening into account, a hardening
rule is needed. A hardening rule gives a description of how the yield criterion is
updated, or modified, by straining beyond the initial yield limit [43].
The yield function, F , can be expressed as [43] :
F = F {σ}, {α}, WP (2.25)
where
{σ} state of stress
{α} vector describing kinematic hardening
WP scalar describing the isotropic hardening
{α} and WP , which take material hardening into account, give a description of how
a "yield surface" in multidimensional stress space is altered due to plastic strain-
ing. The "yield surface" can change both location and size, depending on how the
hardening rule is implemented [43].
Depending on which value the yield function, F , takes on, the following can be
stated [43]:
F < 0 elastic domain
F = 0 yielding
F > 0 physically impossible
It should be noted that the statements given above are only valid for elastic-plastic
materials. In theory of viscoplasticity, the yield function, F , is allowed to take values
greater than zero. Theory of viscoplasticity basically describes the rate-dependent
inelastic behaviour of solids [37]. The Johnson-Cook material model is an example
in which viscoplasic behaviour is allowed.
26
CHAPTER 2. THEORY
Strain Measure
There exist a large number of measures for strain and strain rate in nonlinear con-
tinuum mechanics. Two of the most widely used measures are [42]:
1. The Green strain, E.
2. The rate-of-deformation tensor, D.
A strain measure must meet several requirements, e.g. vanish for any rigid body
motions, especially for rigid body rotations. If not, it will predict nonzero strains
for rigid body motions, which again lead to nonzero stresses. It can be shown that
E and D vanish for rigid body motions [42].
27
CHAPTER 2. THEORY
When elastic strains are small compared to the plastic strains, a hypoelastic-plastic
model is typically used [42]. Therefore, the rate-of-deformation tensor, D, will be
selected as the measure of strain rate.
The rate-of-deformation, D, is defined as the symmetric part of the velocity gradient,
L, which is defined as [42]:
∂v
L= = (∇v)T = (grad v)T (2.26)
∂x
where the symbol ∇ specifies the spatial gradient. Further, the velocity gradient,
L, can be decomposed into a symmetric and a skew-symmetric part:
1 1
L = (L + LT ) + (L − LT ) (2.27)
2 2
The first term on the RHS in Equation (2.27), is defined as the rate-of-deformation,
D, while the second term is referred to as the spin, W . Thus, the following can be
stated:
1
D = (L + LT ) (2.28)
2
1
W = (L − LT ) (2.29)
2
Use of an additive decomposition gives the opportunity of separating instantaneous
deformation (symmetric part) form rotations (skew-symmetric part) [46].
Stress Measure
28
CHAPTER 2. THEORY
the stress can be expressed in any of the aforementioned stress measures [42]. The
deformation gradient, F , is also discussed in Chapter 2.8.
Having selected a strain rate measure, one should be careful in the choice of a stress
measure. An important aspect is that the measure of strain and the measure of stress
should be energetically conjugate, or conjugate in power [42]. It can be shown that
the Cauchy stress, σ, and the rate-of-deformation, D, is energetically conjugate. As
a consequence of this, the Cauchy stress, σ, is taken as the measure of stress.
In nonlinear problems, large rotations may be present. This can complicate the
handling of structural elements and anisotropic materials. To simplify this, a coro-
tational approach may be introduced. In short, a coordinate system is then attached
to the material or the element, which will cause it to rotate with the material. The
rate-of-deformation and the Cauchy stress are then expressed in terms of the corota-
tional coordinate system, and are thus referred to as the corotational Cauchy stress,
σ̂, and the corotational rate-of-deformation, D̂. The corotated formulations can be
expressed as [42]:
σ̂ = RT · σ · R or σ̂ij = Rik
T
σkl Rlj (2.30)
D̂ = RT · D · R or D̂ij = Rik
T
Dkl Rlj (2.31)
Hypoelastic Law
Dσ
σ ∇J = − W · σ − σ · WT (2.32)
Dt
where W is the spin tensor, defined in Equation (2.29). The superscripts ∇ and J
refer to objective stress rate and use of Jaumann rate, respectively.
29
CHAPTER 2. THEORY
Denoting the Jaumann rate of the Cauchy stress as an objective rate, indicates that
it is independent of the frame of reference. This is an important property, since
the mechanical response should not depend on the frame of reference [47]. One
appropriate hypoelastic constitutive equation is given by [42]:
σ ∇J = C σJ : D (2.33)
Dσ
= σ ∇J + W · σ + σ · W T = |C σJ{z: D} + |W · σ +{zσ · W T} (2.34)
Dt
material rotation
Equation (2.34) gives the material response in terms of the Jaumann objective stress
rate. As indicated, the material response consists of two parts: the rate of change
due to material response and the change of stress due to rotation [42]. Thus, a valid
rate-constitutive relation has been established.
D = De + Dt + Dp (2.35)
1 + ν ∇J ν
De = σ − tr σ ∇J I (2.36)
E E
where
E Young’s modulus
ν Poisson’s ratio
σ ∇J the Jaumann rate of the Cauchy stress, defined in Equation (2.32)
I 2nd order unit tensor
When it comes to the thermal part, D t , this is defined as [45]:
D t = αṪ I (2.37)
30
CHAPTER 2. THEORY
where Ṫ is the time derivative of the temperature and α is the linear thermal ex-
pansion coefficient.
The last part of the decomposed rate-of-deformation tensor, namely the plastic part,
D p , is defined by the associated flow rule as [45]:
∂f 3 σ0
D p = ṗ = ṗ (2.38)
∂σ 2 σeq
where
ṗ equivalent plastic strain rate
f yield function
σ0 deviatoric stress tensor, i.e. σ 0 = σq− 13 (σ)I
σeq equivalent von Mises stress, σeq = 32 σ 0 : σ 0
Isotropic hardening and thermal softening can be implemented in the von Mises
yield function in the following manner [45]:
where
σ0 = A[1 − (T ∗ )m ] and R = Bpn [1 − (T ∗ )m ]
The material parameters A and B are governing strength, while n and m are gov-
erning strain hardening and thermal softening. T ∗ = TTm−T 0
−T0
is the homologous
temperature, T0 and Tm are the room temperature and melting temperature of the
material, respectively [45]. Having defined the yield function, a multiplicative con-
stitutive relation for the strain-rate dependent material can be defined by:
0
for f ≤ 0
ṗ = h i (2.40)
ṗ0
exp 1
C
σeq (σ)
σ0 (T )+R(p,T )
−1 for f > 0
where ṗ0 is a reference strain rate and C is a material constant governing strain rate
sensitivity. By solving for σeq (σ) in the viscoplastic domain, i.e. for f > 0, the
following expression can be obtained:
ṗ
σeq (σ) = [A + Bpn ] 1 + C ln [1 − (T ∗ )m ] (2.41)
ṗ0
Equation (2.41) gives the equivalent stress as a function of plastic strain, plastic
strain rate and temperature. These characteristics are typical for a Johnson-Cook
type of constitutive equation.
31
CHAPTER 2. THEORY
In fast transient, dynamic problems, e.g. in blast load problems, adiabatic condi-
tions may be assumed. This leads to the following expression for the temperature
evolution [45]:
χ χ
Ṫ = σ : Dp = σeq ṗ (2.42)
ρCp ρCp
ṗ
σeq (σ) = [A + Bpn ] 1 + C ln (2.43)
ṗ0
32
CHAPTER 2. THEORY
The Finite Element Method (FEM), is a numerical method for solving partial differ-
ential equations or integral equations. Conceptually, the method consists of dividing
a domain into many smaller elements. Based on the fact that the characteristics and
behaviour of individual elements are known, the behaviour of the whole domain, or
structure, may be obtained [49]. The computerized version is called Finite Element
Analysis (FEA), and is employed in a wide variety of application areas. FEA in
engineering problems consist of the following three main steps [50]:
1. Pre-processing or modelling of the structure
2. Analysis
3. Post processing
To obtain high-quality results, a thorough understanding of the method is recom-
mended. Thus, the user carries a lot of responsibility.
Two commonly used element formulations are the Lagrangian and the Eulerian for-
mulation. In this section, both the Lagrangian and the Eulerian formulation will be
introduced and described in some detail. In addition, the combined version of these
two formulation, namely the Arbitrary Lagrangian-Eulerian (ALE) formulation, will
be looked into.
The conservation laws, or the balance laws, make up the basis for a group of funda-
mental equations of continuum mechanics. The aim in FEA, related to continuum
mechanics, is to describe a physical system. Therefore, it is of outmost importance
that these fundamental equations are satisfied. For a thermomechanical system, the
following four conservation laws apply [42]:
• Conservation of mass
• Conservation of linear momentum
• Conservation of energy
• Conservation of angular momentum
In short, the conservation laws state that particular measurable properties of an
isolated physical system, do not change as the system evolves [51].
The conservation laws are often expressed by Partial Differential Equations (PDEs),
which are referred to as the strong form. When applying the conservation laws on to
a domain of the body, an integral relation appears, and this results in the so-called
weak form [42]. Due to the importance of these conservation laws, some of them
will be presented in the upcoming sections.
33
CHAPTER 2. THEORY
Kinematic Description
To be able to provide models for the behaviour of fluids, solids and structures,
kinematic relations have to be established. In other words, this means descriptions of
the movements of a body in space and time, which again forms a basis for calculations
of strains and stresses.
The domain of a body at an initial state, at time t = 0, is denoted by Ω0 . This is
said to be the initial configuration. Furthermore, at time t > 0, the body is said to
occupy the domain Ω, and this is often referred to as the current configuration [42].
A graphical description is given in Figure 2.10.
y, Y
φ(X, t)
u Ω
Γ
Ω0 x
X
Γ0
x, X
34
CHAPTER 2. THEORY
Essential Definitions
x = φ(X, t) (2.44)
where x is the position of the material point, X, at time t. The function φ(X, t)
is often referred to as a mapping function, since it maps the reference configuration
onto the current configuration at time t. Equation (2.44) gives the spatial position
of the body as a function of the time, t, and spatial, or Eulerian coordinates, x [42].
The displacement of a material point is defined as the difference between the current
and the original position. Thus, it can be expressed as [42]:
The velocity of a material point is defined as the rate of change of the belonging
position vector. Thus, it is the time derivative with X held constant. The velocity
can be given as [42]:
∂φ(X, t) ∂u(X, t)
v(X, t) = = ≡ u̇ (2.46)
∂t ∂t
Acceleration is defined as the rate of change of the velocity of a material point, and
is thus defined by:
∂v(X, t) ∂ 2 u(X, t)
a(X, t) = = ≡ v̇ (2.47)
∂t ∂t2
Another important variable in the process of calculating deformations and strains
in nonlinear continuum mechanics, is the deformation gradient. The deformation
gradient is defined by [42]:
∂φ ∂x
F = ≡ ≡ (∇0 φ)T (2.48)
∂X ∂X
The expression for the deformation gradient, F , is on the form of a Jacobian matrix,
which is a matrix of all first order derivatives of a vector valued function. Thus, the
deformation gradient can be seen as the Jacobian matrix of the motion φ(X, t) [42].
An interesting relation related to the deformation gradient comes from its determi-
nant, J=det(F ). By finding the determinant of F , a relation between integrals in
the current and the reference configurations can be established [42]:
35
CHAPTER 2. THEORY
Z Z
f (x, t) dΩ = f (φ(X, t), t)J dΩ0 (2.49)
Ω Ω0
36
CHAPTER 2. THEORY
φ(X, t)
Y Y
X X
Material point
Grid node
37
CHAPTER 2. THEORY
Energy conservation
∂wint (X, t) T
ρ0 ẇint = ρ = Ḟ : P − ∇0 · q̃ + ρ0 s (2.53)
∂t
where
ρ0 , ρ original and current density
w hyperelastic potential on reference configuration
b force per unit mass
J determinant of Jacobian between spatial and material coordinates,
J=det[∂xi /∂Xj ]
F deformation gradient, Fij = ∂xi /∂Xj
S second Piola-Kirchhoff (PK2) stress
P nominal stress (transpose of first Piola-Kirchhoff stress)
q heat flux per unit area
s specific heat source term
The corotational formulation, on the other hand, is regarded as an additional for-
mulation, that can be employed together with other formulations, e.g. a total La-
grangian formulation [52]. This type of formulation is useful in structural elements
such as bars, beams and shells [42].
38
CHAPTER 2. THEORY
39
CHAPTER 2. THEORY
φ(X, t)
Y Y
X X
Material point
Grid node
Dρ Dρ
+ ρdiv(v) = 0 or + ρvi,j = 0 or ρ̇ + ρvi,j = 0 (2.54)
Dt Dt
Linear momentum conservation
40
CHAPTER 2. THEORY
Dv Dvi ∂vi
ρ = ∇ · σ + ρb ≡ div(σ) + ρb or ρ = + ρbi (2.55)
Dt Dt ∂xj
Energy conservation
Dwint
ρ = D : σ − ∇ · q + ρs (2.57)
Dt
where
v velocity field
D the rate-of-deformation or velocity strain, D = sym(∇v)
σ cauchy stress tensor
The difficulties with an Eulerian formulation are related to the tracking of moving
material boundaries and interfaces. This require use of different tracking methods
in order to keep track of the material deformations. In contrast to a Lagrangian
formulation, where the mesh is defined by the material domain, an Eulerian mesh has
to exceed the material boundaries to account for deformed configurations. Otherwise
instabilities may occur, leading to inaccurate results [42].
41
CHAPTER 2. THEORY
In ALE formulations, Eulerian and Lagrangian formulations are combined. The aim
of these hybrid techniques are to capture the advantages of the two aforementioned
formulations, while minimizing the disadvantages [42].
The basic concept of this type of formulation is to optimize the meshing process, and
thus reduce the chance of large distortions of elements, which may give numerical
instabilities. In order to establish an ALE formulation, another reference coordinate
system, known as a referential system, is needed [42]. Also, for this more general
framework, velocities and accelerations for both the material and the mesh have
to be defined. Since the ALE is based on Eulerian and Lagrangian formulations,
several attributes from these two formulations can be recognized, e.g. the motion of
the material [42]. Thus, this is still given by Equation (2.44).
The reference domain in the ALE formulation is called the referential domain, and
is denoted by Ω̂. Having defined the referential domain, the motion of the mesh
can be described independently of the motion of the material. Referential, or ALE
coordinates, are denoted by χ. Now, the motion of the mesh can be expressed by
[42]:
x = φ̂(χ, t) (2.58)
From Equation (2.58), the mapping function, φ̂, simply maps points, χ, from the
ALE domain, Ω̂, to points, x, in the spatial domain, Ω. Thus, Equation (2.58)
clearly plays a crucial role in the ALE finite element formulation [42].
In order to define displacement, velocity and acceleration of the mesh in an ALE
formulation, a relation between the ALE coordinates and the material coordinates
must be established. This relation can be expressed by a composition of functions
[42]:
A graphical description of the relation between different domains and the mapping
expressions are given in Figure 2.13.
42
CHAPTER 2. THEORY
Spatial domain, Ω
x
φ̂(χ,t) = φ ◦ ψ −1 φ(X,t) = φ̂ ◦ ψ
Figure 2.13: Maps between Lagrangian, Eulerian and ALE domains [42].
By combining the motion of the mesh and the ψ map, the material motion can be
expressed by the following composition [42]:
From the foregoing equations, expressions for the mesh displacement, mesh velocity
and the mesh acceleration can be obtained. The mesh displacement, û, is defined
by [42]:
∂ φ̂(χ, t) ∂ φ̂
v̂(χ, t) = ≡ ≡ φ̂,t[χ] (2.62)
∂t ∂t χ
Finally, an expression for the mesh acceleration, â, can be given as:
∂ v̂(χ, t) ∂ 2 û(χ, t)
â(χ, t) = = = û,tt[χ] (2.63)
∂t ∂t2
Unless the ALE mesh is Lagrangian, neither the mesh acceleration, nor the mesh
velocity, have any physical meaning. If the mesh is Lagrangian, mesh acceleration
and mesh velocity correspond to material acceleration and material velocity [42].
In order to obtain the conservation laws for the ALE formulation, expressions for the
convective velocity, c, and the material time derivative should be stated beforehand.
43
CHAPTER 2. THEORY
The convective velocity, c, which is the difference between the material velocity and
the mesh velocity, is defines by [42]:
ci = vi − v̂i (2.64)
∂f xj
= f,t[χ] +f,j wi = f,t[χ] +f,j cj (2.65)
∂t ∂χi
Note that the comma, followed by an index, represents the spatial derivative with
respect to an Eulerian coordinate.
Having defined the convective velocity, c, and the material derivative, expressions
for the conservation laws can be presented.
Conservation of mass
σ = σT (2.68)
Conservation of energy
where
θ temperature
k thermal conductivity of the material
44
CHAPTER 2. THEORY
Kinetic molecular theory is the study of gas molecules and their interaction on a mi-
croscopic level. This study has brought forth a relation on macroscopic level, namely
the ideal gas law. A set of assumptions regarding the behaviour and characteristics
of the gas molecules have been made, and these are listed below [54].
• The average distance between the molecules are large compared to their size.
• There is a thermodynamical equilibrium, i.e. the molecules are in random
motion.
• The molecules obey Newton’s law of motion. Relativistic quantum mechanical
effects are negligible.
• The collisions between molecule-molecule and molecule-structure interactions
are perfectly elastic.
The idea behind the kinetic molecular theory dates back to 1738, when D. Bernoulli
suggested that air pressure against a piston could be represented by discrete molec-
ular collisions. In the light of this, J. C. Maxwell derived an expression for the
molecular velocity distribution at thermal equilibrium. From Maxwell’s statistical
descriptions, quantities, such as mean free path and frequency of collision, can be
derived [54].
45
CHAPTER 2. THEORY
By employing the kinetic molecular theory and the belonging assumptions, an ex-
pression for the pressure can be obtained. To illustrate the concept, a rectangular
box with side lengths Lx , Ly and Lz is introduced. Inside this box, a number of N
particles with mass mi and velocity vi = [vx,i , vy,i , vz,i ] are placed, as can be seen in
Figure 2.14.
vi
mi
Lz
Lx
Ly
At thermal equilibrium, the kinetic energy will be evenly distributed between the
different Cartesian directions. Therefore, the pressure in the different directions can
be set equal and expressed as:
2Wk 2wk
Px = Py = P z = P = = (2.70)
3V 3
where
Wk total translational kinetic energy of the particles
V the volume of the box, V = Lx · Ly · Lz
wk the specific translational kinetic energy
From Equation (2.70), it can be seen that at thermodynamic equilibrium, a direct
relation between the pressure and the specific translational kinetic energy can be
established [54].
When it comes to implementation of the DPM, e.g. in LS-DYNA, some assumptions
are necessary in order to sustain computational efficiency. Olovsson et al. [15] have
discussed these type of assumptions for modeling in LS-DYNA. First of all, they
46
CHAPTER 2. THEORY
stated that all molecules for a detonation product cannot possibly be modeled.
Therefore, each particle will represent a large amount of molecules. Particles are
also assumed to be rigid with a spherical shape, since this will speed up the contact
treatment. For a more thoroughly discussion of the DPM, literature by Olovsson
[54] and Olovsson et al. [15] should be revised.
47
CHAPTER 2. THEORY
where
ξ damping ratio in mode j
ωi eigenfrequency of mode j
In order to evaluate the critical time step for a given finite element mesh, the highest
frequency, ωmax , for the unsupported elements should be revealed. In can be shown
that for an undamped material, the critical time step can be written as:
2 L
∆tcr ≤ = (2.72)
ωmax cd
48
CHAPTER 2. THEORY
where
ωmax maximum eigenfrequency of the FE mesh
L element length
cd dilatational wave speed
Equation (2.72) is often referred to as the CFL condition, after Courant, Friedrichs
and Lewy [43]. The equation gives a limit for which information will not propagate
more than the distance between adjacent nodes during a single time step [56]. If the
time step is within this limit, the solution will be bounded.
When it comes to solving of contact problems, like buckling and material failure,
explicit methods are also the most suitable choice [57]. In sum, this makes explicit
methods preferable for blast and explosion problems.
where
{Rine (t)} the inertia force vector, which may be written as [M ]{D̈(t)}
{Rdmp (t)} the damping force vector, which may be written as [C]{Ḋ(t)}
{Rint (t)} the internal force vector
{Rext (t)} the external force vector
{Ḋ(t)} and {D̈(t)} represent the nodal point velocities and the nodal point accel-
erations, respectively.
Approximations of velocity, u̇n (t), and acceleration, ün (t), can be obtained by per-
forming a Taylor series expansion of the displacements un+1 and un−1 , about time
tn . For simplicity, the equations are given for a Single-Degree-of-Freedom (SDOF)
system. This enable use of a graphical description, as can be seen in Figure 2.15.
This will in turn lead to the following expressions for velocity and acceleration [56]:
un+1 − un−1
u̇n = (2.74)
2∆t
49
CHAPTER 2. THEORY
Displacement u(t)
Approximate
Exact
un+1
un−1 un
tn−1 tn tn+1
Time t
∆t ∆t
The above approximations to velocity, u̇n , and acceleration, ün , are referred to as
the conventional central difference equations. Note that these expressions are valid
for SDOF systems. However, expressions for MDOF systems may be obtained in
a similar manner. By substituting these approximations into Equation (2.73), the
following expression may be obtained for a nonlinear MDOF system [56]:
where
1 1
[K eff ] = [M ] + [C] (2.77)
∆t 2 2∆t
2 1 1
{R }n = {R
eff ext
}n −{R }+ 2 [M ]{D}n − 2 [M ]−
int
[C]{D}n−1 (2.78)
∆t ∆t 2∆t
Unless the mass matrix, [M ], and the damping matrix, [C], are diagonal in Equation
(2.76), the effective stiffness matrix, [K eff ], has to be established and factorized in
order to obtain displacements, {D}n+1 . A computational efficient way of finding
displacements, {D}n+1 , arises if only lumped mass matrices are employed. In this
case, damping is assumed either set to zero, or to be given as mass-proportional
damping, [C] = α[M ]. However, the dynamic response will contain high-frequency
50
CHAPTER 2. THEORY
1 1
u̇n− 1 = (un − un−1 ) and u̇n+ 1 = (un+1 − un ) (2.79)
2 ∆t 2 ∆t
This leads to the following expression for the acceleration:
c c
(un−1 − un+1 ) ≈ (un−1 − un ) (2.81)
2∆t ∆t
the preferred form of the central difference method, for a MDOF system, can be
written as [56]:
where
1
{Reff }n = {Rext }n − {Rint }+ [M ](2{D}n + {D}n−1 )
∆t2
1
− [C]({D}n − {D}n−1 ) (2.83)
∆t
The drawback of Equation (2.83), is that only first-order accuracy can be guaranteed,
due to the approximations of the viscous forces. However, for practical structures
that have light damping, and require very small time steps, the two methods will
yield results with insignificant differences when it comes to accuracy [43].
51
CHAPTER 2. THEORY
2.11 EUROPLEXUS
In this section, a brief introduction to the computer code EUROPLEXUS will be
given.
Descriptions below are given according to the EUROPLEXUS user’s manual [46].
EUROPLEXUS is a computer code made for numerical simulations. The code has
been developed by the French Commissariat à l’Energie Atomique (CEA Saclay)
and the Joint Research Centre of the European Commission (JRC Ispra).
EUROPLEXUS is used for analysis of 1D, 2D and 3D domains composed of solids
(continua, shells or beams) and fluids. Fluid-structure interactions (FSI) are also
taken into account, and it is in fact regarded as the state-of-the-art solver with
respect to FSI [58].
EUROPLEXUS fits perfectly well for problems involving rapid dynamic phenomena
(fast transient dynamics) such as explosions, impacts, crashes etc. This is due to the
fact that the code is based on an explicit, central-difference algorithm. In addition,
both geometric and material nonlinearities are fully taken into account.
When it comes to spatial discretization, the Finite Element (FE) and Finite Vol-
ume (FV) are the two most commonly employed formulations. Other discretization
formulations may also be used, such as Smoothed Particle Hydrodynamics (SPH).
Depending on the type of problem to be modeled, the following element formulations
can be used in EUROPLEXUS:
• Lagrangian
• Eulerian
• ALE
Discussions regarding these three element formulations have been given in Chapters
2.8.1, 2.8.2 and 2.8.3. In general, Lagrangian descriptions are used for structural
domains, while Eulerian descriptions are used for fluids.ALE formulations, on the
other hand, are often used for fluids which are interacting with structural domains.
52
CHAPTER 2. THEORY
Structural Domain
where
M lumped mass matrix
a nodal acceleration vector
F ext external forces
F int internal forces
Fext
S1
x
S2
Figure 2.16: Discretization of structural subdomain S1 [59].
To be more precise, the vector containing internal forces, F int , can be given as:
53
CHAPTER 2. THEORY
N els Z
F int = B T σ dV (2.85)
X
i=1 Vn
where
Nels number of elements
Vn element volume
B matrix of shape functions derivatives
σ the Cauchy stress
Equation (2.84) is solved explicitly, since a lumped mass matrix is employed. Time
integration is achieved via the central differences scheme [59].
Fluid Domain
Lagrangian
v
Eulerian w
ALE
54
CHAPTER 2. THEORY
expressed in integral form. For more details, see Chapter 2.8.2. In addition, a
suitable equation of state is needed, in order to model the fluid domain in EURO-
PLEXUS. Often, this implies use of the ideal gas law. Further descriptions of how
these equations are treated in EUROPLEXUS are given by Casadei et al. [59].
When it comes to discretization of the fluid domain in EUROPLEXUS, the following
three techniques may be employed [58]:
• finite elements (FE)
• finite volumes (FV)
• smoothed particle hydrodynamics (SPH)
Regarding modeling of the structural domain, different kinds of finite elements are
generally employed [58].
Figure 2.18 shows three different approaches for discretizing the fluid domain. These
are related only to the finite elements and finite volumes formulations.
FE NCFV CCFV
Figure 2.18: Three different ways of discretizing the fluid domain: finite elements
(FE), Node-Centered Finite Volumes (NCFV) and Cell-Centered Finite Volumes
(CCFV) [59].
If the fluid domain is built up by the use of finite elements, the kinetic variables are
discretized at the nodes, while state variables are discretized at Gauss points. The
internal pressure is an example of a state variable, and is thus collected in Gauss
points. When using a NCFV discretization technique, all variables are discretized
at the nodes, while all variables are discretized at the volume center for CCFV
discretization [58].
55
CHAPTER 2. THEORY
In Table 2.3, two different enforcement strategies are listed, namely the "weak" and
"strong" formulation. The differences are related to how constraints are imposed.
For the weak formulation, the pressure from the fluid is applied as an external force
on the structure. In a strong approach, on the other hand, the constraints are
imposed on the particle velocity and at the mesh velocity in the interface between
the two domains. Constraints on velocity are enforced exactly by use of Lagrange
multipliers. In this way, essential conditions set forth by the user, e.g. restraints
on movement, contact and symmetries, may be coupled automatically with FSI
conditions. Regarding descretizations, finite elements (FE) are traditionally used
in strong FSI algorithms, while finite volumes (FV) are most common in weak
algorithms [58].
Depending on whether an algorithm is able to handle structural failure or not,
it may be classified as either basic or embedded, as indicated in Table 2.3. A basic
formulation is able to handle large deformations, except from structural failure, given
a suitable mesh rezoning technique. A basic formulation can further be divided into
conforming and non-conforming discretizations. If a finite element mesh is classified
as nodally conforming, this simply means that for each node on the Lagrangian
structure, there is a corresponding ALE fluid node, and the other way around. In
the case of a non-conforming finite element mesh, the nodes of the structure and
56
CHAPTER 2. THEORY
the fluid nodes do not need to occupy the same positions in space. However, some
of the nodes will correspond [58].
Embedded FSI algorithms are able to handle structural failure. The concept is to
embed the structure in a fixed Eulerian mesh, and the fluid-structure interactions in
the vicinity of the structure are handled. In this way, the work related to meshing
is reduced compared to a basic algorithm, since it is no longer necessary to make
sure that nodes conform, or that nodes merge. For an embedded algorithm, there
is no possibility of distortion in the fluid mesh, which implies that use of ALE
formulation becomes redundant [58]. Further discussions regarding embedded FSI
algorithms will be taken in Chapter 2.11.4.
In Table 2.4 a descriptive summary of the different FSI algorithms in EUROPLEXUS
is presented.
The embedded approach is, as mentioned, able to handle structural failure. This
is in contrast to algorithms based on an ALE formulation, where the structure
is Lagrangian and the fluid ALE. Approaches based on these formulations, will
experience problems with keeping a suitable form of the mesh if the deformations
are too severe. In particular, this apply if rotations are involved. Also, difficulties
related to merging of the fluid domains after structural failure are regarded as a
drawback [58]. Figure 2.19 shows a discretized structure in a completely independent
Eulerian fluid mesh.
57
CHAPTER 2. THEORY
In order to include FSI in the embedded approach, information about which nodes
that are located close to the structure is needed. The structure’s so-called influence
domain, defines these nodes. Due to deformations of the structure, the influence
domain must be obtained for each time step. Thus, a fast search algorithm is
required. In short, a user defined circle (or sphere in 3D) with a given radius is
defined, and centered at every structural node. By connecting nodes inside these
circles, an area (or a volume in 3D) appears. This will then constitute the influence
domain, consisting of coupled nodes [58].
58
CHAPTER 2. THEORY
p(t)
Below, three different approaches for simulating the explosive charge in EURO-
PLEXUS are given.
• Pressure-time function (AIRB)
• Compressed balloon (GAZP, BUBB)
• Solid TNT (JWLS)
The AIRB approach is a quick and simple way of simulating air blast problems.
Here, a pressure-time function is used for describing the behaviour of the explosive
and the intermediate air. The calculation time is relatively low, since only the struc-
ture is modeled. Consequently, this is an inexpensive and efficient approach. The
pressure-time function is implemented through the modified Friedlander equation,
which has been discussed in Chapter 2.2.4. However, the negative phase cannot be
described by this equation, and is instead represented by a bilinear curve. This is
described in detail by Larcher [27]. The AIRB approach can handle structural fail-
ure and fragmentation [59]. However, some drawbacks are related to this approach.
Reflections, shadowing, channeling effects and other more complex phenomena are
not taken into account. The AIRB approach is therefore limited to simple problems
[59].
Another approach, which is given by the GAZP command in EUROPLEXUS, is
often referred to as the compressed balloon method. The idea behind this method
is to create a "balloon" of compressed gas to represent the explosive. In this way,
structures can be loaded without having to model the explosive [46]. Regarding
the calculation time, this approach is somewhat more expensive than the AIRB
approach.
Initial conditions of the compressed gas are evaluated on the basis of the TNT
equivalent mass, WE , defined in Equation (2.12), and the initial volume of the
compressed gas [59]. In the paper by Larcher et al. [36] the following procedure is
given for obtaining the initial pressure, pbal , for given values of the TNT equivalent
mass, WE , and balloon volume, Vbal .
59
CHAPTER 2. THEORY
ET N T
pT N T = pBrode − p0 = (γ − 1) (2.87)
Vbal
where γ is the heat capacity ratio of the gas in the balloon, p0 is the atmospheric
pressure and
ET N T (γ − 1)
pBrode = + p0 (2.88)
Vbal
3. In order to obtain correct values for impulses, a dimensional scaling factor, αbal ,
is introduced. This factor is defined by
4. Then, an expression for the initial pressure, pbal , in the balloon can be given as
pT N T
pbal = + p0 (2.90)
αbal
5. In order to distribute the scaling factor, αbal , between the internal energy and the
density of the gas inside the balloon, a square root scaling factor, fbal , is introduced
s
pbal
fbal = (2.91)
p0
6. In the end, values for the internal specific energy, eint,bal , and the density of the
gas inside the balloon, ρbal , are obtained in the following manner
where ρ0 is the density of the uncompressed air and eint,0 is the corresponding
internal specific energy.
In the paper by Larcher et al. [36], a more detailed description of the balloon
approach can be found. In addition, discussions regarding the introduction of the
60
CHAPTER 2. THEORY
correction factor, αbal , and the scaling factor, fbal , are given. However, the basic
concept of the balloon approach is that a pressure-time function, resulting from a
compressed balloon, can be employed to represent the corresponding curve for an
air blast wave [36].
EUROPLEXUS also has a command called BUBB. The BUBB command automat-
ically calculates the input values for the balloon method with the help of similar
equation as explained for GAZP. Because of this, only the TNT equivalent mass and
the assigned elements have to be specified.
In the solid TNT approach, the explosive material is described by use of a material
law. This can for instance be the Jones-Wilkins-Lee (JWL) equation of state (EoS)
for explosives, which is described in Chapter 2.3. Thus, an initially solid explosive
material is the starting point in this way of simulating an air blast problem.
In the paper by Larcher and Casadei [36], it is emphasized that elements close
to the explosive should be smaller than the elements which make up the explosive
material. If not, sufficient accuracy may not be obtained, leading to awkward results.
Since small elements increase the computational time, this method is regarded as
expensive with respect to computational cost. In order to reduce the computational
time, so-called spatial partitioning may be employed. This method reduces the
computational time for models with large variation in element sizes [36]. For more
details about this method see the paper by Casadei and Halleux [61].
In Figure 2.21, a graphical comparison between the three approaches is given. Start-
ing from the pressure-time curve in the AIRB approach, an increase in realism,
complexity and computational cost can be observed when going from the AIRB ap-
proach to the GAZP or BUBB approach. The same tendency can be observed when
moving on to the JWLS approach.
p(t)
Compressed
balloon TNT
61
CHAPTER 2. THEORY
Use of DIC renders the possibility of validating different numerical models, which
is of outmost importance in a scientific context. In addition, DIC rarely requires a
comprehensive experimental setup, and this is regarded as one of the main reasons
for the increasing popularity of DIC within mechanical testing. Further, the con-
siderable development of digital cameras in the later years, e.g. within resolution,
recording rate and data storage, have also contributed to the increasing interest for
this post-processing technique.
62
CHAPTER 2. THEORY
Two main approaches exist when it comes to formulation of DIC. One of them is a
subset-based formulation, while the other is based on a finite element (FE) formu-
lation. A graphical description of the FE formulation is given in Figure 2.23. Both
have advantages and disadvantages regarding practical applications in mechanical
testing. For the FE formulation, which is based on conservation of optical flow,
similar formulations in experimental measurements (DIC) and in numerical mod-
eling (FEM) may be employed. This possibility, of combining DIC measurements
with numerical modeling, makes the FE formulation of DIC favourable. However,
regarding the accuracy of measurements, it is not clear from the literature which of
the two approaches that yields the most accurate results.
X u1 u3
1 4
v1 v3
Y
Q4 element
2 u2 3 u4
v2 v4
63
CHAPTER 2. THEORY
64
Chapter 3
Preliminary Studies
65
CHAPTER 3. PRELIMINARY STUDIES
In the paper by Nurick and Martin [6], a method for comparison of results of de-
formed plates with similar geometries, boundary and loading conditions is desired.
All other variables that impact the deformation are therefore put into dimensionless
groups. The basis for their approach is a dimensionless damage number defined by
W. Johnson. Johnson’s damage number, α, is defined as [6]:
ρv 2
α= (3.1)
σd
where ρ is the density of the material, v is the impact velocity and σd is the damage
stress. For metals in impact situations, this number can be used as a guide for
assessing the material behaviour. However, it should be noted that this measure
considers neither structural geometry nor loading conditions. Johnson’s damage
number can be written in terms of the impulse as [6]:
I2 I02
α0 = = (3.2)
A20 t2 ρσd t2 ρσd
Here, I is the total impulse, A0 is the area over which the impulse is imparted, I0
is the impulse per area and t represents the plate thickness. Nurick and Martin
[6] further introduced a geometrical number,β , in the case of quadrangular plates,
defined as:
L
β= (3.3)
B
where L is the plate length and B is the breadth of the plate. Further, the geomet-
rical damage number, ψ, can be written as:
2 1
A0
2
ψ = βα0 (3.4)
A
where A is the area of the plate. In addition, an aspect ratio was introduced, giving
the relation between the plate thickness and the distance from the center of the
plate to the nearest boundary. For quadrangular plates the aspect ratio, λ, is given
as:
B
λ= (3.5)
2t
66
CHAPTER 3. PRELIMINARY STUDIES
By combining the aforementioned equations, Nurick and Martin [6] proposed a mod-
ified damage number, which take both dimensions and loading into account. The
modified damage number, φ, is defined as:
φ = ψλζ (3.6)
For quadrangular plates, the expression for the damage number can be given as:
I
φq = 1 (3.7)
2t2 (BLρσ0 ) 2
δ
= 0.480φq + 0.277 (3.8)
t
where δ is the mid-point deflection. Similar expressions for circular plates are also
given by Nurick and Martin [6]. Jacob et al. [19] performed modifications of Equa-
tion (3.7), in order to include a loading parameter for quadrangular plates subjected
to circular blast loading. Jones [64] presented a dimensionless numbers for both
uniformly loaded circular and quadrangular plates. However, only the expression
obtained by Nurick and Martin [6] be will used for a preliminary evaluation of the
mid-point deflection.
Input parameters
67
CHAPTER 3. PRELIMINARY STUDIES
From Table 3.2, it can be seen that the predicted mid-point deflections increase for
decreasing stand-off distances.
It should be noted that CONWEP is less accurate for close-range blast loading, as
discussed in Chapter 2.2.7. This means that there is a greater source of error as
the stand-off distance decreases. In addition, the total impulse is introduced by
assuming that the peak-impulse value is constant throughout the entire plate. In a
design way of thinking this would have been regarded as a conservative approach,
since the pressure in fact vary considerably over the plate area. Also, the empirical
formulas that Equation (3.8) is based on, do only consider material density and yield
stress. In sum this indicates that the analytical approach suggested by Nurick and
Martin [6], should be considered only as a quick assessment in estimating values for
mid-point deflections. In order to get additional insight in expected values for the
mid-point deflection, some preliminary numerical analysis were performed.
68
CHAPTER 3. PRELIMINARY STUDIES
3.2.1 Procedure
Loading was applied on the modeled plate by use of the AIRB command in EURO-
PLEXUS. This command is described in Chapter 2.11.5. In this approach, pressure
curves are generated on the form of the modified Friedlander equation. These are
generated based on the location of the charge with respect to the structure, as well
as the TNT equivalent mass, and updated at each time step. Since only the plate
is modeled, and the load is applied onto the structure as pressure, no reflections or
additional effects due to the explosion are included [59].
In order to apply the pressure load onto the structure in EUROPLEXUS, two over-
lapping meshes have to be defined. One mesh consist of elements of the type CL3D.
This mesh serves as an absorbing boundary for transferring and distributing the
pressure load generated by the AIRB approach onto the structure. Since this mesh
only functions as an absorbing boundary, no physical properties are assigned.
The plate was modeled by the use of Q4GS elements. A Q4GS element is a 4-node
shell element that can be used in 3D models. The element has 4 integration points
69
CHAPTER 3. PRELIMINARY STUDIES
in the plane and 5 integration points through the thickness [46]. Shell elements can
represent bending properly, and were therefore chosen due to their computational
efficiency compared to solid elements [13].
The quadratic elements in the model were assigned a thickness of 0.8 mm. In Figure
3.1, a mesh with an element size of 10 mm can be seen. The charge was placed in
the middle of the plate with a stand-off distance of 500 mm. The mass of the C4
was set to 30 g, which corresponds to a TNT equivalent mass of 40 g.
Apart from the size of the elements, all parameters were kept constant throughout
the testing. The element sizes that were tested can be seen in Table 3.4. The total
time interval simulated, was set to 5 ms for all analysis. One of the EUROPLEXUS
input-files used in these analyses, can be found in Appendix E.2.
3.2.2 Results
Subsequent to the simulations, the midpoint displacements for all the analyses were
plotted against time. This is shown in Figure 3.2. In Table 3.5, the computational
70
CHAPTER 3. PRELIMINARY STUDIES
time for the different analyses are given. Keep in mind that although all analyses
were performed on the same computer, the computer was also used for other work
simultanously. This means that the available memory can vary between the analyses.
As a result, the listed times are not directly representative for the different element
sizes. However, they can be used for revealing possible tendencies.
40
35
Midpoint-displacement [mm]
30
25
20
15
10 25 mm elements
15 mm elements
5
10 mm elements
5 mm elements
2.5 mm elements
1.5 mm elements
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Time [ms]
From Figure 3.2 it can be seen that the displacement approaches convergence for
the smaller element sizes. The computational time however, increases dramatically.
By comparing the results from the analyses where 2.5 mm and 5 mm meshes were
used, it can be seen that the displacement increases with about 1 mm. At max
displacement, this corresponds to an increase of about 2-3%. Based on these results,
a mesh with an element size of 5 mm is ideal, as the increase in displacement is fairly
low for further reduction in element size. In addition, the computational time is still
relatively low.
71
CHAPTER 3. PRELIMINARY STUDIES
The computational times listed in Table 3.5 are based on a Lagrangian approach,
which is regarded as inexpensive when it comes to computational costs. When
performing e.g. an FSI-analysis, modeling of the air is required. The size of the
elements for the modeled air have to be less than, or equal, to the size of the
elements in the structural domain. As a result of this the computational time will
increase exponentially when reducing the element size. This is because the model
will contain more elements, and since the elements become smaller, the critical time
step will also be reduced. Taking these considerations into account, an element size
of 10 mm might be the most reasonable choice as a basis for comparison between
the different approaches.
72
CHAPTER 3. PRELIMINARY STUDIES
3.3.1 Procedure
When evaluating the stand-off distance, the total time interval simulated was in-
creased to 10 ms. This was done to ensure that the displacements for all stand-off
distances were stabilized when approaching the end of the analyses. Apart from
this, the analyses were set up in the same way as described in Chapter 3.2.
The plate was modeled in accordance with the mesh displayed in Figure 3.1. In
order to get an initial impression of expected results, the possibility of performing
analyses at a lot of different stand-off distances within a short period of time was
regarded as important. In addition, it was expected that future analyses taking FSI
into account, would require the elements of the plate to be 10 mm or larger due to
high computational costs. By selecting 10 mm elements, future analyses would be
more suited for comparison.
The AIRB command in EUROPLEXUS does not work as expected for stand-off
distances smaller than 250 mm. Because of this, only distances larger than 250 mm
were examined. Table 3.6 lists the tested stand-off distances.
3.3.2 Results
The mid-point displacements for all the analyses were plotted against time. This
can be seen in Figure 3.3. To investigate the relation between stand-off distance and
displacement, the maximum positive displacement was identified for each stand-off
distance. Positive direction is defined away from the C4 charge. Figure 3.4 shows
the maximum displacement plotted against the stand-off distance.
73
CHAPTER 3. PRELIMINARY STUDIES
100
250 mm
375 mm
500 mm
80
625 mm
750 mm
875 mm
1000 mm
Displacement [mm]
60
40
20
−20
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
100
90
Maximum displacement [mm]
80
70
60
50
40
30
20
10
200 300 400 500 600 700 800 900 1000
Stand-off distance [mm]
From Figure 3.3 it can be seen that all stand-off distances give significant defor-
mations. As expected, shorter stand-off distances result in larger displacements
than longer stand-off distances. From Figure 3.4 it can be seen that the maximum
displacement increases exponentially when the stand-off distance is reduced.
74
CHAPTER 3. PRELIMINARY STUDIES
Most of the deformation is plastic, as all the analysis stabilizes at a deformed state.
The oscillations around the residual plastic deformation can be attributed to the so-
called springback effect, reported by Neuberger et al. [11]. This effect describes the
difference between the peak transient deflection and the residual plastic deflection.
An interesting observation is that negative deformations can be observed for the
analyses performed for stand-off distances of 875 mm and 1000 mm. This phe-
nomenon may be referred to as counter-intuitive behaviour, and have been reported
by many authors, e.g. by Flores-Johnson and Li [16]. They investigated numerically
the counter-intuitive behaviour of square plates subjected to blast loading.
75
CHAPTER 3. PRELIMINARY STUDIES
76
Chapter 4
Experimental Work
77
CHAPTER 4. EXPERIMENTAL WORK
Si Fe Cu Mn Mg Zn Ti Al
Min 99.50
Max 0.25 0.40 0.05 0.05 0.05 0.07 0.05
Used in experiments 0.03 0.36 0.001 0.002 0.000 0.003 0.010 99.61
The following material properties are assumed known and are therefore not inves-
tigated further. Young’s modulus, E, and the density, ρ, are assumed equal to the
values presented by Aalco Metals Limited [65]. Poison’s ratio, ν, is assumed equal
to the one used by Spranghers et al. [13]. Table 4.2 lists the selected values.
78
CHAPTER 4. EXPERIMENTAL WORK
4.2.1 Procedure
In order to get the best possible results when simulating the airblast, it was decided
to do some material tests to establish a material model. This was done through
material tensile tests. The experiments were performed in a Zwick/Roell Z030 30
kN test machine. A total of 9 specimens were loaded with a deformation rate of
2.1 mm/min. This corresponds to a strain rate of 5 · 10−4 . During loading, the
machine were logging the displacements as well as the applied force. In addition,
a high-speed camera was installed to take pictures during the deformation. These
pictures were taken to be used in a DIC analysis. Figure 4.1a shows a picture of the
test setup. Figure 4.1b shows a close-up view of the specimen installed in the test
machine.
Specimen
Camera
The test specimens were made of the aluminum alloy EN AW-1050A-H14. All
specimens were cut from 0.8 mm thick, rolled aluminum plates. For rolled plates,
the material behaviour can vary depending on the orientation of the material relative
to the direction of the rolling. Because of this, three specimens were made from each
orientation of 0, 45 and 90 degrees, relative to the rolling direction. The nominal
dimensions of the specimens can be seen in Figure 4.2.
79
CHAPTER 4. EXPERIMENTAL WORK
25 25.05 70 65
t=0.8
40
⊘
12.
7
12.5
15
=
R
200
The thickness, as well as the width, of the specimens were measured in three places.
At the midpoint, as well as at each end of the gauge length. The measurements were
taken by the help of a sliding caliper and a micrometer. Average values, as well as
the standard deviations, are listed for all 3 orientations in Table 4.3.
80
CHAPTER 4. EXPERIMENTAL WORK
4.2.2 Results
After testing it was apparent that the necking initiated close to the shoulder of the
specimen. This can be seen in Figure 4.4. This is not ideal as the behavior can be
affected by the change in geometry close to the shoulder.
Figure 4.5 shows force-displacement curves for all specimens. The test specimens
with the same orientation, relative to the rolling direction, are plotted against each
other. This was done to investigate the variation, as well the impact of the physical
extensometer. Note that the displacement used, were logged by the machine. These
values are inaccurate, but as the source of error is the same for all tests, the curves
can give a good basis for comparison.
Force [kN]
Force [kN]
1 1 1
From the curves in Figure 4.5, it can be seen that the variations between the tests
are small until the specimens approach necking. The difference is larger when ap-
proaching failure, but the general behaviour is similar. Specimen 3, 6 and 9 were
equipped with physical extensometers. From their respective curves it is apparent
that the impact of this is minimal. Figure 4.6 shows force-displacement curves for
the three tests that were equipped with extensometers.
81
CHAPTER 4. EXPERIMENTAL WORK
1.5
1
Force [kN]
0.5
0 degrees orientation
45 degrees orientation
90 degrees orientation
0
0 1 2 3 4 5 6 7
Displacement [mm]
From Figure 4.6 it can be seen that the material behaviour is fairly similar for the
three tests in the elastic region. After necking however, the material shows highly
anisotropic behavior.
As the extensometer did not affect the results, extensometer values were used to
calculate the strains. The extensometer malfunctioned during the test of the 45
degrees specimen. As a result, only curves for the specimens at 0 and 90 degrees
orientation will be shown. Figure 4.7 shows the engineering stress plotted against
engineering strain, and Figure 4.8 shows the true stress plotted against true strain.
As stated in Chapter 2.5.1, the formula for true stress is only valid until necking.
Because of this, all values after necking have been removed from the true stress-true
strain curves.
To be able to calculate the plastic strain, a value for Young’s modulus, E, is needed.
From literature, this value has been assumed to be E = 71 GPa [65]. Figure 4.9
shows true stress plotted against the resulting plastic strain.
82
CHAPTER 4. EXPERIMENTAL WORK
140
120
Engineering stress [MPa]
100
80
60
40
20
0 degrees orientation
90 degrees orientation
0
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
Engineering strain
140
120
100
True stress [MPa]
80
60
40
20
0 degrees orientation
90 degrees orientation
0
0 1 2 3 4 5 6
True strain −3
x 10
83
CHAPTER 4. EXPERIMENTAL WORK
140
120
100
True stress [MPa]
80
60
40
20
0 degrees orientation
90 degrees orientation
0
0 0.5 1 1.5 2 2.5 3 3.5 4
Plastic strain −3
x 10
From Figure 4.9, it can be seen that very little remains of the curves when removing
all points after necking. This means that necking occurs at a really low plastic
strain, which indicate that the material has very little strain hardening.
84
CHAPTER 4. EXPERIMENTAL WORK
4.3.1 Procedure
The light source in Figure 4.10 was later moved to the right side of HS CAM 2. This
was done in order to obtain more favourable light conditions when taking the DIC
analysis into consideration.
85
CHAPTER 4. EXPERIMENTAL WORK
HS CAM 1
Component
Clamping frame
HS CAM 2
3300
Stand-off-distance
Figure 4.12 shows the testing frame with an installed aluminum plate in between
the steel mounting frame and the clamping frame.
86
CHAPTER 4. EXPERIMENTAL WORK
Aluminum plate
Clamping frame
400
HUP50x50mm
HUP50x50mm
1000
400 300
1000
(a) Seen from the side where the charge was placed.
H
U
P5
0 x5
0m
m
1000
Clamping frame
HUP50x50mm
1000
(b) Side view of the testing frame.
The explosive charges, which consisted of Composition C4, were placed as indicated
in Figure 4.11. The mass of the charges were approximately 30 g and the charges
had a diameter close to 34.5 mm. The charges were shaped by hand. Because of this,
it is possible that geometry and mass varied between the tests. In the experimental
87
CHAPTER 4. EXPERIMENTAL WORK
setup, the charges were suspended from the upper girder, as indicated in Figure
4.13.
Explosive charge
Figure 4.13: Pictures showing the placement and shape of an explosive charge.
The vertical position of the charges was adjusted, such that the center of mass was
at the same height as the center point of the tested aluminum plates. The horizontal
distance from the plate was adjusted according to the desired stand-off distances.
Prior to the testing of the aluminum plates, experimental validations were carried
out on a steel plate with a thickness of 15 mm. Pressure transducers were installed
in the steel plate, according to Figure 4.14. A breakwire was installed in order
to determine the time of detonation. Three stand-off distances were selected, and
these are listed in Table 4.4. The steel plate was mounted on the opposite side of
the clamping frame. It should be noted that the stand-off distances given in Table
4.4 are measured with respect to the placement of the steel mounting frame. This
means that stand-off distances for the calibration plate are actually 15 mm longer.
For simplicity, the stand-off distances in Table 4.4 and Table 4.6 are referred to
the same reference point. By conducting these calibration tests, estimates of the
pressure distribution over the plates in the actual experiments were obtained.
88
CHAPTER 4. EXPERIMENTAL WORK
400
5
55
4
40
3
40
2
40
.5
400
6 62 1 7
6 2.5 8
62.5 9
10
Figure 4.14: Calibration plate seen from the side where the charge was placed.
Subsequent to the calibration tests, the aluminum plates were tested. In Figure
4.15, the geometry of the tested plates are given. Thicknesses were measured for
each plate prior to the tests by use of a micrometer. Points measured were located
in the mid-span of each side and were taken close to the boundary. This was done
in order to investigate the variation in the plate thickness. The average value and
the standard deviation from the measurements are presented in Table 4.5.
Thickness [mm]
Average Standard deviation
0.790 0.004
Tested stand-off distances are given in Table 4.6. Stand-off distances used for the
calibration plate and the component tests differ, as can be seen by comparing the
contents in Table 4.4 and Table 4.6. The reason for this is discussed in Chapter
4.3.3.
89
CHAPTER 4. EXPERIMENTAL WORK
Pressure transducers of the type Kistler 603B [66], were used for both the cali-
bration tests and the component tests. However, the arrangement of the pressure
transducers differed. Pressure transducers inside the clamping frame were removed
and pressure transducer 1 was moved down to the lower part of the clamping frame,
according to Figure 4.14. Thus, pressure transducers were installed in the mid-span
of each side of the aluminum plates, corresponding to the holes with diameter of 24
mm in Figure 4.15.
400
13
400
24
During the experiments, two high-speed cameras of the type Phantom v1610 [67],
were taking pictures with a frequency of 21.000 fps. Image size was set to 896x800
pixels. The positioning of the cameras can be seen in Figure 4.11. Synchronized
stereo images taken by the two cameras, together with use of DIC, enabled subse-
quent abstraction of 3D displacement fields from the dynamic response of the plates.
An aluminum plate painted in a speckle-pattern is shown in Figure 4.16. This was
done in order to obtain contrasts, which is advantageous for use in a DIC analysis.
90
CHAPTER 4. EXPERIMENTAL WORK
In order to determine the time of arrival for the shock waves, a blast pencil was
installed in between the two cameras, as can be seen in Figure 4.11. This was
used to define the applicable time window for the DIC analyses, as shock waves
potentially could disturb the camera setup and consequently also the DIC analyses.
91
CHAPTER 4. EXPERIMENTAL WORK
As explained in Chapter 4.3.1, the calibration plate had a total of 10 pressure trans-
ducers installed. Each transducer recorded the pressure during the calibration tests.
This enabled a graphical description of the variation in pressure as a function of
time. Before processing data from the pressure curves, the curves had to be cor-
rected according to the time of detonation. At the time of detonation, a breakwire
transmitted a signal that was plotted against time. The pressure curves could then
be corrected by identifying the time interval between the start of the time series and
the time of detonation.
There were problems with some of the tests. In test R11, pressure transducer 5 gave
strange readings. After being exposed to the blast pressure, the pressure failed to
return to initial level. In test R21 the breakwire malfunctioned, and in test R22
the response of the breakwire was recorded after the pressure waves arrived at the
pressure transducers. For test R31, R32 and R33, transducers 6 through 10, as well
as the breakwire, malfunctioned. Because of this, three more test were performed.
These were labeled R31-2, R32-2 and R33-2. In test R31-2 and R32-2, transducer 9
malfunctioned.
All tests with functional breakwire were corrected according to the time of det-
onation. Figure 4.17 shows the breakwire response from test R13. The time of
detonation is indicated by a red line.
4
Response
−1
0 2 4 6 8 10 12 14 16 18 20
Time [ms]
The pressure transducers were located as shown in Figure 4.14. It was expected
that the blast wave would arrive at the middle pressure transducer first, and at
92
CHAPTER 4. EXPERIMENTAL WORK
the others shortly after. The pressure curves were plotted against each other to
compare differences in pressure. To get a good basis for comparisons, transducer
1 through 4 were plotted together in one plot and 1 plus 8 through 10 in another.
This was done because these transducers were aligned in vertical and diagonal direc-
tions, respectively. In addition, transducer 1 and 4 through 7, were plotted together.
Transducers 4 through 7 were ignored in the other figures due to the fact that the
transducers were not located at the same distance as the plate. The time of arrival
would therefore not be comparable to the other transducers. A set of figures is
presented for one characteristic test at each stand-off distance. Figures for the re-
maining tests can be found in Appendix A.1. Figure 4.18 shows plots of transducers
1 through 4 for test R13. Figure 4.19 shows plots of transducers 1 plus 8 through 10
for test R13. Figure 4.20 shows plots of transducers 1 plus 5 through 7 for test R13.
Figure 4.21, Figure 4.22 and Figure 4.23 show corresponding plots for test R23 while
Figure 4.24, Figure 4.25 and Figure 4.26 show corresponding plots for test R33-2.
12000
Pressure transducer 1
Pressure transducer 2
Pressure transducer 3
10000 Pressure transducer 4
8000
Pressure [kPa]
6000
4000
2000
0
0 0.2 0.4 0.6 0.8 1
Time [ms]
93
CHAPTER 4. EXPERIMENTAL WORK
12000
Pressure transducer 1
Pressure transducer 8
Pressure transducer 9
10000 Pressure transducer 10
8000
Pressure [kPa]
6000
4000
2000
0
0 0.2 0.4 0.6 0.8 1
Time [ms]
12000
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
10000 Pressure transducer 7
8000
Pressure [kPa]
6000
4000
2000
0
0 0.2 0.4 0.6 0.8 1
Time [ms]
Figure 4.20: Transducers in the middle and in the clamping frame for test R13.
94
CHAPTER 4. EXPERIMENTAL WORK
5000
Pressure transducer 1
4500 Pressure transducer 2
Pressure transducer 3
4000
Pressure transducer 4
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
5000
Pressure transducer 1
4500 Pressure transducer 8
Pressure transducer 9
4000
Pressure transducer 10
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
95
CHAPTER 4. EXPERIMENTAL WORK
5000
Pressure transducer 1
4500 Pressure transducer 5
Pressure transducer 6
4000
Pressure transducer 7
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
Figure 4.23: Transducers in the middle and in the clamping frame for test R23.
3000
Pressure transducer 1
Pressure transducer 2
2500
Pressure transducer 3
Pressure transducer 4
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
96
CHAPTER 4. EXPERIMENTAL WORK
3000
Pressure transducer 1
Pressure transducer 8
2500
Pressure transducer 9
Pressure transducer 10
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
3000
Pressure transducer 1
Pressure transducer 5
2500
Pressure transducer 6
Pressure transducer 7
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
Figure 4.26: Transducers in the middle and in the clamping frame for test R33-2.
From Figure 4.18 it can be seen that the pressure arrives at the pressure transducers
in the expected order of 1 to 4. Although transducer 2 indicates a slightly lower peak
pressure than transducer 3, the curves still indicate a reduce in peak pressure when
moving away from the center of the plate. This can also be observed from Figure
4.19. When looking at the time of arrival, it can be seen that it is approximately
the same for transducer 1 and 8. Looking at Figure 4.20 it can be seen that the
97
CHAPTER 4. EXPERIMENTAL WORK
pressure arrives at transducer 6 before transducers 5 and 7. This can indicate that
the load was slightly off-center. As transducer 6 is located at the same half of the
plate as transducer 8, this can also explain why the pressure arrives at the same
time in pressure transducers 1 and 8.
From Figure 4.19, it can be seen that the positive impulse from transducer 9 and
10 are large compared to the other transducers. It is possible that this is a result of
pressure accumulation in the corner of the plate. Similar response can be observed
for the other stand-off distances in Figure 4.25 and Figure 4.22.
Figure 4.21 and Figure 4.24 show the same tendency as Figure 4.18. The pressure
decreases and the time of arrival increases, when moving away from the centre of
the plate. Figure 4.22 and Figure 4.25, both show that the pressure in pressure
transducer 1 arrives after some of the others. Looking at Figure 4.23, it is hard to
conclude whether this could be due to off-centering of the load. The pressure in
transducer 7 arrives before the others. If the load were slightly displaced in this
direction, it would explain the big increase in time of arrival for pressure transducer
10. However, it would not explain why the pressure arrives at transducer 9, before
transducer 1. This might indicate either other influencing factors, or the possibility
that some of the pressure curves may be slightly out of sync. Figure 4.26 indicates
the same as Figure 4.20, that the load might be slightly off in the direction of
transducer 6.
The pressure curves from the transducers were used to calculate the time of arrival,
impulse and the duration of the impulse. Figure 4.27, Figure 4.28 and Figure 4.29
show the pressure in pressure transducer 1 for test R13, R23 and R33-2, respec-
tively. The red lines indicate the time of arrival and the end of the positive impulse.
Corresponding plots for all transducers for all the tests are included in Appendix
A.2.
98
CHAPTER 4. EXPERIMENTAL WORK
12000
10000
8000
Pressure [kPa]
6000
4000
2000
5000
4500
4000
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
99
CHAPTER 4. EXPERIMENTAL WORK
3000
2500
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
By inspecting the pressure curves, it can be seen that the pressure does not return to
ambient pressure. This indicates that the zero level has changed. A weakness with
the Kistler 603B pressure transducers, is that they are designed for temperatures up
to 200°C [66]. It is therefore possible that the zero level has drifted due to the high
temperatures created by the explosion. This means the zero level has changed for
all pressure curves. Because of this, it is very hard to identify the negative phase
of the explosion. In Figure 4.29, the negative phase can be seen fairly well after
the second red line. However, the negative phase is disrupted by a reflection wave
at around 1.5 ms. Because the reflection waves give a significant contribution to
the pressure, as well as the fact that the negative phase is hard to locate due to
the change of the zero level, it was decided not to calculate the negative impulse.
When calculating the positive impulse, it was decided to include any positive peaks
that hit the pressure transducers before the positive phase ended. This is especially
apparent in Figure 4.28.
In addition to the positive impulse, Is+ , and the positive phase duration, t+ , the time
of arrival, ta , and the maximum pressure, pmax , were calculated. These parameters
are explained in Chapter 2.2.4. In order to calculate the time of arrival, it was
important that all tests were corrected according to the time of detonation. The
tests that lacked results from the breakwire, were corrected according to the time
of arrival of the tests with functional breakwires. Test R21 and R22 were then
corrected according to time of arrival for transducer 1 from test R23. Test R31,
R32 and R33 were corrected according to time of arrival for transducer 1 from test
R33-2. This adjustment means that the time of arrival for the corrected tests are
based on the assumption that the time of arrival for transducer 1 is the same as
for the test they were corrected against. Table 4.7 lists the calculated parameters
100
CHAPTER 4. EXPERIMENTAL WORK
from pressure transducer 1 and 5 for tests R13, R23 and R33-2. The complete set
of results is listed in Appendix A.3. One of the MATLAB-scripts used to calculate
the results is included in Appendix C.1.
From the values listed in Table 4.7, it can be seen that the differences in maximum
pressure and the time of arrival decrease as the stand-off distance increase. This
can be attributed to the fact that the difference in stand-off distance between the
pressure transducers will decrease as the stand-off distance increase.
101
CHAPTER 4. EXPERIMENTAL WORK
Before performing the tests presented in Table 4.6, it was decided to perform one
test at stand-off distance 250 mm. This was done in order to investigate whether it
would be possible to perform tests at the same stand-off distances as the calibration
tests, rather than the stand-off distances presented in Table 4.6. After testing, it
was apparent that this was impossible, as the plate was completely cut off at the
supports. Figure 4.30 shows the plate after testing.
Figure 4.30: Plate after testing with stand-off distance 250 mm.
From Figure 4.30a, it can be seen that the plate is heavily deformed. Comparing it
to the failure modes presented in Chapter 2.4, it appears that the failure corresponds
to failure mode II. However, looking at Figure 4.30b, it can be seen that the cut is
very clean along the support, apart from in the corners. This may indicate that the
failure is actually due to shear at the supports. This suggests that the failure might
be of mode III, or possibly something in between mode II and III.
In Figure 4.30b, it can be seen that there is very little deformation around the bolts.
Around the pressure transducers however, the deformation is larger. As there are
no bolts at the midpoint of the edges, it is expected that the friction force is lower.
This means that the force required for sliding on the supports is lower than along
the rest of the frame.
It was decided to go with stand-off distance 375 mm, 500 mm and 625 mm for the
experiments. Stand-off distance 250 mm was omitted due to the complete tearing
along the supports. Figure 4.31 shows one plate after testing for each stand-off
distance. Two plates are shown at stand-off distance 625 mm, in order to highlight
the slight differences in the deformed configuration. The remaining plates can be
seen in Appendix B.2.
102
CHAPTER 4. EXPERIMENTAL WORK
From Figure 4.31a, it can be seen that stand-off distance 375 mm results in fracture
along the support close to the corners. As with stand-off distance 250 mm, the
fracture is very clean. A close-up view of a fracture can be seen in Figure 4.32.
103
CHAPTER 4. EXPERIMENTAL WORK
The fracture in Figure 4.32 is from plate A12. The number of, and length of the
fractures, varies from test to test. The length of the fracture shown in Figure 4.32,
is approximately 45 mm, but the length can be up to around 60 mm, which is the
case for one of the fractures in plate A13.
Figure 4.31b shows that no fracture is observed at stand-off distance 500 mm. This
is also the case for the plates at stand-off distance 675 mm.
The most interesting observation come from Figure 4.31c and Figure 4.31d. Here
it can be seen that the final deformation is in the opposite direction of the other
tests. That means the plate is now closer to the load than before the detona-
tion. This is similar to the counter-intuitive behavior that was investigated by e.g.
Flores-Johnson and Li [16]. This counter-intuitive behavior was also observed in the
simulations performed in Chapter 3.3. However, this was only observed at stand-off
distances of 875 mm or more.
While the plates completely returned to the opposite direction for test A32 and A33,
it can be seen from Figure 4.31c, that the midsection of the plate in test A31, still
had positive displacement relative to the rest of the plate.
Similar to what was observed in Figure 4.30b, it can be seen that there is more
deformation around the pressure transducers than along the rest of the support.
Figure 4.33 shows the deformation around pressure transducer 5 for plate A11.
From Figure 4.33, it can be seen that there is a significant deformation along the
edge. The displacement in the middle is approximately 4 mm relative to the initial
configuration. It can be seen from the figure that the displacement decreases as
it gets closer to the bolt holes. After passing the first bolt, there is no visible
deformation further along the edge. Figure 4.34 shows the hole to the right in
Figure 4.33, as well as the hole further along to the right.
104
CHAPTER 4. EXPERIMENTAL WORK
Figure 4.34: Bolt holes between pressure transducer 5 and upper right corner of
plate A11.
From Figure 4.34a, it can be seen that the hole has been slightly deformed in the
direction of the pressure transducer to the left. From Figure 4.34b, it can be seen that
this is not the case when getting close to the corner. Controlling with a sliding caliper
indicated that there were very little difference in the diameter. It was apparent from
the tests at higher stand-off distances, that the deformation around the pressure
transducers were decreasing as the stand-off distance increased. This also meant
that the deformation of the holes decreased, and already at stand-off 500 mm, very
little deformation was noticeable. This indicates that the friction, produced by the
bolts, is enough to prevent sliding of the plate along the support. Because of this,
it is possible that the boundary conditions will function almost as fixed boundaries
when the load is small. The friction force will be the highest right under the bolts,
and then gradually decrease as the distance to the bolts increases. Due to less
friction forces around the transducers, sliding along the supports will happen at
weaker loading.
The displacement of the plates were measured by the help of a sliding caliper, before
the plates were removed from the clamping frame. In addition, the plates were
measured again after being brought back to Gløshaugen, NTNU. This was done
to account for possible restrained, elastic deformations. The measurements were
taken at the points marked in Figure 4.35. The displacements were not measured
in points 2, 4, 7 and 9 at Østøya. These points were included when measuring
at Gløshaugen, in order to give a more complete representation of the deformation.
Table 4.8 lists the measurements from the plates shown in Figure 4.31. The complete
set of measurements is included in Appendix B.1. The displacements are considered
positive when moving away from the detonation.
105
CHAPTER 4. EXPERIMENTAL WORK
300
Clamping frame boundary
50
75
1
3/8 40
2
35
6 7 9 10
35
4
40
5
75
50
50 75 40 35 35 40 75 50
From Table 4.8, it can be seen that the difference in displacements measured at
Østøya and at Gløshaugen, are fairly small. In most points, the difference is lower
than 2 mm, and in some points the displacements are almost identical. The mea-
surements taken at Østøya were taken rather quickly, and it is likely that the mea-
surements were not taken at the exact same points on Østøya and at Gløshaugen.
Because of this, any variation, or lack thereof, may potentially be due to error in
the measurements. Regardless, the results from Table 4.8 still indicate that the
deformation is pretty similar between the two sets of measurements. This indicates
that any remaining elastic deformations were small, or the clamping frame did not
prevent the plates from releasing the elastic deformations. Comparing the displace-
ments at the midpoint, with the results from the preliminary studies, it can be seen
106
CHAPTER 4. EXPERIMENTAL WORK
that the displacements are way lower than the ones resulting from the numerical
simulations. The analytical calculations, on the other hand, resulted in values lower
than the actual displacements, but still closer than what the numerical simulations
suggested.
In order to investigate the behavior of the plates, they were all prepared for DIC
analysis. On most plates, the experiments were performed before the paint were
completely dry. For plate A23 and A32 however, this were not the case. When
the plates reached maximum displacement, the paint at the center of the plates
detached. This prevented the DIC analyses from giving reliable results after maxi-
mum displacement, as the DIC-mesh was heavily destorted around the affected area.
Figure 4.36a shows where the paint disappeared on plate A23. The affected area is
marked in red. A close-up view of the affected area is shown in Figure 4.36b.
Before performing DIC analyses on the plates, a mesh had to be applied. Figure
4.37 shows the mesh used in the analyses. The same mesh was used for all the
plates. eCorr was used to carry out the DIC analyses.
107
CHAPTER 4. EXPERIMENTAL WORK
When performing the analyses, it was apparent that there were some additional,
minor problems. For plate A11, A12 and A13, the fracture made it possible for light
from the explosion to get through. This affected the DIC-mesh close to the fracture,
but didn not have any negative impact on the remaining mesh. When investigating
the pictures used in the DIC analyses, it was apparent that the fracture at stand-off
distance 375 mm appeared at around 17 mm midpoint displacement. This is lower
than the maximum displacement of the plates at 500 mm stand-off. Because the
plates at 500 mm stand-off, did not get any fracture, this further indicates that the
fractures are a result of shear at the supports.
Some of the pictures used for the DIC analysis of plate A33, were disturbed by light
during the experiment. This may have been due to an unexpected light source, or
due to reflection created by deforming the plate. This light source impacted the
mesh around the affected area. Most of the mesh however, was still unaffected.
During the experiments, the pressure were logged by pressure transducers 1, 5, 6
and 7. The pressure were logged by the pressure transducers at the same frequency
as in the calibration tests. In order to know the pressure associated to each picture,
a separate set of measurements were logged at the frequency of the camera. Be-
cause the cameras took pictures and logged pressure at a lower frequency than the
pressure transducers, the DIC-pictures had to be synchronized with the complete
pressure curves logged by the pressure transducers. This was done by adjusting the
sets of pressure curves so that the pressure in transducer 5, had the same time of
arrival when logged by the transducers, as when logged by the camera. Because
no breakwire was installed, the time of arrival had to be adjusted according to the
calibration tests. The tests at stand-off distance 375 mm were corrected according
to the time of arrival of transducer 5 in test R23, while the tests at stand-off distance
500 mm were corrected according to transducer 5 in test R33-2. As no calibration
tests were performed at stand-off distance 625 mm, these tests had to be corrected
108
CHAPTER 4. EXPERIMENTAL WORK
in another way. At the time of detonation, an impulse were sent out from the deto-
nation mechanism. This impulse was captured by the pressure transducers. Because
the response only appeared in a few of the transducers, it was decided to correct all
the tests according the same time of arrival. The tests were therefore adjusted ac-
cording to transducer 5 from test A32. It should be noted that the potential source
of error is great, as the time of arrival of the pressure created by the detonation
mechanism, is unknown. Still, comparing the time of arrival to the values used for
the other tests, the value seems reasonable. Table 4.9 lists the time of arrival used
to correct the different tests.
After correcting the pressure curves according to the time of arrival, the pressure
transducers were compared. Figure 4.38 shows the pressure in transducers 1, 5, 6
and 7 for test A11. Figure 4.39 and Figure 4.40 show the same plots for test A21
and A31, respectively. Corresponding plots for the remaining tests can be found in
Appendix B.3.
4000
Pressure transducer 1
Pressure transducer 5
3500 Pressure transducer 6
Pressure transducer 7
3000
2500
Pressure [kPa]
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
Figure 4.38: Pressure transducers in the clamping frame for test A11.
109
CHAPTER 4. EXPERIMENTAL WORK
2000
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5
Time [ms]
Figure 4.39: Pressure transducers in the clamping frame for test A21.
1500
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
1000
Pressure [kPa]
500
−500
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Time [ms]
Figure 4.40: Pressure transducers in the clamping frame for test A31.
In all the above figures, it can be observed that the time of arrival of the pressure
transducers are very similar. In addition, the peak pressure is almost the same
between the pressure transducers. This indicates that the load have been centered
pretty well during the experiments. There are slightly more variations for a couple
of the other tests, but the results are still fairly consistent.
110
CHAPTER 4. EXPERIMENTAL WORK
From the DIC analyses, it is possible to extract the displacement of the midpoint.
This can then be plotted against time in order to investigate the behavior of the
plate. Figure 4.41 shows the midpoint displacement plotted against time for the
tests. Test A23 and A32 have not been included in the figure. This was because the
removed paint made it impossible to get accurate values after maximum positive
displacement. Based on the maximum displacements meassured, it can be assumed
that the behaviour of plate A23 was similar to plate A21, and the behavior of plate
A32 similar to A33.
50
A11
A12
40 A13
Displacement at midpoint [mm]
A21
A22
30 A31
A33
20
10
−10
−20
−30
0 2 4 6 8 10 12 14 16 18 20
Time [ms]
From Figure 4.41 it can be seen that the displacement stabilizes fairly quickly at
around 41 mm displacement for test A11. Similar behavior can be observed from all
the plates at 375 mm stand-off, although the maximum displacement vary at values
between 40 to 45 mm. It can be seen that the displacement seems to stabilize at
around 27 mm displacement for test A21. After reaching maximum displacement,
the plate returns back to almost 17 mm displacement before it starts to increase
again. Test A22, on the other hand, stabilizes at around 31 mm. Similar behaviour
as for plate A21 can be observed, although with less change in the displacements.
For test A31 and A33, it can be seen that the displacements start as positive. After
reaching maximum displacement, the plates snap back to a negative displacement.
Comparing the curves at stand-off 500 mm and at stand-off 625 mm, it seems like the
counter-intuitive behaviour seen from test A31 and A33, may result from the same
behaviour observed from plate A21, when the displacements decrease after initally
reaching the peak. It is possible that the plates at 625 mm stand-off reach an
instable configuration, where the plates continue their displacement in the negative
111
CHAPTER 4. EXPERIMENTAL WORK
112
CHAPTER 4. EXPERIMENTAL WORK
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(a) t = 0.39 ms (b) t = 0.39 ms (c) t = 0.39 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(d) t = 0.58 ms (e) t = 0.58 ms (f) t = 0.58 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(g) t = 0.77 ms (h) t = 0.77 ms (i) t = 0.77 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(j) t = 0.96 ms (k) t = 0.96 ms (l) t = 0.96 ms
113
CHAPTER 4. EXPERIMENTAL WORK
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(a) t = 0.59 ms (b) t = 0.59 ms (c) t = 0.59 ms
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(d) t = 0.78 ms (e) t = 0.78 ms (f) t = 0.78 ms
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(g) t = 0.97 ms (h) t = 0.97 ms (i) t = 0.97 ms
114
CHAPTER 4. EXPERIMENTAL WORK
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(j) t = 1.16 ms (k) t = 1.16 ms (l) t = 1.16 ms
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(m) t = 4.92 ms (n) t = 4.92 ms (o) t = 4.92 ms
115
CHAPTER 4. EXPERIMENTAL WORK
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
116
CHAPTER 4. EXPERIMENTAL WORK
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
117
CHAPTER 4. EXPERIMENTAL WORK
From Figure 4.42, Figure 4.43 and Figure 4.44, it can be seen that the positive
impulse ends very early compared with the displacement curves. That means most
of the deformations are caused by the impulse load. It should also be noted that
pressure transducer 5 is located in the clamping frame, above the plate. That means
the pressure actually hits the plate slightly earlier than what is shown in the figures.
Because the impulse is finished early in the displacement progress, it can be assumed
that there is very little fluid-structure interaction during the deformation.
From Figure 4.44, it can be seen that the negative impulse is likely to be finished
as the displacement start to decrease. It is possible that the negative impulse gives
a similar effect as the positive impulse, resulting in deformations after the impulse
has ended. However, looking at the pressure curves, this is unlikely. Due to pressure
from reflection waves, the negative impulse is almost nullified, making it unlikely
that this could be the reason for the negative displacements. In addition, very little
change in the displacements during the negative phase, can be observed. Due to
the duration of the negative phase, it should be expected that any effect would be
noticeable earlier than what can be observed from the figures.
Comparing Figure 4.43m and Figure 4.43o with Figure 4.44m and Figure 4.44o,
it can be seen that there are similarities between the deformation profiles. This
indicate that the behaviour would have been similar if the displacement had reached
an unstable configuration. After the configuration in Figure 4.43o, the plate returned
to a deformed configuration similar to Figure 4.43l.
When investigating the results from the blast pencil, it was apparent that the blast
pressure hits the cameras at around 10 ms. The exact time varies depending on
the stand-off distance. Because the pressure hits the cameras, it is possible that the
cameras are distorted, and that the results therefore will be affected. Looking at the
results from the DIC analyses, it was apparent that the later deformation profiles
were displaced slightly more along the whole cross-section. Consequently, the results
from the DIC has been assumed incorrect after around 10 ms. Because the plates
were still in motion at this point, it was not possible to get accurate deformation
profiles of the final deformation configuration from the DIC.
When comparing all the DIC-results, it can be seen that the behaviour is fairly
similar between the tests at corresponding stand-off distances. There are still slight
variations, and there are multiple possible causes for this. As the C4 loads were
made by hand, it is possible that there were slight variations in mass and geometry.
It is also possible that there are slight differences in the positioning of the loads, or
differences in the plates that were tested. Considering all the potential sources of
error, the experiments still indicate decent repeatability, given approximately the
same setup and conditions.
118
Chapter 5
Material model
119
CHAPTER 5. MATERIAL MODEL
would still be too hard to get an accurate result. Due to this, it is better to use
inverse modeling to determine the material model.
120
CHAPTER 5. MATERIAL MODEL
5.2.1 Mesh
Due to symmetry, it was decided to only model a quarter of the specimen. The mesh
was modeled in a program called SALOME [69]. The mesh created in SALOME,
is not directly compatible with EUROPLEXUS. Because of this, the mesh was
converted by the help of another program called Cast3M [70]. The resulting mesh
can be seen in Figure 5.1. The Cast3M-file used to convert the mesh, is included in
Appendix D.1. The mesh is composed of shell elements with a thickness of 0.8 mm.
The element size is varying between approximately 1.5 mm and 2.0 mm. At the
midpoint of the specimen, the element width is varying. This was done to ensure
that the elements were as square as possible at larger displacements. Q42L elements
were used. This is a 4-node quadrilateral element, that can be used in Lagrangian
calculations with full integration [46].
The displacement values logged by the machine were inaccurate. In addition, the
physical extensometer was placed over the middle 50 mm of the specimen. From
Figure 4.4, it can be seen that the fracture appears close to the shoulder of the
specimens. This means that the physical extensometer is unable to fully record the
121
CHAPTER 5. MATERIAL MODEL
122
CHAPTER 5. MATERIAL MODEL
5.3 Results
When testing different material models, the reference strain rate, ṗ0 , was set equal
to 1.0. C is usually determined by dynamic tensile tests. As no such tests were
performed, it was decided set this parameter equal to a value found in the literature.
C was then assumed equal to 0.014, as used in Chapter 3.
To reduce the computational time, all analyses were run with mass scaling. Initial
tests indicated that EUROPLEXUS were unable to run the analyses when the strain
hardening parameter, n, was set lower than 0.21. After removing the mass scaling,
it turned out that the analyses were unable to run when n were lower than 0.26.
Because of this, all alternatives that were tested, had n equal to 0.26 or larger.
After a lot of trial and error, several alternatives were found. Table 5.1 shows the pa-
rameters used in the different alternatives. Figure 5.4 shows the force-displacement
curves found in EUROPLEXUS, compared with the results from the DIC-analysis.
2000
1800
1600
1400
1200
Force [N]
1000
800
600
400
DIC
200 Spranghers et al.
Alternative 1
Alternative 2
0
0 1 2 3 4 5 6 7 8
Displacement [mm]
From Figure 5.4, it can be seen that both alternatives give pretty good results,
compared with the DIC. The material model used by Spranghers et al. [48], on
123
CHAPTER 5. MATERIAL MODEL
the other hand, appears to be too stiff. Looking at the employed parameter values,
it can be assumed that alternative 1 is the better alternative. A reason for this is
that the strain hardening parameter, n, is unusually large for alternative 2, while
alternative 1 has values in a more reasonable area.
Setting the strain-rate component of the Johnson-Cook model equal to 1, only the
strain dependent part is remaining. By inserting the parameters from Table 5.1 into
Equation (2.43), the true stress - plastic strain curve of the material model can be
plotted. Figure 5.5 shows the true stress generated by the material models. The
results from the material tests have not been included in the figure, as the results
from the material tests only cover very small strains. Plotting them in the same
figure would render them almost invisible, due to the large strains included in the
figure.
260
240
220
True stress [MPa]
200
180
160
140
From Figure 5.5, it can be seen that the behaviour of the models are a lot closer to
an asymptotic behaviour than the material model used by Spranghers et al. [48].
At the same time, the material models still have a very steep slope. Because of this,
it is possible that the material models are still too stiff.
124
Chapter 6
Numerical Simulations
6.1 Lagrange
When comparing experimental results, with results from numerical simulations,
there are a lot of factors that can affect them. This can include sources of error
from the experimental work, e.g. the placement of the load or deviation of the
mass. In addition, there can be differences due to the material model or other pa-
rameters, like the value used for the relative effectiveness factor. To investigate the
effect of some of these parameters, sensitivity studies were performed.
The test setup were the same as previously used in Chapter 3.3. It was decided to
stick with an element size of 10 mm, due to the low computational time, and to
keep the results comparable with the results from the FSI-analyses. In most of the
simulations, the stand-off distance was kept at 500 mm. This was decided, as stand-
off distance 500 mm produced neither fracture, nor counter-intuitive behaviour. The
load is still generated by the help of the AIRB command in EUROPLEXUS.
Procedure
In Chapter 5.2, two potential material models were established. It was concluded
that alternative 1 probably was the most suitable. However, it would still be in-
teresting to see how the other material models affect the results of the numerical
simulations. Because of this, simulations were performed with the established ma-
terial model from Chapter 5, as well as the material model presented by Spranghers
et al. [48]. All the material models, as well as the Johnson-Cook parameters, are
listed in Table 6.1.
125
CHAPTER 6. NUMERICAL SIMULATIONS
Results
Figure 6.1 shows the midpoint displacement, plotted against time for all the analyses.
In addition, the midpoint displacement from plate A22 is included for comparison.
40
35
30
Displacement [mm]
25
20
15
10
Spranghers et al.
5 Alternative 1
Alternative 2
A22
0
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
From Figure 6.1, it can be seen that all material models give larger deformations
than observed in the physical experiments. Alternative 1 and Alternative 2 give
almost identical results, while the material model used by Spranghers et al. [48],
gives slightly smaller deformations. Despite the fact that the material model used
by Spranghers et al. [48] produces displacements closer to the experimental values,
it was decided to stick with the material models established in Chapter 5.3. This
was decided because the inverse modeling indicated that Alternative 1 and 2 had a
behaviour closer to the actual material behaviour. Because both Alternative 1 and
2 gave almost the same results, it was decided to stick with Alternative 1, due to
the reasons mentioned in Chapter 5.3.
126
CHAPTER 6. NUMERICAL SIMULATIONS
Procedure
In Chapter 4.3.3, it was observed that the plates slided along the supports close to the
pressure transducers. In the area around the bolts however, very little deformations
were visible. Based on this, it can be assumed that the boundary conditions of the
plates were somewhere between fixed and simply supported. In order to investigate
the effect of the boundary conditions, it was decided to model the plate, both as
fixed and as simply supported.
When modeling the plate as fixed, the setup was kept the same as in previous
simulations. In the case of modeling a simply supported plate, it was decided to
neglect all effect from the clamping frame, besides preventing displacements out of
the plane. It was decided to only include effects from the bolts, when restricting in-
plane displacements. As the pressure transducers were made of plastic, it is unlikely
that they would be able to prevent much deformation anyway. Because of this, it
was decided to ignore them.
Figure 6.2 shows the mesh that were used when modeling the simply supported
boundary conditions. The middle 300x300mm2 section of the plate, was modeled
as when fixed. That means the mesh is made up of Q4GS shell elements. The
surrounding frame was meshed with triangular shell elements. This was done in
order to model the holes in the plate as circular as possible, without having to use
extremely small elements. The element type employed in EUROPLEXUS, is called
T3GS. This is a 3-node, thick shell element, with one integration point in the plane.
T3GS is considered to be the best suited element, in combination with the Q4GS
element [46].
127
CHAPTER 6. NUMERICAL SIMULATIONS
In Figure 6.2, some of the nodes around the holes have been marked in red. These
nodes have been fixed in all directions. Because the bolts would prevent displacement
of the outer edge of the holes, it has been assumed that the selected nodes do not
change their position. This simplification makes it redundant to model the actual
bolts. The other side of the bolt holes are still free to move in the plane.
All the nodes in the outer frame has been restricted against movement out of the
plane. This was done to represent the restriction made by the clamping frame. This
simplification has one significant weakness. Because the nodes are restricted against
movement out of the plane, it means that the nodes will still be restricted, even
though they move out of the area restricted by the clamping frame. As a result,
the model will not be able to describe the behaviour seen at the shorter stand-off
distances in Chapter 4.3.3. Despite this, the model should still give a decent basis
for comparing fully fixed and simply supported boundary conditions.
A boundary condition where friction restricts the movement, would have been an in-
teresting study. However, no tests were performed in order to determine the friction
coefficient. Because of the uncertainty, it was decided not to test this possibility.
Results
Figure 6.3 shows the midpoint displacement, plotted against time, for the two dif-
ferent boundary conditions. The results from test A22 are included for comparison.
128
CHAPTER 6. NUMERICAL SIMULATIONS
50
45
40
Displacement [mm]
35
30
25
20
15
10
5 Fixed
Simply supported
A22
0
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
From Figure 6.3, it can be seen that the midpoint displacement increases by about
8 mm when using simply supported boundary conditions. Figure 6.4 shows a pro-
jection of the deformed mesh.
From Figure 6.4, it can be seen that the deformation around the pressure transducers
are similar to what was observed in Chapter 4.3.3. The mesh is also deformed around
129
CHAPTER 6. NUMERICAL SIMULATIONS
the bolt holes, in contrast to what was observed in the physical experiments. Because
the nodes were restricted against movement out of the plane, it is possible that the
deformation shown in Figure 6.3, would have been even larger if this have not been
the case. Because both boundary configurations resulted in larger deformations than
in the physical experiments, it was decided to stick with fixed boundaries, as this
behaviour was closest to what was experienced at stand-off distance 500 mm.
Procedure
The charge was suspended by the use of a wire in the experiments. Even though the
placement of the wire was closely measured, it is likely that the placement of the
load had small variations between the tests. This could either be due to variation
in length of the wires, the wires not being completely centered, or the shape of the
charge causing the center of mass to vary. Because of this, it would be interesting
to see how small variations in the load placement would affect the simulations.
It was decided to investigate the effect of placing the charge slightly closer or further
away, from the plate than intended. In addition, it was decided to investigate the
effect if the charge had been offset slightly to the side. Because the AIRB approach
does not account for reflections from the surroundings, it is expected that the same
results would be found if the charge was slightly elevated or lowered. Table 6.2 lists
the investigated configurations. x is defined as the direction of the load, y is along
the horizontal, and z is along the vertical. The origin is defined as the midpoint of
the plate.
Location of charge
Alternatives
x [mm] y [mm] z [mm]
Reference 500 0 0
x + 10 mm 510 0 0
x – 10 mm 490 0 0
y + 10 mm 500 10 0
Results
Figure 6.5 shows the midpoint displacement when varying the location of the charge.
The results from test A22 have been included for comparison.
130
CHAPTER 6. NUMERICAL SIMULATIONS
40
35
30
Displacement [mm]
25
20
15
10
Reference
x + 10 mm
5 x – 10 mm
y + 10 mm
A22
0
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
From Figure 6.5, it can be seen that the deformations change slightly when varying
the placing of the explosive charge along the x-axis. It is unlikely that the variations
during the experiments were as high as 10 mm, but it is still possible that minor
variations have contributed to the differences between the experiments. At the same
time, the contribution is fairly low when the variations are small. Because of this,
it is possible that any effect from varying the load placement may be neglected
compared to other factors.
By varying y and keeping x constant, it can be seen that the displacement curve
almost completely overlaps the reference curve. In order to investigate this further,
the 2D cross-section were plotted at different time steps up to maximum displace-
ment. Figure 6.6 shows the 2D cross-section along the y-axis at different time steps.
The cross-sections from alternative "y + 10 mm" have been mirrored in order to
visualize the corresponding curves for a theoretical case of "y – 10 mm".
131
CHAPTER 6. NUMERICAL SIMULATIONS
40
t = 0.65 ms y + 10 mm
y – 10 mm
35
30
Displacement [mm]
25
20
t = 0.30 ms
15
10
t = 0.15 ms
5
0
−150 −100 −50 0 50 100 150
Cross section [mm]
Figure 6.6: Cross section with off-center load at different time steps.
From Figure 6.6, it can be seen that there are variations in the cross-sectional
displacement during the deformation. As the displacement gets closer to maximum,
the differences disappear. The reason for this seems to be that the deformations
follow the yield lines. That means that as long as the load is large enough to
produce yield lines all the way to the midpoint, the final deformations should be
fairly symmetric, even though the load is slightly off-centered.
Procedure
Literature reports different values for the relative effectiveness factor for C4. The
values used can vary a lot. CONWEP uses 1.190 or 1.370 depending on whether it
is the impulse or the peak pressure that is of interest [13]. The most common value
however, is 4/3. The relative effectiveness factor has been discussed in Chapter
2.2.2. In order to investigate the effect of the relative effectiveness factor, it was
decided to perform simulations with the most common value, as well as the values
used by CONWEP. Table 6.3 lists the alternatives that were investigated.
132
CHAPTER 6. NUMERICAL SIMULATIONS
Results
Figure 6.7 shows the midpoint displacement, for the different relative effectiveness
factors. The results from test A22 have been included for comparison.
40
35
30
Displacement [mm]
25
20
15
10
Standard = 4/3
5 Min = 1.190
Max = 1.370
A22
0
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
In Figure 6.7, it can be seen that the displacement increases, as the mass increases.
The increase when using a relative effectiveness factor of 1.370 is only about 1 mm.
Using 1.190 however, decreases the displacement quite a lot.
The curves in Figure 6.7 can also give an indication of what would happen if there
were slight variations is mass. As the loads in the experiments were created by hand,
133
CHAPTER 6. NUMERICAL SIMULATIONS
it is likely that this was the case. Looking at the standard curve, compared to the
max curve, it can be seen that an increase of about 1 g in mass, would increase
the displacement with about 1 mm. Based on this, it is possible that variations in
mass have been a contributing factor to the variations in the displacements from
the experiments in Chapter 4.3.3.
To investigate the pressure generated by the different relative effectiveness factors,
the pressure at the midpoint was plotted against time. This can be seen in Figure
6.8. The pressure from pressure transducer 1 in test R33-2, has been included for
comparison.
2500
Standard = 4/3
Min = 1.190
Max = 1.370
2000 R33-2
Pressure [kPa]
1500
1000
500
Figure 6.8: Pressure in the middle with different relative effectiveness factors.
From Figure 6.8, it can be seen that the pressure curves generated from the three
different relative effectiveness factors, are fairly similar. The curve from test R33-2
declines faster than the curves from the simulations, but the second peak, hitting
the plate at around 0.3 ms, makes up for some of the difference. In order to get a
closer look at the difference in peak pressure, it is necessary to zoom in on a smaller
area. Figure 6.9 shows a close-up view of the peaks of the pressure curves.
134
CHAPTER 6. NUMERICAL SIMULATIONS
Standard = 4/3
2000 Min = 1.190
Max = 1.370
R33-2
1800
Pressure [kPa]
1600
1400
1200
1000
800
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
Time [ms]
From Figure 6.9 it can be seen that the peak pressure is lower than in the exper-
iment for all the simulations. Note that the peak pressure is a little higher in the
simulations than shown in the figure. Due to the size of the time steps, the peak
value is not registered in the plots. The peak pressures shown in the figure, does
not explain why the displacement from the simulations are larger than from the
experiments. To investigate further, the positive impulse were calculated from the
curves. The resulting values are listed in Table 6.4.
In Table 6.4, it can be seen that the positive impulse is higher in the numerical
simulations than in the physical experiments. This is the case, even when using the
lowest value for the relative effectiveness factor, and including additional peaks as
part of the positive impulse for the physical experiments. Based on this, it is possible
that EUROPLEXUS overestimates the deformations because the AIRB approach
generates an impulse that is too high compared to the physical experiments. In
order to say more about this, it is necessary to investigate the pressure distribution
over the whole plate.
135
CHAPTER 6. NUMERICAL SIMULATIONS
6.2.1 Procedure
When comparing results between the simulations and the experiments, it was de-
cided to compare the pressure transducers in the plate, and ignore the ones in the
frame. The elements shown in Figure 6.10 were selected in order to extract the
pressure.
136
CHAPTER 6. NUMERICAL SIMULATIONS
6.2.2 Results
Figure 6.11, Figure 6.12 and Figure 6.13, show the pressure plotted against time for
stand-off 265 mm, 390 mm and 515 mm respectively.
9000
Pressure transducer 1
Pressure transducer 2
8000 Pressure transducer 3
Pressure transducer 4
7000
6000
Pressure [kPa]
5000
4000
3000
2000
1000
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time [ms]
(a) Pressure transducers along the vertical.
9000
Pressure transducer 1
Pressure transducer 8
8000 Pressure transducer 9
Pressure transducer 10
7000
6000
Pressure [kPa]
5000
4000
3000
2000
1000
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time [ms]
(b) Pressure transducers along the diagonal.
Figure 6.11: Pressure curves from simulation at stand-off distance 265 mm.
137
CHAPTER 6. NUMERICAL SIMULATIONS
3500
Pressure transducer 1
Pressure transducer 2
3000
Pressure transducer 3
Pressure transducer 4
2500
Pressure [kPa]
2000
1500
1000
500
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time [ms]
(a) Pressure transducers along the vertical.
3500
Pressure transducer 1
Pressure transducer 8
3000
Pressure transducer 9
Pressure transducer 10
2500
Pressure [kPa]
2000
1500
1000
500
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time [ms]
(b) Pressure transducers along the diagonal.
Figure 6.12: Pressure curves from simulation at stand-off distance 390 mm.
138
CHAPTER 6. NUMERICAL SIMULATIONS
1500
Pressure transducer 1
Pressure transducer 2
Pressure transducer 3
Pressure transducer 4
1000
Pressure [kPa]
500
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time [ms]
(a) Pressure transducers along the vertical.
1500
Pressure transducer 1
Pressure transducer 8
Pressure transducer 9
Pressure transducer 10
1000
Pressure [kPa]
500
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time [ms]
(b) Pressure transducers along the diagonal.
Figure 6.13: Pressure curves from simulation at stand-off distance 515 mm.
From the figures, it can be seen that the time of arrival increases and the peak
pressure decreases when moving away from the center of the plate. Note that the
peak pressures of pressure transducer 1 did not register due to the time step when
logging the results. As a result, it may seem from the figures that the peak pressure
is lower in pressure transducer 1 than in pressure transducers 2 and 8. This is not
the case though, but only a miss caused by the logging frequency.
139
CHAPTER 6. NUMERICAL SIMULATIONS
In order to get an impression of the load, compared to the calibration tests, the
positive impulse was calculated. Table 6.6 lists the positive impulse for pressure
transducer 1, 2, 3, 4, 8, 9 and 10 for the three simulations, as well as for test R13,
R23 and R33-2.
From Table 6.6, it can be seen that the differences, between experiments and simu-
lations, get smaller as the stand-off distance increases. The results from the exper-
iments show an increase in positive impulse from pressure transducer 8 to pressure
transducers 9 and 10. This is possibly due to accumulation of pressure in the cor-
ner, and is therefore not representative of the load under ideal conditions. Based
on this, the AIRB approach seems to generate pressure that corresponds well with
actual experiments when the stand-off distance increases past 500 mm. At distances
shorter than 500 mm however, the load is overestimated a lot.
140
CHAPTER 6. NUMERICAL SIMULATIONS
6.3 FSI
In Chapter 4.3.3 it was observed that most of the deformations of the plate hap-
pened after the positive load was finished. This means that the impulse from the
blast only initiates the deformations, while the inertia keeps the deformation going.
Based on this, it can be assumed that there is little fluid-structure interaction in the
physical experiments. In Chapter 6.1 it could be seen that the simulations with the
AIRB approach in EUROPLEXUS overestimated the displacements of the plate. It
is possible that the other methods explained in Chapter 2.11.5 can give a better
representation of the loading. FSI-simulations are therefore performed as both the
GAZP approach and JWL approach require FSI in order to represent the loading.
6.3.1 Mesh-Sensitivity
FSI-simulations are mesh-sensitive, and smaller element size will usually give more
accurate results [59]. As cubic elements are used, the number of elements will
increase cubically when decreasing the element size. In addition, the critical time
steps will decrease with the element size. Combined, this can result in a considerable
increase in the computational time. In order to find the element size best suited for
FSI-simulations, a mesh-sensitivity study were performed.
Procedure
Before running the simulations, a proper model has to be created. In the following
section, the x-axis will be assumed out of the plate, the y-axis will represent the
horizontal axis and the z-axis the vertical axis. The origin is assumed at the midpoint
of the plate.
By ignoring all the surrounding area, the setup presented in Chapter 4.3.1 can be
assumed symmetric over the xz-plane. Because of this it could have been possible
to model the setup with symmetry conditions over this plane. It would not have
been possible to model with symmetry conditions over the xy-plane. The reason for
this is that the floor will work as a reflecting boundary while the pressure can pass
freely over the top of the frame. Despite the possibility of symmetry, it was decided
to model the whole plate. This was done in order to make it possible to investigate
the effect of non-symmetric loading or other non-symmetric influences.
The model was created by the help of Cast3M. The Cast3M-file used to generate
the mesh is included in Appendix D.2.
The plate was modeled the same way as in Chapter 6.1, with Q4GS elements. The
frame was modeled as seen in Figure 6.14. Around the frame, 200 mm of air was
included. 300 mm of air was modeled at the back of the plate in order to include
possible suction effects from the back. The amount of air included in the front of the
141
CHAPTER 6. NUMERICAL SIMULATIONS
plate was dependent on the stand-off distance used, and how the explosive charge
was modeled. When modeling the fluid with 25 mm elements, and the charge as
shown in Figure 6.15d, 625 mm of air was included in the front of the frame. The
fluid mesh was modeled as finite volumes. All elements used for both the fluid and
the explosives, were assigned the element type CUVF in EUROPLEXUS. This is a
3D cubic finite volume element, that is defined as cell centered [46]. The frame and
the floor were modeled as reflecting boundaries, while the remaining sides of the
fluid mesh were modeled as absorbing boundaries. The absorbing boundaries were
assigned the CL3D element type.
200
Air
Frame
500
Plate
300
(0,0,0)
500
300
Charge
500
Stand-off distance
xair 300
(b) Seen from the side.
142
CHAPTER 6. NUMERICAL SIMULATIONS
Figure 6.15 shows the finished mesh when using element size 25 mm. In Figure
6.15a, the plate is marked in red. In Figure 6.15c, the plate is marked by a thick
black line, while the explosive charge is marked in red. Figure 6.15d shows a close-up
view of the explosive charge. The charge was modeled with the same dimensions
as used in Chapter 4.3.3. This means that the explosive charge was modeled as
a sphere with diameter 34.5 mm. In order to model the charge properly, smaller
elements were used for the charge than for the rest of the model. A transition zone
was made around the sphere, eventually leading out to the regular fluid mesh. As
the elements in the center of the charge were really small, the critical time step
would be very low. In order to prevent this, a EUROPLEXUS command called
PART, was activated. This command makes the analysis pay more attention to the
smaller elements without carrying out useless computation on the big ones [46].
The loading was simulated by the use of the GAZP command in EUROPLEXUS.
This method is explained in Chapter 2.11.5. The input parameters had to be calcu-
lated by the help of the equations explained there. The value used for γ was set to
143
CHAPTER 6. NUMERICAL SIMULATIONS
1.4 [59]. The values used for p0 and ρ0 are the same as used by Larcher et al. [36].
Due to the way the explosive charge is modeled, the volume of the balloon, Vbal , is
the volume used in the experimental work. Table 6.7 lists the parameters used in
the model, as well as parameters necessary to compute them.
Parameter Value
ET N T [J] 180800
p0 [Pa] 105
ρ0 [kg/m3 ] 1.3
Vbal [m3 ] 2.15 · 10−5
γ 1.4
pbal [Pa] 2.05 · 1010
ρbal [kg/m3 ] 588.1
From Chapter 4.3.3, it could be seen that the transition from positive to negative
displacement during the counter-intuitive behaviour, ended at around 5 ms. In order
to make sure that the displacement has stabilized, the simulations should be run a
little longer. Because of this, it was decided to use 10.5 ms as the total time interval
simulated.
The simulations performed in Chapter 6.1, used a mesh with element size 10 mm.
This was done as this was the element size expected to be used in the FSI-simulations.
Initial tests discovered that the computational time with 10 mm elements was very
high, which made it necessary to investigate the use of larger elements. Because the
loads had stand-off distances of 375 mm, 500 mm and 625 mm, it was desired that
the stand-off distance was dividable by the element size. This meant that 25 mm
elements and 5 mm elements were the best possibilities. 10 mm would also work
well, but would make it harder to model identical loads at the three stand-off dis-
tances. Table 6.8 lists the element sizes that were tested in the end. It was planned
to test element size 5 mm as well, with the plate using both 5 mm elements and 10
mm elements. However, this was not done due to lack of time.
The simulations were run with the embedded FSI approach. This means that the
plate can deform freely within the fluid mesh. Embedded FSI approach is explained
in detail in Chapter 2.11.4. The EUROPLEXUS input-file used for this analysis is
included in Appendix E.3.
144
CHAPTER 6. NUMERICAL SIMULATIONS
Results
When testing, it was apparent that the computational time increased more than
expected for the mesh with element size of 10 mm. Due to lack of time, it was
decided to stop the calculation when the displacement passed the maximum point.
Because of this, the total time interval simulated was only 1 ms.
Figure 6.16 shows the midpoint displacement from the two simulations. The results
from test A22 have been included for comparison. Note that the curve from the
experiment has been adjusted so that the displacement starts in the origin. This
was done as the displacements from the simulations start a lot earlier than for the
experiment when not corrected.
35
30
Displacement [mm]
25
20
15
10
5
25 mm elements
10 mm elements
A22
0
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
From Figure 6.16, it can be seen that maximum displacement increases by about
1.8 mm when using 10 mm elements compared with 25 mm elements. Both element
sizes still show lower displacements than the experiments. This is different from the
Lagrangian simulations performed in Chapter 6.1, as the simulations there produced
larger displacements than the experiments. If the displacement continue to increase
when reducing the element size further, it is possible that the displacement would
be fairly similar to the experiments.
The simulation with element size 10 mm was run for a shorter time interval. Table
6.9 lists the computational times for the two analyses.
145
CHAPTER 6. NUMERICAL SIMULATIONS
It should be noted that the computational times listed in Table 6.9, are not directly
comparable. The simulation with element size 25 mm, was run at a computer which
was used for other work at the same time. As the computer was not powerful enough
to run the analysis with element size 10 mm, this analysis had to be run at a more
powerful computer.
The computational time with element size 10 mm would be almost 70 times as large,
if the time continued to increase at the same rate. The reason for this is most likely a
combination of an increased number of elements, as well as a lower critical time step
due to smaller element size. The elements at the center of the explosive charge gets
particularly small. Another way of modeling the explosive charge could therefore
give much improvement to the computational time.
The version of EUROPLEXUS used for the simulations, did not support Massively
Parallel Processing (MPP). If this option had been available, it would have been
possible to decrease the computational time substantially. With the current ver-
sion, even using a model employing symmetry, the computational time would be
unreasonably high at element size 10 mm or lower. Because of this, it was decided
to stick to an element size of 25 mm for the FSI-simulations.
Modeling the explosive charge with exact dimensions, results in very small elements
in the center of the explosive charge. This increases the computational time, and it
is therefore interesting to see whether it is possible to model the explosive in a more
efficient way.
Procedure
When testing different load shapes, it was decided to keep the stand-off distance
constant at 500 mm. In addition, it was decided to stick with 25 mm elements in
order to reduce the computational time needed for each simulation. The volume of
the load, was kept as close to the actual volume as possible. Figure 6.17 shows the
different load shapes that were investigated. The figures presented are vertical cross
sections.
146
CHAPTER 6. NUMERICAL SIMULATIONS
(d) Detail
The alternative presented in Figure 6.17a, was made up of 8 elements with element
size of 25 mm. This was the lowest possible number of elements to use, if the load
was going to be symmetrical.
Figure 6.17b is based on the mesh in Figure 6.17d. In order to model a box with
correct volume, the transition zone, as well as the cubic mesh in the center of the
charge, were expanded. The cubic mesh was expanded such that it was possible to
model the explosive charge as different shapes, with correct volume. In order to get
the correct volume, the small box presented in Figure 6.17b, was composed of 512
elements with element size 3.4758 mm.
147
CHAPTER 6. NUMERICAL SIMULATIONS
The mesh presented in Figure 6.17c, is similar to the one in Figure 6.17b. Instead of
modeling the charge as a box, the elements were selected by the help of the COND
SPHERE command in EUROPLEXUS. This command selects all elements within
a sphere defined by a radius and a center point. By using the correct radius in the
same mesh, as for Figure 6.17b, 480 elements were selected. In order to get the
correct volume, the element size was increased to 3.5514 mm.
The fluid mesh were modeled such that there was one row of 25 mm elements along
the absorbing boundary, in the back of the model. That means the length, xair ,
shown in Figure 6.14b, is equal to 550 mm for the the alternative in Figure 6.17a,
675 mm for the alternative in Figure 6.17b and Figure 6.17c, and 625 mm for the
alternative shown in Figure 6.17d.
Table 6.10 shows the volumes of the charges, as well as the GAZP parameters used in
the simulations. The GAZP parameters were calculated by the equations presented
in Chapter 2.11.5.
Results
Figure 6.18 shows the displacement of the midpoint, plotted against time, for all the
simulations. The results from test A22 have been included for comparison.
148
CHAPTER 6. NUMERICAL SIMULATIONS
35
30
Displacement [mm]
25
20
15
10
Big box
5 Small box with transition zone
Small sphere with transition zone
Detail
A22
0
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
In Figure 6.18, it can be seen that the big box model gives displacements slightly
larger than the physical experiments, while the other simulations give lower dis-
placements. While the detail results in larger displacements than the small box,
the small sphere results in even larger ones. This indicates that both the volume
and the shape of the charge are important for the pressure generated by the blast.
In order to get a closer look at the reason for the differences in displacement, the
pressure at the midpoint was plotted against time. This is presented in Figure 6.19.
It should be noted that the pressure curves were extracted from one of the fluid
elements at the middle, right next to the plate. As the embedded FSI approach
was used, the selected fluid element did not move during the displacement of the
plate. This means that the pressure is not necessarily representative for the pressure
hitting the plate.
149
CHAPTER 6. NUMERICAL SIMULATIONS
2000
Big box
1800 Small box with transition zone
Small sphere with transition zone
1600 Detail
Relative pressure [kPa]
1400
1200
1000
800
600
400
200
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
From Figure 6.19, it can be seen that the variations between the tests, correspond
to the variations seen in Figure 6.18. Because the pressure is extracted from a fluid
element that does not follow the plate, it is hard to define the end of the positive
impulse affecting the plate. Based on the results from the calibration tests, as well as
the displacement curves in Figure 6.18, it can be assumed that the positive impulse
ends at around 0.5 ms for the three simulations. It is possible that the second peak
is a result of pressure reflected from the plate. By assuming the first peak as the
positive impulse, the values listed in Table 6.11 were calculated. The result from
test R33-2 is included for comparison.
From Table 6.11, it can be seen that the positive impulse is larger for the big box than
for test R33-2, and lower for the rest. This corresponds well with the displacements.
The impulse from the big box is a lot higher than from test R33-2. Comparing the
150
CHAPTER 6. NUMERICAL SIMULATIONS
displacements, the difference is not too large. The difference is a lot bigger for the
other simulations, when comparing the difference in impulse and in displacements.
Based on the observations from Chapter 3.2 and Chapter 6.3.1, it is possible that
the differences had been more reasonable with a smaller element size. Considering
the assumptions made when calculating the positive impulse, it is difficult to draw
any clear conclusions.
Table 6.12 lists the computational time of the different analyses. It should be noted
that the workload on the computer were varying when running the simulations.
Because of this, the time can only be used as an estimate.
Computational time
Alternatives
[hh:mm:ss]
Big box 02:15:22
Small box with
02:07:08
transition zone
Small sphere with
02:26:56
transition zone
Detail 05:12:31
From Table 6.12, it can be seen that the time increases a lot for the simulation with
detailed model of the charge load. This is most likely due to the very small element
size in the center of the charge. The time when simulating the charge as a big box,
is longer than the time for the small box and sphere. This is unexpected, due to
an increased amount of elements, as well as smaller elements in the small box and
small sphere models. Because of the small differences, it is likely that the big box
model was run at a time when there were less memory available, than for the other
simulations.
In EUROPLEXUS, there are several different ways to simulate the explosive load.
In addition to the AIRB and GAZP commands, it is also possible to use BUBB or
JWL. BUBB works in the same way as GAZP, but calculates the input parameters
automatically based on the mass of the load and the volume of the elements. JWL
models the explosive charge with the Jones-Wilkins-Lee equations of state. The
available methods are explained further in Chapter 2.11.5.
151
CHAPTER 6. NUMERICAL SIMULATIONS
Procedure
When investigating the different approaches, it was decided to use detail model
from Chapter 6.3.2. When using the JWL approach, the input parameters listed in
Chapter 2.3 were used. In order to get the correct representation of the load, the
volume had to be exact. The reason for this, is that the equation of state does not
take a theoretical volume and TNT-equivalent mass into consideration, but rather
assigns an equation of state specific to C4 and its abilities. The input parameters
for the JWL equation is listed in Table 2.2.
The GAZP model uses the same input parameters as in Chapter 6.3.2. The BUBB
model calculates the input parameters automatically based on a TNT-equivalent
mass. Because of this, the TNT-equivalent mass was the only input. Table 6.13
lists the alternatives investigated.
Method
GAZP
BUBB
JWL
Results
Figure 6.20 shows the midpoint displacement plotted against time for the different
simulation methods. The results from test A22 have been included for comparison.
35
30
25
Displacement [mm]
20
15
10
−5
GAZP
−10 BUBB
JWL
A22
−15
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
152
CHAPTER 6. NUMERICAL SIMULATIONS
In Figure 6.20, it can be seen that the displacements from the GAZP model and the
BUBB model are very similar. The BUBB model produces a little more displacement
than the GAZP model. It is possible that the BUBB model uses a different value
for the heat capacity ratio, γ, as well as a different function in order to calculate the
initial balloon pressure, pbal .
The JWL model results in a lot smaller displacements than the GAZP and BUBB
models. Larcher et al. [36] states that a fine mesh is essential when using the
JWL model. Ideally, the element size around the charge should be a lot smaller
than the ones used to model the explosive charge. This means that this particular
approach is especially mesh sensitive, which could explain the lower deformations.
An interesting observation from the JWL curve, is that the JWL model is able to
simulate the counter-intuitive behavior, observed for the aluminum plates at stand-
off distance 625 mm.
Figure 6.21 shows the pressure at the midpoint, plotted against time, for the different
approaches. As in Chapter 6.3.2, the pressure is extracted from an element in
the fluid mesh, next to the center of the plate. As the model uses the embedded
approach, this element does not move together with the plate.
2000
GAZP
1800 BUBB
JWL
1600
Relative pressure [kPa]
1400
1200
1000
800
600
400
200
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
From Figure 6.21, it can be seen that the pressure from the GAZP model and the
BUBB model are almost the same. Due to the time steps of the simulations, the
peak pressure has not been registered on the curves. The pressure of the JWL curve,
is much lower than the other models. It can also be seen that the JWL curve has a
fairly long negative phase. It is therefore possible that the displacements in Figure
6.20, may have been influenced by this negative phase.
153
CHAPTER 6. NUMERICAL SIMULATIONS
By making the same assumptions as in Chapter 6.3.2, the positive impulse can be
calculated. Table 6.14, lists the positive impulse calculated for the three methods.
The value from test R33-2 has been included for comparison.
From Table 6.14, it can be seen that the impulse is slightly higher for the BUBB
model than for the GAZP model. This explains the slightly larger displacements
produced by the BUBB model compared to the GAZP model.
Table 6.15 lists the computational times for the different models. The times are not
directly comparable as the workload on the computer varied between the simulations.
Computational time
Alternatives
[hh:mm:ss]
GAZP 05:12:31
BUBB 06:04:44
JWL 03:06:27
From the times listed in Table 6.15, it can be seen that the GAZP model and
BUBB model are pretty close. This is expected, considering the fact that both
models are based on the same equation of state. The time from the JWL model,
on the other hand, is much lower. It should be noted that an ideal model of JWL
probably would have a lot smaller elements. This means that the computational
time for a fully functional JWL model, probably would be a lot longer than for a
fully functional balloon model.
In Chapter 2.11.3 and Chapter 2.11.4, the difference between an embedded FSI
approach, and ALE is explained. It was decided to run simulations with both
approaches in order to investigate the difference.
154
CHAPTER 6. NUMERICAL SIMULATIONS
Procedure
The difference between an embedded FSI approach and ALE, will be most important
when reaching fracture in the model. In addition, embedded will usually be better
when dealing with large deformations. In order to get the largest deformations
possible, it was decided to model the experiments performed at stand-off distance
375 mm. In addition, it was decided to model the explosive charge as the "big box"
used in Chapter 6.3.2.
In order to perform an ALE simulation, the fluid elements around the plate have
to be allowed to deform. The elements selected are shown in Figure 6.22. All the
selected elements are within a box with dimension 200x300x300 mm3 . The center
of the box is in the origin of the plate. Selected elements are colored blue in the
figure. The nodes in the plate are marked in red. The EUROPLEXUS input-file for
the ALE-simulation, is included in Appendix E.4.
Results
Figure 6.23 shows the midpoint displacement, plotted against time, for the two
FSI-approaches. The results from test A22 have been included for comparison.
155
CHAPTER 6. NUMERICAL SIMULATIONS
50
45
40
Displacement [mm]
35
30
25
20
15
10
5 Embedded
ALE
A22
0
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
From Figure 6.23, it can be seen that both simulations result in larger displacements
than the experiments. This is expected due to the choices that were made when
running the simulations. It is interesting however, that the ALE-simulations result
in smaller displacements than the embedded FSI-approach.
Figure 6.24 shows the fluid mesh around the plate for both the ALE-simulation and
the embedded FSI-approach. The nodes in the plate are marked in red. Both of the
figures present the displacements at the end of the time series.
156
CHAPTER 6. NUMERICAL SIMULATIONS
From Figure 6.24b, it can be seen that the mesh is deformed around the plate. The
elements on one side of the plate, have been stretched out, while the ones on the
other side have been compressed. The displacements present, are not large enough
to cause any major problems. However, if the displacements are even larger, it may
eventually become a problem. Looking at the mesh in Figure 6.24a, it can be seen
that the plate has deformed without affecting the surrounding mesh. As the fluid
mesh in ALE follows the structure when deforming, it causes problems when fracture
occurs. This is because the nodes are splitting and the fluid mesh then has trouble
following the nodes. In problems where large deformations occur, or fracture might
be a problem, an embedded FSI-approach is the preferred method.
157
CHAPTER 6. NUMERICAL SIMULATIONS
6.4.1 Procedure
In Chapter 6.1 it was seen that the Lagrange simulations, generally overestimated
the displacements. On the other hand, it was observed in Chapter 6.3, that the
FSI-simulations resulted in too small deformations. Because of this, it was decided
to run simulations of the three stand-off distances used in Chapter 4.3.3, with both
the Lagrangian approach as well as with en embedded FSI-approach.
The model for the Lagrangian approach was the same as used in Chapter 6.1.3.
This means that the plate was modeled with 10 mm elements and fixed boundary
conditions. Alternative 1, from the inverse modeling, was used as material model,
and 4/3 was used as the relative effectiveness factor of the charge load. Table 6.16
lists the important choices made for the analyses.
Plate
Element size 10 mm
Boundary conditions Fixed
Material constants
Young’s modulus, E 71 GPa
Density, ρ 2710 kg/m3
Poisson’s ratio, ν 0.33
Johnson-Cook parameters
A 115 MPa
B 57 MPa
n 0.26
C 0.014
Loading
EUROPLEXUS approach AIRB
TNT equivalent mass 40 g
Stand-off distances 375 mm
500 mm
625 mm
The model used for the embedded FSI-approach, was the same as used in Chapter
6.3.1. Element size of 25 mm were used, for both the plate and the fluid. As in the
Lagrangian simulation, the plate was modeled with fixed boundaries. The charge
158
CHAPTER 6. NUMERICAL SIMULATIONS
load was modeled with exact dimensions and with the use of the GAZP approach.
Table 6.17 lists the important choices made for the analyses.
Plate
Element size 25 mm
Boundary conditions Fixed
Fluid
Element size 25 mm
Material constants
Young’s modulus, E 71 GPa
Density, ρ 2710 kg/m 3
Johnson-Cook parameters
A 115 MPa
B 57 MPa
n 0.26
C 0.014
Loading
EUROPLEXUS approach GAZP
TNT equivalent mass 40 g
Stand-off distances 375 mm
500 mm
625 mm
6.4.2 Results
Figure 6.25 shows the results from the simulations at 375 mm stand-off, compared
with the results from test A11. Figure 6.26 shows the results from the simulations at
500 mm stand-off, compared with the results from test A22. Figure 6.27 shows the
results from the simulations at 625 mm stand-off, compared with the results from
test A31. Note that the deformation curves from the component tests have been
adjusted so that the deformation starts at t = 0 ms. This was done in order to give
159
CHAPTER 6. NUMERICAL SIMULATIONS
the best basis for comparison, as the displacements from the physical experiments,
start fairly late, compared with the results from the simulations.
70
60
Displacement [mm]
50
40
30
20
10
Lagrange
FSI
A11
0
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
40
35
30
Displacement [mm]
25
20
15
10
5 Lagrange
FSI
A22
0
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
160
CHAPTER 6. NUMERICAL SIMULATIONS
30
25
20
Displacement [mm]
15
10
−5
−10
−15 Lagrange
FSI
A31
−20
0 1 2 3 4 5 6 7 8 9 10
Time [ms]
In Figure 6.25, it can be seen that the displacements from the Lagrangian analysis
is way too large compared with the experiments. Looking at Figure 6.26, it can be
seen that the displacements are getting closer to the experiments, and in Figure 6.27,
the maximum positive displacement is very similar for the Lagrangian analysis and
the experiment. However, the Lagrangian analysis does not manage to represent
the counter-intuitive behaviour present at this stand-off distance.
Looking at the displacements resulting from the FSI-simulations, it can be seen that
the displacements in general are lower than from the experiments. It is possible that
this is due to the element size, and would therefore improve with a finer mesh. An
interesting observation from the figures, is that the maximum displacements from
the FSI-simulations actually vary fairly little between the stand-off distances.
161
CHAPTER 6. NUMERICAL SIMULATIONS
162
Chapter 7
Discussion
Calibration tests
The calibration tests in general provided fairly good results, as seen in Chapter
4.3.2. By comparing pressure measured in transducers, decreasing peak pressures
and increasing values for the time of arrivals could be observed when moving away
from the center of the plate. This was the case for most of the tests.
In some of the tests, one or two pressure transducers provided results that deviated
from the rest, e.g. unusually low peak pressure or that the pressure arrived at the
pressure transducers in unexpected order. It can be seen that the time of arrival is
different in the pressure transducers in the clamping frame, for some of the tests.
As all three pressure transducers have the same stand-off distance, it is expected
that the pressure should arrive at the same time.
Comparing the different tests at the same stand-off, it can be seen that there are
slight variation in peak pressure and positive impulse. For most pressure transduc-
ers, the values are within range of each other, and it can clearly be seen which results
are from which stand-off distance. The time of arrival is pretty similar between the
tests, but the peak pressure and positive impulse have a lot more variation.
By inspecting the pressure curves, it can be seen that most of the curves have
additional pressure peaks right after the initial peak. It is possible that these peaks
are results of reflection pressure from the floor or other parts of the surroundings.
The second peak however, occurs very shortly after the first one. It is hard to say
what surface would be able to reflect the pressure this early. It is therefore possible
that this has another explanation. Spranghers et al. [63] used a very similar test
setup, but did not experience this behavior.
It can be seen that the pressure and positive impulse in the pressure transducers
close to the edge of the plate, tend to have higher values than the ones a little closer
to the center. It is likely that this happens due to an accumulation of pressure,
created when the pressure moves out to the side and is stopped by the thickness of
163
CHAPTER 7. DISCUSSION
the frame.
There are multiple possible sources of error that can have caused the deviations in
the results. The primary source is probably related to the shape and location of the
blast load. As the charge was created by hand, it is possible that both the mass and
geometry vary between the tests. To ensure that the mass and shape were the same,
it would have been better if the charges had been casted. Another likely source of
error, is the location of the explosive charge. The charge was elevated by a wire that
was not constant between the tests. Because of this it is possible that there have
been variations, both in terms of the centering of the load, as well as the stand-off
distance. A variation in the centering, would explain the variations in time of arrival
of the pressure transducers in the frame. Slight variations in the stand-off could also
explain the variations in peak pressure between the tests.
As explained in Chapter 4.3.2, the pressure transducers had drifting of the zero level
due to temperature. Based on the pressure curves, it can be seen that the amount
of drifting varies between the stand-off distance. The change is biggest at the lowest
stand-off distance. This is most likely because the temperature is higher when the
stand-off distance is low. As the heat created by the explosion moves slower than
the pressure, it is assumed that the drifting did not have much effect on the positive
phase, but it is likely that this is a problem for the negative phase. As the zero
level has changed, it is very hard to determine the end of the negative phase, and
therefore the negative impulse.
Component tests
In Chapter 4.3.3, it can be seen that the component tests provided good results,
but there were still minor differences in the displacement curves when comparing
the tests at the same stand-off distances. For the tests at stand-off distance 375
mm, the difference was only a couple of millimeters between the tests. Fractures
occurred several places along the supports, for the plates at stand-off distance 375
mm. The DIC-pictures indicate that this happens due to shear at the support. As
the fractures occur before the plate reaches maximum displacement, it is possible
that the fractures have contributed to an increase in the midpoint displacement.
The DIC paint was destroyed at one of the plates at stand-off distance 500 mm, and
one plate at stand-off distance 625 mm. This made it hard to get a good impression
of the displacement development of the two plates. The permanent deformations
measured still indicate that the two plates behaved similar to the other plates at
the same stand-off distance.
For the plates at stand-off distance 500 mm, it can be seen that plate A21 and A22
had around 4 mm difference in the final displacement of the midpoint. In addition,
it can be seen that the displacement development is different.
The plates at stand-off distance 625 mm ended with negative displacement. Similar
164
CHAPTER 7. DISCUSSION
behavior was investigated by e.g. Flores-Johnson et al. [16]. The reason for this
may be a combination of many factors. The negative impulse may have some effect.
In addition, it is possible that some change in displacement is caused by the release
of elastic deformations in the plate. Similar behavior can be seen in test A21, but in
this case the displacements return back to a positive. It is possible that the counter-
intuitive behaviour happens if the deformation reaches an unstable state, where the
inertia of the plate causes the plate to displace in the opposite direction. A lot is
unknown about the phenomenon, and it is therefore difficult to determine the exact
cause.
There are multiple possible causes for the slight differences in the displacement of
the plates. One possibility is that there might be minor differences between the
plates in terms of material quality, or the geometry of the plates. When measuring
the thickness of the plates, the differences were pretty small. If it is assumed that
the quality of the plates were similar, differences in load is the most likely cause
of the variations. This can be due to slight changes in mass, or the placement of
the load. From the simulations in Chapter 6.1.3, it can be seen that there are little
difference in deformation if the load is off-center. If the explosive charge is slightly
closer or further away, the displacement will change. Similar could be seen from the
simulations in Chapter 6.1.4. The results indicate that slight changes in mass can
cause changes in the displacement.
Analytical methods
In Chapter 4.3.3 it can be seen that the midpoint displacement of the plates are a
lot larger than the displacements calculated in Chapter 3.1. This is despite the fact
that maximum impulse was assumed distributed over the whole plate. Looking at
the impulses calculated from the experiments in Chapter 4.3.2, it can be seen that
the positive impulse at the middle of the plate are larger than the impulses used
in Chapter 3.1. This indicate that CONWEP underestimates the positive impulse
from the blast.
Figure 7.1 shows the results from the experiments compared with the empirical
relation established by Nurick and Martin [6]. The impulse has been taken from
the calibration tests, and the displacement has been taken from the component test.
The impulse used is not directly representative for the displacement. This is due
to the fact that the stand-off in the calibration tests were slightly larger due to
the thickness of the frame. In addition, the value used is not representative for
the variations between the tests. In order to keep the sources of error as small as
possible, it was decided to use the positive impulse from the experiments with the
median value. Because of this, the impulse from test R21 and R31-2 were used.
The impulse at the middle of the plate has been assumed constant and distributed
across the whole plate. Only stand-off distances 375 mm and 500 mm are included,
as these were the only stand-off distances tested in both the calibration tests and
the component tests.
165
CHAPTER 7. DISCUSSION
60
50
Displacement/thickness ratio
40
30
20
10
Empirical relation
375 mm stand-off
500 mm stand-off
0
0 10 20 30 40 50 60 70 80 90 100
Non-dimensional impulse φq
From Figure 7.1 it can be seen that the results from the tests at stand-off distance
375 mm are missing the empirical relation by a lot. It is possible that this is due
to increased deformation beacuse of the fractures along the support. The tests
at stand-off distance 500 mm however, fit fairly well. When calculating the non-
dimensional impulse, total impulse over the plate should be used. In Figure 7.1, the
total impulse is assumed to be the same as the impulse at the middle multiplied by
the area. This is most likely higher than the total impulse actually present. This
means that the non-dimensional impulse in reality should have a lower value. Based
on the difference in the values between the stand-off distances, it is still possible
that the empirical relation would give a better estimate at smaller displacements.
In order to investigate this, a more accurate value has to be measured for the total
impulse over the plate.
Numerical simulations
From the results in Chapter 3.3 and Chapter 6.1, it can be seen that the Lagrangian
simulations overestimated the displacements. Based on the result from Chapter 6.2,
it can be assumed that the main reason for this, is that the AIRB command that
generates the load, overestimates the load applied to the structure. In Chapter 6.2
and Chapter 6.4, it can be seen that the pressure and displacement gets closer to
the experiments when the stand-off distance increases.
The AIRB command in EUROPLEXUS is based on the empirical equations estab-
lished by Kingery and Bulmash [33]. This is explained further in Chapter 2.2.7.
166
CHAPTER 7. DISCUSSION
167
CHAPTER 7. DISCUSSION
168
Chapter 8
Concluding Remarks
In this thesis, experiments and numerical simulations have been performed in order
to investigate the response of thin aluminium plates exposed to blast loading. The
goal was to investigate what numerical or analytical methods could provide the best
results when simulating the response.
Material tensile tests were performed in order to determine a material model. Neck-
ing occurred at low plastic strains. As a result, the material model had to be
determined through inverse modeling. Inverse modeling was performed with the
computer code EUROPLEXUS. Due to limitations in the code, the material prop-
erties were identified and fitted to the Johnson-Cook material model.
Experiments were performed with free air explosions generated by using 30 g of
composition C4. The load was applied in calibration tests with a stiff steel plate,
as well as in component tests of thin plates made of the commercial aluminum alloy
EN AW-1050A-H14.
The calibration plate was outfitted with pressure transducers. These were used to
measure the pressure at different points on the surface of the plate. Examining the
results revealed that the time of arrival was consistent between the experiments,
while the peak pressure and positive impulse varied a lot more between the tests.
It is probable that some of the variation is a result of variation in mass and shape
of the explosive charges, and some is likely due to variations in the location of the
load.
The component tests showed that a stand-off distance of 375 mm produced fracture
along the support. It is possible that the fractures are a result of shear failure. The
plates at 625 mm stand-off showed counter-intuitive behavior, as the final deformed
configuration had deformation towards the explosive charge, instead of away from
the charge as observed at the other stand-off distances. All the stand-off distances
indicated a good degree of repeatability, as there were only minor differences in
the displacements. Numerical simulations performed, indicates that the differences
can be a result of variation in mass, or deviations in the stand-off distance. The
169
CHAPTER 8. CONCLUDING REMARKS
deformation appears to be following the yield lines, which causes the deformation to
be similar, even if the explosive charge should be somewhat off-center. By comparing
the pressure curves with the displacement curves, it is apparent that most of the
deformation happens after the positive impulse is finished. Based on this, it can be
assumed that the deformations happens as a result of impulsive load, and that very
little fluid-structure interaction is present during the deformation.
The analytical approach investigated, appears to underestimate the displacements
from the experiments. Results indicate that this may be due to the large displace-
ments present in the experiments. From the tests at stand-off distance 500 mm, it
appears that the analytical approach might give better indication of the response
when the expected displacements are smaller.
Numerical simulations were performed in EUROPLEXUS. The results show that
neither the Lagrangian, nor the embedded FSI approach give convincing results
when simulating the component tests. The AIRB approach, used to simulate the
load in the Lagrangian simulations, seem to overestimate the load on the plate. As
a result, the deformations resulting from the simulations are too large. When the
stand-off distance increases, the load gets closer to what was measured during the
experiments. It is therefore possible that the approach will be able to represent the
loading better, if the stand-off distance is increased even further.
The embedded FSI approach resulted in deformations that were too low, compared
to the experiments. Comparing the pressure and impulse, indicated that the GAZP
approach was able to give a pretty good representation of the load. Mesh-sensitivity
studies indicated that the element size was important when determining the dis-
placements from the simulations. It is therefore possible that the embedded FSI
approach would be able to give a good representation of the experiments using a
smaller element size.
170
Chapter 9
Further Work
During the work with this thesis, there have been revealed several aspects that
could be interesting to investigate further. This includes aspects related to the ex-
perimental work, as well as to the analytical and numerical computational methods.
Following is a list of suggestions for further work:
• More experimental work should be performed in order to get a larger set of
results to use for comparison with computational methods. This includes both
calibration tests and component tests, in order to get results for both pressure
and displacement. It should be prioritized to do experiments at larger stand-
off distances, but experiments with variation in mass could be interesting as
well. The explosive charges should be casted in order reduce the chance of
large geometric deviations.
• Additional tests should be performed in order to investigate the causes of the
additional positive peaks in the pressure curves. This particularly applies for
the second pressure peak, as the time between this peak and the initial positive
peak is very low. It is therefore possible that the peak is a result of something
else than reflection waves from the surroundings.
• The drifting of the zero level of the pressure transducers, should be investigated
in order to see if it is possible to predict when the zero level changes, and
by how much. By doing this, it would be possible to get a more accurate
understanding of the load conditions.
• The counter-intuitive behaviour observed at stand-off distance 625 mm, should
be investigated further, in order to identify possible causes.
• The boundary conditions of the experimental setup should be investigated
in order to make it possible to do numerical simulations with more accurate
boundary conditions.
• Numerical simulations that take fracture into consideration, should be per-
formed in order to investigate the effect of the fractures in the plates at stand-
171
CHAPTER 9. FURTHER WORK
172
Chapter 10
Bibliography
[1] Applied Research Associates Inc, “Showcase Project: Finite Element Simula-
tion for Blast Response of Structures.” URL http://www.ara.com/Projects/
SVO/blast_response.html. Accessed: May 12, 2014.
[2] S. B. Menkes and H. J. Opat, “Broken beams - tearing and shear failures
in explosively loaded clamped beams,” Experimental Mechanics, vol. 13(11),
pp. 480–486, 1973.
[3] R. G. Teeling-Smith and G. N. Nurick, “The deformation and tearing of thin
circular plates subjected to impulsive loads,” International Journal of Impact
Engineering, vol. 11(1), pp. 77–91, 1991.
[4] G. N. Nurick and G. C. Shave, “The deformation and tearing of thin square
plates subjected to impulsive loads - an experimental study,” International
Journal of Impact Engineering, vol. 18(1), pp. 99–116, 1996.
[5] G. N. Nurick and J. B. Martin, “Deformation of thin plates subjected to im-
pulsive loading - a review. Part I: Theoretical considerations,” International
Journal of Impact Engineering, vol. 8(2), pp. 159–170, 1989.
[6] G. N. Nurick and J. B. Martin, “Deformation of thin plates subjected to impul-
sive loading - a review. Part II: Experimental Studies,” International Journal
of Impact Engineering, vol. 8(2), pp. 171–186, 1989.
[7] G. N. Nurick, M. E. Gelman, and N. S. Marshall, “Tearing of blast loaded
plates with clamped boundary conditions,” International Journal of Impact
Engineering, vol. 18(7-8), pp. 803–827, 1996.
[8] T. Wierzbicki and G. N. Nurick, “Large deformation of thin plates under
localised impulsive loading,” International Journal of Impact Engineering,
vol. 18(7-8), pp. 899–918, 1996.
173
CHAPTER 10. BIBLIOGRAPHY
174
CHAPTER 10. BIBLIOGRAPHY
175
CHAPTER 10. BIBLIOGRAPHY
176
CHAPTER 10. BIBLIOGRAPHY
Lecture_Notes_files/FEA_Lectures/Chapter_1_%20Introduction.pdf.
Accessed: April 27, 2014.
[51] Wikipedia, “Conservation law.” URL http://en.wikipedia.org/wiki/
Conservation_law. Accessed: March 10, 2014.
[52] K. M. Hsiao, “Corotational total lagrangian formulation for three-dimensional
beam element,” AIAA Journal, vol. 30(3), pp. 797–804, 1992.
[53] Wikipedia, “Newtonian fluid.” URL http://en.wikipedia.org/wiki/
Newtonian_fluid. Accessed: March 26, 2014.
[54] L. Olovsson, Corpuscular method for airbag deployment simulations. 6th Euro-
pean LS-DYNA User’s Conference.
[55] Wikipedia, “Numerical integration.” URL http://en.wikipedia.org/wiki/
Numerical_integration. Accessed: May 31, 2014.
[56] K. M. Mathisen, “Lecture 7: Solution of the dynamic equilibrium equations
by explicit direct integration.” Lecture Notes in TKT4197 Nonlinear Finite
Element Analysis, NTNU (Norwegian University of Science and Technology),
2013.
[57] K. M. Mathisen, “Lecture 9: Solution of the nonlinear dynamic equilibrium
equations.” Lecture Notes in TKT4197 Nonlinear Finite Element Analysis,
NTNU (Norwegian University of Science and Technology), 2013.
[58] M. Kristoffersen, “Fluid-structure interaction simulations of pipeline impact
using europlexus.” NTNU (Norwegian University of Science and Technology),
2013.
[59] F. Casadei, G. Solomos, and G. Valsamos, “Numerical Simulation of Fast Tran-
sient Dynamic Phenomena in Fluid-Structure Systems,” tech. rep., European
Commission, Joint Research Centre, Institute for the Protection and Security
of the Citizen, 2013.
[60] Wikipedia, “Fluid-structure interaction.” URL http://en.wikipedia.org/
wiki/Fluid%E2%80%93structure_interaction. Accessed: April 24, 2014.
[61] F. Casadei and J.-P. Halleux, “Binary spatial partitioning of the central-
difference time integration scheme for explicit fast transient dynamics,” Inter-
national Journal for Numerical Methods in Engineering, vol. 78(12), pp. 1436–
1473, 2009.
[62] E. Fagerholt, “Field Measurements in Mechanical Testing Using Close-Range
Photogrammetry and Digital Image Analysis,” Doctoral thesis, NTNU (Nor-
wegian University of Science and Technology), 2012.
[63] K. Spranghers, I. Vasilakos, D. Lecompte, H. Sol, and J. Vantomme, “Full-
Field Deformation Measurements of Aluminum Plates Under Free Air Blast
Loading,” Experimental Mechanics, vol. 52, pp. 1371–1384, 2012.
177
CHAPTER 10. BIBLIOGRAPHY
[64] N. Jones, Structural Impact. 32 Avenue of the Americas, New York, NY 10013-
2473, USA: Cambridge University Press, 2nd ed., 2012.
[65] Aalco Metals Limited, “Aluminium Alloy - Commercial Alloy - 1050A
H14 Sheet.” URL http://www.aalco.co.uk/datasheets/Aluminium-Alloy_
1050-H14_57.ashx. Accessed: April 26, 2014.
[66] Kistler, “Type 603B - Piezoelectric Pressure Sensor (Miniature).” URL http://
www.kistler.com/ca/en/product/pressure/603B. Accessed: June 9, 2014.
[67] Phantom Camera Products, “Phantom v1610.” URL http://www.
visionresearch.com/Products/High-Speed-Cameras/v1610/. Accessed:
June 9, 2014.
[68] Livermore Software Technology Corporation, “LS-OPT - Graphical Optimaza-
tion Tool.” URL http://www.lstc.com/products/ls-opt. Accessed: June 9,
2014.
[69] Open Cascade EDF CEA, “SALOME - The Open Source Integration Plat-
form for Numerical Simulations.” URL http://www.salome-platform.org/.
Accessed: June 8, 2014.
[70] Commissariat à l’énergie atomique et aux énergies alternatives (France),
“Cast3M.” URL http://www-cast3m.cea.fr/index.php. Accessed: June 8,
2014.
178
Appendix A
Calibration Tests
14000
Pressure transducer 1
Pressure transducer 2
12000 Pressure transducer 3
Pressure transducer 4
10000
Pressure [kPa]
8000
6000
4000
2000
A1
APPENDIX A. CALIBRATION TESTS
14000
Pressure transducer 1
Pressure transducer 8
12000 Pressure transducer 9
Pressure transducer 10
10000
Pressure [kPa]
8000
6000
4000
2000
10000
Pressure [kPa]
8000
6000
4000
2000
A2
APPENDIX A. CALIBRATION TESTS
4
x 10
2.2
Pressure transducer 1
2 Pressure transducer 2
Pressure transducer 3
1.8
Pressure transducer 4
1.6
Pressure [kPa]
1.4
1.2
0.8
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
Time [ms]
(a) Pressure transducers along the vertical.
4
x 10
2.2
Pressure transducer 1
2 Pressure transducer 8
Pressure transducer 9
1.8
Pressure transducer 10
1.6
Pressure [kPa]
1.4
1.2
0.8
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
Time [ms]
(b) Pressure transducers along the diagonal.
A3
APPENDIX A. CALIBRATION TESTS
4
x 10
2.2
Pressure transducer 1
2 Pressure transducer 5
Pressure transducer 6
1.8
Pressure transducer 7
1.6
Pressure [kPa]
1.4
1.2
0.8
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
Time [ms]
(c) Pressure transducers in the middle and in the clamping frame.
A4
APPENDIX A. CALIBRATION TESTS
12000
Pressure transducer 1
Pressure transducer 2
Pressure transducer 3
10000 Pressure transducer 4
8000
Pressure [kPa]
6000
4000
2000
0
0 0.2 0.4 0.6 0.8 1
Time [ms]
(a) Pressure transducers along the vertical.
12000
Pressure transducer 1
Pressure transducer 8
Pressure transducer 9
10000 Pressure transducer 10
8000
Pressure [kPa]
6000
4000
2000
0
0 0.2 0.4 0.6 0.8 1
Time [ms]
(b) Pressure transducers along the diagonal.
A5
APPENDIX A. CALIBRATION TESTS
12000
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
10000 Pressure transducer 7
8000
Pressure [kPa]
6000
4000
2000
0
0 0.2 0.4 0.6 0.8 1
Time [ms]
(c) Pressure transducers in the middle and in the clamping frame.
A6
APPENDIX A. CALIBRATION TESTS
5000
Pressure transducer 1
4500 Pressure transducer 2
Pressure transducer 3
4000
Pressure transducer 4
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
(a) Pressure transducers along the vertical.
5000
Pressure transducer 1
4500 Pressure transducer 8
Pressure transducer 9
4000
Pressure transducer 10
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
(b) Pressure transducers along the diagonal.
A7
APPENDIX A. CALIBRATION TESTS
5000
Pressure transducer 1
4500 Pressure transducer 5
Pressure transducer 6
4000
Pressure transducer 7
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
(c) Pressure transducers in the middle and in the clamping frame.
A8
APPENDIX A. CALIBRATION TESTS
6000
Pressure transducer 1
Pressure transducer 2
Pressure transducer 3
5000 Pressure transducer 4
4000
Pressure [kPa]
3000
2000
1000
4000
Pressure [kPa]
3000
2000
1000
A9
APPENDIX A. CALIBRATION TESTS
6000
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
5000 Pressure transducer 7
4000
Pressure [kPa]
3000
2000
1000
A10
APPENDIX A. CALIBRATION TESTS
5000
Pressure transducer 1
4500 Pressure transducer 2
Pressure transducer 3
4000
Pressure transducer 4
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
(a) Pressure transducers along the vertical.
5000
Pressure transducer 1
4500 Pressure transducer 8
Pressure transducer 9
4000
Pressure transducer 10
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
(b) Pressure transducers along the diagonal.
A11
APPENDIX A. CALIBRATION TESTS
5000
Pressure transducer 1
4500 Pressure transducer 5
Pressure transducer 6
4000
Pressure transducer 7
3500
Pressure [kPa]
3000
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
(c) Pressure transducers in the middle and in the clamping frame.
A12
APPENDIX A. CALIBRATION TESTS
3000
Pressure transducer 1
Pressure transducer 2
2500
Pressure transducer 3
Pressure transducer 4
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(a) Pressure transducers along the vertical.
3000
Pressure transducer 1
Pressure transducer 5
2500
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(b) Pressure transducers in the middle and in the clamping frame.
A13
APPENDIX A. CALIBRATION TESTS
2000
Pressure transducer 1
Pressure transducer 2
Pressure transducer 3
Pressure transducer 4
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(a) Pressure transducers along the vertical.
2000
Pressure transducer 1
Pressure transducer 5
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(b) Pressure transducers in the middle and in the clamping frame.
A14
APPENDIX A. CALIBRATION TESTS
3000
Pressure transducer 1
Pressure transducer 2
2500
Pressure transducer 3
Pressure transducer 4
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(a) Pressure transducers along the vertical.
3000
Pressure transducer 1
Pressure transducer 5
2500
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(b) Pressure transducers in the middle and in the clamping frame.
A15
APPENDIX A. CALIBRATION TESTS
2000
Pressure transducer 1
Pressure transducer 2
Pressure transducer 3
Pressure transducer 4
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(a) Pressure transducers along the vertical.
2000
Pressure transducer 1
Pressure transducer 8
Pressure transducer 10
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(b) Pressure transducers along the diagonal.
A16
APPENDIX A. CALIBRATION TESTS
2000
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(c) Pressure transducers in the middle and in the clamping frame.
A17
APPENDIX A. CALIBRATION TESTS
3000
Pressure transducer 1
Pressure transducer 2
2500
Pressure transducer 3
Pressure transducer 4
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(a) Pressure transducers along the vertical.
3000
Pressure transducer 1
Pressure transducer 8
2500
Pressure transducer 10
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(b) Pressure transducers along the diagonal.
A18
APPENDIX A. CALIBRATION TESTS
3000
Pressure transducer 1
Pressure transducer 5
2500
Pressure transducer 6
Pressure transducer 7
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(c) Pressure transducers in the middle and in the clamping frame.
A19
APPENDIX A. CALIBRATION TESTS
3000
Pressure transducer 1
Pressure transducer 2
2500
Pressure transducer 3
Pressure transducer 4
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(a) Pressure transducers along the vertical.
3000
Pressure transducer 1
Pressure transducer 8
2500
Pressure transducer 9
Pressure transducer 10
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(b) Pressure transducers along the diagonal.
A20
APPENDIX A. CALIBRATION TESTS
3000
Pressure transducer 1
Pressure transducer 5
2500
Pressure transducer 6
Pressure transducer 7
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
(c) Pressure transducers in the middle and in the clamping frame.
A21
APPENDIX A. CALIBRATION TESTS
14000 14000
12000 12000
10000 10000
Pressure [kPa]
Pressure [kPa]
8000 8000
6000 6000
4000 4000
2000 2000
0 0
12000 12000
10000 10000
Pressure [kPa]
Pressure [kPa]
8000 8000
6000 6000
4000 4000
2000 2000
0 0
12000 12000
10000 10000
Pressure [kPa]
Pressure [kPa]
8000 8000
6000 6000
4000 4000
2000 2000
0 0
A22
APPENDIX A. CALIBRATION TESTS
14000 14000
12000 12000
10000 10000
Pressure [kPa]
Pressure [kPa]
8000 8000
6000 6000
4000 4000
2000 2000
0 0
12000 12000
10000 10000
Pressure [kPa]
Pressure [kPa]
8000 8000
6000 6000
4000 4000
2000 2000
0 0
A23
APPENDIX A. CALIBRATION TESTS
4 4
x 10 x 10
2.2 2.2
2 2
1.8 1.8
1.6 1.6
Pressure [kPa]
Pressure [kPa]
1.4 1.4
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time [ms] Time [ms]
2 2
1.8 1.8
1.6 1.6
Pressure [kPa]
Pressure [kPa]
1.4 1.4
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time [ms] Time [ms]
2 2
1.8 1.8
1.6 1.6
Pressure [kPa]
Pressure [kPa]
1.4 1.4
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time [ms] Time [ms]
A24
APPENDIX A. CALIBRATION TESTS
4 4
x 10 x 10
2.2 2.2
2 2
1.8 1.8
1.6 1.6
Pressure [kPa]
Pressure [kPa]
1.4 1.4
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time [ms] Time [ms]
2 2
1.8 1.8
1.6 1.6
Pressure [kPa]
Pressure [kPa]
1.4 1.4
1.2 1.2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time [ms] Time [ms]
A25
APPENDIX A. CALIBRATION TESTS
12000 12000
10000 10000
8000 8000
Pressure [kPa]
Pressure [kPa]
6000 6000
4000 4000
2000 2000
0 0
10000 10000
8000 8000
Pressure [kPa]
Pressure [kPa]
6000 6000
4000 4000
2000 2000
0 0
10000 10000
8000 8000
Pressure [kPa]
Pressure [kPa]
6000 6000
4000 4000
2000 2000
0 0
A26
APPENDIX A. CALIBRATION TESTS
12000 12000
10000 10000
8000 8000
Pressure [kPa]
Pressure [kPa]
6000 6000
4000 4000
2000 2000
0 0
10000 10000
8000 8000
Pressure [kPa]
Pressure [kPa]
6000 6000
4000 4000
2000 2000
0 0
A27
APPENDIX A. CALIBRATION TESTS
5000 5000
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
A28
APPENDIX A. CALIBRATION TESTS
5000 5000
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
A29
APPENDIX A. CALIBRATION TESTS
6000 6000
5000 5000
4000 4000
Pressure [kPa]
Pressure [kPa]
3000 3000
2000 2000
1000 1000
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
5000 5000
4000 4000
Pressure [kPa]
Pressure [kPa]
3000 3000
2000 2000
1000 1000
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
5000 5000
4000 4000
Pressure [kPa]
Pressure [kPa]
3000 3000
2000 2000
1000 1000
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
A30
APPENDIX A. CALIBRATION TESTS
6000 6000
5000 5000
4000 4000
Pressure [kPa]
Pressure [kPa]
3000 3000
2000 2000
1000 1000
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
5000 5000
4000 4000
Pressure [kPa]
Pressure [kPa]
3000 3000
2000 2000
1000 1000
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
A31
APPENDIX A. CALIBRATION TESTS
5000 5000
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
A32
APPENDIX A. CALIBRATION TESTS
5000 5000
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
4500 4500
4000 4000
3500 3500
Pressure [kPa]
Pressure [kPa]
3000 3000
2500 2500
2000 2000
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms] Time [ms]
A33
APPENDIX A. CALIBRATION TESTS
3000 3000
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
A34
APPENDIX A. CALIBRATION TESTS
2000 2000
1500 1500
Pressure [kPa]
Pressure [kPa]
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
1500 1500
Pressure [kPa]
Pressure [kPa]
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
A35
APPENDIX A. CALIBRATION TESTS
3000 3000
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
A36
APPENDIX A. CALIBRATION TESTS
2000 2000
1500 1500
Pressure [kPa]
Pressure [kPa]
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
1500 1500
Pressure [kPa]
Pressure [kPa]
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
1500 1500
Pressure [kPa]
Pressure [kPa]
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
A37
APPENDIX A. CALIBRATION TESTS
2000 2000
1500 1500
Pressure [kPa]
Pressure [kPa]
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
A38
APPENDIX A. CALIBRATION TESTS
3000 3000
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
A39
APPENDIX A. CALIBRATION TESTS
3000 3000
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500
2000
Pressure [kPa]
1500
1000
500
−500
0 0.5 1 1.5 2 2.5 3
Time [ms]
A40
APPENDIX A. CALIBRATION TESTS
3000 3000
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
A41
APPENDIX A. CALIBRATION TESTS
3000 3000
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
2500 2500
2000 2000
Pressure [kPa]
Pressure [kPa]
1500 1500
1000 1000
500 500
0 0
−500 −500
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Time [ms] Time [ms]
A42
APPENDIX A. CALIBRATION TESTS
A43
APPENDIX A. CALIBRATION TESTS
A44
APPENDIX A. CALIBRATION TESTS
A45
APPENDIX A. CALIBRATION TESTS
A46
Appendix B
Component Tests
B1
APPENDIX B. COMPONENT TESTS
B2
APPENDIX B. COMPONENT TESTS
(i) A33
B3
APPENDIX B. COMPONENT TESTS
4000
Pressure transducer 1
Pressure transducer 5
3500 Pressure transducer 6
Pressure transducer 7
3000
2500
Pressure [kPa]
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
Figure B.2: Pressure transducers in the clamping frame for test A11.
4000
Pressure transducer 1
Pressure transducer 5
3500 Pressure transducer 6
Pressure transducer 7
3000
2500
Pressure [kPa]
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
Figure B.3: Pressure transducers in the clamping frame for test A12.
B4
APPENDIX B. COMPONENT TESTS
4500
Pressure transducer 1
4000
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
3500
3000
Pressure [kPa]
2500
2000
1500
1000
500
−500
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time [ms]
Figure B.4: Pressure transducers in the clamping frame for test A13.
2000
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5
Time [ms]
Figure B.5: Pressure transducers in the clamping frame for test A21.
B5
APPENDIX B. COMPONENT TESTS
2000
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5
Time [ms]
Figure B.6: Pressure transducers in the clamping frame for test A22.
2000
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
1500
Pressure [kPa]
1000
500
−500
0 0.5 1 1.5 2 2.5
Time [ms]
Figure B.7: Pressure transducers in the clamping frame for test A23.
B6
APPENDIX B. COMPONENT TESTS
1500
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
1000
Pressure [kPa]
500
−500
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Time [ms]
Figure B.8: Pressure transducers in the clamping frame for test A31.
1500
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
1000
Pressure [kPa]
500
−500
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Time [ms]
Figure B.9: Pressure transducers in the clamping frame for test A32.
B7
APPENDIX B. COMPONENT TESTS
1500
Pressure transducer 1
Pressure transducer 5
Pressure transducer 6
Pressure transducer 7
1000
Pressure [kPa]
500
−500
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Time [ms]
Figure B.10: Pressure transducers in the clamping frame for test A33.
B8
APPENDIX B. COMPONENT TESTS
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(a) t = 0.39 ms (b) t = 0.39 ms (c) t = 0.39 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(d) t = 0.58 ms (e) t = 0.58 ms (f) t = 0.58 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(g) t = 0.77 ms (h) t = 0.77 ms (i) t = 0.77 ms
B9
APPENDIX B. COMPONENT TESTS
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(j) t = 0.96 ms (k) t = 0.96 ms (l) t = 0.96 ms
B10
APPENDIX B. COMPONENT TESTS
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(a) t = 0.39 ms (b) t = 0.39 ms (c) t = 0.39 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(d) t = 0.58 ms (e) t = 0.58 ms (f) t = 0.58 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(g) t = 0.77 ms (h) t = 0.77 ms (i) t = 0.77 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(j) t = 0.96 ms (k) t = 0.96 ms (l) t = 0.96 ms
B11
APPENDIX B. COMPONENT TESTS
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(a) t = 0.39 ms (b) t = 0.39 ms (c) t = 0.39 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(d) t = 0.58 ms (e) t = 0.58 ms (f) t = 0.58 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(g) t = 0.77 ms (h) t = 0.77 ms (i) t = 0.77 ms
50 5000 50
Displacement [mm]
Displacement [mm]
Pressure [kPa]
40 4000 40
30 3000 30
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 0 100 200 300
Time [ms] Cross section [mm]
(j) t = 0.96 ms (k) t = 0.96 ms (l) t = 0.96 ms
B12
APPENDIX B. COMPONENT TESTS
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(a) t = 0.59 ms (b) t = 0.59 ms (c) t = 0.59 ms
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(d) t = 0.78 ms (e) t = 0.78 ms (f) t = 0.78 ms
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(g) t = 0.97 ms (h) t = 0.97 ms (i) t = 0.97 ms
B13
APPENDIX B. COMPONENT TESTS
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(j) t = 1.16 ms (k) t = 1.16 ms (l) t = 1.16 ms
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000 20
15 1500
10
10 1000
5 500 0
0 0
−5 −500 −10
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(m) t = 4.92 ms (n) t = 4.92 ms (o) t = 4.92 ms
B14
APPENDIX B. COMPONENT TESTS
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000
20
15 1500
10 1000
5 500 10
0 0
−5 −500 0
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(a) t = 0.54 ms (b) t = 0.54 ms (c) t = 0.54 ms
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000
20
15 1500
10 1000
5 500 10
0 0
−5 −500 0
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(d) t = 0.73 ms (e) t = 0.73 ms (f) t = 0.73 ms
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000
20
15 1500
10 1000
5 500 10
0 0
−5 −500 0
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(g) t = 0.92 ms (h) t = 0.92 ms (i) t = 0.92 ms
B15
APPENDIX B. COMPONENT TESTS
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000
20
15 1500
10 1000
5 500 10
0 0
−5 −500 0
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(j) t = 1.11 ms (k) t = 1.11 ms (l) t = 1.11 ms
35 3500
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
25 2500
20 2000
20
15 1500
10 1000
5 500 10
0 0
−5 −500 0
0 1 2 3 4 5 6 7 8 9 10 0 100 200 300
Time [ms] Cross section [mm]
(m) t = 4.40 ms (n) t = 4.40 ms (o) t = 4.40 ms
B16
APPENDIX B. COMPONENT TESTS
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 2.5 0 100 200 300
Time [ms] Cross section [mm]
(a) t = 0.54 ms (b) t = 0.54 ms (c) t = 0.54 ms
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 2.5 0 100 200 300
Time [ms] Cross section [mm]
(d) t = 0.73 ms (e) t = 0.73 ms Displacement [mm] (f) t = 0.73 ms
Displacement [mm]
30 3000 30
Pressure [kPa]
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 2.5 0 100 200 300
Time [ms] Cross section [mm]
(g) t = 0.92 ms (h) t = 0.92 ms (i) t = 0.92 ms
Displacement [mm]
Displacement [mm]
30 3000 30
Pressure [kPa]
20 2000
20
10 1000
10
0 0
0
0 0.5 1 1.5 2 2.5 0 100 200 300
Time [ms] Cross section [mm]
(j) t = 1.11 ms (k) t = 1.11 ms (l) t = 1.11 ms
B17
APPENDIX B. COMPONENT TESTS
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
B18
APPENDIX B. COMPONENT TESTS
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
B19
APPENDIX B. COMPONENT TESTS
30
Displacement [mm]
40 1000
Displacement [mm]
Pressure [kPa]
30 750
20
20 500
10 250
10
0 0
−10 −250 0
0 0.5 1 1.5 2 2.5 3 0 100 200 300
Time [ms] Cross section [mm]
(a) t = 0.74 ms (b) t = 0.74 ms (c) t = 0.74 ms
30
Displacement [mm]
40 1000
Displacement [mm]
Pressure [kPa]
30 750
20
20 500
10 250
10
0 0
−10 −250 0
0 0.5 1 1.5 2 2.5 3 0 100 200 300
Time [ms] Cross section [mm]
(d) t = 0.93 ms (e) t = 0.93 ms (f) t = 0.93 ms
30
Displacement [mm]
40 1000
Displacement [mm]
Pressure [kPa]
30 750
20
20 500
10 250
10
0 0
−10 −250 0
0 0.5 1 1.5 2 2.5 3 0 100 200 300
Time [ms] Cross section [mm]
(g) t = 1.12 ms (h) t = 1.12 ms (i) t = 1.12 ms
30
Displacement [mm]
40 1000
Displacement [mm]
Pressure [kPa]
30 750
20
20 500
10 250
10
0 0
−10 −250 0
0 0.5 1 1.5 2 2.5 3 0 100 200 300
Time [ms] Cross section [mm]
(j) t = 1.31 ms (k) t = 1.31 ms (l) t = 1.31 ms
B20
APPENDIX B. COMPONENT TESTS
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
B21
APPENDIX B. COMPONENT TESTS
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
40 1000
Displacement [mm]
Displacement [mm]
30
Pressure [kPa]
30 750
20
20 500
10
10 250
0 0 0
B22
Appendix C
MATLAB-files
clear all
close all
clc
testdata=importdata('R11.csv');
t=testdata.data(:,1);
p=testdata.data(:,2:11)*10^2; % Converting from bar to kPa
bw=testdata.data(:,12);
max_pressure=max(p);
t=t*10^3;
C1
APPENDIX C. MATLAB-FILES
top=bw(1);
count=1;
while top<=0.5
count=count+1;
top=bw(count);
end
figure(11);
plot(t,bw,'linewidth',1.6);
hold on;
grid on;
tbw=t(count);
bwline=[tbw tbw];
line=[-1000 1000]*10^2;
plot(bwline,line,'r','linewidth',1.6);
axis([0 20 -1 6]);
xlabel('Time [ms]','fontsize',26,'interpreter','latex');
ylabel('Response','fontsize',26,'interpreter','latex');
set(gca,'fontsize',20);
pause(0.01);
jFrame=get(handle(gcf),'JavaFrame');
set(jFrame,'maximized',true);
set(gcf,'PaperPositionMode','auto');
pause(0.001);
export_fig R11-bw.eps -eps -transparent
t=t(count:end);
t=t-t(1);
p=p(count:end,:);
ta=zeros(1,10);
t1=zeros(1,10);
t2=zeros(1,10);
impulse_pos=zeros(1,10);
impulse_neg=zeros(1,10);
for i=1:10
% Find time of arrival
top=p(1,i);
count=1;
while top<=300
count=count+1;
top=p(count,i);
end
C2
APPENDIX C. MATLAB-FILES
ta(i)=t(count);
tacount=count;
while t(count)<=t1pick(i)
count=count+1;
end
t1count=count;
t1(i)=t(t1count);
while t(count)<=t2pick(i)
count=count+1;
end
t2count=count;
t2(i)=t(t2count);
tpos=t(tacount:t1count);
ppos=p(tacount:t1count,i);
tneg=t(t1count:t2count);
pneg=p(t1count:t2count,i);
impulse_pos(i)=trapz(tpos,ppos);
impulse_neg(i)=trapz(tneg,pneg);
taline=[ta(i) ta(i)];
t1line=[t1(i) t1(i)];
t2line=[t2(i) t2(i)];
figure(i);
plot(t,p(:,i),'linewidth',1.6);
hold on;
grid on;
plot(taline,line,'r','linewidth',1.6);
plot(t1line,line,'r','linewidth',1.6);
plot(t2line,line,'r','linewidth',1.6);
axis([0 1.2 -10*10^2 140*10^2]);
xlabel('Time [ms]','fontsize',26,'interpreter','latex');
ylabel('Pressure [kPa]','fontsize',26,'interpreter','latex');
set(gca,'fontsize',20);
pause(0.01);
jFrame=get(handle(gcf),'JavaFrame');
set(jFrame,'maximized',true);
set(gcf,'PaperPositionMode','auto');
end
C3
APPENDIX C. MATLAB-FILES
tplus=t1-ta;
tminus=t2-t1;
figure(12);
plot(t,p(:,1),t,p(:,2),t,p(:,3),t,p(:,4),'linewidth',1.6);
grid on;
axis([0 1.2 -10*10^2 140*10^2]);
xlabel('Time [ms]','fontsize',26,'interpreter','latex');
ylabel('Pressure [kPa]','fontsize',26,'interpreter','latex');
set(gca,'fontsize',20);
leg=legend('Pressure transducer 1','Pressure transducer 2', ...
'Pressure transducer 3','Pressure transducer 4');
set(leg,'fontsize',26,'interpreter','latex');
pause(0.01);
jFrame=get(handle(gcf),'JavaFrame');
set(jFrame,'maximized',true);
set(gcf,'PaperPositionMode','auto');
pause(0.001);
export_fig R11-1-4.eps -eps -transparent
figure(13);
plot(t,p(:,1),t,p(:,8),t,p(:,9),t,p(:,10),'linewidth',1.6);
grid on;
axis([0 1.2 -10*10^2 140*10^2]);
xlabel('Time [ms]','fontsize',26,'interpreter','latex');
ylabel('Pressure [kPa]','fontsize',26,'interpreter','latex');
set(gca,'fontsize',20);
leg=legend('Pressure transducer 1','Pressure transducer 8', ...
'Pressure transducer 9','Pressure transducer 10');
set(leg,'fontsize',26,'interpreter','latex');
pause(0.01);
jFrame=get(handle(gcf),'JavaFrame');
set(jFrame,'maximized',true);
set(gcf,'PaperPositionMode','auto');
pause(0.001);
export_fig R11-1+8-10.eps -eps -transparent
figure(14);
plot(t,p(:,1),t,p(:,5),t,p(:,6),t,p(:,7),'linewidth',1.6);
grid on;
axis([0 1.2 -10*10^2 140*10^2]);
xlabel('Time [ms]','fontsize',26,'interpreter','latex');
ylabel('Pressure [kPa]','fontsize',26,'interpreter','latex');
set(gca,'fontsize',20);
leg=legend('Pressure transducer 1','Pressure transducer 5', ...
'Pressure transducer 6','Pressure transducer 7');
set(leg,'fontsize',26,'interpreter','latex');
C4
APPENDIX C. MATLAB-FILES
pause(0.1);
jFrame=get(handle(gcf),'JavaFrame');
set(jFrame,'maximized',true);
set(gcf,'PaperPositionMode','auto');
pause(0.01);
export_fig R11-1+5-7.eps -eps -transparent
figure(1);
pause(0.001);
export_fig R11-p1.eps -eps -transparent
figure(2);
pause(0.001);
export_fig R11-p2.eps -eps -transparent
figure(3);
pause(0.001);
export_fig R11-p3.eps -eps -transparent
figure(4);
pause(0.001);
export_fig R11-p4.eps -eps -transparent
figure(5);
pause(0.001);
export_fig R11-p5.eps -eps -transparent
figure(6);
pause(0.001);
export_fig R11-p6.eps -eps -transparent
figure(7);
pause(0.001);
export_fig R11-p7.eps -eps -transparent
figure(8);
pause(0.001);
export_fig R11-p8.eps -eps -transparent
figure(9);
pause(0.001);
export_fig R11-p9.eps -eps -transparent
figure(10);
pause(0.001);
export_fig R11-p10.eps -eps -transparent
C5
APPENDIX C. MATLAB-FILES
clear all
close all
clc
testdata=importdata('A11.csv');
figure(1);
plot(t,p);
grid on;
xlabel('Time [ms]');
ylabel('Pressure [kPa]');
position=defodata.data(:,1:2:2*nframes-1);
defoprofile=defodata.data(:,2:2:2*nframes)*-1;
for i=1:nframes
position(:,i)=position(:,i)+(300-position(end,i))*7/17;
end
for i=1:nframes
figure(framedefo(i));
plot(position(:,i),defoprofile(:,i),'k','linewidth',1.6);
grid on;
xlabel('Cross section [mm]','fontsize',45,'interpreter','latex');
ylabel('Displacement [mm]','fontsize',45,'interpreter','latex');
set(gca,'fontsize',32);
axis([0 300 0 50]);
set(gca,'Position',[0.2 .18 .6 .79]);
set(gcf,'units','normalized','outerposition',[0.2 0.1 0.6 0.8])
C6
APPENDIX C. MATLAB-FILES
set(gcf,'PaperPositionMode','auto');
end
logdata=importdata('A11 - log.txt');
frame=logdata.data(:,1);
tlog=logdata.data(:,2);
plog=logdata.data(:,3:6)*10^2;
p5log=plog(:,2);
bp=logdata.data(:,7);
middefo=logdata.data(:,8)*-1;
figure(3);
plot(tlog,plog);
grid on;
xlabel('Time [ms]');
ylabel('Pressure [kPa]');
figure(4);
plot(tlog,p5log);
hold on;
scatter(tlog(framedefo),p5log(framedefo),'r');
grid on;
xlabel('Time [ms]');
ylabel('Pressure [kPa]');
while top<=1.5*10^2
count=count+1;
top=plog(count,2);
end
talog=tlog(count-1);
talogcount=count-1;
while top<=0.5*10^2
count=count+1;
top=p(count,2);
end
C7
APPENDIX C. MATLAB-FILES
tafull=t(count);
tafullcount=count;
fullpos=zeros(1,length(framedefo));
for i=1:length(framedefo)
while t(count)-tafull<=tlog(framedefo(i))-talog
count=count+1;
end
fullpos(i)=count;
end
tfullcorr=t-t(tafullcount)+ta;
tlogcorr=tlog-tlog(talogcount)+ta;
figure(5);
plot(tfullcorr,p5);
hold on;
scatter(tfullcorr(fullpos),p5(fullpos),'r');
grid on;
xlabel('Time [ms]');
ylabel('Pressure [kPa]');
figure(6);
plot(tfullcorr,p,'linewidth',1.6);
grid on;
xlabel('Time [ms]','fontsize',26,'interpreter','latex');
ylabel('Pressure [kPa]','fontsize',26,'interpreter','latex');
leg=legend('Pressure transducer 1','Pressure transducer 5', ...
'Pressure transducer 6','Pressure transducer 7');
set(leg,'fontsize',26,'interpreter','latex');
set(gca,'fontsize',20);
axis([0 2 -500 4000]);
pause(0.01);
jFrame=get(handle(gcf),'JavaFrame');
set(jFrame,'maximized',true);
set(gcf,'PaperPositionMode','auto');
figure(7);
plot(tlogcorr,middefo,'linewidth',1.6);
hold on;
grid on;
xlabel('Time [ms]','fontsize',26','interpreter','latex');
C8
APPENDIX C. MATLAB-FILES
pause(0.01);
jFrame=get(handle(gcf),'JavaFrame');
set(jFrame,'maximized',true);
set(gcf,'PaperPositionMode','auto');
line=[-1000 1000];
for i=1:nframes-1
figure(100+framedefo(i));
[haxes,hline1,hline2]=plotyy(tlogcorr,middefo,tfullcorr,p5);
ylabel(haxes(1),'Displacement [mm]','fontsize',45, ...
'interpreter','latex');
ylabel(haxes(2),'Pressure [kPa]','fontsize',45,'interpreter','latex');
xlabel(haxes(2),'Time [ms]','fontsize',45,'interpreter','latex');
set(haxes,'fontsize',32);
set(hline1,'color','k','linewidth',1.6);
set(hline2,'color','b','linewidth',1.6);
set(haxes,{'ycolor'},{'k';'b'});
hold on;
scatter(tlogcorr(framedefo(i)),middefo(framedefo(i)),'r', ...
'linewidth',1.6);
scatter(tfullcorr(fullpos(i)),p5(fullpos(i))*10^-2,'r', ...
'linewidth',1.6);
grid(haxes(1),'on');
axis([0 2 -5 50]);
set(haxes(1),'ytick',[-10:10:50]);
set(haxes(2),'ytick',[-1000:1000:5000]);
set(haxes,'xtick',[0:0.5:2]);
ylim(haxes(2),[-500 5000]);
xlim(haxes(2),[0 2]);
tpickline=[tfullcorr(fullpos(i)) tfullcorr(fullpos(i))];
plot(tpickline,line,'r','linewidth',1.6);
set(gca,'Position',[0.2 .18 .6 .79]);
set(gcf,'units','normalized','outerposition',[0.2 0.1 0.6 0.8]);
set(gcf,'PaperPositionMode','auto');
end
figure(45);
export_fig A11-2d-45.eps -eps -transparent
figure(49);
export_fig A11-2d-49.eps -eps -transparent
C9
APPENDIX C. MATLAB-FILES
figure(53);
export_fig A11-2d-53.eps -eps -transparent
figure(57);
export_fig A11-2d-57.eps -eps -transparent
figure(145);
export_fig A11-plot-45.eps -eps -transparent
figure(149);
export_fig A11-plot-49.eps -eps -transparent
figure(153);
export_fig A11-plot-53.eps -eps -transparent
figure(157);
export_fig A11-plot-57.eps -eps -transparent
C10
APPENDIX C. MATLAB-FILES
clear all
close all
clc
name = '';
if(id < 10)
name = ['000' num2str(id)];
elseif(id < 100)
name = ['00' num2str(id)];
elseif(id < 1000)
name = ['0' num2str(id)];
elseif(id < 10000)
name = [num2str(id)];
end
dicmesh = ReadMeshBinary(fpath);
dicmesh = RigidBodyTransform(dicmesh, -0.60, -14.29, 0, 0, 0, -44.0);
corrX = 60;
corrY = 187;
if(~isempty(dicmesh.nlocX))
[elmGrid, nodeGrid] = Mesh2ElementQ4Grid(dicmesh);
ndX = dicmesh.nlocX + dicmesh.ndefX;
X = zeros(size(nodeGrid)); Y = X; Z = X; dZ = X;
for i = 1:length(nodeGrid(:,1))
for j=1:length(nodeGrid(1,:))
X(i,j) = ndX(nodeGrid(i,j), 1);
Y(i,j) = ndX(nodeGrid(i,j), 2);
Z(i,j) = ndX(nodeGrid(i,j), 3);
dZ(i,j) = dicmesh.ndefX(nodeGrid(i,j), 3);
end
end
% Plot 3D-surface
C11
APPENDIX C. MATLAB-FILES
figure(1)
surf(X,Y,-Z,-dZ);
camlight left;
lighting phong
colormap gray
material dull
shading interp
h = camlight;
set(h,'Position',[0 1200, 100])
scale = 20;
iLen = length(X(:,1))*scale;
jLen = length(X(1,:))*scale;
Xi = zeros(iLen,jLen);
Yi = zeros(iLen,jLen);
for i= 1:iLen,
for j= 1:jLen,
Xi(i,j) = j / scale;
Yi(i,j) = i / scale;
end
end
ZZ = interp2(Z,Xi,Yi);
XX = interp2(X,Xi,Yi);
YY = interp2(Y,Xi,Yi);
mini = 1;
for i = 1:length(ZZ(:,1))
if(length(find(isnan(ZZ(i,:)))) < length(ZZ(i,:)))
mini = i;
break;
end
end
minj = 1;
for i = 1:length(ZZ(1,:))
if(length(find(isnan(ZZ(:,i)))) < length(ZZ(:,1)))
minj = i;
break;
end
end
XX = XX(mini:length(XX(:,1)), minj:length(XX(1,:)));
YY = YY(mini:length(XX(:,1)), minj:length(XX(1,:)));
ZZ = ZZ(mini:length(XX(:,1)), minj:length(XX(1,:)));
C12
APPENDIX C. MATLAB-FILES
planeimg = -ZZ;
imgzposition = -60;
colorimg = mat2im(scaledimg,jet(256));
hold on;
% Plot image
surf([min_x max_x],[min_y max_y],repmat(imgzposition, [2 2]),...
colorimg,'facecolor','texture');
view(-45,30);
f3 = 22;
set(gca,'XTick',0:75:300,'YTick',0:75:300,'ZTick',[],...
'FontSize',f3);
grid off;
xlabel('x [mm]','fontsize',28,'interpreter','latex');
ylabel('y [mm]','fontsize',28,'interpreter','latex');
colormap(jet(256));
cbfreeze(cb);
colormap(gray);
set(gca,'Position',[0.1 .35 .8 .55]);
set(gcf,'units','normalized','outerposition',[0.3 0 0.4 1])
set(gcf,'PaperPositionMode','auto');
end
pause(0.1);
print -depsc A11-3d-57.eps
C13
APPENDIX C. MATLAB-FILES
C14
Appendix D
Cast3M-files
OPTI ECHO 1;
OPTI REST FORM 'tt2d.sauv';
REST FORM;
LIST;
OPTI SAUV FORM 'tt2d.msh';
LIST (nbel allele);
mesh = allele;
TASS allele NOOP;
SAUV FORM allele;
FIN;
D1
APPENDIX D. CAST3M-FILES
OPTI ECHO 1;
OPTI DIME 3 ELEM QUA4;
OPTI TITR 'PLATE EXPLOSION FSI';
OPTI TRAC PSC FTRA 'expfsi11_mesh.ps';
opti SAUV FORM 'expfsi11.msh';
* Defining plate
den = 0.025;
n1 = 12;
n2 = 25;
DENS den;
p1 = 0 -0.15 -0.15;
p2 = 0 0.15 -0.15;
p3 = 0 0.15 0.15;
p4 = 0 -0.15 0.15;
cont1 = p1 d p2 d p3 d p4 d;
* Defining frame
p5 = 0 -0.5 -0.5;
p6 = 0 0.5 -0.5;
p7 = 0 0.5 0.5;
p8 = 0 -0.5 0.5;
cont2 = p5 d p6 d p7 d p8 d;
D2
APPENDIX D. CAST3M-FILES
* Defining air
p9 = 0 -0.7 -0.5;
p10 = 0 -0.7 0.7;
p11 = 0 0.7 0.7;
p12 = 0 0.7 -0.5;
* Defining boundaries
* Top
* Back
* Front
* Left
* Right
D3
APPENDIX D. CAST3M-FILES
* Bottom
dt1 = -0.3 0 0;
dt2 = xend 0 0;
dtlb = 0.2;
* dtlb is 2 times lbsize
* dtlb is the length when extruding the large box
D4
APPENDIX D. CAST3M-FILES
nsc = 8;
* nsc is the number of elements along the extruded edge
* nsc is usually nsc = dtlb/den
sizeex = 0.01725;
* sizeex is the radius of the circle
sizeai = 0.1;
* sizeai is half the length of the sides in the cube
width = 0.4;
height = 2.0;
dini = 3.141*sizeex/(4.*nel0);
dfin = 3.141*sizeai/(4.*nel0);
* Reference
o0 = 0. 0. 0.;
x0 = (sizeex) 0. 0.;
x01 = ((-1)*sizeex) 0. 0.;
y0 = 0. (sizeex) 0.;
z0 = 0. 0. (sizeex);
z01 = 0. 0. ((-1)*sizeex);
x1 = (sizeai) 0. 0.;
x11 = ((-1)*sizeai) 0. 0.;
y1 = 0. (sizeai) 0.;
z1 = 0. 0. (sizeai);
z11 = 0. 0. ((-1)*sizeai);
symp1 = (sizeex/2.) (sizeex/2.) (sizeex/2.);
symp11 = ((-1)*sizeex/2.) (sizeex/2.) (sizeex/2.);
symp12 = ((-1)*sizeex/2.) (sizeex/2.) ((-1)*sizeex/2.);
D5
APPENDIX D. CAST3M-FILES
* Building surfaces
D6
APPENDIX D. CAST3M-FILES
D7
APPENDIX D. CAST3M-FILES
* Finishing touches
FIN;
D8
Appendix E
EUROPLEXUS-files
TENSION TEST
ECHO
!CONV win
CAST mesh
CPLA LAGR
GEOM Q42L mesh TERM
COMP EPAI 0.8e-3 LECT mesh TERM
COUL TURQ LECT mesh TERM
NGRO 1 'disp' LECT mesh TERM
COND BOX X0 0.043 Y0 -0.0001 Z0 -0.0001
DX 0.002 DY 0.0002 DZ 0.0002
FONC 1 TABL 22
0 0
0.0025 0.0000129
0.0050 0.0000514
0.0075 0.0001144
0.0100 0.0002005
0.0125 0.0003075
0.0150 0.0004328
0.0175 0.0005733
0.0200 0.0007255
0.0225 0.0008857
0.0250 0.0010500
0.0275 0.0012143
0.0300 0.0013745
0.0325 0.0015267
0.0350 0.0016672
0.0375 0.0017925
0.0400 0.0018995
0.0425 0.0019856
0.0450 0.0020486
0.0475 0.0020871
E1
APPENDIX E. EUROPLEXUS-FILES
0.0500 0.0021000
100.00 0.0021000
REGI 'force' RESU POIN LECT boundx TERM
MATE VMJC RO 2710 YOUN 7.1E10 NU 0.33 COA1 1.15E08
COA2 0.57E08 CLB1 0.014 CLB2 0.26 SRRF 1.0
LECT mesh TERM
LINK COUP
VITE 1 1.0 FONC 1 LECT load TERM
CONT SPLA NX -1 NY 0 LECT boundx TERM
CONT SPLA NX 0 NY -1 LECT boundy TERM
ECRI DEPL VITE CONT ECRO TFREQ 0.04E-1
FICH SPLI ALIC TFRE 0.04E-1
OPTI NOTEST
LOG 1
CALCUL TINI 0.0 NMAX 100000000 TEND 2.0
FIN
E2
APPENDIX E. EUROPLEXUS-FILES
E3
APPENDIX E. EUROPLEXUS-FILES
EXPFSI
ECHO
!CONV win
CAST MESH
TRID ALE
DIME NALE 1 NBLE 1 TERM
GEOM
CUVF fluid
CUVF tnt tnt2 air air2
Q4GS plate
CL3D absorb
TERM
COMP GROU 2 'half' LECT fluid TERM COND YB GT 0.0
'vcut' LECT fluid plate TERM COND YB GT 0
COUL ROUG LECT tnt tnt2 TERM
TURQ LECT fluid air air2 TERM
ROSE LECT absorb TERM
VERT LECT plate TERM
EPAI 0.8e-3 LECT plate TERM
GRIL LAGR LECT plate TERM
EULE LECT fluid TERM
MATE
GAZP RO 588.1309 GAMM 1.4 PINI 2.0467E10
PREF 1.E5 LECT tnt tnt2 TERM
GAZP RO 1 GAMM 1.4 PINI 1.E5
PREF 1.E5 LECT fluid air air2 TERM
CLVF ABSO RO 1.
LECT absorb TERM
VMJC RO 2710 YOUN 7.1E10 NU 0.33 COA1 1.15E08
COA2 0.57E08 CLB1 0.014 CLB2 0.26 SRRF 1.0
LECT plate TERM
LINK COUP
BLOQ 123456 LECT bound TERM
LINK DECO
FLSW STRU LECT plate TERM
FLUI LECT fluid TERM
R 0.0219
HGRI 0.0256
DGRI
FACE
BFLU 1
FSCP 1
ECRI VITE ECRO TFRE 10.5E-4
ELEM LECT 1 TERM
POIN LECT 1 TERM
FICH SPLI ALIC TFRE 10.5E-6
OPTI PART NOTEST
CSTA 0.125
LOG 1
VFCC FCON 6
E4
APPENDIX E. EUROPLEXUS-FILES
ORDR 2
OTPS 1
RECO 0
FLSW 1
CALCUL TINI 0 TEND 10.5E-3
FIN
E5
APPENDIX E. EUROPLEXUS-FILES
EXPFSI
ECHO
!CONV win
CAST MESH
TRID ALE
DIME NALE 1 NBLE 1352 TERM
GEOM
CUVF fluid
Q4GS plate
CL3D absorb
TERM
COMP GROU 4 'explo' LECT fluid TERM
COND BOX X0 0.35775 Y0 -0.01725 Z0 -0.01725
DX 0.0345 DY 0.0345 DZ 0.0345
'half' LECT fluid TERM COND YB GT 0.0
'vcut' LECT fluid plate TERM COND YB GT 0
'rezo' LECT fluid TERM
COND BOX X0 -0.1 Y0 -0.151 Z0 -0.151
DX 0.2 DY 0.302 DZ 0.302
COUL ROUG LECT explo TERM
TURQ LECT fluid DIFF explo TERM
ROSE LECT absorb TERM
VERT LECT plate TERM
BLEU LECT rezo TERM
EPAI 0.8e-3 LECT plate TERM
GRIL LAGR LECT plate TERM
EULE LECT eule TERM
AUTO NOEU LECT rezo DIFF eule plate TERM
MATE
GAZP RO 171.3846 GAMM 1.4 PINI 1.7380E09
PREF 1.E5 LECT explo TERM
GAZP RO 1 GAMM 1.4 PINI 1.E5
PREF 1.E5 LECT fluid DIFF explo TERM
CLVF ABSO RO 1.
LECT absorb TERM
VMJC RO 2710 YOUN 7.1E10 NU 0.33 COA1 1.15E08
COA2 0.57E08 CLB1 0.014 CLB2 0.26 SRRF 1.0
LECT plate TERM
LINK COUP
BLOQ 123456 LECT bound TERM
ECRI VITE ECRO TFRE 10.5E-4
ELEM LECT 1 TERM
POIN LECT 1 TERM
FICH SPLI ALIC TFRE 10.5E-6
OPTI NOTEST
CSTA 0.25
LOG 1
REZO LIAI MVRE MODU GAM0 0.5
VFCC FCON 6
ORDR 2
E6
APPENDIX E. EUROPLEXUS-FILES
OTPS 1
RECO 0
CALCUL TINI 0 TEND 10.5E-3
FIN
E7