Simulation of blast waves by using

mapping technology in EUROPLEXUS

Martin Larcher
Folco Casadei
George Solomos

2014
Report EUR 26735 EN
JRC91102
European Commission
Joint Research Centre
Institute for the Protection and Security of the Citizen

Contact information
Martin Larcher
Address: Joint Research Centre, Via Enrico Fermi 2749, TP 480, 21027 Ispra (VA), Italy
E-mail: [email protected]
Tel.: +39 0332 78 9563
Fax: +39 0332 78 9049

http://ipsc.jrc.ec.europa.eu/
http://www.jrc.ec.europa.eu/

This publication is a Technical Report by the Joint Research Centre of the European Commission.

Legal Notice
This publication is a Technical Report by the Joint Research Centre, the European Commission’s in-house science service.
It aims to provide evidence-based scientific support to the European policy-making process. The scientific output expressed
does not imply a policy position of the European Commission. Neither the European Commission nor any person
acting on behalf of the Commission is responsible for the use which might be made of this publication.

JRC91102
EUR 26735 EN
ISBN 978-92-79-39249-8
ISSN 1831-9424
doi: 10.2788/98310

Luxembourg: Publications Office of the European Union, 2014

© European Union, 2014

Reproduction is authorised provided the source is acknowledged.

Printed in Italy
 CONTENTS

1 Introduction........................................................................................................................................ 2
2 1D calculation .................................................................................................................................... 4
2.1 Reference model ........................................................................................................................ 4
2.2 General procedure ...................................................................................................................... 5
2.3 Writing map files ....................................................................................................................... 6
2.4 Reading the map file .................................................................................................................. 7
3 Examples and results ......................................................................................................................... 8
3.1 Mapping files ............................................................................................................................. 8
3.2 Comparison with Kingery curves ............................................................................................ 10
3.3 Alia tests data ........................................................................................................................... 18
3.3.1 1D simulations to create map file ..................................................................................... 19
3.3.2 3D simulations .................................................................................................................. 21
3.4 Simulation of an explosion inside a train ................................................................................ 23
4 How to use mapping ........................................................................................................................ 26
5 Conclusions...................................................................................................................................... 27
6 References........................................................................................................................................ 28
7 Appendix.......................................................................................................................................... 29
7.1 Creation of a map file .............................................................................................................. 29
7.2 map3d test cases ...................................................................................................................... 29
7.3 alia tests ................................................................................................................................... 30
7.3.1 Parameter test for creation of mapping file ...................................................................... 30
7.3.2 Alia simulation with small model .................................................................................... 31
7.3.3 Alia simulation with full model........................................................................................ 32
7.3.4 Explosion inside a train .................................................................................................... 33

1
1 Introduction
Finite element or finite volume simulations of the propagation of blast waves by using a model for the
explosion of the explosive solid itself require very fine meshes in the explosive material and in the
immediately zone around it. Structures may be at a long distance from the source of the explosive and
this leads often to very big meshes with many elements. On the other hand, the option of meshing
coarsely the explosive leads to results that are not accurate. There are several possibilities to deal with
this problem (Figure 1, Larcher, Casadei [10]).

Figure 1: Simulation models for blast waves

The application of a pressure-time function on the structures reduces computational costs but it is
often not applicable especially when wave reflections, shadowing etc. have to be considered. The
balloon model (Larcher, Casadei [8]) is sometimes a good compromise. It uses a compressed balloon
under high pressure, which is released and forms a pressure wave that is somehow similar to a blast
wave. Large 3D calculations with a solid TNT model using a JWL-equation can be used but they are
more effective when the results of one finer mesh are mapped onto a coarser mesh after some
calculation steps. When the blast wave reaches a certain distance to the charge, the small elements
inside the charge are not needed any more since the pressure ratio decreases strongly. These small
elements result in very small time steps for the full model and the calculation is inefficient in terms of
CPU.

2
In addition, mapping can be seen as a much more general methodology. A fine solid TNT simulation
in 1D can be mapped in 2D when the first reflection surface is reached and can be mapped in 3D when
the second reflection surface is reached. But also a mapping from 3D to 3D could have advantages in
cases with very complex geometries where results of an inner zone can be mapped to a wider zone
after reaching the end of the first model. Mapping from one material law to another should also be
possible. This could be a possibility in order to eliminate the very sophisticated and time consuming
JWL-material which is not needed after a certain time.

Another possibility is layer mapping. If fluid-structure interaction is needed in a case where the wave
is nearly plane, i.e. a structure far away from the charge, the blast could be applied on a thin layer of
one side of the fluid mesh.

It is very important that not only the pressure and the density (or two other equivalent parameters of
the ideal gas equation) are mapped into the other mesh, but also the particle velocity because this plays
an important role in a blast wave.

The commercial programs AUTODYN and LS-DYNA both offer in a certain way such mapping
technologies. These codes could be used as a reference.

This report presents first the general idea and the implementation of the mapping algorithm in
EUROPLEXUS. After that several blast examples computed using the mapping approach, are
compared with analytical or experimental results in order to show the capabilities of the method. The
results are also compared with the balloon model in order to highlight the advantages and
disadvantages of both methods. The report also shows how the mapping can be used in the code, and
terminated with some conclusions. The appendix contains all relevant input files of the examples.

3
2 1D calculation
2.1 Reference model
1D models are seldom used in EUROPLEXUS. Their capability to describe explosion phenomena is
limited, and none of the 1D elements was until now usable together with the JWL material law, which
is fundamental to describe detonation effects. All developed models are first validated using a
spherical detonation of 1kg TNT whose density is 1.654 kg/dm3. This corresponds to a spherical
charge with a radius of about 5.2 cm. The detonation wave is compared when it arrives at a distance of
1.0 m from the centre of the charge. No reflections are considered since this will be done later in the
3D calculation.

The form of a detonation wave follows in general the curve given in Figure 2.

Figure 2: Theoretical pressure-time history for an air blast detonation

This form can be represented very well for the positive part by the modified Friedlander equation. The
measured pressure at a given point can be described by the modified Friedlander equation and depends
on the time t from the arrival of the pressure wave at the point under consideration ( t  t0  ta ), where

t0 is the absolute time

 t  
bt

p(t )  p0  pmax 1   e td (1)

 td 

The other parameters involved are the atmospheric pressure p0, the maximum overpressure pmax, the
duration of the positive pressure td and the decay parameter b. These parameters can be taken from
several analytical-empirical equations like Kingery [5] and Kinney [6]. The parameter b describes the
decay of the curve. Some formulas to calculate this parameter are given in Goel [4].

4
Table 1 shows the analytical incident wave parameters obtained by the AIRB command in
EUROPLEXUS for the given geometry [7]. These parameters can be compared with the numerical
ones to validate the models.

Table 1: Pressure wave for 1 kg TNT @ 1 m according to Kingery [5] and Kinney [6]

pmax [Pa] tstart [s] td [s] pneg [Pa] tneg [s] b [-] Ipos

Kingery 0.9349E6 0.5318E-3 0.1795E-2 0.7943E5 0.1201E-1 8.416 0.1745E+03

Kinney 0.9956E6 n.a. 0.5202E-3 n.a. n.a. n.a. n.a.

2.2 General procedure

The following variables are needed to map a blast wave from one model to another:

 Pressure (scalar)

 Density (scalar)

 Particle velocity (vector)

Pressure and density are related to the element and are stored in the ECRO variables. The particle
velocity is calculated in the elements in case of cell centred finite volumes and on the nodes in case of
finite elements and node centred finite volumes. The idea is to write the three variables in an output
file that can be read during a subsequent calculation and mapped onto another mesh.

Figure 3: General mapping procedure; 1D (left) to 3D (right)

Figure 3 shows the general procedure. A 1D calculation with JWLS material law for the explosive and
the air is performed. The boundary condition is chosen in such a way that the model represents a
spherical configuration. Tube (or bar) elements with an increasing radius are used to model this

5
behaviour. When the pressure wave reaches a certain distance to the centre the mapping file is written.
This information is mapped into the 3D configuration (Figure 3, right hand side).

Figure 4: Interpolation of the blast wave parameters for a certain element / integration point

Several procedures are possible to calculate the blast wave parameters in a certain element. An
averaged value is more difficult to implement. Therefore, the value of the parameters at a certain point
(integration point) is interpolated between from the parameter-distance curves of the mapping file
(Figure 4). The same procedure is used for the particle velocity at the nodes.

In case of mapping from 2D or 3D a more sophisticated interpolation procedure must be used.

As it will be seen later on, the use of cell centred finite volumes with suitable parameters is
recommended since finite elements results in more inaccurate pressure waves in case of detonations.

