Simulation of the Effects of an Air Blast Wave

Martin Larcher

a) t=2e-5

c) t=1e-4

b) t=6e-5

d) t=1.4e-4

PUBSY JRC41337 - 2007

The Institute for the Protection and Security of the Citizen provides researchbased, systemsoriented support to EU policies so as to protect the citizen against economic and technological risk. The Institute maintains and develops its expertise and networks in information, communication, space and engineering technologies in support of its mission. The strong crossfertilisation between its nuclear and non-nuclear activities strengthens the expertise it can bring to the benefit of customers in both domains.

European Commission Joint Research Centre Institute for the Protection and Security of the Citizen Contact information Address: Martin Larcher, T.P. 480, Joint Research Centre, I-21020 Ispra, ITALY E-mail: [email protected] Tel.: +390332789004 Fax: +390332789049 http://ipsc.jrc.ec.europa.eu http://www.jrc.ec.europa.eu Legal Notice 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. 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/ JRC 41337 ISSN 1018-5593 Luxembourg: Office for Official Publications of the European Communities © European Communities, 2007 Reproduction is authorised provided the source is acknowledged

Printed in Italy

Distribution List Lechner S. Anthoine A. Casadei F. Dyngeland T. Géradin M. Giannopoulos G. Gutierrez E. Halleux J.P. Larcher M. Paffumi E. (JRC Petten) Pegon P. Solomos G. DG Tran External: Mr. Bung H. (CEA) Faucher V. (CEA) Galon P. (CEA) Kill N. (Samtech) Cheruet A. (Samtech) Potapov S. (EDF) S. Lechner

The information contained in this document may not be disseminated, copied or utilized without the written authorization of the Commission. The Commission reserves specifically its rights to apply for patents or to obtain other protection for the matter open to intellectual or industrial protection.

CONTENTS 1 2

Introduction ...................................................................................................................................5 Air Blast Waves.............................................................................................................................6 2.1 Introduction ..........................................................................................................................6 2.1.1 Detonations ......................................................................................................................6 2.1.2 Air Blast Waves ...............................................................................................................7 2.2 Literature Data .....................................................................................................................8 2.2.1 Pressure-Time Distribution ..............................................................................................8 2.2.2 Maximum / Minimum Pressure .....................................................................................10 2.2.3 Impulse...........................................................................................................................11 2.2.4 Negative Phase ...............................................................................................................11 2.2.5 Wave Form Parameter ...................................................................................................13 2.2.6 Shock Front Velocity .....................................................................................................16 2.2.7 Specific Heat Ratio ........................................................................................................16 3 Numerical Loading of a Structure with Air Blast Waves............................................................18 4 Investigations with Explosives as a Charge (Solid TNT) ...........................................................20 4.1 Modelling of the explosive ................................................................................................20 4.2 Behaviour in the Explosive ................................................................................................24 4.3 Cone with two Symmetry Axes .........................................................................................29 4.4 Cubic Charge with two Symmetry Axes............................................................................38 4.5 Spherical Charge with two Symmetry Axes ......................................................................42 4.6 Comparison between the Different Models .......................................................................46 4.6.1 Maximum Pressure.........................................................................................................46 4.6.2 Impulse...........................................................................................................................47 4.6.3 Arrival Time...................................................................................................................49 4.6.4 Positive Phase Duration .................................................................................................49 4.6.5 Comparison with results of other authors ......................................................................50 4.7 Influence of several parameters .........................................................................................50 4.7.1 Specific heat ratio (CON15) ..........................................................................................50 4.7.2 Values for γ, E0, ρ .........................................................................................................51 4.7.3 Parameters for the explosive ..........................................................................................51 4.7.4 Burn mass fraction .........................................................................................................51 5 Bubble model...............................................................................................................................52 6 Control Volume ...........................................................................................................................58 6.1 Flux between the CL3D and the fluid element ..................................................................58 6.2 Several models ...................................................................................................................58 7 Implementation of an Air Blast Loading Function .....................................................................62 7.1 Motivation ..........................................................................................................................62 7.2 Used Function ....................................................................................................................62 7.3 Implementation ..................................................................................................................62 7.4 Verification with Examples................................................................................................63 8 Mesh generation for LS-DYNA ..................................................................................................64 9 References ...................................................................................................................................65 10 Apendix ..................................................................................................................................67 10.1 EUROPLEXUS Code ........................................................................................................67 10.2 Miscellaneous code ............................................................................................................71 10.3 Sample input files...............................................................................................................74

4

1 Introduction This work is being conducted in the framework of the project RAILPROTECT, which deals with the security and safety of rail transport against terrorist attacks. The bombing threat is only considered, and focus is placed on predicting the effects of explosions in railway and metro stations and rolling stock and on assessing the vulnerability of such structures. The project is based on numerical simulations, which are carried out with the explicit Finite Element Code EUROPLEXUS that is written for the calculation of fast dynamic fluid-structure interactions. This program has been developed in a collaboration of the French Commissariat à l'Energie Atomique (CEA Saclay) and the Joint Research Centre of the European Commission (JRC Ispra).

As the aim of this project is to calculate the behaviour of structures under a loading produced by air blast waves, an indispensable starting point in this study is the ability to simulate the generation of such waves from a given quantity of explosive, and to follow their propagation through 3D spaces as they finally impinge onto the structures under consideration. The results of such numerical tests of free air blasts are presented in this report and are compared to experimental data available in the literature. In the absence of such data the EUROPLEXUS results are compared to the results of other codes, in particular to LS-DYNA, which is run in collaboration with the University of Karlsruhe. This analysis is preceded by an exposition of some basic concepts on blast wave characteristics, explosives, and a description of the equation of state adopted herein for the modelling of the detonation of a solid explosive.

5

2 Air Blast Waves 2.1 Introduction 2.1.1 Detonations Explosions can be distinguished in detonations and deflagrations. The difference between detonations and deflagrations is the velocity of the reaction zone in the explosive. Deflagrations have a slower reaction zone than the sound speed. Examples for deflagrations are the burning of gas-airmixtures and slow explosives like gun powder. Detonations have a faster reaction zone than the sound speed. The most common explosives react with detonations. To compare different explosives the TNT equivalent can be used. The TNT equivalent is a method for quantifying the energy released in the detonation of an explosive substance, by comparing it to that of an equal quantity of TNT. It is known that 1 kg TNT releases the energy of 4.520x106 J. The TNT equivalent is available for standard explosives and for some of them it is summarized in Table 1. Explosive

Mass Specific

TNT

energy [kJ/kg]

Equivalent

TNT

4520

1

Torpex

7540

1.667

Semtex 1A

4980

1.102

C4

6057

1.34

Table 1: TNT equivalent for different explosives

The effects of an explosion can be distinguished in three ranges: •

Contact detonation: The explosive is in contact with the loaded material. The load-time function depends on the loaded material, which, in most cases, is destroyed. Occurrences are the blasting of concrete (demolition etc.) or terrorist attacks where the explosive is located directly on the structure.

•

Near zone of the explosion: In most cases he material is also directly damaged like in the contact zone.

6

•

Far zone. The blast wave resulting from the detonation dominates the effects on humans and structures.

The size of all these zones depends on the quantity of the explosive charge. Additional parameters for a detonation, depending on the size of the explosive, can be defined. For example, the radius in which debris from the explosion (not from the blast wave) are possible is given by Kinney [15] as 1

r = 45 W 3

(1)

where, r is expressed in m and W is the TNT equivalent of the explosive in kg.

2.1.2 Air Blast Waves The pressure that arrives at a certain point depends on the distance and on the size of the explosive. p

pmax

pmin t ta

td

p0

tn

Figure 1: Pressure-time curve for a free air blast wave

The main characteristics of the development of this pressure wave are the following: -

The arrival time ta of the shock wave to the point under consideration. This includes the time of the detonation wave to propagate through the explosive charge.

-

The peak overpressure pmax . The pressure attains its maximum very fast (extremely short rise-time), and then starts decreasing until it reaches the reference pressure po (in most cases the normal atmospheric pressure).

-

The positive phase duration td, which is the time for reaching the reference pressure. After this point the pressure drops below the reference pressure until the maximum negative

7

pressure pmin. The duration of the negative phase is denoted as tn. -

The incident overpressure impulse, which is the integral of the overpressure curve over the positive phase td.

The idealised (free air blast) form of the pressure wave of Figure 1 can be greatly altered by the morphology of the medium encountered along its propagation. For instance, peak pressure can be increased up to 8 times if the wave is reflected on a rigid obstacle. The effects of the reflection depend on the geometry, the size and the angle of incidence. By setting γ = 1.4 (ratio of specific heats of air), it can be shown that the reflected overpressure pr is ⎡ 7 p + 4 pmax ⎤ pr = 2 pmax ⎢ 0 ⎥ ⎣ 7 p0 + pmax ⎦

(2)

All parameters of the pressure time curve are normally written in terms of a scaled distance

Z=

3

d W

(3)

where W is the mass of the explosive charge and d the distance to the centre of the charge.

2.2 Literature Data 2.2.1 Pressure-Time Distribution There are available in the literature several pressure-time-curves for different kinds of explosions. The effects of nuclear explosions here should be disregarded. The pressure at a known point can be described by the modified Friedlander equation (from Baker [2]) and depends on the time t from the arrival of the pressure wave at this point ( t = t0 − ta ) ⎛ t ⎞ p (t ) = p0 + pmax ⎜ 1 − ⎟ ⎝ td ⎠

−

bt td

(4)

The other parameters involved are the atmospheric pressure p0, the maximum overpressure pmax and the duration of the positive pressure td. The parameter b describes the decay of the curve. It can be calculated with a known minimum pressure after the positive phase. Alternatively, the parameter b can be calculated with the knowledge of the impulse. This will be done in chapter 2.2.5. All parameters for the pressure-time curve can be taken from different diagrams and equations (Baker [2], Kinney [15], Kingery [14], see e.g. Figure 2).

8

Figure 2: Model of Kingery [14] with scaled distances

9

2.2.2 Maximum / Minimum Pressure Kingery [14] developed in 1984 curves for the description of the different air blast parameters by using a rich body of experimental data, which had been properly homogenised. The parameters are presented in double logarithmic diagrams with the scaled distance Z as abscissa, but are also available as polynomial equations. These diagrams and equations enjoy the greatest overall acceptance and are widely used as reference by most researchers. The parameters are also implemented in different computer programs that can be used for the calculation of air blast wave values. e.g. they are implemented in ConWep – a program developed from the US-Army that calculates conventional weapons effects. The same curves are also used for an easy air blast load model (*LOAD_BLAST) in LS-DYNA. Also in [14] curves are provided for reflection effects (surface burst of hemispherical charges) and free air conditions (spherical charge). Another equation has been proposed by Kinney [15], in which the overpressure-distance relation for chemical explosions can be written as

pmax = p0

⎡ ⎛ Z ⎞2 ⎤ 808 ⎢1 + ⎜ ⎟ ⎥ ⎣⎢ ⎝ 4.5 ⎠ ⎦⎥ ⎛ Z ⎞ 1+ ⎜ ⎟ ⎝ 0.048 ⎠

2

⎛ Z ⎞ 1+ ⎜ ⎟ ⎝ 0.32 ⎠

2

⎛ Z ⎞ 1+ ⎜ ⎟ ⎝ 1.35 ⎠

2

Figure 3 shows the small differences between the two models.

Figure 3: Difference of the model of Kingery and the model of Kinney with 1kg TNT

10

(5)

2.2.3 Impulse The impulse of the air blast wave has a big influence on the response of the structures. The impulse is defined here as the area under the pressure time curve with the unit of pressure*sec. The impulse can be calculated with (Kinney [15]) I=

0.067 1 + ( Z / 0.23) 4

(6)

Z 2 3 1 + ( Z /1.55) 4

Another possibility is the polynomial equation of Kingery [14]. The comparison of the impulse resulting from both equations shows that the equation of Kinney simplifies the curve of the impulse between a scaled distance of 0.5 and 1.5 m/kg1/3. 800

impulse [Pa sec]

600 Kinney Kingery 400

200

0 0

0.5

1

1.5

2

2.5

3

1/3

scaled distance [m/kg ]

Figure 4: Different equations for the impulse (Kinney [15] and Kingery [14])

2.2.4 Negative Phase Detonations produce an overpressure peak, and afterwards the pressure decreases and drops below the reference pressure (generally the atmospheric pressure). The influence of the so-called negative phase depends on the scaled distance. For scaled distances Z larger than 20 and especially for Z larger than 50 the influence of the negative phase can not always be neglected. The size of the positive impulse and of the negative impulse is then nearly the same. If the structure can react

11

successfully to the positive pressure but is more sensitive to a negative pressure, failure of parts of the structure can result from this negative pressure phase (see Krauthammer [16]). However, in several cases the negative phase is neglected e.g. in the air blast function of the CONWEP-Code. Smith [22] presents the following equation to calculate the value of the negative pressure

pmin =

0.35 105 Pa Z

for Z > 1.6

(7)

The duration time of the negative pressure pmin can be calculated with

tn = 0.00125W 1/ 3 [sec]

(8)

Another possibility to get these parameters is a diagram (see Figure 5) in Krauthammer [16]. By using this diagram the limitation of equation (7) can be overcome by assuming 0.35 5 10 Pa Z = 104 Pa

pmin =

for Z > 3.5

pmin

for Z < 3.5

(9)

The duration of the negative phase in the diagram of Krauthammer can be described with the following function tneg = 0.0104 ⋅ W 1/ 3 [sec] for Z < 0.3 1/ 3 tneg = ( 0.003125 ⋅ log( Z ) + 0.01201) ⋅ W [sec] for Z ≥ 0.3 ∧ Z ≤ 1.9 1/ 3 tneg = 0.0139 ⋅ W [sec] for Z > 1.9

(10)

The difference between the duration of the negative phase of Smith and Krauthammer is observable. Smith quotes a relatively old paper from Brode [5] that uses one of the first computational machines to calculate the behaviour in the negative phase. This paper from Brode doesn’t describe the used units. Therefore the value for the duration of the negative phase should be used with caution. The value for the duration of the negative phase from Smith seems too small. In contrast the value from Krauthammer seems to be sometimes too big. In the absence of other results in the literature the values of Krauthammer will be used in this work.

12

Figure 5: Different parameters for the negative phase (see Krauthammer [16])

2.2.5 Wave Form Parameter The decay or form parameter b in the Friedlander equation (4) describes the decay of the pressuretime curve. The Friedlander equation has the parameters pmax, td and b. pmax and td can be readily found as explained before. There are several possibilities to calculate the decay parameter b by using another known value of the pressure-time curve: 1. Using the minimal pressure in the negative phase. Then, as it will be shown, the impulse of the positive phase is not accurate. 2. Using the impulse of the positive phase. Then, as it will be shown, the minimal pressure in the negative phase is not accurate. An additional equation for the negative phase should be used to avoid a smaller underpressure than the atmospheric pressure. Kinney [14] and Baker [3] calculate the parameter b by using the impulse of the positive phase. They use different equations for the pressure, for the impulse and for the duration of the positive phase. Therefore, the results for the parameter b differ (see Figure 6). Both methods, described above, for the calculation of the decay parameter b should be used here to see the difference between the results. The Friedlander equation is too complex to solve

13

analytically, and a program written in C++ can be used for the approximation. The listing is shown in the appendix. At the first step the negative pressure with the values from Kingery are used. The results for b differ from the function of Kinney [13] and Baker [3] (see Figure 6). The comparison of the resulted impulses (see Figure 7) shows that the parameter b calculated with the minimal pressure in the negative phase gets a too small positive impulse and should not be used. 10 Kinney Baker

8

Parameter b by using pmin Parameter b by using the impulse of the postive phase

b [-]

6

4

2

0 0

5

10

15

20

25

30

1/3

Z [m/kg ]

Figure 6: Decay parameter b – different methods

14

35

40

250 Impulse from Kingery Impulse with b from Kinney Impulse with b from Baker Impulse with b by using pmin

Impulse [Pa sec]

200

150

100

50

0 0

2

4

6

8

10

12

14

16

18

20

1/3

Scaled distance Z [m/kg ]

Figure 7: Decay parameter b – resulting impulse in comparison with the impulse from Kingery

Therefore, the parameter b is next calculated using the impulse of the positive phase. Then, the resulting curve of b is similar to the curves of Kinney and Baker. The exponential trend line given by Excel has the following equation

b = 5.2777 ⋅ Z −1.1975

(11)

The pressure time curve that is built with this b doesn’t fulfil the minimal pressure in the negative phase. Sometimes the pressure p is smaller than the atmospheric pressure. This results in an impossible state of the air. Therefore, the approximation of the negative phase is done with a bilinear curve shown in equation (12) and Figure 8 by using the values of the negative phase shown in section 2.2.4.

⎛ t ⎞ p = p0 + pmax ⎜1 − ⎟ ⎝ td ⎠ 2p p = p0 − n (t − td ) tn p = p0 − p = p0

−

bt td

2 pn (td + tn − t ) tn

for t < td for t > td for t > td +

∧ t < td +

tn ∧ t < td + tn 2

for t > td + tn

15

tn 2

(12)

Figure 8: Pressure-time curve for a free air blast wave – approximation of the negative phase

2.2.6 Shock Front Velocity The arrival time of the shock front at different points can be used to calculate the velocity of the shock front. With the knowledge of this velocity the pressure can be obtained with the RankineHugoniot relationship. Kingery [14] calculates also the shock front velocity depending on the pressure as 1/ 2

⎛ γ + 1 pmax ⎞ u = c0 ⎜ 1 + ⎟ 2γ p0 ⎠ ⎝

(13)

The parameter γ (ratio of specific heats of air) depends also on the overpressure and can be taken from a table in [13]; c0 is the sound velocity in air (331 m/sec); pmax is the peak overpressure and p0 is the atmospheric pressure (101.3 kPa).

2.2.7 Specific Heat Ratio The specific heat ratio γ is defined as

γ=

cp cv

(14)

with cp being the specific heat at constant pressure and cv the specific heat at constant volume. Both

16

the specific heat ratio and the speed of sound depend on the temperature, the pressure, the humidity, and the CO2 concentration. Kingery [14] defines the variation of the specific heat ratio with a range of 1.402 to 1.176.

17

3 Numerical Loading of a Structure with Air Blast Waves There are several ways of numerical modelling in order to load a structure with an air blast wave. These methods differ in the number of used elements and with them in the calculation time. •

Model with the mechanical modelling of the explosive (JWL-equation (15)). A fine mesh is essential to get realistic results. The size of the element in the range around the explosive should be approximately 1 mm. These calculations are very expensive. To reduce the computation time partitioning can be used. This method reduces the calculation time for models with a large variation of element sizes.

•