2.3 Writing map files

The map data is written in a separate mapping file. This file contains in the first header line the
following information

Dimension Up to now only 1D is possible

Time Arrival time at the point where it is written

Points Number of points n for which the output is written

The following lines are repeated n times and contain the following information separated by blanks:

Distance Distance from the centre to the point where output is written

Pressure Pressure at this point in the units of the calculation

Density Density at this point in the units of the calculation

Particle velocity at this point in the units of the calculation. Three directions are given. For
1D the second and the third dimension of the particle velocity are 0.

6
The format of these files is kept as general as possible in order to allow further developments for other
cases.

The output of the file is done by defining an additional ECRI command in the input file. This ECRI
command defines the distance from the charge at which the pressure should be monitored. In addition
a critical value of the pressure is given. When this pressure value is reached at the defined point, the
map file is written. The FREQ parameter of ECRI should be defined to 1 in order to check the pressure
at each time step. The following input line could for example be used:

ECRI … FICH MAPB DIPR 1.0 PCHE 1.5e5 FREQ 1

Internally, the values of pressure, density and internal energy are calculated so that mapping can be
used both for finite elements and finite volumes. Mapping from one type of fluid formulation to
another one (finite elements, finite volumes) is possible. However, it is strongly recommended to use
finite volumes for the calculation of the map files since the results are much more accurate.

2.4 Reading the map file

The map file described in the previous section is read by using the keyword INIT MAPB. The number
of points where the information about the blast wave is stored is in most cases much bigger than the
one used in the 3D file. The map file is therefore first read and stored. Three subroutines are used in
order to compute the values for the 3D model. MAP_BLAST_ELEMENT calculates the pressure and
the density for finite elements as well as for finite volumes. MAP_BLAST_ELEMENT_V calculates
the fluid velocity for finite volumes. MAP_BLAST_NODE calculates the velocities of the fluid at the
nodes for finite elements. All these subroutines use a generic subroutine GET_MAP_VALUE that
interpolates the needed value from the map file.

7
3 Examples and results
3.1 Mapping files
In a first step the blast wave due to the detonation of a charge of 1 kg TNT is mapped into different
geometries. The map file is checked against the Kingery-data to ensure that the mapping is done with
reasonable input data from the map file.

The map file is created by using a 1D model made of finite volume elements (TUVF). This element is
used with a linearly increasing diameter in order to get a cone-type geometry. Finite Volumes have at
their borders fixed (reflecting) boundary conditions. Therefore, this model represents a cone with
symmetry conditions in all directions (spherical symmetry). The JWLS material is used for the
explosive as well as for the air with the parameters given in Table 2. This is done for simplification
since multi-material simulations are much more complicated to perform. The JWL equation can also
describe the behaviour of the air. Probably this is numerically not effective but it is more robust. Since
the calculation time of the 1D model is very small, efficiency is not important at this stage.

Table 2: Parameters for the JWL equation for TNT

Parameter Description Parameters for explosive Parameters for the air

A (Pa) 3.738e11 3.738e11
B (Pa) 3.747e9 3.747e9
R1 4.15 4.15
R2 0.90 0.90
beta 0.25 0.25
ro density 1630 1.3
eint current internal 3.68e6 2.1978E5
energy per unit mass
omeg specific heat ratio 0.35 0.35
d detonation speed 6930 -

Figure 5 shows the result of the 1D calculation (file map1d.epx in the appendix). It can be seen that the
mapping file can represent the pressure-distance curve quite well. The comparison is done in such a
way that the wave at a similar distance is compared (i.e. the time of the pressure distance curve is
taken from the Kingery data). The accordance decreases when the distance to the charge is getting
larger. The negative phase is not represented by the Kingery data. Therefore, this part cannot be
compared. Several formulas for the negative phase are given in Larcher [7]. The values are shown in

8
Table 3. The numerical results from the 1D calculation show a relatively good agreement with the data
from the literature.

Table 3: Pressure wave for 1 kg TNT @ 1 m according to Kingery [5] and Kinney [6] and Larcher [7]

pmax [Pa] tstart [s] td [s] pneg [Pa] tneg [s] b [-] Ipos

Kingery 0.9349E6 0.5318E-3 0.1795E-2 0.7943E5 0.1201E-1 8.416 0.1745E+03

Kinney 0.9956E6 n.a. 0.5202E-3 n.a. n.a. n.a. n.a.

Map file 0.897E6 0.5615E-3 0.688E-3 0.64827E5 0.189E-2 - -

As shown in Cullis [3] a shock wave is reflected also at the charge centre. The part after the peak drops
at a certain point nearly to zero. This wave travels towards the centre of the charge and results in a
large pressure that can even reach the same magnitude as the peak pressure of the first blast wave at
the same time.

This effect can be very well represented with the numerical simulation as it can be seen in Figure 5 for
the situation at a propagated distance of 1.8 m. Kingery data can also not represent this behaviour since
this behaviour is not included in the modified Friedlander equation.

It has also been checked that the angle of the cone of the 1D model that is used for this type of
simulations has no significant influence on the pressure-distance curve.

Figure 5: Comparison of the mapping file for 1 kg TNT

9
3.2 Comparison with Kingery curves
The mapping files are used in order to map the 1D pressure-distance curve into 3D cubical cells of air
and several modelling approaches are tried. It is very early found that the case with finite elements
(map3d) is not stable and terminates at a certain point by reaching a too small time step. This was
observed in several test cases where the pressure became too high in one particular element and
therefore the time step size decreased too much. Since the results were much more successful with
finite volumes, no further investigations with finite elements were done. It is generally not
recommended to use the mapping algorithm with finite elements.

Table 4 shows the models used for comparison. For all models three symmetry planes are used in order
to reduce the numerical costs, thus only 1/8 of the reality is built up in the numerical model. The
symmetry conditions must be considered in the two numerical models (finite elements, finite volumes)
in a different way. While for the finite volumes the boundary condition is reflecting, the boundary
condition for finite elements is absorbing.

The mapping file from the previous section is used where a charge of 1 kg TNT is modelled. The
element size in the 1D calculation is 1 mm, i.e. the charge consists from 30 elements. The total length
of the 1D model is 3 m. The mapping file is written when the blast wave reaches 1 m distance from the
charge. This is the case 0.571 ms after the start of the detonation.

The mapping builds up the charge in the 1D calculation with spherical symmetry. Further attention in
case of the mapping algorithm is not needed since the pressure is mapped and not the charge. For the
bubble model the total charge must be divided by 8 in order to consider the right amount of explosive.

10
 Table 4: Models for comparison of mapping algorithm with analytical data from Kingery

Model Length Element Element Number of Calculation time

type size [m] elements

Map3d 5m FE (CUBE) 0.10 125000 termination

Map3dvfcc 5m VF (CUVF) 0.10 125000 401 s (0.01 s)

Map3dvdcc_fine 5m VF (CUVF) 0.05 1000000 5509 (0.01 s)

Map3dvfcc_fine_fine 2m VF (CUVF) 0.02 1000000 3464 (0.003 s)

Map3dvf_bub 5m VF (CUVF), 0.10 125000 338 s (0.01 s)

bubble
model

Map3dvf_bub_fine 5m VF (CUVF)), 0.05 1000000 6610 s (0.01 s)

bubble
model

Map3dvf_bub_fine_fine 2 m VF (CUVF)), 0.02 1000000 5235 (0.003 s)

bubble
model

For the finite volume modelling Map3dvfcc_fine_fine Figure 6 shows the situation at the beginning of
the 3D calculation (time is reset to zero) when the 1D blast wave is mapped to the 3D fluid volume. It
can be seen that the resolution of the pressure wave is very fine. The thickness of the shock front with
the maximum pressure is small as is actually observed in reality. Inside the sphere of the high pressure
wave the pressure falls below the atmospheric pressure (underpressure). This negative phase results in
a reflection at the central point of the charge as it can be seen in Figure 7. The resolution of the
positive pressure wave at this time (t=1 ms) is still very good.

11
Figure 6: map3dvf_fine_fine immediately after the mapping (t=0 ms + 0.57 ms mapping). Blast front has
reached a distance from the charge of 1 m.

Figure 7: map3dvf_fine_fine at t=1 ms +0.57 ms mapping

The third possibility, which has been tried out, to impose in a fluid a blast pressure wave is the so-
called compressed balloon or bubble model [8]. A sphere with a radius much higher than the one of the

12
explosive if filled with compressed gas that is then released. The pressure wave resulting from that
pressure release can describe a blast wave quite well especially when the pressure of the compressed
air is fitted by comparing the resulting pressure-time curves with the ones of experimental data.