The method proposed by Clutter [9] is also a solid TNT model and uses only one element for the explosive. This is possible by different not specified methods in combination with the Becker-Kistiakowsky-Wilson EOS for the explosive.

•

1D to 3D. This modelling is proposed in [4] and is also a solid TNT model. A 1D calculation is used until the wave reaches a surface. Then the values of the density, energy, velocity and pressure are mapped into the 3D mesh. Rose [20] maps the 1D model to 2D when the wave arrives the first surface and maps the 2D model to 3D when the wave arrives a second surface with another direction. EUROPLEXUS allows the implementation of this method. The method should also be validated with a calculation of the first model. The model is a mixture of the first and the third model. The calculation time should be larger than for the second model.

•

Model with a compressed bubble. The pressure-time function resulting from a compressed bubble can not easily match the curve of an air blast wave. The size of the compression can be calibrated with the maximum pressure or the impulse. The calculation time is smaller than for the first model.

•

Control volume. A volume around the explosive is loaded by a pressure-time curve. This pressure time curve can be calculated with a model based on the modelling of the explosive. Alternatively, the well known curves from Kingery [14] can be used. This method should be validated through comparisons with calculations of the first model. The computation time is in the range of the second model.

•

Load-time function. This is only usable for an estimation of the behaviour of a structure loaded by an air blast wave. The structure is loaded by a load-time function built with the pressure-time function e.g. from Kingery [14]. This function is implemented in EUROPLEXUS (see chapter 7). The calculation is relatively inexpensive. Alternatively, the pressure-time function can be 18

determined with a fluid pre-calculation with fixed boundaries for the structure. The structure is then loaded by the pressures resulting from this fluid calculation. The choice among these methods depends on the scope of the analysis and on further investigations about their advantages and shortcomings. Figure 9 shows different models for the simulation of an air blast wave.

Compressed bubble

Solid TNT

Load on a control volume

Load-time function

Figure 9: Several models for air blast wave simulations

19

4 Investigations with Explosives as a Charge (Solid TNT) The aim of the RAILPROTECT project is to contribute to alleviating the vulnerability of Europe's passenger transport infrastructures. The effects of a terrorist attack should be simulated numerically, and for a numerical investigation the knowledge of the loading of the structures is necessary. There are different approaches and possibilities for the calculation of a detonation inside buildings, as discussed in chapter 3. A detonation releases a large amount of energy in a very short time. This results in an air shockwave which is spread outwards from the charge. Then, the air blast wave reaches the structure, which, depending on the size of the charge and on the distance, will respond to this wave loading. A calculation of the behaviour of the air blast wave requires the knowledge of the behaviour of the explosive and of the air around the explosive. The results of the numerical investigation can be compared for the validity of the calculations with existing experimental-analytical data. As will be shown in this chapter, the experimental-analytical results of Kingery [14] will constitute the basis for these comparisons.

4.1 Modelling of the explosive The explosive for the numerical investigation can be built up e.g. with the Jones–Wilkins–Lee (JWL)-equation. This equation of state (EOS) is widely used because of its simplicity and due to the fact that most high explosives are well modelled by this equation. According to it, the value of pressure is given as ⎛ ⎛ ω ⎞ − R1V ω ⎞ − R2V E pEOS = A ⎜1 − + B ⎜1 − +ω ⎟e ⎟e V ⎝ R1V ⎠ ⎝ R2V ⎠

(15)

In this equation A, B, R1, R2 and ω are the model parameters, V is the ratio ρsol/ρ, where ρ=current density and ρsol =density of solid explosive, and E is the internal energy per unit volume of the explosive. It is noted that E=ρsol eint, where eint is the current internal energy per unit mass. The parameters of this equation for most explosives are shown in Dobratz [10]. Different authors use slightly differing parameter values for this equation, as shown in Table 2. Note that the equation will be reduced only to its last term if the solid explosive is exhausted and the resulting gases fully expanded. The last term of equation (15) is the EOS of an ideal gas that can be used e.g. for the air. pEOS = ω

20

E V

(16)

From this asymptotic form it can also be concluded that ω=γ-1 (γ=ratio of specific heats). Parameter

Description

ref.[6]

AUTODYN ref.[21]

parameters used

- manual

for the air

A (Pa)

3.738e11

3.7377e11

3.712e11

B (Pa)

3.747e9

3.7471e9

3.21e9

R1

4.15

4.15

4.15

R2

0.90

0.90

0.95

1630

1630

1630

1.3

3.68e6

3.68e6

4.29e6

2.1978E5 1.30

ρsol (kg/m3)

density

eint (J/kg)

current internal energy per unit mass

γ

specific heat ratio

1.35

1.35

1.30

vdet (m/sec)

detonation speed

6930

6930

6930

Table 2: Parameters for the JWL equation for TNT

For the air the same EOS will be used without a detonation and different starting density and internal energy. By ignoring the explosion, the last part of the JWL-equation will prevail and therefore, an ideal gas will be used. Different FE-Codes smear the detonation front over different time steps. This procedure is called burn fraction and its motivation is to control the release of the chemical energy for the simulation. The effects of the combustion on the pressure can be considered with this formula p = pEOS min (1, max ( F1 , F2 ) )

(17)

The burn mass functions F1 and F2 are computed by (see LS-DYNA and [18]) ⎧ 2 ( t − t1 ) d ⋅ Ae ,max ⎪ F1 = ⎨ 3ν e ⎪ 0 ⎩ F2 =

1−V 1 − VCJ

if t > t1

(18)

if t ≤ t1 (19)

where t1 is the ignition time of the observed element (calculated with the detonation velocity d), Ae,max is the maximum surface area and νe the volume of the element. V is the actual specific volume and VCJ the specific volume at the Chapman-Jouguet-pressure. The Chapman-Jouguet

21

pressure is reached if the sonic velocity of the reaction gases reaches the detonation velocity. The volume at the Chapman-Jouguet-point is VCJ = 1 −

PCJ ρ0 d 2

(20)

The term F1 of the burn mass function intends to spread the burn front over several elements. The second term should control the releasing of the energy. Interestingly, MSC-Dytran uses only the term (19), whereas ABAQUS uses only one burn mass function ⎧ ( t − t1 ) d ⎪ Fb = ⎨ BS ⋅ le ⎪ 0 ⎩

if t > t1

(21)

if t ≤ t1

In this formula BS is the constant that controls the width of the burn wave (set to a value of 2.5) and le is the element length. This function is very similar to (18). A 3D model of explosive is used to show the influence of the different burn mass fractions (Figure 12, explosive4), whose different components can be compared in Figure 10 and Figure 11. The influence of the term (21) is visible. The slope of the pressure peak decreases with Fb, and the arrival of the pressure peak is later. The calculations with EUROPLEXUS do not show an influence of the second term (F2). Both curves are almost identical. LS-DYNA

uses

the

functions

(18)

and

(19)

for

calculations

with

the

MAT_HIGH_EXPLOSIVE_BURN model. The input syntax allows calculations with both functions (beta=0) or only with the function (19) (beta=1). The model in LS-DYNA was built with hexahedral with an element size of 0.001 m. The results of a 3D model show that the burn mass fraction in LS-DYNA reduces also the slope of the pressure peak. Nevertheless, the difference between the calculations with beta=0 and beta=1 is negligible. The differences of the calculations between EUROPLEXUS and LS-DYNA are the smaller peak in LS-DYNA and the higher pressure values behind the peak (here for a distance less than 0.04 m). The impulses are nearly the same (LSDYNA 2.07 108 Pa sec, EUROPLEXUS 2.00 108 Pa sec). Therefore, the burn mass fraction is for the future work implemented only with the equation (21).

22

2.0E+10 without burn mass fraction burn mass fraction Fb burn mass fraction Fb and F2 LS-DYNA beta=0 LS-DYNA beta=1

Pressure [Pa]

1.6E+10

1.2E+10

8.0E+09

4.0E+09

0.0E+00 0

0.01

0.02

0.03

0.04

0.05

0.06

Distance [m]

Figure 10: Influence of the burn mass fraction in EUROPLEXUS (t=7 10-6 sec)

2.0E+10 EUROPLEXUS bmf=0

1.6E+10

EUROPLEXUS bmf=2.5 LS-DYNA beta=0

Pressure [Pa]

LS-DYNA beta=1

1.2E+10

8.0E+09

4.0E+09

0.0E+00 0.04

0.05

0.06

0.07

0.08

0.09

0.1

Distance [m]

Figure 11: Influence of the burn mass fraction in EUROPLEXUS (t=1.4 10-5 sec)

Another EOS for explosive is the “Ignition and Growth Reactive Model” based on Tarver et al. [23]. This model can also be used for the burning of propellant (deflagrations). 23

4.2 Behaviour in the Explosive The blast behaviour in the air is affected by the development of the pressure in the explosive. Therefore, it should be investigated, whether the numerical simulation can sufficiently represent the behaviour of the explosive. The numerical model for the explosive calculates the pressure with the JWL-equation (15). An accurate model for the development of the detonation front is used (See [1]). The detonation starts at the initiation point, and the detonation front is moved with the given velocity of the detonation. An element detonates if the detonation front reaches this element. From this time the JWL-equation will be used for this element. A spherical TNT charge of volume of 8000 cm3, i.e. with a radius of 0.124 m is considered. To control the behaviour of the pressure in the explosive a conical model is used (similar to the model of chapter 4.3, see Figure 12).

length

opening x

d pyra d ex,in

d ex,end

Figure 12: 3D simplified model for the behaviour in the explosive

The models use hexahedrons and a pyramid for the top of the model. The FSR-condition is used for all surfaces of the model. The difference between the meshes, listed in Table 3, is the refinement.

24

Case

opening

dpyra

dex,in

dex,end

Number of elements

explosive1

Eulerian

2e-2

0.001

0.001

0.01

33

explosive2

Eulerian

1e-2

0.0005

0.0005

0.005

65

explosive3

Eulerian

4e-3

0.0005

0.0002

0.002

159

explosive4

Eulerian

2e-3

0.0002

0.0001

0.001

318

explosive5

Eulerian

2e-3

0.0002

0.0001

0.0005

499

explosive6

Eulerian

2e-3

0.0002

0.0001

0.0002

859

explosive7

ALE

2e-3

0.0002

0.0001

0.0005

499

Table 3: Comparison of different models for the explosive

The pressures are here analysed in intervals of 1µsec. The procedure to get the values in space (SCOURBE command in EUROPLEXUS) uses the averaged elemental pressure values for the nodal values. This results in a half value of the pressure at the node between an element inside and an element outside of the detonation zone (see Figure 13). The vertical lines in this figure correspond to the distances of the elements that are detonated at a certain time step (can be get from the listing). The detonation front can be identified by the steep increasing of the pressure. Therefore, it is important to consider the value of the last element of the detonation front and not the value of the last node.

Figure 13: Distribution of pressures over distance (explosive1)

The maximum pressure of the detonation is increasing with the distance from the initial detonation point. The Chapman-Jouguet-pressure is the maximum experimental resulted pressure. The 25

parameters of the JWL-equation should represent this limitation (see Shin [21]). Shin shows the influence of the discretisation with a 1D model, where a finer mesh reproduces larger pressures. The results with the models explosive1 to explosive4 show the same dependency of the element sizes (Figure 14). The Chapman-Jouguet-pressure seems to be the limit in the convergence study. 2.40E+10

Pressure [Pa]

2.00E+10

1.60E+10 t=7e-6 t=1.4e-6 Shin t=7e-6 Shin t=1.4e-5 Chapman-Jouguet-pressure

1.20E+10

8.00E+09

4.00E+09

0.00E+00 0

200

400

600

800

1000

1200

1400

1600

1800

2000

Number of elements

Figure 14: Maximum pressures in the explosive depending on the number of elements

The rather unexpected behaviour of the models with 499 and 859 elements has to be clarified. The reason could lie in the location of the detonation point in conjunction with the location of the integration points. Figure 15 shows the pressure depending on the distance from the initiation point. The curve does not differ for models with fine discretisation. It is observed that the pressures behind the pressure peak in the model explosive6 are definitely smaller than the results of Shin. Therefore, the area under the pressure-distance curve is also smaller. This area reaches only 67 % of the area of Shin. The reason could be the missing burn fraction (see chapter 4.1).

26

2.4E+10

explosive1, t=7e-6 explosive1, t=1.4e-5 explosive4, t=7e-6 explosive4, t=1.4e-5 explosive6, t=1.4e-5 Shin t=1.4e-5 Shin t=7e-6

2.0E+10

Pressure [Pa]

1.6E+10

Area Shin 5.54e8 Area explosive6 3.71e8 (between 0.051 and 0.97)

1.2E+10

8.0E+09

4.0E+09

0.0E+00 0

0.02

0.04

0.06

0.08

0.1

0.12

Distance [m]

Figure 15: Pressure distance curve in the explosive – comparison with Shin [21]

The model explosive7 uses an ALE mesh instead of an Eulerian mesh. The results are almost identical for the maximum pressure as well as for the bending of the curve. A calculation with a Lagrangian mesh fails. The next question is how the pressures are developed at the border of the explosive. At this time experimental results for this region are not available. Here, the numerical results of EUROPLEXUS could be compared with numerical results of LS-DYNA. For this calculation a conical model is used (see chapter 4.3, CON20, Table 4). In LS-DYNA the boundary condition like the FSR condition is relatively complex (BOUNDARY_SPC). To avoid an influence of the usage of this boundary condition the model in LS-DYNA is built as a full 3D model. The LS-DYNA model has a length of 15 cm, a height of 5 cm and a width of 5 cm. The boundaries are built as fixed. Therefore, after the reflection of the wave at the boundaries the pressure patterns are no longer spherical. The element size in LS-DYNA is chosen as 1 mm with an ALE multi-material formulation (see Figure 16).

27

Figure 16: Model in LS-DYNA

The calculations show that the numerical results of both programs are nearly the same inside the explosive as well as in the air near the explosive (see Figure 17). The calculation with EUROPLEXUS is done without the burn mass fraction. 2.0E+10 EUROPLEXUS x=0.105 EUROPLEXUS x=0.125 EUROPLEXUS x=0.136 LS-DYNA x=0.105 LS-DYNA x=0.125 LS-DYNA x=0.136

Pressure [Pa]

1.5E+10

1.0E+10

5.0E+09

Reflection Reflection

0.0E+00 1.4E-5

1.6E-5

1.8E-5

2.0E-5

2.2E-5

2.4E-5

2.6E-5

2.8E-5

3.0E-5

Time [sec]

Figure 17: Pressure time curve in the explosive and in the air near the explosive; length of the explosive = 0.124 m

28

EUROPLEXUS x=0.105 EUROPLEXUS x=0.125 EUROPLEXUS x=0.136 LS-DYNA x=0.105 LS-DYNA x=0.125 LS-DYNA x=0.136

5.0E+09

Pressure [Pa]

4.0E+09

Reflection

3.0E+09

2.0E+09 Reflection 1.0E+09

0.0E+00 1.4E-5

Reflection

1.6E-5

1.8E-5

2.0E-5

2.2E-5

2.4E-5

2.6E-5

2.8E-5

3.0E-5

Time [sec]

Figure 18: Detail of the pressure time curve

4.3 Cone with two Symmetry Axes A two dimensional model does not consider the behaviour of the fluid in the third direction. Therefore, a three dimensional model for the simulation of the detonation should be used. The calculation costs are minimized by using a conical model (pyramid, see Figure 19). This model can be built with one hexahedral element in the radial direction and symmetry axes at the edges of the elements. The element on the top is then a pyramid. To get elements with similar geometry the length of the elements should be increased with an increasing distance from the centre. The size of the opening angle depends also on the size of the elements. This angle should be chosen so that the aspect ratio of the elements is not too large. Normally, the time step size for this model is relatively small because of the very small element on the top of the cone.

29

length Explosive opening x

d pyra d ex,in d ex,end d air,in

d air,end

Figure 19: 3D simplified model with hexahedrons or tetrahedrons

Alternatively four tetrahedral elements can be used instead of one hexahedral element. This results in a higher number of elements. The meshing can be done with scripts that convert the hexahedrons in tetrahedrons. These scripts are presented in [7]. The script pxhex2te converts the hexahedrons in tetrahedrons; the script pxqua2tr converts the quadrangles to triangles (necessary for the surfaces). For the simulation a model is used which has in most cases a TNT charge of a volume of 8000 cm3, accordingly a cube of 20 x 20 x 20 cm or a sphere with radius of 12.4 cm. The mass of this charge is 12.8 kg. An ALE-calculation is used for all further investigations. The explosive is build as an Eulerian mesh, the air is build as an ALE mesh. The symmetry at the surfaces can be considered by defining symmetry planes with the CONT SPLA AUTO command. This command defines symmetry conditions orthogonal to the surface of the defined nodes. Another possibility is the definition with the FSR-command as a sliding surface. Both methods give the same results. The calculations performed with the conical model are summarized in Table 4 and shortly described hereafter.

30

Case

d pyr/d ex,in /dex,end/d air,in/d air,end

opening/

Charge

Element type

tend

length CON1

0.05/0.05/0.08/0.2/0.2

0.2/4.0

12.8

FL34 (NF34)

4.0e-3

CON2

0.01/0.02/0.02/0.02/0.05

0.2/4.0

12.8

FL34 (NF34)

4.0e-3

CON3

0.008/0.01/0.01/0.01/0.05

0.1/4.0

12.8

FL34 (NF34)

4.0e-3

CON4

0.008/0.005/0.005/0.005/0.02

0.03/1.3

12.8

FL34 (NF34)

4.7e-4

CON5

0.008/0.005/0.005/0.005/0.02

0.03/1.3

12.8

FL35/FL38

4.7e-4

CON6

0.002/0.002/0.002/0.002/0.01

0.01/1.3

12.8

FL35/FL38

4.7e-4

CON7

0.001/0.001/0.001/0.001/0.005

0.01/1.3

12.8

FL35/FL38

4.7e-4

CON8

0.001/0.001/0.001/0.001/0.01

0.05/3.0

12.8

FL35/FL38

1.5e-3

CON9

0.001/0.001/0.001/0.001/0.01

0.05/3.0

1.0

FL35/FL38

1.5e-3

CON10

0.001/0.001/0.001/0.001/0.005

0.05/3.0

1.0

FL35/FL38

1.5e-3

CON11

0.0005/0.0005/0.0005/0.0005/0.01

0.025/1.5

1.0

FL35/FL38

1.5e-3

CON12

0.005/0.001/0.001/0.001/0.01

0.1/3.0

1.0

FL35/FL38

1.5e-3

CON13

0.005/0.001/0.001/0.001/0.04

0.2/3.0

1.0

FL35/FL38

1.5e-3

CON14

0.005/0.001/0.001/0.001/0.02

0.2/1.5

1.0

FL35/FL38