Figure 8 shows the bubble model after filling the bubble with the compressed air. It is obvious that the
main idea is similar to the mapping algorithm but the bubble model uses a constant pressure instead of
the physical pressure-distance relation that is derived by a numerical 1D pre-calculation.

The comparison between Figure 7 and Figure 9 shows the advantage of the mapping algorithm. The
pressure wave at the time t=1 ms is much more pronounced for the mapping algorithm than for the
bubble model. This corresponds much better with the physical pressure-time curve as it can be seen in
Figure 10 and Figure 11.

The triggering of the two models is different. The arrival time of the blast wave in case of mapping
algorithm (3D model) can be determined by adding the time that the wave needs in the 1D model. In
case of the balloon model the arrival time can be estimated by using the arrival time at the radius of the
balloon. This time can be calculated using Kingery formulas (AIRB).

Figure 8: map3dvf_bub_fine_fine bubble filled with compressed gas (t=0 s + 0.1 ms bubble)

13
 Figure 9: map3dvf_bub_fine_fine at t=1 ms + 0.1 ms bubble

In Figure 12 and Figure 13 the maximum pressures of all models are compared against the empirical
values from Kingery. As expected the numerical simulations underestimate the pressure waves since
numerical simulations are a kind of filter or damping for high frequencies. This means that the blast
wave is smoothed in each step of the calculation. The results of the bubble model are slightly better in
comparison to the Kingery values. This is also clear since the compressed pressure in the balloon has
initially been calculated in such a way that the pressure wave can best represent the empirical data. For
the mapping the blast wave is filtered much longer and therefore the pressure wave is slightly smaller
than the bubble one. Nevertheless, the mapping algorithm can represent also the maximum pressure
over distance function for that charge size. The finer the mesh is (for bubble or mapping model), the
better the results are. This is also evident since the high frequency filtering is then smaller.

The pressure-time curves shown in Figure 10 and Figure 11 indicate also some slight differences
between the two models. The peak pressure is slightly smaller for the mapping algorithm due to the
fact that the bubble model uses a correction factor to get better results. This correction is obtained by
fitting the numerical results to Kingery’s empirical data. Kingery’s empirical data are very smooth
since only an exponential equation is used to describe the decay of the pressure (modified Friedlander
equation (1)). In reality, the physical process is much more complicated e.g. includes a negative phase
and reflections of the negative wave at the centre. This process cannot be described neither with
Kingery’s formulation nor with the bubble model. A comparison with experimental values would most
probably be more similar to the mapping results, especially the pressure history after reaching the
peak.

14
The curves in Figure 10 and Figure 11 are moved in such a way that the blast wave arrives always at
t = 0 ms to allow a comparison between all models. The numerical arrival times can also be compared
here more in detail. For the bubble model, the arrival time at the maximum radius of the bubble is
about 0.09 ms (from AIRB). The time for the wave to reach the distance of 1.2 m is about 0.75 ms
resulting in a total arrival time of 0.84 ms. For the mapping model, the time of the wave to reach the
map point (1 m distance from the charge) is about 0.57 ms. The time of the wave in the 3D model to
reach the 1.2 m distance is 0.2 ms corresponding in total to an arrival time of 0.77 ms. The empirical
arrival time that can be calculated for a distance of 1.2 m with AIRB is about 0.75 ms and corresponds
quite well with the one of the mapping model but slightly worse with the one of the balloon model.

Figure 10: Pressure history for the comparison between mapping and bubble model at a distance of
1.2 m. The starting time is set to the arrival of the wave in order to make it comparable.

15
Figure 11: Pressure history for the comparison between mapping and bubble model at a distance of
1.2 m, detail

Figure 12: Comparison of the pressure-distance curve for a charge of 1 kg TNT equivalent

16
 Figure 13: Comparison of the detailed pressure-distance curve

Also the impulse-distance curve is an indicator of how well a numerical model can represent the
development of a blast wave. In Figure 14 the impulse-distance curves of the mapping method are
compared with the bubble method and with the empirical values also from Kingery. It can be seen that
the impulse is better represented by the bubble method. The reason is that the pressure wave is
represented in a much sharper way in case of the mapping algorithm which together with the lower
peak pressure produce a smaller impulse value. In case of the bubble model the numerical impulse is
fit to the Kingery curves by using a correction parameter. Therefore, the bubble model can represent
the impulse better, although the correction factor does not have any physical background.

Figure 14: Comparison of the 1 kg TNT: impulse-distance curve

17
The calculation time needed is (except for the very coarse mesh) slightly smaller for the mapping
technology. This was not foreseeable since the bubble model has at least at the beginning a more
uniform pressure field than the mapped model. Nevertheless, this fact could be an additional advantage
of the mapping technology if this difference can be confirmed by other comparative calculations.

3.3 Alia tests data

The next comparisons are done with data from a real experiment. As in a previous report [8], the
experiment of Alia [1] has been used since it provides blast data that show also a reflection on a
surface (Figure 15).

The experiment uses a charge of 227 g C4, which corresponds to a charge of 268 g TNT equivalent.
Three symmetry planes can be considered in order to reduce the computational costs, i.e. 1/8 of the
reality is modelled.

The previous report shows results from solid TNT and bubble models testing several parameters. The
comparison should be completed here by adding the mapping model results.

Figure 15: Experimental setup of the experiment of Alia [1]

18
3.3.1 1D simulations to create map file
In the first step, the mapping file is created. The first simulations of that problem get about 15 % lower
peak pressures in comparison to Kingery’s empirical values. To understand this difference several
parameters are changed in order to identify substantial influences on the pressure values.

Figure 16 shows the shape and dimensions of the 1D model. To get a spherical symmetry a cone is
built using 1D elements. This means that the cross section of the elements is increased linearly with the
distance to the centre. Since a cross section of 0.0 for the inner element is not possible in the recent
EUROPLEXUS code, a small enough value should be used, which however is not too small to cause
the premature termination of the calculations. In addition to the element size the minimum and the
maximum cross sections are changed for the parameter study. The parameters for the different runs are
given in Table 5.

Figure 16: Dimension of the model that is used to create the mesh file

One of the runs was not completed properly. The element size was too big and the calculation stopped
due to negative energy in one element. The peak pressures of all the other simulations were quite
similar except for the model map_para6. In this model the maximum cross section was chosen to high.
This resulted in a configuration that cannot be seen as a spherical condition and is therefore not
considered.

A very small influence of the element size on the peak pressure is found. map_para4 produces the
highest and sharpest peak pressure. The needed calculation time for all models is very small. Even the
very fine model with Δel=0.1 mm and 20000 elements uses only 20 minutes of CPU time.

In comparison with the empirical data from Kingery the current results for the peak pressure are
smaller (Figure 17). There are several possible reasons for this difference. The used 1D finite volume
implementation can until now be used only under first order formulation in space and time, and this
can result in smaller peak pressure values. In addition, the values of Kingery have been mainly
obtained using very large explosions and theoretically scaling them down to small events. It is not
clear if these values fit also very well the small charges. Bogosian et al. [2] presented a comparison of

19
Kingery’s data with data from several other sources and experiments. The comparison showed that
these data systematically underestimated Kingery’s data values.

Table 5: Parameter study to create the map file

Model Element Truncated cone: Truncated cone: Vfcc options Result

length el min dimension (amin) max dimension
(amax)

map_para1 0.5 mm 25e-6 0.05 m 6

map_para2 1.0 mm 25e-6 0.05 m 6

map_para3 2.0 mm 25e-6 0.05 m 6 Tilt

map_para4 0.1 mm 25e-6 0.05 m 6

map_para5 0.5 mm 25e-6 0.1 m 6

map_para6 0.5 mm 25e-6 0.01 m 6

map_para7 0.5 mm 25e-6 0.05 m 3

map_para8 0.5 mm 5e-6 0.05 mm 6

Figure 17 shows that all pressure results are slightly lower than Kingery’s data but they are all very
similar. Therefore, the map file created by the model map_para1 is used for the mapping.

Figure 17: Pressure-distance curve of map files (Table 5)

20
3.3.2 3D simulations
The results of the 1D simulations were next mapped into a 3D model. This model has an (hexahedral)
element size of 2 cm. Two different models are used: a small model (Figure 15) and a full model with
44400 and 284400 elements, respectively. To see the influence of the region where the mapping is
applied, the simulations are performed with four different configurations each with a different distance
of the blast wave to the centre when the blast wave is mapped. This distance varies from 0.2 m to 1.1
m. The models can be compared with the results from the previous report [8] where a bubble (balloon)
and a solid TNT model have been used. All model configurations are presented in Table 6.