1.5e-3

CON20

0.0002/0.0001/0.001/0.001/0.005;

6.5e-3/0.4

1.0

FL35/FL38

7.0e-5

Table 4: Calculations with a conical mesh

CON1

This calculation uses the modified tetrahedral elements (see [7]) also for the tip of the model. The mesh is relatively coarse. The development of the air pressure wave is presented in Figure 20. The pressure-time curve for a distance of 1 m shows the increase of the pressure until a value of 1.83 106 Pa (time step size for evaluation 2 10-6 sec) which is equivalent to 18.3 times the atmospheric pressure. The resulting pressure does not depend on the time step size for the evaluation. By choosing every calculated time step for the evaluation (1.75 10-7 sec) the maximum pressure is also 1.830958 106 Pa. After the pressure peak the pressure is decreasing up to 232.2 Pa. This value of the “negative” pressure is relatively high. The pressure-time curve for other distances follows also these trends.

31

2.0E+06

1.5E+06

1.0E+06

1

1

5.0E+05 2

3

3

2

1

0.0

3 2

2

1

1

1

Time [s] -5.0E+05

0.0 5.0E-04 FLUID_3D_SF7 -1- Distance = 1m

1.0E-03

1.5E-03

-2- Distance = 2m

2.0E-03

2.5E-03

3.0E-03

3.5E-03

4.0E-03

EUROPLEXUS 14 JUNE 2007 DRAWING 1

-3- Distance = 3m

Figure 20: Development of the pressure depending on the distance to the charge

CON2

This model uses smaller tetrahedral elements than the mesh CON1. The finer mesh results in a steeper air blast wave at a distance of 1 m as well as in a distance of 2 m. The steeper wave causes also a higher pressure (See Figure 21). Figure 22 shows the maximum pressure depending on the distance to the charge. These curves can be obtained by using a CAST3M macro that selects all nodes lying near a line. The macro pxpdroi1 (See [7]) requires a tolerance in which the nodes are lying. Then, with the macro pxordpoi the selected nodes are ordered. The line which is used to select the nodes is one of the edges of the CON models. With EUROPLEXUS it is possible to get out the data depending on the space variable. This can be done at a certain time step. From the results of different time steps the maximal values, the impulse and

the

positive

phase

duration

can

be

extracted

with

the

FORTRAN-executable

AIRBLASTRESULT (see Appendix). A comparison of the different discretisations (see Figure 22) with the results of Kingery [14] shows that a finer tetrahedral mesh produces smaller maximum pressures (disregarding the model CON1). Except for distances less than 25 cm, these pressures are smaller than the experimental pressures. The experience shows that Eulerian meshes normally react smoother than the reality.

32

Figure 21: Comparison of a coarse with a fine mesh – pressure time curve

CON3 and CON4

These calculations use finer meshes. To reduce the computation time the length of model CON4 (see Figure 19) is limited to 1.3 m. For the results of the models CON1 to CON4 it is important to note that the definition of the pressure at a certain point is arguable. There are several tetrahedrons at the same distance with different orientations, and these elements have different pressures. Thus, the results with the tetrahedrons is a field that requires further work, as probably, tetrahedrons will be used for the calculation of the air inside the structures. For the employment of the tetrahedrons in a productive framework, these elements should show safer and more reliable results. CON5

This calculation uses hexahedral elements instead of the tetrahedrons. The mesh is relatively coarse. In comparison to the tetrahedrons the maximum pressures versus the distance show a relatively good correlation with the analytical results. CON6 and CON7

Both models use finer meshes than the model CON5. The behaviour of models with hexahedral elements regarding the dependency on the element size shows the reverse trends to those of the tetrahedral elements (See also Figure 22). A finer mesh results in a higher maximum pressure. This should be the normal behaviour of a refinement of the elements. The results of the model CON7 show the best correlation, even though the differences 33

between the experimental and the numerical results are still quite big. Further work has to be done to check this discrepancy. 1E+8

CON1 ▲ CON2 ▲ CON3 ▲ CON4 ▲ CON5 ■ CON6 ■ CON7 ■ Kingery

Max. pressure [Pa]

8E+7

6E+7

4E+7

2E+7

0E+0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Distance [m]

Figure 22: Comparison of max. pressure – distance relationships from different conical models.

CON8

This calculation is done with a length of 3.0 m. The element sizes are nearly the same as in model CON7. In comparison to the results of model CON7 the difference is small (see Figure 23). 1E+7

Kingery CON7 CON8 ■ CON9 ■ 1kg CON10 CON11 CON12

Max. pressure [Pa]

8E+6

6E+6

4E+6

2E+6

0E+0 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Scaled distance [m/kg1/3]

Figure 23: Comparison of peak pressures – scaled distance relationships from different conical models.

34

CON9

This model has approximately the element sizes of the mesh CON8 but has a charge of only 1 kg. The resulting pressures in a scaled size are smaller than those of model CON8. CON10

This is a finer mesh for the air of the model CON9. However, all values are nearly the same. CON11

CON11 uses a finer mesh in the explosive and near the explosive. The pressures near the explosive are higher than in the previous models, but the difference in a larger distance is small. CON12

Figure 24 shows the pressure versus the distance at a time of 1.5 10-3 sec for the model CON9. The pressures in the first elements are approximately 10 104 times higher than expected. Controlling the aspect ratio of the elements at the tip (pyramid and hexahedrons), it is realised that the elements are not conformant. Therefore, model CON12 is tried, which uses a larger angle of the cone than the model CON9 and has also a bigger pyramid at the tip. The largest aspect ratio of the hexahedrons is now 3.0 and of the pyramid 15. The aspect ratio of the pyramid depends on the aspect ratio of the whole model ( =

length ). The pressures in the pyramid and in the three first hexahedrons are 2 ⋅ opening

larger than expected (blue rhombus points in Figure 24). To calculate the aspect ratio of an element the procedure aspectra can be used (see appendix).

Figure 24: Pressure at the top of the conical model

35

CON13

The angle of this model is increased. The aspect ratio of the hexahedrons is 1.5, of the pyramid 7.5. The results for the elements near the tip show that the pressure value in the hexahedrons is represented better if the aspect ratio is smaller. CON14

The pressures in the pyramid are decreasing with an aspect ratio of 3.75. The calculations do not show better results. CON22, CON23, CON24

These models use only cubic elements. Therefore, the element size near the explosive is similar to the other models. The size of the elements increases with the distance to the explosive. Then, the only parameter for this model is the opening of the cone (or the angle of the cone). The number of elements is relatively small but the time step size depends on the smallest element.

Case

opening/ length

Number of

Charge

elements CON22

0.1/4.0

386

12.8

CON23

0.05/4.0

711

12.8

CON24

0.02/4.0

1595

12.8

Table 5: Calculations with a conical mesh

The peak pressures show a big dependency on the element size. (see Figure 25). The influence of the element size on the impulse is much smaller (see Figure 26). Therefore, for further calculations it has to be checked if the element size is small enough.

36

4E+7

Kingery 12.8 kg SPHE9 CON8 ■ CON22 CON23 CON24

Max. pressure [Pa]

3E+7

2E+7

1E+7

0E+0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

1/3

Scaled distance [m/kg ]

Figure 25: Peak Pressure – Models CON22, CON23, CON24

2.5E+3

SPHE9 CON8 ■ CON22 CON23 CON24 Kingery 12.8kg

Impulse [Pa sec]

2.0E+3

1.5E+3

1.0E+3

5.0E+2

0.0E+0 0.1

0.2

0.2

0.3

0.3

Scaled distance [m/kg1/3]

Figure 26: Impulse – Models CON22, CON23, CON24

37

0.4

0.4

Summary

The pressures resulting from a calculation with tetrahedrons are too small. The results with hexahedrons are also small, but the difference to the analytical-experimental results is less. The difference depends on the size of the elements. Smaller elements result in higher pressure and a better correlation with the experimental results. The smallest element sizes are of the order of 1 mm. It is not possible to calculate a fluid-structure interaction problem of realistic dimensions with a mesh of this size. So, the question which kind of simulation of the loading is most appropriate to be used is important (see chapter 3).

4.4 Cubic Charge with two Symmetry Axes To check if the boundary condition of the cone model represents a full 3D model two different models are used. Both models are one eighth of the geometry with FSR conditions at the symmetry faces. Another advantage with such a model is the possibility to compare the calculations with other finite element programs like LS-DYNA. In the first model the charge and the air are built as cubes. A regular mesh with hexahedrons is applicable with the same size for all elements. The second model uses a spherical charge (see chapter 4.5). An interaction with a structure is not of interest for these calculations. Therefore, the explosive and the air are built with fluid meshes. An ALE mesh is used for the air, an Eulerian mesh is used for the explosive (see Figure 27).

Figure 27: Model CUB1 – explosive as Eulerian (red), air as ALE (grey)

The calculations performed are summarized in Table 6 and shortly described hereafter.

38

Case CUB1

Size of the model

Description

0.5 X 0.5 X 0.5 m

Charge 12.8 kg TNT,

tend

Steps

CPU [s]

Elements

4.7e-4

415

181.5

15625

4.7e-4

270

790.6

125000

4.7e-4

911

2275.9

125000

6.0e-4

291

1367.7

216000

element size 0.02 m CUB2

1.0 X 1.0 X 1.0 m

Charge 12.8 kg TNT, element size 0.02 m

CUB3

0.5 X 0.5 X 0.5 m

Charge 12.8 kg TNT, element size 0.01 m

CUB4

1.0 X 1.0 X 1.0 m

Charge 1.0 kg TNT, element size 0.0167 m

Table 6: Calculations with cubic charges

CUB1

In this calculation a small model with a coarse mesh is used. The charge has a mass of 12.8 kg. The development of the pressure can be shown in Figure 28. At t=2 10-5 the explosive has burned down except for the explosive in the corners of the cube. The pressure in the air is the atmospheric pressure of 105 Pa. After the explosion (e.g. t=6 10-5) the pressure of the gas produced from the explosive is decreasing and the pressure wave is running into the air. At t=10-4 there can be observed a discontinuous pressure along the three axes at the same distance from the detonation initiation point. The reason could be the cubical shape of the charge. At t=1.4 10-4 the pressure wave is reaching the surface (FSR condition) and is reflected.

39

a) t=2e-5

b) t=6e-5

c) t=1e-4

d) t=1.4e-4

Figure 28: CUB1 – Pressure at different time steps

CUB2

This calculation uses a bigger model but also with the coarse mesh. Figure 29 shows the maximum pressure versus the distance from the charge. The maximum pressures resulting from model CUB1 and from model CUB2 are the same in the range of the size of model CUB1. The pressures are overestimated up to a distance of 0.46 m; in relation to the size of the explosive (0.13 m) this range of overestimation is negligible. The resulting maximum pressures after the range of 0.4 m are definitely smaller than the experimental results. CUB3

The size of the model is the same as that used in model CUB1 but with a finer mesh. The pressures resulting from this calculation are bigger than in model CUB1. For CUB1 and CUB2 only the pressures along the orthogonal to the surface of the explosive are considered. The values are smaller if the pressures along the diagonal of the cube are used. The difference between both locations is very big. Therefore to compare the results in a region near the explosive with experimental results a cubic charge is not suitable. The differences between experimental and numerical results can be compared better with a spherical model.

40

CUB4

This model uses a finer mesh than model CUB2. The ending time of the calculation is longer so that the wave can be also observed in regions with a bigger distance from the charge. In comparison to the model CUB2 the pressures at CUB4 are smaller because of the smaller charge. A better method to compare is the use of the scaled distance, and in fact the values of the pressure are of the same order of magnitude in this case. There is also a big difference between the pressures along the diagonal and the pressures along the orthogonal to the charge surface. In comparison to the results with LS-DYNA (cub2) the results of EUROPLEXUS for the pressure are higher and represent the experimental values better but still not in a very satisfactory manner.

Figure 29: Maximum pressures versus distance – cubical models

41

Figure 30: Maximum pressures versus scaled distance – cubical models

4.5 Spherical Charge with two Symmetry Axes The weaknesses of using a cubical charge can be avoided by using a spherical charge. There are no problems with a spherical charge from the different ending times of the detonation at the corners and with effects resulting from the different surfaces. As mentioned initially, these first investigations disregard the interaction with the structure. So it is possible to use a spherical volume for the surrounding air, too. This allows an easier modelling of the mesh, which can be done in several ways. At this point, the first approach of a spherically symmetric regular mesh with the centre of the charge is in CAST3M neither possible in an easy way nor is it necessary in considering the behaviour of the explosive and the air. The method used instead converts a cubical model via the INCL operator to a spherical model. The resulting model seems to be relatively regular. The method can be done with the following steps:

1. Modelling of a cubic surface for the outer charge (half volume) (Figure 31a, green). 2. Modelling of a cubic surface for the inner charge. This part will be meshed with regular hexahedrons (half volume) (Figure 31a, blue). 3. Modelling of a cubic surface for the air volume around of the charge (half volume) (Figure 31a, red). 4. Projection of the outer charge surface and of the air surface with the PROJ operator to spheres (Figure 31b). 5. Filling the volumes between the surfaces with hexahedrons (Figure 31c). 42

6. Complete the volumes by adding the turned model. 7. Defining of bounding box for an one-eighth model. 8. Choosing the elements inside of the bounding box with the INCL operator (Figure 31d). 9. Defining the bounding surfaces with the POIN operator.

a) Cubic surface for outer charge (green) for inner charge (blue) and for air (red)

b) Projection to a sphere

c) Filling the volume between the surfaces

d) Elements inside the bounding box

Figure 31: Modelling of the spherical model

Alternatively an eighth of the whole spherical model can be used form the beginning. This saves the required memory. Another possibility is the modelling with two macros (pxhex2te and pxqua2tr, see [7]). With them it is possible to rotate first a line around a point to get a part of a circle. A second macro is available to turn this part of a circle around a line to get a part of a sphere. However, the resulting model seems to be more irregular then the model built with the method described before. The calculation cases performed are summarized in Table 7 and shortly described hereafter, together with the modifications introduced in the code to achieve this type of calculations. The regions with different element densities are shown in Figure 32. The charge of all models is 12.8 kg TNT, the end of the calculation depends on the size of the model.

43

Figure 32: Modelling of the spherical model – number of elements in the different models

Case

Size of the model

Length

Distance with

Number of

nel0/nel1/nel2

[m]

aspect ratio = 10

elements

SPHE1

10/10/10

1

(5.55)

1625

SPHE2

10/20/70

1

0.79

6875

SPHE3

10/20/90

1

0.62

8375

SPHE4

10/20/150

1

0.37

12875

SPHE5

10/20/200

1

0.28

16625

SPHE6

10/40/200

1

0.28

18125

SPHE7

20/40/200

1

0.56

73000

SPHE8

10/10/70

2

1.70

6125

SPHE9

20/20/140

2

1.70

12125

Table 7: Calculations with spherical charges

SPHE1

In this calculation a small model with a coarse mesh is used. The results show a relatively low maximum pressure (See Figure 33).

44

Figure 33: Maximum pressures versus distance – spherical models

SPHE2, SPHE3, SHPE4 and SPHE5

These calculations use finer meshes. The maximum pressures are increasing. The models SPE2 and SPE4 give realistic results. The models SPHE4 and SPHE5 give larger maximum pressure than the experimental values. An Eulerian calculation should normally result in smaller values than the experimental values (in contrast to a Lagrangian calculation which reacts stiffer with smaller element sizes). SPHE6

The larger values for the maximum pressure in model SPHE5 start near the charge. Therefore, the model SPHE6 tests if a smaller element size in the outer charge (See d) red range) can result in more realistic values. The overestimating of the pressure begins in this model earlier than in the model SPHE5. The size of the elements in the explosive should not be the reason for the overestimation. SPHE7

This calculation uses a finer mesh in the inner charge (See Figure 31 d) green range). With the smaller elements in the inner charge also the number of the elements in tangential direction is increased. The results are quite better. For every model with small elements there is a limit from which the cal45

culation fails. This limit depends on the radial and the tangential size of the elements. It seems that the elements are not usable when they are too thin. The pressures are overestimated if the aspect ratio is about 10. When exceeding this ratio, the pressure is also not symmetric. It appears to be higher near the edges and smaller in the diagonal of the model. This difference between the edge and the diagonal begins approximately at a distance of 0.45. The aspect ratio is there 8.0. It seems that the boundary conditions can not act with this exceptional element sizes. Respecting this, also the values for the calculation SPHE3 would not be applicable over a distance of 0.60 m. This has to be considered for further calculations with the conical model and the spherical model. The aspect ratio reaches the critical value 10 at a certain distance. This distance is calculated in Table 7. SPHE8

This model is a coarse spherical model with a size of 2.0 m. The results follow the results of the models with a size of 1 m. The pressure is approximately half the experimental value. SPHE9

This calculation uses a finer mesh than model SHPE8. The pressures are a little bit larger than the pressures in model SHPE8.

4.6 Comparison between the Different Models 4.6.1 Maximum Pressure The comparison will be done with the fine meshes of the different models summarized in Table 8. Case

Element size

Charge

Size of the model

[kg TNT]

[m]

Steps

CPU [s]

Elements

CON4 (Tet)

0.02-0.005

12.8

1.3

8956 9

914.3

1582

CON8 (Hex)

0.001-0.01

12.8

3.0

2205 13

3625

860

CON9 (Hex)

0.001-0.01

1.0

3.0

2205 13

3482

805

SPHE3

0.01

12.8

1.0

698

98

8375

SPHE9

0.013

12.8

2.0

917

691

49000

Table 8: Comparison of different models

The summery of the different curves (Figure 34) shows that the conical model with the tetrahedral elements results not in a realistic maximum pressure-distance curve. The closest agreement with the

46

experimental values is obtained with the finest conical mesh. This mesh has a smallest element size of 1 mm! The calculation time is very high because of the small time step sizes that are necessary with the small elements. The difference between the spherical and the conical model shows that the element size has a big influence for the maximum pressures. 1E+7 Kingery SPHE3 SPHE9 CON4 ▲ CON8 ■ CON9 ■ 1kg

Max. pressure [Pa]

8E+6

6E+6

4E+6

2E+6

0E+0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1/3

Scaled distance [m/kg ]

Figure 34: Maximum pressures versus scaled distance