Table 6: Parameter study for mapping of Alia experiment

Model Mapping of Small/full Element End of calculation Calculation

blast wave at model size (3D) time (3D)

alia7_map1 0.4 m small 2 cm 3.00 ms 99 s

alia7_map2 1.0 m small 2 cm 3.00 ms 87 s

alia7_map3 0.2 m small 2 cm 3.00 ms 107 s

alia7_map4 1.1 m small 2 cm 4.00 ms 101 s

alia7_map1_full 0.4 m full 2 cm 5.00 ms 1074 s

alia7_map2_full 1.0 m full 2 cm 5.00 ms 853 s

alia7_map3_full 0.2 m full 2 cm 5.00 ms 1051 s

alia7_map4_full 1.1 m full 2 cm 4.00 ms 1031 s

alia_bubb4 bubble small 2 cm 3.00 ms 156 s

alia_bubb7 bubble full 2 cm 6.00 ms 1524 s

alia7 Solid TNT small 2 cm 3.00 ms 2252 s

alia20 Solid TNT small 1 cm 3.00 ms 36334 s

Figure 18 indicates that the differences in the pressure-time histories for the different starting points of
the mapping (distance of the initial blast wave from the centre) are very small. This means first, that
the meshes are sufficiently small and represent the blast wave development accurately enough and

21
second, that the mapping can also be done shortly before the point where the first reflection takes
place. The calculation time of the 3D model is then smaller. This distance is here 1.2 m.

While the first peak pressure is represented quite well, the reflections “from the stand” and from the
wall cannot be described accurately enough in the full models. Both reflections reach the investigated
point (pressure gage) earlier and in the case of the reflection from the wall the pressure value is
substantially smaller than the experimental one. Nevertheless, the fundamental behaviour is well
described. Clearly the small model can capture only the initial features of the pressure record.

Also Alia [1] shows a difference of the arrival time of the blast wave between the experiment and the
numerical simulations. While the experiment shows the arrival of the blast wave at about 1.5 ms the
arrival time in the numerical simulations are about 0.5 ms later. To allow a better comparison of the
numerical and the experimental results, all curves are moved to a similar arrival time of 1.5 ms.

Figure 18: Pressure history for Alia experiment: comparison of mapping parameters (Table 6)

When the results of the mapping models are compared with the bubble and solid TNT models (Figure
19) it can be seen that the peak-pressure arrives earlier for the bubble model and is slightly higher.
This results most probably from the effects already mentioned for the calculation of the mapping file.
The peak pressure in the solid TNT model is slightly smaller, most probably since the element size
inside the explosive was not small enough. The results of the small model are all quite similar except
the already mentioned difference for the reflections for the mapping model.

22
Since the results are quite similar the calculation time can be compared in order to choose the
appropriate simulation model. The solid TNT models have a very high computational cost but they are
only used for calculating the small model. The calculation time for the full model is too long. In
comparison, the mapping simulation runs much faster. The calculation time of the bubble model is
slightly higher than the one of the mapping model. If the calculation time for creating the map file is
added, the computational costs are most probably of the same order of magnitude. The two models
(mapping and bubble) will be compared with their advantages and drawbacks in the conclusions
section.

Figure 19: Pressure history for Alia experiment: comparison with previous results (Table 6)

3.4 Simulation of an explosion inside a train

In order to test the mapping algorithm also for larger structures the simulations of explosions inside
trains (previous report [9] and publication [11]) were repeated and instead of the bubble model the
mapping technology is used. The model consists of a regular fluid mesh with a mesh size of 10 cm.
The structure is built from stiff beam elements and weaker shell elements for the structural surface.
The same geometry as reported in the previous work [9] is used.

The charge is detonated inside the train. A charge with 9 kg TNT equivalent is used which correlates
to a rucksack bomb.

23
The charge is first calculated in a 1D mesh and then mapped into the 3D case. The 1D calculation uses
9 kg TNT equivalent (radius of the spherical charge of 11 cm) and a total length of the 1D model of 1
m. The element size in that model is 1 mm. In that way, the charge is represented by 110 elements. The
mapping file is written when the blast wave arrives a distance of 1 m from the charge centre. This
corresponds to the distance of the train structure to the charge.

In comparison to the calculation with the compressed balloon the 3D calculation is quite unstable. At
some points the time step size drops to 1e-8 s while for the beginning and for several periods in
between the parts with small time step sizes the time step size is in the rage of 1e-6 s. This very small
time step occurs only in one or some elements. There could be two reasons:

1. Mapping is an exact physical model. That means that the pressure wave is reflected inside the
charge after some time. That leads to very high pressures and with them the time step size
decreases a lot.

2. Some of the elements with very small time step sizes are along the boundary condition. It could
be that the boundary condition has some problems with very high pressures. This phenomenon
should be checked with additional models.

Nevertheless, the calculation results in a deformed and partly failed train carriage, as it can be seen in
Figure 20 (bottom). The comparison between balloon and mapping method shows that the
displacements in case of the mapping method are smaller. It cannot be stated which one of the
approaches provides a better solution due to the lack of possibilities to compare these results with real
experimental data.

24
Figure 20: Comparison of the displacements and the failure of the train carriages using balloon (top)
and mapping (bottom) method

25
4 How to use mapping
In order to facilitate the use of the mapping technology the most important steps are presented here:

 Create a 1D mesh for the calculation of the data to be stored in the mapping file. This
calculation should be done using 1D finite volume elements (TUVF). The element size can be
very small such the calculation time of the 1D problem is typically very short. A discretization
of about 50 elements in the charge could be an indicative value.

 The cone-type geometry for the 1D model should be used, as shown in the examples in the
appendix. The mapping file is written by using the command MAPB in the ECRI command.
The propagated distance of the mapped blast wave from the centre must also be given. This
should be chosen in such a way that no reflections have occurred yet. The distance should also
be chosen as long as possible in order to reduce the computational cost of the 3D model.

 The mapping file must be read by using INIT MAPB. The centre of the spherical blast wave
can be indicated by using POS. The 3D model should also be modelled with finite volumes
since they appear to provide the best description of the blast wave development. Finite
elements often reduce the blast wave too much. FCONV 6 (HLLC-solver) is recommended for
the 1D as well as for the 3D model.

 The simulation time of the 3D model starts again with zero timing. In order to get the time
from the initiation of the explosion the time indicated in the mapping file must be added to the
current one.

26
5 Conclusions
The report presents the mapping method that can be used to simulate the development of blast waves
in the air very effectively. The idea is to do a 1D simulation with a very fine mesh for the explosive
and the air until the blast wave reaches the first obstacle. Just before the wave reaches this obstacle the
wave is mapped into a 3D mesh (in the future also 2D meshes would be possible).

This procedure provides a physically based approach to speed up the simulation of blast wave
development. In comparison to other methods (bubble method, solid TNT method, AIRB) several
advantages and drawbacks can be summarized.

Mapping is physical while the bubble model is only an approximation. In both cases a pressure field is
applied on the fluid elements. While for the bubble model the pressure is constant and no particle
velocity is considered the mapping considers the pressure-distance curve and also the particle velocity.
Nevertheless, the bubble model is easier to use and easy to adapt to other charges and distances. For
the mapping algorithm two runs are needed. It is slightly more complicated to perform and to adapt to
other parameters (e.g. charge size). The description of the wave after the peak pressure is better in case
of the mapping model but is shows slightly smaller peak pressures most probably due to first order
approach in 1D.

Solid TNT models are very expensive. To reach appropriate values very small elements are needed at
least in the explosive and around it. Also if bigger elements are used in zones with a smaller pressure
ratio (far field of the explosion) the time step is very small. The results are not better than the one of
the mapping method. Full solid TNT models use mainly the JWLS material law for the explosive as
well as for the air. The additional terms for the calculation of the explosive are only needed during the
explosion effect. Afterwards this terms results in unnecessary numerical costs.

In summary, the mapping method is an additional method in EUROPLEXUS that can be efficiently
used to simulate the development of air blast waves.

27
6 References
[1] A. Alia and M. Souli. High explosive simulation using multi-material formulations. Applied
Thermal Engineering, 26:1032–1042, 2006.

[2] D. Bogosian, J. Ferritto, and Y. Shi. Measuring uncertainty and conservatism in simplified
blast models. In 30th Explosives Safety Seminar, Atlanta, Georgia, 2002.

[3] I.G. Cullis. Blast waves and how they interact with structures. Journal of the Royal Army
Medical Corps, 147:16–26, 2001.