4.6.2 Impulse The impulse is the parameter which has a capital importance for the loading of a structure. The numerical results are shown in Figure 36. The impulse has been calculated with an integration of the pressure- time curve over the time. The impulse is larger for the models with an explosive of 12.8 kg and is smaller for an explosive of 1.0 kg. This depends on the positive phase duration which has also a big difference in the scaling. After the detonation the compressed combustion gas needs time to expand. The dashed line for the model CON9 (Z2=LIBRE, NO REFLECHIE,KINGERY, >3=LIBRE, REFLECHIE,BAKER, >4=DEMI-ESPACE, REFLECHIE,KINGERY, >5=DEMI-ESPACE, NO REFLECHIE,KINGERY)'/ CELSE CHARACTER(32), PARAMETER :: NOM='IMPEDANCE (AIR BLAST)' DATA GET_FMT(:) / >'X-COOR OF THE CHARGE', >'Y-COOR OF THE CHARGE', >'Z-COOR OF THE CHARGE', >'MASS OF THE CHARGE', >'INITIAL TIME OF THE EXPLOSION', >'CONFINEMENT (1=UNCONFINED, REFLECTED,KINGERY, >2=UNCONFINED, NOT REFLECTED,KINGERY, >3=UNCONFINED, REFLECTED,BAKER, >4=HALF-SPACE, REFLECTED,KINGERY, >5=HALF-SPACE, NOT REFLECTED,KINGERY)'/ CENDIF * CALL CREATE_MATERIAL (LENX, LENI, 0) NEWMAT%NAME = NOM NEWMAT%TYPE = 16 !! IMPEDANCE NEWMAT%MATENT(1) = 26 !! NUMERO DE L'IMPEDANCE NEWMAT%NUMLDC = INUMLDC NEWMAT%LGECR = LGECR(NEWMAT%TYPE) * *----- lecture des parametres CALL LIRVAL(NMOT,MOT,NEWMAT%MATVAL,KOPT) * *----- donnees completes ? (1, 2, 4 indispensables) IF((KOPT(1)+KOPT(2)+KOPT(4)) /= 3 ) THEN WRITE(BLABLA,1001) CALL ERRMSS('MAT_AIRB',BLABLA) STOP 'MAT_AIRB 3' ENDIF * * default values IF(KOPT(3) == 0) NEWMAT%MATVAL(3)=0.D0 !Z IF(KOPT(5) == 0) NEWMAT%MATVAL(5)=HUGE(0.D0) ! TINT IF(KOPT(6) == 0) NEWMAT%MATVAL(6)=1.D0 ! CONF * *----- impressions WRITE(BLABLA,1000) NEWMAT%NUMLDC,NEWMAT%NAME CALL MECTSG(BLABLA) IF(IMPRIM(8)) THEN CALL MECVAL(GET_FMT(1),NEWMAT%MATVAL(1)) CALL MECVAL(GET_FMT(2),NEWMAT%MATVAL(2)) CALL MECVAL(GET_FMT(3),NEWMAT%MATVAL(3)) CALL MECVAL(GET_FMT(4),NEWMAT%MATVAL(4)) CALL MECVAL(GET_FMT(5),NEWMAT%MATVAL(5))

GO TO 200 * 200 CONTINUE * 210 CONTINUE * CALL CLOFIC(49) * END

m_material_i_airb.ff MODULE M_MATERIAL_I_AIRB * * material of type "air blast" (16_26) * USE M_MATERIALS USE M_MATERIALS_ARRAY * IMPLICIT NONE SAVE * PRIVATE PUBLIC :: MI_AIRB, CL_AIRB * CONTAINS *============================================ SUBROUTINE MI_AIRB(INUMLDC) * -----------------------------------------------------------------* impedance : air blast m.larcher 04-07 * -----------------------------------------------------------------* iop=26 (airb) * xmut(1) = x-coordinate of explosive charge * xmut(2) = y-coordinate of explosive charge * xmut(3) = z-coordinate of explosive charge * xmut(4) = mass of explosive charge (in kg) * xmut(5) = initial time of the explosion * xmut(6) = choose of different explosion models * 1 = unconfined (full space), reflected (kingery) * 2 = unconfined (full space), not reflected (kingery) * 3 = unconfined (full space), not reflected (baker) * 4 = half-confined (half space), reflected (kingery) * 5 = half-confined (half space), not reflected (kingery) * INCLUDE 'NONE.INC' * INCLUDE 'CARMA.INC' INCLUDE 'POUBTX.INC' * *----- variables globales : INTEGER, INTENT(IN) :: INUMLDC * *----- variables locales : INTEGER, PARAMETER :: NMOT=6 , LENX=6 , LENI = 1 CHARACTER*4 :: MOT(NMOT) INTEGER :: KOPT(NMOT) LOGICAL :: IMPRIM * DATA MOT/'X ','Y ','Z ','MASS','TINT','CONF'/ * INTEGER, PARAMETER :: N_MSG = 6 , LG_FMT = 350 CHARACTER(LG_FMT) :: GET_FMT(N_MSG) *IF FRANCAIS CHARACTER(32), PARAMETER :: NOM='IMPEDANCE (EXPLOSION EN AIR)' DATA GET_FMT(:) /

68

CALL MECVAL(GET_FMT(6),NEWMAT%MATVAL(6)) ENDIF * *IF FRANCAIS 1000 FORMAT('LOI NUMERO',I5,' : ',A) 1001 FORMAT('LA DIRECTIVE "AIRB" EST INCOMPLETE') CELSE 1000 FORMAT('LAW NUMBER',I5,' : ',A) 1001 FORMAT('THE DIRECTIVE "AIRB" IS INCOMPLETE') CENDIF * END SUBROUTINE MI_AIRB *============================================ SUBROUTINE CL_AIRB (MAT_CUR, ECR, D, P, DTAIRB) * -----------------------------------------------------------------* cond. aux limites air blast m.larcher 04-07 * -----------------------------------------------------------------* mat_cur : current material * matval(1) = x-coordinate of explosive charge * matval(2) = y-coordinate of explosive charge * matval(3) = z-coordinate of explosive charge * matval(4) = mass of explosive charge (in kg) * matval(5) = initial time of the explosion * matval(6) = choose of different explosion models * 1 = unconfined (full space), reflected (kingery) * 2 = unconfined (full space), not reflected (kingery) * 3 = unconfined (full space), not reflected (baker) * 4 = half-confined (half space), reflected (kingery) * 5 = half-confined (half space), not reflected (kingery) * d : distance between charge and clxx element centroid * (already computed in the clxx element) * p : blast pressure (output) * ecr(1) : pressure (output) * dtairb : maximum time step (output) * IMPLICIT NONE * INCLUDE 'TEMPX.INC' ! T = CURRENT TIME INCLUDE 'TEMPS1.INC' ! TINIZI = INITIAL TIME * *--- variables globales TYPE(MATERIAL), INTENT(INOUT) :: MAT_CUR REAL(8), INTENT(IN) :: D REAL(8), INTENT(OUT) :: P REAL(8), INTENT(OUT) :: DTAIRB REAL(8), INTENT(INOUT) :: ECR(*) *--- variables locales REAL(8) :: T_START ,T_D , T_CURR, T_NEG REAL(8) :: PARAM1, PARAM2, PARAM3,PARAM4 REAL(8) :: PARAM5, PARAM6, PARAM7,PARAM8 REAL(8) :: P_MAX, B_BLAST , Z_BLAST, U_T_D, Y_CONWEP, > P_MAX1, T_MAX1, T_D1, U_P_MAX, U_T_START, P_NEG REAL(8), DIMENSION(9) :: POLY_T_D REAL(8), DIMENSION(10) :: POLY_T_START REAL(8), DIMENSION(12) :: POLY_P_MAX INTEGER :: I * P = 0.D0 ! initialization DTAIRB = 1e99 ! initialisation - only in the positive phase the time step size will be changed