[4] Manmohan Goel, Vasant Matsagar, Anil Gupta, and Steffen Marburg. An abridged review of
blast wave parameters. Defence Science Journal, 62(5):300–306, 2012.

[5] Charles N. Kingery and Gerald Bulmash. Airblast parameters from TNT spherical air burst and
hemispherical surface burst. Technical report, Defence Technical Information Center, Ballistic
Research Laboratory, Aberdeen Proving Ground, Maryland, 1984.

[6] Gilbert F. Kinney and Kenneth J. Graham. Explosive Shocks in Air. Springer, Berlin;
Heidelberg; New York; Tokyo, 1985.

[7] Martin Larcher. Pressure-time functions for the description of air blast waves. Technical Report
JRC46829, Joint Research Centre, Ispra, Italy, 2008.

[8] Martin Larcher. Simulation of air blast waves in urban environment. Technical Report
JRC43400, Joint Research Centre, Ispra, Italy, 2008.

[9] Martin Larcher. Simulations of a metro carriage exposed to an internal detonation. Technical
Report JRC50327, 2009.

[10] Martin Larcher and Folco Casadei. Explosions in complex geometries - a comparison of
several approaches. International Journal of Protective Structures, 1(2):169–195, 2010.

[11] Martin Larcher, Folco Casadei, Georgios Giannopoulos, George Solomos, Jean-Luc Planchet,
and Anne Rochefrette:. Determination of the risk due to explosions in railway systems. Journal of Rail
and Rapid Transit, in press:1–10, 2011.

28
7 Appendix
Cour 5 'VF 500'ecro comp 1 elem lect 500 term
Cour 6 'VF 1000' ecro comp 1 elem lect 1000
7.1 Creation of a map file term
Cour 7 'VF 1500' ecro comp 1 elem lect 1500
1 kg TNT equivalent in 1 m distance term
TRAC 1 axes 1.e-5 'pressure (bar)'
TRAC 2 axes 1.e-5 'pressure (bar)'
bm_flu_map1.dgibi TRAC 3 axes 1.e-5 'pressure (bar)'
TRAC 4 axes 1.e-5 'pressure (bar)'
TITRE 'raccord TUBM et VF' ; TRAC 5 axes 1.e-5 'pressure (bar)'
OPTION DIME 3 ELEM CUB8 ; TRAC 6 axes 1.e-5 'pressure (bar)'
option sort 'map1d.msh' ; TRAC 7 axes 1.e-5 'pressure (bar)'
*
tl_tnt = 0.052;
* about 1 kg tnt (0.96 kg) 7.2 map3d test cases
tl_air = 1.948;
origine = 0. 0. 0. ; map3d-fine-fine.dgibi
pt0 = origine plus (0 0 0) ;
pt1 = pt0 plus (tl_tnt 0 0) ; OPTI echo 1;
pt2 = pt1 plus (tl_air 0 0) ; OPTI dime 3 elem cub8;
tu_tnt = pt0 DROIT 52 pt1 ; DENS 0.02;
tu_air = pt1 DROIT 1948 pt2 ; sizx = 2.0;
cl_2 = MANU POI1 pt2 ; sizy = 2.0;
* sizz = 2.0;
tube = tu_air et tu_tnt; p0 = 0 0 0;
mesh = tube et cl_2; x0 = (sizx) 0. 0.;
* y0 = 0. (sizy) 0.;
sortier mesh; z0 = 0. 0. (sizz);
fin; *volume of the air
vol0 = p0 droi x0 tran z0 volu tran y0 coul
bm_flu_map1.epx bleu;
*bubble
Bench to test writing of the map file ages = vol0 enve;
ECHO
CAST FORM mesh a_abso1 = p0 droi x0 tran z0 plus y0;
TRID euler a_abso2 = p0 droi y0 tran z0 plus x0;
GEOM a_abso3 = p0 droi x0 tran y0 plus z0;
TUVF tube a_abso = a_abso1 et a_abso2 et a_abso3;
TERM
COMP a_symm1 = p0 droi x0 tran z0;
DIAM CONE D1 0.000025 D2 0.05 a_symm2 = p0 droi y0 tran z0;
ORIG LECT pt0 TERM LIST LECT tube TERM a_symm3 = p0 droi x0 tran y0;
MATE a_fsr = a_symm1 et a_symm2 et a_symm3;
*jwl air
JWLS a 3.738e11 b 3.747e9 r1 4.15 r2 0.90 elim (ages et a_abso et a_fsr);
omeg 0.35 ros 1630 BETA 0.25
ro 1.3 pini 1e5 eint 0.21978e6 pref 0 geom_new = (vol0 et a_abso et a_fsr);
LECT tu_air TERM
* TASS geom_new;
JWLS a 3.738e11 b 3.747e9 r1 4.15 r2 0.90 OPTI sauv form 'map3d-fine-fine.msh';
omeg 0.35 d 6930 BETA 0.25 sauv form geom_new;
ro 1630 pini 1e5 eint 3.68e6 pref 0 fin;
xdet 0. ydet 0. zdet 0.
LECT tu_tnt TERM map3dvfcc-fine-fine.epx
LINK BLOQ 1 LECT pt0 pt2 TERM
ECRI Titel
FICH PVTK FORM TFRE 1.0e-4 VARI ECRO VCVI ECHO
FICH MAPB MSPA DIPR 1.0 PCHE 1.5e5 FREQ 1 CAST 'map3d-fine-fine.msh' vol0
FICH ALIC tfre 1e-4 TRID EULER RISK
FICH ALIC TEMP tfre 1e-6 OPTI TOLC 1e-1
poin lect 1 term $
elem lect 1 26 45 50 52 60 500 1000 1500 GEOM CUVF vol0 TERM
term COMP
OPTI notest noprint GROU 1 'x1' LECT vol0 TERM COND NEAR POIN
cstab 0.9 LOG 1 0.0 0.0 0.5
VFCC FCONV 6 ORDR 2 NGRO 1 'ns1' LECT vol0 TERM COND LINE X1
CALC tini 0 nmax 30000 tfin 2.0e-3 0.1 Y1 0.1 Z1 0
* X2 0.1 Y2 0.1 Z2 5.0 TOL 0.01
SUITE *
Post-treatment MATE
RESU ALIC TEMP SORT GRAP $ air
AXTE 1.0 'Time [s]' GAZP RO 1.3 PINI 1E5 GAMMA 1.35 PREF 1e5
Cour 1 'VF 1 ' ecro comp 1 elem lect 1 term LECT vol0 TERM
Cour 2 'VF 26 ' ecro comp 1 elem lect 26 term INIT MAPB 'map1d.map' POS 0 0.0 0.0
Cour 3 'VC 45 ' ecro comp 1 elem lect 45 term ECRI FICH ALIC TEMP TFREQ 2e-6
Cour 4 'VF 60 ' ecro comp 1 elem lect 60 term ELEM LECT 60 x1 TERM

29
 POIN LECT ns1 TERM MATE
FICH PVTK TFREQ 1.0e-4 VARI ECRO CONT VCVI $ air
RISK GAZP RO 1.3 PINI 1E5 GAMMA 1.35
$ PREF 1e5
OPTI NOTE LOG 1 STEP IO LECT flui TERM
CSTA 0.50 VFCC BUBB MASS 0.125 ! must be 1/8 - symmetry
VFCC LECT bubb TERM
FCONV 6 ECRI FICH ALIC TEMP TFREQ 2e-6
ORDRE 2 OTPS 2 ELEM LECT 60 x1 x2 x3 x4 x5 x6 x7 x8
RECONS 1 x9 s1 ns1 TERM
LMAS 3 LQDM 3 LENE 3 POIN LECT s1 ns1 TERM
KMAS 0.75 KQDM 0.75 KENE 0.75 FICH PVTK TFREQ 1.0e-4 VARI ECRO CONT VCVI
CENER RISK
$ $
CALC TINI 0 TEND 3e-3 PAS1 1.E-8 OPTI NOTE LOG 1 STEP IO
*============================================= CSTA 0.50 VFCC
SUIT VFCC
POSTTREATMENT FROM ALICE TEMPS FILE FCONV 6
$ ORDRE 2 OTPS 2
ECHO RECONS 1
$ LMAS 3 LQDM 3 LENE 3
RESU ALIC TEMP GARD PSCR KMAS 0.75 KQDM 0.75 KENE 0.75
* CENER
OPTI PRIN $
SORT GRAP CALC TINI 0 TEND 3e-3 PAS1 1.E-8
* *=============================================
AXTE 1.0 'Time [s]' POSTTREATMENT FROM ALICE TEMPS FILE
* ECHO
COUR 1 'x1' ECRO COMP 1 ELEM LECT 60 TERM RESU ALIC TEMP 'map3dvfcc-fine-fine_bub.alt'
* GARD PSCR
trac 1 axes 1.0 'PRESS. [PA]' *
list 1 axes 1.0 'PRESS. [PA]' OPTI PRIN
FIN SORT GRAP
AXTE 1.0 'Time [s]'
*
map3dvfcc-fine-fine-bub.epx COUR 1 'x1' ECRO COMP 1 ELEM LECT 60 TERM
trac 1 axes 1.0 'PRESS. [PA]'
Titel list 1 axes 1.0 'PRESS. [PA]'
ECHO FIN
! VERI
* CONV WIN
CAST 'map3d-fine-fine.msh' vol0 7.3 alia tests
TRID EULER RISK
OPTI TOLC 1e-1
$ 7.3.1 Parameter test for creation of
GEOM
CUVF vol0 mapping file
TERM
COMP map_para1.dgibi
GROU 12 'x1' LECT vol0 TERM COND NEAR POIN
0.0 0.0 0.5 TITRE 'raccord TUBM et VF' ;
'x2' LECT vol0 TERM COND NEAR POIN OPTION DIME 3 ELEM CUB8 ;
0.0 0.0 1.0 option sort 'map_para1.msh' ;
'x3' LECT vol0 TERM COND NEAR POIN *
0.0 0.0 1.5 tl_tnt = 0.034;
'x4' LECT vol0 TERM COND NEAR POIN * about 268 g tnt
0.0 0.0 2.0 tl_air = 1.966;
'x5' LECT vol0 TERM COND NEAR POIN origine = 0. 0. 0. ;
0.0 0.0 2.5 pt0 = origine plus (0 0 0) ;
'x6' LECT vol0 TERM COND NEAR POIN pt1 = pt0 plus (tl_tnt 0 0) ;
0.0 0.0 3.0 pt2 = pt1 plus (tl_air 0 0) ;
'x7' LECT vol0 TERM COND NEAR POIN tu_tnt = pt0 DROIT 68 pt1 ;
0.0 0.0 3.5 tu_air = pt1 DROIT 3932 pt2 ;
'x8' LECT vol0 TERM COND NEAR POIN cl_2 = MANU POI1 pt2 ;
0.0 0.0 4.0 *
'x9' LECT vol0 TERM COND NEAR POIN tube = tu_air et tu_tnt;
0.0 0.0 4.5 mesh = tube et cl_2;
's1' LECT vol0 TERM COND BOX X0 0 Y0 sortier mesh;
0 Z0 0 fin;
DX 0.1
DY 0.1 DZ 5.0 map_para1.epx
'bubb' LECT vol0 TERM
COND SPHE R 0.35 XC 0.0 YC VF
0.0 ZC 0.0 ECHO
'flui' LECT vol0 DIFF bubb TERM CAST FORM mesh
NGRO 1 'ns1' LECT vol0 TERM COND LINE X1 TRID euler
0.1 Y1 0.1 Z1 0 GEOM
X2 0.1 Y2 0.1 tuvf tube
Z2 5.0 TOL 0.01 TERM
* COMP

30
 DIAM CONE D1 0.000025 D2 0.05 SCOUR 101 '0.2ms' ECRO COMP 1 T 0.2e-3 SAXE 1.0
ORIG LECT pt0 TERM LIST LECT tube TERM 'ABSC' LECT tube TERM
MATE SCOUR 102 '0.4ms' ECRO COMP 1 T 0.4e-3 SAXE 1.0
*jwl air 'ABSC' LECT tube TERM
jwls a 3.738e11 b 3.747e9 r1 4.15 SCOUR 103 '0.5ms' ECRO COMP 1 T 0.5e-3 SAXE 1.0
r2 0.90 'ABSC' LECT tube TERM
omeg 0.35 ros 1630 BETA 0.25 SCOUR 104 '1.0ms' ECRO COMP 1 T 1.0e-3 SAXE 1.0
ro 1.3 pini 1e5 eint 'ABSC' LECT tube TERM
0.21978e6 pref 0 TRACE 99 100 101 102 103 104 TEXT axes 1.0
LECT tu_air TERM 'Pression (Pa)'
*
** Le TNT : on donne directement ro = ros SCOUR 200 '0.1ms' VCVI COMP 1 T 0.1e-3 SAXE 1.0
* avec ignition au point P7 'ABSC' LECT tube TERM
* la vitesse de detonation est celle de SCOUR 201 '0.2ms' VCVI COMP 1 T 0.2e-3 SAXE 1.0
Chapman-Jouguet 'ABSC' LECT tube TERM
* SCOUR 202 '0.4ms' VCVI COMP 1 T 0.4e-3 SAXE 1.0
jwls a 3.738e11 b 3.747e9 r1 4.15 'ABSC' LECT tube TERM
r2 0.90 SCOUR 203 '0.5ms' VCVI COMP 1 T 0.5e-3 SAXE 1.0
omeg 0.35 d 6930 BETA 0.25 'ABSC' LECT tube TERM
ro 1630 pini 1e5 eint SCOUR 204 '1.0ms' VCVI COMP 1 T 1.0e-3 SAXE 1.0
3.68e6 pref 0 'ABSC' LECT tube TERM
xdet 0. ydet 0. zdet 0. TRACE 200 201 202 203 204 TEXT axes 1.0
LECT tu_tnt TERM 'Particle velocity (m/s)'
LINK BLOQ 1 LECT pt0 pt2 TERM
ECRI TRACE 100 axes 1.0 'Pression (Pa)'
FICH MAPB MSPA DIPR 1.0 PCHE 1.5e5 FREQ 1 TRACE 101 axes 1.0 'Pression (Pa)'
FICH ALICE tfreq 1e-4 TRACE 102 axes 1.0 'Pression (Pa)'
FICH ALICE TEMPS tfreq 1e-6 TRACE 103 axes 1.0 'Pression (Pa)'
poin lect 1 term TRACE 104 axes 1.0 'Pression (Pa)'
elem lect 1 26 45 50 52 60 500 1000 1500
term FIN
*
OPTI notest noprint
cstab 0.9 7.3.2 Alia simulation with small
*
VFCC ! DUMP
model
FCONV 6 ORDR 2
alia_map.dgibi
CALC tini 0 nmax 30000 tfin 2.0e-3
* *Construction d'une sphere a partir d'un cube
SUITE opti dime 3 elem cub8;
Post-treatment 1 *Nombre de bissections
RESU ALICE TEMPS nel0 = 7;
SORT GRAP nel1 = 7;
AXTE 1.0 'Time [s]' *Cote du cube intermediaire
r0 = .25;
Courbe 1 'VFCC 1 ' ecrou comp 1 elem lect 1 sizeex = 0.0323;
term sizeai = 0.14;
Courbe 2 'VFCC 26 ' ecrou comp 1 elem lect width = 0.4;
26 term height = 2.0;
Courbe 3 'VFCC 45 ' ecrou comp 1 elem lect dini = 3.141*sizeex/(4.*nel0);
45 term dfin = 3.141*sizeai/(4.*nel0);
* *Reference
Courbe 31 'VFCC 60 ' ecrou comp 1 elem lect o0 = 0. 0. 0.;
60 term x0 = (sizeex) 0. 0.;
Courbe 4 'VFCC 500 ' ecrou comp 1 elem lect xa0 = 0 (sizeex) (sizeex);
500 term xb0 = (sizeex) (sizeex) (sizeex);
Courbe 5 'VFCC 1000 ' ecrou comp 1 elem lect xc0 = (sizeex) 0 (sizeex);
1000 term xd0 = 0 0 (sizeex);
Courbe 6 'VFCC 1500 ' ecrou comp 1 elem lect y0 = 0. (sizeex) 0.;
1500 term z0 = 0. 0. (sizeex);
TRACE 1 axes 1.e-5 'Pression (bar)' x1 = (sizeai) 0. 0.;
TRACE 2 axes 1.e-5 'Pression (bar)' y1 = 0. (sizeai) 0.;
TRACE 3 axes 1.e-5 'Pression (bar)' z1 = 0. 0. (sizeai);
TRACE 31 axes 1.e-5 'Pression (bar)' c0 = x0 plus y0 plus z0 / 2.;
TRACE 4 axes 1.e-5 'Pression (bar)' c1 = x1 plus y1 plus z1 / 2.;
TRACE 5 axes 1.e-5 'Pression (bar)' symp1 = (sizeex/2.) (sizeex/2.) (sizeex/2.);
TRACE 6 axes 1.e-5 'Pression (bar)' symp2 = (sizeai/2.) (sizeai/2.) (sizeai/2.);
* *Cube intermediaire (centre=o0 et arete=r0)
SUITE vol0 = o0 droi nel0 x0 tran nel0 z0 volu tran
Post-treatment 2 nel0 y0
RESU ALICE coul bleu homo o0 r0;
SORT GRAP vol10 = o0 droi nel0 x0 tran nel0 z0 volu tran
AXTE 1.0 'Time [s]' nel0 y0
* coul bleu homo o0 0.2;
SCOUR 99 '0.05ms' ECRO COMP 1 T 0.05e-3 SAXE cub0 = (o0 droi nel0 y0 tran nel0 x0)
1.0 'ABSC' LECT tube TERM et (o0 droi nel0 z0 tran nel0 y0) et
SCOUR 100 '0.1ms' ECRO COMP 1 T 0.1e-3 SAXE 1.0 (o0 droi nel0 x0 tran nel0 z0)
'ABSC' LECT tube TERM syme 'POINT' symp1 homo o0 r0;
cub1 = (o0 droi nel0 y0 tran nel0 x0)