*----- air blast pressure *----- parameters of detonation for all models Z_BLAST = D/MAT_CUR%MATVAL(4)**0.333333D0 B_BLAST = 5.2777D0*(Z_BLAST**(-1.1975D0)) P_NEG = (0.35/Z_BLAST)*1.E5 IF(P_NEG MAT_CUR%MATVAL(4)**(1./3.)*1E-3 END IF *----- time of explosion------------------------------------------IF (MAT_CUR%MATVAL(5) -0.06017700, 0.0696360, 0.0215297, -0.01616589, 0.0023253, > 0.00147752/) U_T_START = -0.2024257+1.37784*LOG10(Z_BLAST) END IF Y_CONWEP = 0.0 DO I = 1, 10 Y_CONWEP = Y_CONWEP + POLY_T_START(I)*U_T_START**(I-1) END DO T_START = 1D-3 * 10 ** Y_CONWEP *----- time in the pressure curve T_CURR = T - MAT_CUR%MATVAL(5) - T_START IF (T_CURR>0) THEN *---------------------------------------------------------------------*begin different models P_MAX=0 T_D=0 P_MAX1=0 T_D1=0 *-------------------------------------------------------------------*baker incident IF (MAT_CUR%MATVAL(6)==3) THEN PARAM1 = 808.D0*(1.D0+(Z_BLAST/4.5D0)**2.0) PARAM2 = SQRT(1.D0+(Z_BLAST/0.048D0)**2.0) PARAM3 = SQRT(1.D0+(Z_BLAST/0.32D0)**2.0) PARAM4 = SQRT(1.D0+(Z_BLAST/1.35D0)**2.0) *------- max pressure for this distance P_MAX = 1D5 * PARAM1/(PARAM2*PARAM3*PARAM4) PARAM5 = 980.D0*(1.D0+(Z_BLAST/0.54D0)**10.D0) PARAM6 = 1.D0+(Z_BLAST/0.02D0)**3.D0 PARAM7 = 1.D0+(Z_BLAST/0.74D0)**6.D0 PARAM8 = SQRT(1.D0+(Z_BLAST/6.9)**2.0) T_D = 1D-3 * MAT_CUR%MATVAL(4)**(1.D0/3.D0) * > PARAM5/(PARAM6*PARAM7*PARAM8) GOTO 1000

69

> 0.336743114941, -0.00516226351334, -0.0809228619888, > -0.00478507266747, 0.00793030472242, 0.0007684469735, > 0.0, 0.0, 0.0/) *-------------------------------------------------------------------* kingery hemispherical, reflected ELSE IF(MAT_CUR%MATVAL(6)==4) THEN U_P_MAX = 0.240657322658+1.36637719229*LOG10(Z_BLAST) POLY_P_MAX = (/3.4028321, -2.2103087, -0.218536586, 0.89531958, > 0.24989, -0.569249, -0.1179168, 0.2241311, 0.0245620, 0.0455116, > -0.001909307, 0.003614711/) *-------------------------------------------------------------------* kingery hemispherical, incident ELSE IF (MAT_CUR%MATVAL(6)==5) THEN U_P_MAX = 0.214362789151+1.35034249993*LOG10(Z_BLAST) POLY_P_MAX = (/2.780769, -1.6958988, -0.1541937, 0.514050, > 0.0988534, -0.2939126, -0.02681123, 0.109097,0.001628467, > -0.0214631, 0.0001456723, 0.001678477/) END IF *-------------------------------------------------------------------* kingery calculation DO I = 1, 12 P_MAX1 = P_MAX1 + POLY_P_MAX(I)*U_P_MAX**(I-1) END DO P_MAX= 1D3 * 10 ** P_MAX1 DO I = 1, 9 T_D1 = T_D1 + POLY_T_D(I)*U_T_D**(I-1) END DO T_D= 1D-3 * 10 ** T_D1 *-------------------------------------------------------------------*end different models 1000 P = P_MAX*(1.D0-T_CURR/T_D)*EXP(B_BLAST*T_CURR/T_D) IF (T_CURRT_D . AND . T_CURRT_D+T_NEG/2. . AND . T_CURRT_D+T_NEG) THEN P=0 ! AFTER LOADING END IF END IF IF(P 0.00269739758043, -0.00361976502798, 0.00100926577934/) ELSE IF(Z_BLAST>2.28) THEN U_T_D = -3.130058+3.152472*LOG10(Z_BLAST) POLY_T_D = (/0.62103, 0.096703, -0.00801302, 0.00482705, >0.00187587, -0.002467385, -0.000841116668, 0.00061932910, 0.0 /) ELSE U_T_D = 1.33361206714+9.2996288611*LOG10(Z_BLAST) POLY_T_D = (/ 0.23031841078, -0.0297944268969, 0.0306329542941, > 0.0183405574074, -0.0173964666286, 0.00106321963576, > 0.0056206003128, 0.0001618217499, -0.0006860188944 /) END IF *-------------------------------------------------------------------* kingery hemispherical ELSE IF(Z_BLAST-0.00475933, -0.00428144, 0.0, 0.0, 0.0/) ELSE IF(Z_BLAST>2.78) THEN U_T_D = -3.53626+3.463497*LOG10(Z_BLAST) POLY_T_D = (/0.6869066, 0.09330353, -0.00058494, 0.0022688499, > -0.000295908, 0.0014802986, 0.0, 0.0, 0.0 /) ELSE U_T_D = -2.124925+9.2996288*LOG10(Z_BLAST) POLY_T_D = (/0.31540924, -0.0297944, 0.0306329, 0.018340557, >-0.0173964, -0.00106321, 0.0056206, 0.000161821, 0.00068601889/) END IF END IF *-------------------------------------------------------------------* kingery spherical, reflected IF(MAT_CUR%MATVAL(6)==1) THEN U_P_MAX = 0.214362789151+1.35034249993*LOG10(Z_BLAST) POLY_P_MAX = (/3.22958031387, -2.21400538997, 0.35119031446, > 0.657599992109, 0.0141818951887, -0.243076636231, > -0.00158699803158, 0.0492741184234, 0.00227639644004, > -0.00397126276058, 0.0 , 0.0/) *-------------------------------------------------------------------* kingery spherical, incident ELSE IF (MAT_CUR%MATVAL(6)==2) THEN U_P_MAX = 0.214362789151+1.35034249993*LOG10(Z_BLAST) POLY_P_MAX = (/2.611368669, -1.69012801396, 0.00804973591951,

* END SUBROUTINE CL_AIRB *============================================= END MODULE M_MATERIAL_I_AIRB

70

*----------------------------------------*loop over time *----------------------------------------do yi=1,n if(y(xi,yi)y0) then tend=t(yi) if(endimp.eq.0) impulse = impulse + tdist(yi)*(y(xi,yi)-1e5) if(t00.49).and.(x(xi,1)1.001e5)then impulse=-impulse td=-td end if write (6,102) x(xi,1), ymax(xi), td, impulse, ta(xi) 102 FORMAT(5E11.3) end do END PROGRAM airblastresult

10.2 Miscellaneous code airblastresult.f * Reads a file which is build with europlexus LIST-operator. * Writes as an output a file with the peak pressure, the impulse, * the arrival time and the duration of the positive phase. * airblastresult dyna7d.txt * arrival time - arrival of the first pressure PROGRAM airblastresult * x(location_number,time)=location, y(location_number,time)=pressure real*8 :: x(100000,1000), y(100000,1000), ymax (1000) real*8 :: t(1000), tdist(1000), flagmax(1000), ta(1000) real*8 :: t0, tend, impulse,td integer ijk(4,100000) integer :: n, ival, i, xi, yi, icom, flag, endimp character(len=30) :: text1 n=0 do n=n+1 xmax=0 read(5,100,END=999) text1, ival, text2, icom 100 FORMAT(A8, I8, A11 , I8) if(ival1) tdist(n)=t(n-1) tdist(n)=t(n)-tdist(n) *----------------------------------------*Reading of the values *----------------------------------------DO i = 1, ival read(5,*,END=999) x(i,n), y(i,n) 2001 FORMAT(2E14.5) if((y(i,n)>1e5).and.(ta(i)==0))ta(i)=t(n) if(y(i,n)>ymax(i)) then if(flagmax(i)1).and.(y(i,n) maxar); maxar=ar; FINS; FIN I0; MESS maxar;

calc_b.cpp // calc_b.cpp : calculation of the b-factor for the Friedlander // equation and the impulses resulting of a given b #include "stdafx.h" #include #include #include #include #include

lis1 = PROG di1 di2 di3 di4 di5 di6 di7 di8 di9 di10 di11 di12; min2 = mini(lis1); max2 = maxi(lis1); asra=0; SI (min2 > 0.0); asra = max2/min2; FINS; finproc asra;

using namespace std; /*Calculation of the maximum pressure*/ double p0(double z,int flag){ double pp; if(flag==1){ //Kinney double zaehler=808.*(1.+pow(z/4.5,2.)); double nenner1=sqrt(1.+pow(z/0.048,2.)); double nenner2=sqrt(1.+pow(z/0.32,2.)); double nenner3=sqrt(1.+pow(z/1.35,2.)); pp=zaehler/(nenner1*nenner2*nenner3); } else{ //Kingery double t=log10(z); double u=-0.214362789151+1.35034249993*t; double y=2.6113686691.69012801396*u+0.00804973591951*pow(u,2.)+0.336743114 941*pow(u,3.)-0.00516226351334*pow(u,4.) -0.0809228619888*pow(u,5.)0.00478507266747*pow(u,6.)+0.00793030472242*pow(u,7.)+0 .0007684469735*pow(u,8.); pp=pow(10.,y); } return pp; } /*Calculation of the length of the positive phase*/ double td(double w,double z,int flag){ double tdd; if(flag==1){ //Kinney double zaehler=980.*(1.+pow(z/0.54,10.)); double nenner1=1.+pow(z/0.02,3.); double nenner2=1.+pow(z/0.74,6.); double nenner3=sqrt(1.+pow(z/6.9,2.)); tdd=pow(w,1./3.)*zaehler/(nenner1*nenner2*nenner3); }

aspect_ratio.dgibi *Example for the calculation of the aspect ratio *Construction d'une sphere a partir d'un cube opti donn 'aspectratio.procedur'; opti dime 3 elem cub8; *Nombre de bissections nel0a =10; nel0 = 10; r0 = .25; sizeex = 0.5; *Reference o0 = 0. 0. 0.; x0 = (sizeex) 0. 0.; xa0 = 0 (sizeex) (sizeex); xb0 = (sizeex) (sizeex) (sizeex); xc0 = (sizeex) 0 (sizeex); xd0 = 0 0 (sizeex); y0 = 0. (sizeex) 0.; z0 = 0. 0. (sizeex); symp1 = (sizeex/2.) (sizeex/2.) (sizeex/2.); *Cube intermediaire (centre=o0 et arete=r0) vol0 = o0 droi nel0a x0 tran nel0a z0 volu tran nel0a y0 coul bleu homo o0 r0; cub0 = (o0 droi nel0 y0 tran nel0 x0) et (o0 droi nel0 z0 tran nel0 y0) et (o0 droi nel0 x0 tran nel0 z0) syme 'POINT' symp1 homo o0 r0; cub1 = (o0 droi nel0 y0 tran nel0 x0)

72

/*Calculation of the pressure at a certain time step*/ double p(double w,double z,double t,double b,int flag){ return 1.+p0(z,flag)*(1.-t/td(w,z,flag))*exp(-b*t/td(w,z,flag)); }

else{ //Kingery double t=log10(z); double y,u; if(z