31
 et (o0 droi nel0 z0 tran nel0 y0) et
(o0 droi nel0 x0 tran nel0 z0)
alia7_map1.epx
syme 'POINT' symp1; alia1 - finite volume with mapping
cub2 = (o0 droi nel0 y1 tran nel0 x1) ECHO
et (o0 droi nel0 z1 tran nel0 y1) et CAST 'alia_map.msh' mesh
(o0 droi nel0 x1 tran nel0 z1) TRID EULER
syme 'POINT' symp2; OPTI TOLC 1e-1
*Pojection sur la sphere de rayon unitaire *
spe1 = cub0 proj 'CONI' o0 'SPHE' o0 x0; GEOM
spe2 = cub2 proj 'CONI' o0 'SPHE' o0 x1; CUVF air
CL3D a_abso
*Remplissage TERM
vol1 = cub0 volu nel1 spe1 coul vert; *
vol1 = vol1 et vol0; MATE
vol2 = spe1 volu 'DINI' dini 'DFIN' dfin cub2 GAZP RO 1.3 PINI 1E5 GAMMA 1.35
coul bleu; PREF 1e5
LECT air TERM
CLVF ABSO RO 1.3 lect a_abso term
x2 = (width) 0. 0.; INIT MAPB 'map1d_a1.map' POS 0 0.0 0.0
y2 = 0. (width) 0.;
z2 = 0. 0. (height); ECRI FICH ALIC TEMP tfreq 3e-6
dw = 20; ELEM LECT 3404 fe1 TERM
dh = 100; FICH ALIC TFRE 5E-5
air = o0 d dw x2 tran dh z2 volu FICH PVTK TFREQ 1.0e-4 VARI ECRO CONT VCVI
tran dw y2 coul bleu; RISK
*
vges = (air); OPTION NOTEST PASAUTO
elim vges 1e-6; *
afull = enve vges; CSTA 0.50 VFCC
VFCC
kod1 = faux; FCONV 6
kod2 = faux; ORDRE 2 OTPS 2
kod3 = faux; RECONS 1
kod4 = faux; LMAS 3 LQDM 3 LENE 3
REPE I0 (NBEL afull); KMAS 0.75 KQDM 0.75 KENE 0.75
symm = faux; CENER
xx yy zz = coor (bary (afull ELEM QUA4 &I0)); *
SI ( xx < 0.0001 ); OPTI NOTE LOG 1
si (kod1);a_symmx=(a_symmx et (afull ELEM CALC TINI 0 TEND 3e-3
QUA4 &I0)); *=============================================
sinon; a_symmx = afull ELEM QUA4 &I0;kod1 = SUIT
vrai; Post-treatment (time curves from alt file)
fins; ECHO
symm = vrai; *
FINS; RESU ALIC TEMP GARD PSCR
SI ( yy < 0.0001); SORT GRAP
si (kod2);a_symmy=(a_symmy et (afull ELEM AXTE 1.0 'Time [s]'
QUA4 &I0)); *
sinon; a_symmy = afull ELEM QUA4 &I0;kod2 = COUR 4 'fe1' ECRO COMP 1 ELEM LECT 3404 TERM
vrai; trac 4 TEXT axes 1.0 'Pressure [Pa]'
fins; list 4 TEXT axes 1.0 'Pressure [Pa]'
symm = vrai; fin
FINS;
SI ( zz < 0.0001);
si (kod3);a_symmz=(a_symmz et (afull ELEM 7.3.3 Alia simulation with full model
QUA4 &I0));
sinon; a_symmz = afull ELEM QUA4 &I0;kod3 =
vrai;
alia_full.dgibi
fins; *ALIA model with the air upto the wall.
symm = vrai; opti dime 3 elem cub8;
FINS; *Nombre de bissections
SI (symm NEG vrai); sizeai = 0.1;
si (kod4);a_abso=(a_abso et (afull ELEM nbel_bub = 5;
QUA4 &I0)); widthx = 1.2;
sinon; a_abso = afull ELEM QUA4 &I0;kod4 = widthy = 0.8;
vrai; height = 2.0;
fins; *Reference
FINS; o0 = 0. 0. 0.;
FIN I0;
*Remplissage
fp1 = vol0 poin proche (0.0 0.0 1.52); x2 = (widthx) 0. 0.;
fe1 = vol0 elem contenant fp1; y2 = 0. (widthy) 0.;
z2 = 0. 0. (height);
mesh = (vges et a_abso et a_symmx et a_symmy et dwx = 60;
a_symmz et fe1); dwy = 40;
dh = 100;
TASS mesh; air = o0 d dwx x2 tran dh z2 volu
OPTI sauv form 'alia_map.msh'; tran dwy y2 coul bleu;
sauv form mesh; vges = air;
elim vges 1e-6;

32
afull = enve vges; GAZP RO 1.3 PINI 1E5 GAMMA 1.35
PREF 1e5
kod1 = faux; LECT air TERM
kod2 = faux; CLVF ABSO RO 1.3 lect a_abso term
kod3 = faux; INIT MAPB 'map1d_a1.map' POS 0 0.0 0.0
kod4 = faux;
kod5 = faux; ECRI FICH ALIC TEMP tfreq 3e-6
REPE I0 (NBEL afull); ELEM LECT 3404 fe1 TERM
symm = faux; FICH ALIC TFRE 5E-5
xx yy zz = coor (bary (afull ELEM QUA4 &I0)); FICH PVTK TFREQ 1.0e-4 VARI ECRO CONT VCVI
SI ( xx < 0.0001 ); RISK
si (kod1);a_symmx=(a_symmx et (afull ELEM GROU 1 OBJE 'air' LECT air TERM
QUA4 &I0)); *
sinon; a_symmx = afull ELEM QUA4 &I0;kod1 = OPTION NOTEST PASAUTO
vrai; *
fins; CSTA 0.50 VFCC
symm = vrai; *
FINS; VFCC
SI ( yy < 0.0001); FCONV 6
si (kod2);a_symmy=(a_symmy et (afull ELEM ORDRE 2 OTPS 2
QUA4 &I0)); RECONS 1
sinon; a_symmy = afull ELEM QUA4 &I0;kod2 = LMAS 3 LQDM 3 LENE 3
vrai; KMAS 0.75 KQDM 0.75 KENE 0.75
fins; CENER
symm = vrai; *
FINS; OPTI NOTE LOG 1
SI ( zz < 0.0001); CALC TINI 0 TEND 5e-3
si (kod3);a_symmz=(a_symmz et (afull ELEM *
QUA4 &I0)); *=============================================
sinon; a_symmz = afull ELEM QUA4 &I0;kod3 = SUIT
vrai; Post-treatment (time curves from alice temps
fins; file)
symm = vrai; ECHO
FINS; *
SI ( xx > (widthx-0.0001)); RESU ALIC TEMP GARD PSCR
si (kod4);a_fix=(a_fix et (afull ELEM QUA4 SORT GRAP
&I0)); AXTE 1.0 'Time [s]'
sinon; a_fix = afull ELEM QUA4 &I0;kod4 = *
vrai; COUR 4 'fe1' ECRO COMP 1 ELEM LECT 3404 TERM
fins; trac 4 TEXT axes 1.0 'Pressure [Pa]'
symm = vrai; list 4 TEXT axes 1.0 'Pressure [Pa]'
FINS; fin
SI (symm NEG vrai);
si (kod5);a_abso=(a_abso et (afull ELEM
QUA4 &I0));
sinon; a_abso = afull ELEM QUA4 &I0;kod5 = 7.3.4 Explosion inside a train
vrai;
fins; map_9kg_1m.dgibi
FINS;
FIN I0;
TITRE 'raccord TUBM et VF' ;
fp1 = vges poin proche (0.0 0.0 1.52); OPTION DIME 3 ELEM CUB8 ;
fe1 = vges elem contenant fp1; option sort 'map_9kg_1m.msh' ;
*
mesh = (vges et a_abso et a_symmx et a_symmy et tl_tnt = 0.11;
a_symmz * about 9 kg tnt ( kg)
et a_fix et fe1); tl_air = 1.89;
origine = 0. 0. 0. ;
TASS mesh; pt0 = origine plus (0 0 0) ;
OPTI sauv form 'alia_full.msh'; pt1 = pt0 plus (tl_tnt 0 0) ;
sauv form mesh; pt2 = pt1 plus (tl_air 0 0) ;
mess(nbel(mesh)); tu_tnt = pt0 DROIT 110 pt1 ;
mess(nbno(mesh)); tu_air = pt1 DROIT 1890 pt2 ;
mess(mesu(bub)); cl_2 = MANU POI1 pt2 ;
*
alia_map1_full.epx tube = tu_air et tu_tnt;
mesh = tube et cl_2;
alia1 - finite volume with mapping *
$ sortier mesh;
ECHO fin;
CAST 'alia_full.msh' mesh
TRID EULER map_9kg_1m.epx
OPTI TOLC 1e-1
*
GEOM Calculation of the blast wave in 1 m distance
CUVF air with a charge of 9 kg TNT
CL3D a_abso *
TERM ECHO
* CAST FORM mesh
MATERIAU TRID euler

33
GEOM TUVF tube VM23 RO 7800. YOUNG 2.1E11 NU 0.333 ELAS
TERM 1.2E8
* FAIL VMIS LIMI 1.6E8
COMP TRAC 2 1.2E8 1.1E-3 1.6E8 0.13
DIAM CONE D1 0.000025 D2 0.05 LECT sall s_tri seats TERM
ORIG LECT pt0 TERM LIST LECT tube TERM VMIS ISOT RO 2700 YOUNG 7E10 NU 0.3 ELAS
MATE 200E6 FAIL 1 LIMI 2.1e8
*jwl air TRAC 2 200e6 2.85e-3
JWLS a 3.738e11 b 3.747e9 r1 4.15 r2 0.90 250e6 0.23
omeg 0.35 ros 1630 BETA 0.25 LECT segm segm_top TERM
ro 1.3 pini 1e5 eint 0.21978e6 pref 0 GLAS RO 2500 YOUN 7E10 NU 0.23 CORR 16.0
LECT tu_air TERM FAIL PSAR LIMI 159.6e6
* LECT wall doors TERM
JWLS a 3.738e11 b 3.747e9 r1 4.15 r2 0.90 VM23 RO 1400.0 YOUN 4E7 NU 0.3
omeg 0.35 d 6930 BETA 0.25 ELAS 10E7
ro 1630 pini 1e5 eint 3.68e6 pref 0 FAIL PEPS LIMI 0.2
xdet 0. ydet 0. zdet 0. TRAC 1 10e7 0.25
LECT tu_tnt TERM LECT hb1 TERM
LINK BLOQ 1 LECT pt0 pt2 TERM CLVF ABSO RO 1.4 LECT ages TERM
ECRI *--------------------------Boundary conditions
FICH PVTK FORM TFRE 1.0e-4 VARI ECRO VCVI LINK DECO FLSW STRU LECT sall s_tri wall doors
FICH MAPB DIPR 1.0 PCHE 1.5e5 FREQ 1 seats hb1 TERM
FICH ALIC tfre 1e-4 FLUI LECT vol0 TERM
FICH ALIC TEMP tfre 1e-6 FACE
poin lect 1 term HGRI 0.14
elem lect 1 26 45 50 52 60 500 1000 term R 0.10
* BFLU 2
OPTI notest noprint FSCP 0
cstab 0.9 LOG 1 BLOQ 123 LECT a_fix TERM
VFCC FCONV 6 ORDR 2 INIT MAPB 'map_9kg_1m.map' POS 0.6 0.75 10.2
CALC tini 0 nmax 30000 tfin 2.0e-3 *--------------------------------------Outputs
FIN ECRI
FICH ALIC TEMP FREQ 1
elem lect 244397 term
FICH PVTK tFREQ 1e-3 VARI ECRO CONT DEPL
train_madrid_21_vf.epx FAIL RISK
GROU 6
train OBJE LECT sall s_tri TERM
*---------------------------------------------- OBJE LECT wall doors TERM
ECHO ! comment OBJE LECT segm segm_top TERM
CAST 'train_madrid21_vf.msh' mesh OBJE LECT vol0 TERM
EROS 0.5 OBJE LECT seats TERM
*---------------------------------Problem type OBJE LECT hb1 TERM
TRID ALE RISK *--------------------------------------Options
*---------------------------------Dimensioning OPTI NOTE
DIME NF34
NALE 1 NBLE 1 csta 0.25
TERM LOG 1
*------------------------------------Geometry VFCC
GEOM FCONV 6
DKT3 s_tri ORDRE 2 OTPS 2
Q4GS sall wall seats hb1 doors RECONS 1
POUT segm segm_top LMAS 3 LQDM 3 LENE 3
CUVF vol0 KMAS 0.75 KQDM 0.75 KENE 0.75
CL3D ages CENER
TERM *-------------------------Transient calculation
*-------------------------Geometric complements CALC TINI 0 TFAI 1e-9 TEND 40.E-3
COMP EPAI 3.E-3 LECT sall s_tri TERM FIN
4.E-3 LECT seats TERM
6.E-3 LECT wall TERM
1.E-2 LECT doors TERM
1.E-1 LECT hb1 TERM
TERM
GEOP QUEL VX 0 VY 0 VZ 1
AIRE 7.6E-4 IY 80.1E-8 IZ 8.49E-8 HY
0.023 HZ 0.04 R 3.41e-2
LECT segm TERM
RECT VX 1 VY 0 VZ 0
AY 0.06 AZ 0.03
LECT segm_top TERM
EROS 0.0 LECT wall doors TERM
*-------------------------------------ALE
GRIL LAGR LECT sall s_tri wall doors seats hb1
TERM
EULE LECT vol0 TERM
*--------------------------------Material data
MATE GAZP RO 1.3 PINI 1E5 GAMMA 1.35
PREF 1e5
LECT vol0 TERM

34
Europe Direct is a service to help you find answers to your questions about the European Union
Freephone number (*): 00 800 6 7 8 9 10 11
(*) Certain mobile telephone operators do not allow access to 00 800 numbers or these calls may be billed.

A great deal of additional information on the European Union is available on the Internet.
It can be accessed through the Europa server http://europa.eu/.

How to obtain EU publications

Our priced publications are available from EU Bookshop (http://bookshop.europa.eu),

where you can place an order with the sales agent of your choice.

The Publications Office has a worldwide network of sales agents.

You can obtain their contact details by sending a fax to (352) 29 29-42758.

European Commission
EUR 26735 EN – Joint Research Centre – Institute for the Protection and Security of the Citizen

Title: Simulation of blast waves by using mapping technology in EUROPLEXUS

Authors: Martin Larcher, Folco Casadei, George Solomos

Luxembourg: Publications Office of the European Union

2014 – 32 pp. – 21.0 x 29.7 cm

EUR – Scientific and Technical Research series – ISSN 1831-9424

ISBN 978-92-79-39249-8
doi:10.2788/98310

Abstract

Finite element or finite volume simulations for the development of blast waves by using a model for the explosion of the solid
itself need very fine meshes in the explosive and in the zone around the explosive. Structures may have a long distance to the
source of the explosive. This leads often to very big meshes with many elements. The explosive is meshed often only coarse and
therefore the results are not very accurate. There are several possibilities to deal with this problem.
Large 3D calculations with a solid TNT model using a JWL-equation can be used but they are more effective when the results of
one finer mesh could be mapped in a coarser mesh after some calculation steps. When the blast wave reaches a certain
distance to the charge, the small elements inside the charge are not needed any more since the pressure ratio is decreased
strongly. These small elements results in very small time steps for the full model. The report shows the implementation of the
mapping algorithm in EUROPLEXUS and several validation tests of the method.
 LB-NA-26735-EN-N

ISBN 978-92-79-39249-8
doi:10.2788/98310