LANL NumerEx MACH2 Explosive Model 1993

A Dynamic Explosive Model for MACH2 With Applications to Magnetic Flux Compression Generators

John J. Watrous and Michael H. Frese March 2, 1993

Report 93-04

NumerEx

1400 Central SE, Suite 2000 Albuquerque, New Mexico 87106-481 1

Los Alamos National Laboratory Contract 9-XTbQ4340-1

This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, make any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein t o any specific commercial product, process, or trademark, manufacturer, or service by trade name, otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

DISC LA I M ER

DISCLAIMER

Portions of this document may be illegible in electronic image products. Images are produced from the best available original document.

… List of Figures … … … … … … … … … … … … … . .

1 Introduction … … … … … … … … … … … … … . . 1 2 The Dynamic Explosive Model … … … … … … … … … .

3 One-Dimensional Demonstration … … … … … … … … . . 4 4 Two-Dimensional Demonstration … … … … … … … … .

5 Conclusions… … … … … … … … … … … … … .

Listing of One-Dimensional Demonstration Input File … . . 19 Appendix A Listing of Two-Dimensional Demonstration Input File … . . 22 Appendix B Listing of MACH2 Modifications File … … … … … . 27 Appendix C

Contents

List of Figures

Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Figure 11 Figure 12 Figure 13 Figure 14 Figure 15 Figure 16 Figure 17 Figure 18 Figure 19 Figure 20 Figure 21 Figure 22 Figure 23

Mass density at 10 ps … … … … … … … … … .

Temperature at 10 ps … … … … … … … … … .

Pressure at 10 ps … … … … … … … … … … .

Detonated material tag at 10 ps … … … … … … … 6 Flow speed at 10 ps … … … … … … … … … . .

Mass density at 22.5 p s … … … … … … … … . .

Temperature at 22.5 ps … … … … … … … … …

Pressure at 22.5 ps … … … … … … … … … . .

Detonated material tag at 22.5 ps … … … … … … . .

Flow speed at 22.5 ps … … … … … … … … … .

Initial simulation grid for two-dimensional demonstration … 11 Simulation grid at 10 p s … … … … … … … … . 12 Fluid velocity at 10 ps … … … … … … … … . .

Azimuthal magnetic field at 10 ps … … … … … … . 13 Internal energy density at 10 ps … … … … … … . .

Simulation grid at 20 ps … … … … … … … … .

Fluid velocity at 20 ps … … … … … … … … . .

Azimuthal magnetic field at 20 ps … … … … … … . 15 Internal energy density at 20 ps … … … … … … . .

Simulation grid at 30 ps … … … … … … … … .

Fluid velocity at 30 ps … … … … … … … … . .

Azimuthal magnetic field at 30 ps … … … … … … . 17 Internal energy density at 30 ps … … … … … … . .

…

A n explosive model has recently been designed and added to MACH2 to enable that code to be used as a tool for studying explosive magnetic flux compression generators. This report describes this model and gives examples of its use in both one- and two-dimensional simulations. Section 2 will provide a description of the model. One-dimensional simulations will be discussed in Section 3. Section 4 will show examples of two-dimensional simulations. Appendices contain input decks for the one- and two-dimensional simulations and a listing of the modifications made to MACH2 for this purpose.

Introduction
The Dynamic Explosive Model

Models useful for computational modeling of detonations have been under con- tinuous use and development for decades. They have achieved a high degree of so- phistication, and are capable of predicting experimental results with quite adequate accuracy. However, such models have not yet found wide use in multidimensional MHD codes. As the behaviors of primary interest in an explosive magnetic flux compression generator are those of the detonating driver and of the compressing magnetic field, an MHD code with an explosive model is potentially of great value for designing such devices.

The intent of this effort was not to produce a new explosive model, nor to endow MACH2 with the most sophisticated state-of-the-art model available. Rather, the intent was to outfit MACH2 in the most expedient means possible with a simple model that would allow it t o simulate the basic features of a detonating explosive. Thus, the primary concerns guiding the design reported here were (1) that the model be capable of reproducing the essential qualitative behavior of an explosive and provide some degree of control over the quantitative behavior, and (2) that the model be compatible with the existing structure of MACH2. The qualitative features which were considered essential were that the explosive exhibit a threshold behavior i.e. that it not detonate unless sufficiently perturbed, that the process of detonation release a set amount of energy which would then be available to drive on the detonation, and that once material had detonated, it would not be able t o detonate again. Quantitative features deemed important t o match were the pressure in the newly detonated material and the propagation speed of the detonation wave. With respect to the structure of MACHB, the goal was to produce a modular package that could fit into the existing multiblock, arbitrary Langragian-Eulerian MHD algorithms without requiring that they be modified extensively.

The model used here is an adaptation of a technique proffered by Dr. G. McCall of Los Alamos National Laboratory. The essence of the technique is to define a set of conditions under which undetonated material will detonate, then to

monitor the conditions in each undetonated cell until the chosen conditions are met. Once that occurs, the value of the internal energy of that cell is incremented by a certain amount, and the material is marked as detonated. The detonation conditions suggested by McCall are (1) the artificial viscosity in the given cell is decreasing and (2) the internal energy in the cell exceed a critical value. These conditions will now be described in greater detail.

Artificial viscosity, denoted here as q, is a numerical technique used predomi- nantly in Lagrangian calculations to enhance the numerical stability of the hydro- dynamics algorithm in the presence of a shock wave. Without artificial viscosity, hydrodynamic quantities tend to oscillate wildly behind a shock transition. The stabilizing nature of q is effected by modifying the momentum equation to read

d . + —2, = - - V ( p + q ) , dt P where the hydrodynamic pressure has been augmented by the inclusion of the artificial viscosity, q. Recipes for calculating q abound, but most share the feature that q is zero unless the cell is under compression, in which case q is set to a positive value tied t o the local gradient of the fluid velocity. MACH2 uses the following recipe:

if V - v ’ < _ O

where p is a dimensionless input constant, typically near unity in value, p is the mass density, and dl is the cell size. Strong shocks are characterized by a value of q equal in magnitude to p . Numerically, the density ahead of a shock increases for several cells, telegraphing the arrival of the shock. The size of this precursor zone can be controlled by the parameter p , but the point is that in this region, since it is undergoing some degree of compression, q has a nonzero value. In fact, q increases over several orders of magnitude in this transition region, until it is comparable t o the peak pressure. Behind this point, that is, on the upstream side of the shock, g decreases over a few cells until the shock compression is over. Thus, the point at which q begins to decrease is a useful indicator for the peak of the shock, and motivates its use as part of the detonation criteria.

A decreasing value of g is not alone sufficient to provide a useful detonation cri- terion. Density gradients can diffuse numerically at the grid speed, cg = % leading to very small compressions and hydrodynamically insignificant nonzero values of q. Thus, keying on decreasing values of Q alone would lead to a detonation wave whose speed would be determined by the local grid speed, and would occur in response to the most minute of perturbations. However, imposing the simultaneous requirement that the cell’s internal energy exceed a critical value has the effect of filtering out these uninteresting and unimportant compressions. The precise value of the critical energy density depends on the explosive being modeled. Likewise, the value of the

internal energy to which to boost a detonating cell depends on the material being modeled. We found that a critical energy density of 5x105 J/m3 and a detonation energy density of 5 x lo6 J/m3 allowed detonation speeds and pressures comparable to those reported for many high explosives t o be obtained in the simulations.

Once the material in a cell has detonated, it is crucial to tag it in some way so that it cannot detonate again. MACH2 has a convenient method for doing this, but its use requires that the following restrictions be followed. A block containing the explosive must consist entirely of the explosive. Only one type of explosive is allowed. Each explosive block must be entirely Lagrangian, and the boundaries of such blocks must also be entirely Lagrangian. Finally, as the designation of explosive material as detonated or not detonated makes use of MACH2’s multimaterial capability, this capability is effectively disabled in explosives calculations. Once the material in a cell has detonated, its value for the variable con2 is changed from zero to one. This instructs the equation of state routine to use a different equation of state for the detonated material than for the undetonated material. As the calculation is required to be entirely Lagrangian, the contents of any given cell do not mix with those of its neighbors. This insures that once a material has detonated, it will not flow into a region of undetonated material and contribute more than once to the ongoing explosion.

The equation of state used here for the detonated material is the Gruneisen equation of state. This is an analytic formulation wherein the pressure is determined by the expression

where E is the specific internal energy, p is the mass density, I’ is the Gruneisen coefficient, po is the reference density, and C, is the reference sound speed. Numerical tests show that the parameters I’ and C , have some degree of influence over the speed of the propagating detonation wave. In tamped 1-d tests, with a reference density of 1.894 x lo3 kg/m3, an initial internal energy density of 0.267 x lo6 J/m3, and a detonation energy release of 5.0 x lo6 J/m3, the speed of the detonation wave for c, = 0.141 cm/ps was observed to increase from 0.527 cm/ps to 0.638 c d p s as r increased from 2.0 to 2.37. For r fixed at 2.0, the detonation wave speed increased from 0.492 c d p s to 0.550 c d p s as c, varied from 0.01 c d p s to 0.282 cm/ps. Thus, these two parameters allow the detonation speed to be tuned, giving the user some degree of control over the detonation characteristics.

The effects of different detonation initiation schemes have not yet been studied here in any great detail. The initiation method used to date has been a simple, but somewhat artificial, method wherein a chosen group of computational cells is

assumed to be undergoing detonation just as the computation starts. This is done by giving to these cells just enough internal energy to insure that their internal energy density exceeds the critical value. Their artificial viscosity is also adjusted by setting the old values of q to 2 and the new values t o 1. This insures that in these chosen cells, q appears to be decreasing, and thus allows the detonation criterion t o be met.

One-Dimensional Demonstration

One-dimensional tests afford a quick and convenient means of surveying the characteristics of the detonation model. The Id demonstration reported here used 20 cells to span a distance of 20 cm. The initial conditions were a uniform distribution of mass with density 1.894 x lo3 kg/m3, temperature of 0.025 eV (room temperature), and zero flow velocity. The Gruneisen parameters I’ and c, had values 2.37 and 0.141 c r d p s , respectively. The detonation was initiated at the right and propagates to the left. Boundary conditions at both the right and left represented impenetrable, immobile walls so that this was a tamped detonation.

Figures 1,2,3,4, and 5 show the mass density, the temperature, the pressure, the detonated material tag, and the flow speed, respectively at t = 10 ps. The detonated material tag has a value of 0 for undetonated material and a value of 1 for detonated material. Plots of these same quantities at t = 22.5 ps are shown in Figures 6, 7, 8, 9, and 10.

CXPmSIVE YODLL: lD D u l o

1.0002?.-05

0.r - 1.OOC-07

CYCLC -

0.05

0 . 1 0 YC (M)

Figure 1 Mass density at 10 ps.

Figure 2 Temperature at 10 ps.

0 . 1 0 YC ( M I

Y o W

I 0.05

0 ’

W I - *

0.15

0 . 1 5

CYCLE - I M - DT - 1.00L-07

1.0002E-05

T M - 1.0002L-05 !n - 1.ooc-07

llzIpL

0.05

O C rn

o t

Figure 3 Pressure at 10 ps.

Figure 4 Detonated material tag at 10 ps.

YC ( M I

0.05

0.10 YC (M)

0.10

0.15

ZxPu)SIVE HODEL: lD DIM0

TI3lZ - 1.00022-05

DT = 1 . 0 0 2 - 0 7

CYCm =

0.15

1 1 4

1-2XeL

rl I

0.05

0.00

0 . 1 0 fM)

Figure 5 Flow speed at 10 ps.

Figure 6 Mass density at 22.5 p s .

0 . 1 0 YC ( M I

0.05

In rl

l*mL

0.15

LXPLOSIa YODEL: 1D DWK)

C Y C m - T D a - 2.25022-05 DT - 1.OOL-07

K X P L O S ~ KODZL: 1D D

D-T I 1 . 0 0 E - 0 7

CYCLE -

TIMI *

0.10 YC ( M I

l U X P L

2.25023-05

0.05

Figure 7 Temperature at 22.5 ,US.

Figure 8 Pressure at 22.5 ps.

I 0.05

0 I

l=KX!?L

239

C I C I Z L

TWIK - 2.2502L-05 DT - 1.OOE-01

LXPLOSIVZ MDKL: 1 D DKMO

0.10 ( M I

0.15

WBLOSIVL MODLL: 1D D M

T T M L 2.2502L-05

m = 1.00L-07

CYCLE -

CYCLE L

1*1(PL

0.05

N o z

0 . 1 0 YC ( M I

Figure 9 Detonated material tag at 22.5 p,s.

Figure 10 Flow speed at 22.5 ps.

$ Y W >

0 -

0.10 ( M )

h V w .

10.00

0.05

n I

ID D M

1 - E x e L

239

rr)

0.15

0 . 1 5

LXPLOSTVL YODLL:

T m L 2 . 2 5 0 2 L - 0 5

Dl = 1.00s-07

4. Two-Dimensional Demonstration

Figure 11 shows the grid used for the two-dimensional demonstration. The geometry is cylindrical with radius increasing from left to right and axial distance increasing from bottom t o top. The region from the left hand boundary at r = 10 cm to r = 15 cm is a void region. The thin layer of cells from r = 15 cm to r = 16 cm represents a n aluminum armature. The SESAME tables are used for the aluminum equation of state. The region from r = 16 cm to r = 20 cm contains the explosive material. The parameters describing the explosive are identical to those used in the one-dimensional demonstration reported above. The detonation is initiated in the row of explosive cells at the top of the &d plot. An azimuthal magnetic field is present in the void region. The field has a maximum value of 1 T at the inner radius of the void region and decreases as l/r. A return current flows on the inner surface of the armature; there is initially no magnetic field in either the armature nor the explosive.

Figures 12, 13, 14, and 15 show grid, fluid velocity, magnetic field, and internal energy density at t = 10 ps. At this point, the detonation is well under way; the portion of the armature nearest the initiation region has been accelerated to a speed of approximately 0.167 cdps. Figures 16, 17, 18, and 19 show the same quantities plotted at t = 20 ps, figures 20,21,22, and 23 at t = 30 ps. Averaged over the 30 ps it has taken the detonation wave t o cross the simulation domain, the detonation speed is 0.67 c d p s . On this time scale, the armature acts as a conductor, sweeping up the magnetic flux contained between it and the inner conductor. Some magnetic field does penetrate into the armature, but essentially none diffuses completely through it. The void region has been reduced in area by roughly a factor of two at t = 30 ps. This corresponds to the approximate doubling observed in the magnetic field intensity.

c1IcuLWJ!ION I a S E

MoDEL:?An NLW-aL

1 - 1.0001-11 mLz I IS? x - 1.001-01 x INC - 2.001-02 1ST ‘I - 0.001+00

Y INC = 5.001-02

8.8.1

Figure 11 Initial simulation grid for two-dimensional demonstration.

8.8.1

CALCUUTION YCSE

M0LlLL:TAB NLW-AL

Figure 12 Simulation grid at 10 p s .

T - 1.0001-05 CYCLE = 1ST x - 1.002-01 x INC - 2.001-02 1ST Y - 0.001+00 Y INC - 5.001-02

Figure 13 Fluid velocity at 10 ps.

MOD1L:TAB ?lKl?-AL T = 1.0002-05 CYCLE = v?xLcITY

kUx - 2.2822+03

8 . 8 . 1

1 1 4

w0Du:TAB NCW-AL

2 - 1.000c-0s CYCLZ - 114

T W I D A L UAGNJETIC TmLD

8.8.1

-= 1.4L-07 E= 1.9L-01 D= 3.9C-01 ,r 5.8I-01 a- 1.m-01 +- 9.7c-01

. I s

Figure 14 Azimuthal magnetic field at 10 ps.

Figure 15 Internal energy density at 10 ps.

2.6lAQ5 B- 5.1W05 D- 1.0C+06 2.0Z+06 E- 3.8t+O6 +- 7.4C+06

NEW-AL

8.8.1

MODU:TAB

U T I 0 . 0 0 1 + 0 0

CALCULATION XCSE

Figure 16 Simulation grid at 20 ps.

T - 2.0001-05 CYCLE - 214 1ST X - 1.001-01 x INC - 2.001-02 I INC - 5.001-02

T - 2.0001-05 CICTZ - 214 IUX - 3.2381+03

Figure 17 Fluid velocity at 20 ps.

M W L L : T ~ NEW-AL

VLL€U=I.lY

8.8.1

X0DLL:TAB mPAL

? = 2.000L-05 C Y C U - 214 TOROIDAL IUGNCTIC ETLLD — 1 . a - 0 5 B- 2.3L-01 DE 4.5L-01 r= 6.8L-01 8- 9.1L-01 +- 1.13+00

Figure 18 Azimuthal magnetic field at 20 ps,

Figure 19 Internal energy density at 20 p s .

8.8.1

M0DCL:TAB N E P A L

V I 2.OOOE-05 CYCCTI: - 214

s a c . 1m. LmR6Y

-I 2.6t+05 8- 5.11+05 D- 1.0L+06 r- 2.0LM6 8- 4.1E+06 +- 8.11+06

8.8.1

M0DLL:TAB NCW-AL

Y I X I 5.OOC-02

CALCUUTION KCSB

Figure 20 Simulation grid at 30 ps.

T - 3.0002-05 CYCLE * 1ST X - 1.001-01 x INC - 2.00c-02 1ST Y - O.OOC+OO

T - 3.000L-05 CYCLE - YU - 8.223C+03

Figure 21 Fluid velocity at 30 ps.

M0DCL:TAB NW-AL

VELOCITY

8.8.1

nwu:ma NW-u

T I 3.0001-05 C Y C I Z - 404

TaRoXDAL -TIC

8.8.1

?RIB

Figure 22 Azimuthal magnetic field at 30 ps.

Figure 23 Internal energy density at 30 p s .

8.8.1

nwr.L:nn IPW-AL

f - 3.OOm-05 CYCLE = SPEC. XM. L a R G Y — 2.62+05 Br 5.31+05 D- 1.1L+06 I= 2.21+06 61 4.52+06 +- 9.21+06

I04

5. Conclusions

The dynamic explosive model described here appears to be a useful device for allowing MACH2 to simulate detonations. It has the positive features of being quite easy to institute, use, and modify, and of causing the detonation t o act under many of the same physical effects that influence a real detonation. The On the negative side, it is most likely not as accurate as the more sophisticated detonation models that are available, nor does it provide the user with the abundance of controls that more advanced explosive equations-of-state provide.

The capability for designing explosive magnetic flux compression generators that this explosive model gives to MACH2 has been demonstrated by a simple two- dimensional simulation. There are still several areas of work that would make this a more useful and more easily used tool. The most significant of these areas concerns the impact of the armature on the inner conductor. In the simulation presented above, the cells in the void region become very highly compressed as the armature approaches the inner conductor. Research into the relevant physical processes and into an acceptable numerical treatment would be of great benefit in treating the collision. Another area of work that would lead to significant improvements in MACHZ’s capacity t o design such devices is in the way the code handles problems Its present method is adequate for the containing several different materials. simple sort of problem present in the two-dimensional demonstration, but far more interesting problems could be treated if MACH2 were outfitted with an improved multimaterials capability.

Explosive model: Id sensitivity characterization

Appendix A Listing of One-Dimensional Demonstration Input File

t = 1.Oe-11, twfn = 28.5e-6, imns = 60, dt = 1.0e-9, dtmax = 1.e-07,

hydron = .true.,

meshon = .true. ,

radiate = .false.,

radsplit = .false., radflxlt = .false.,

! !

Scont rl

nsmooth = 4, wrelax = 0.25,

volratm = 0.8, courmax = 1.0, rmvolrm = 0.2, itopt = 20, mu = 5.6,

thmldif = .false., tdtol = 1 .e-4,

bdiff = .false.,

rdtol = 1.e-4, aresfdg = 0.05,

ciron = .false. ,

con2on = .true.,

rofvac = 1.e-1, rofjoule = 1.e-1, rof = 1.e-5, rofsiecp = 1.e-4, siecap = 1. e9,

make the thing slabindrical cy1 = 0,

donormn = l., t h e b = 1. , e p s = 1 . e - 3 , conserv = 0 . ,

mglmax = 1,

dtm = l . e l 0 , d t r s t = 1O.Oe9, d t o = 1000.0e-9, d t p = 2.5e-6,

d t s l i c = 2.5e-6, i b d y s l i c = 1, l b l k s l i c = 1, i j s l i c = 2,

i n t t y = ’ e d i t s , l O ’ ,

i n t b o u n d = . f a l s e . , k c o n ( 1 ) = 11, c o n t y p ( 1 ) = ’ l o g ’ p l o t ( 8 ) = ’ n u m v i s ’ , p l o t ( 9 ) = p l o t (11) = sie’ , p l o t ( l 2 ) = ’ d i r k e

$e n d$ o u t p u t

Send $e z geom

1, 2 , 3 , 4,

j o u l h e a t ’ ,

n p n t s = 4 ,

p o i n t x (1) = 0 . 0 0 e - 2 , p o i n t y ( 1 ) = 20.00e-2,

p o i n t x ( 2 ) = 1 . 0 0 e - 2 , p o i n t y ( 2 ) = 2O.OOe-2,

p o i n t x (3) = 1.OOe-2, p o i n t y ( 3 ) = 0.OOe-2,

p o i n t x ( 4 ) = 0 . 0 0 e - 2 , p o i n t y ( 4 ) = 0.OOe-2,

n b l k = 1, corners (1,l) =

icellsg = 2, jcellsg = 20, eosmodlg = “grun”,

roig = 1.894e3, densityg = 1.894e03, csqOg = 16.e6, tempig = 2.5e-2, sieig = 2.67e05, gdvlg = 0.95, gdvlg = 1.0, donorg = .true. , radmodlg = “none”,

ang = 6., awg = 12.,

Send $ inme s h

Send Sezphys

! Send ! Sinmesh

$end ! Smodtim

! $end

! !

name(5) = ‘l=expl’, nigen = 0, niter = 3 , eqvol = 2500., eosmodl(1) = “explosiv” , pretend the material is pbx-9502

siei (1) = 0.2674e06, roi(1) = 1.894e03, an(l) = 5.598417, aw(l) = 10.980013, tempi (1) = 0.025, gmlO(1) = 1.37, csqO(1) = 2.e06, tfusi(1) = 5.e6, hfusi(1) = l.eO2, tvap(1) = 6.e6, hvap(1) = l.eO2, tdiss(1) = 7.e6, hdiss(1) = l.eO2, tionize(1) = 8.e6, hion(1) = l.eO2,

tflow (4,l) = 2.5e-02, velbc(4,l) = ‘freeslip’, probc(4,l) = ‘wall’,

! tmod = 100.e-9,

Appendix B Listing of Two-Dimensional Demonstration Input File

E x p l o s i v e mode1:TAB new-AL

t = 1.e-11, twfn = 50.0e-6, imns = 6 0 , d t = 1 . 0 e - 9 , dtmax = 1 . e - 0 7 ,

hydron = . t r u e . ,

meshon = . t r u e . ,

r a d i a t e = . f a l s e . ,

r a d s p l i t = . f a l s e . , r a d f l x l t = . f a l s e . ,

S c o n t r l

! !

c y 1 = 1,

t h m l d i f = . f a l s e . , t d t o l = 1 . e - 4 ,

b d i f f = . t r u e . ,

r d t o l = 1 . e - 4 , a r e s f d g = 0 . 0 5 ,

c i r o n = . f a l s e . ,

con2on = . t r u e . ,

gdvlmod = . t r u e . ,

rofvac = l., r o f j o u l e = l., r o f = 1. e-5, r o f s i e c p = 1 . e - 4 ,

s i e c a p = 1 . e 9 ,

nsmooth = 4 , w r e l a x = 0 . 2 5 ,

! make t h e t h i n g c y l i n d r i c a l

mglmax = 1,

$e n d$ output

v o l r a t m = 0 . 8 , courmax = 1 . 0 , rmvolrm = 0 . 2 , i t o p t = 20, mu = 5 . 6 , donormn = l., t h e b = l., e p s = 1 . e - 3 , conserv = O . ,

$end Sezgeom

n p n t s = 8 ,

n c y c h i s t = 1,

h i s t y ( 1 ) = 1 5 . e - 0 2 ,

h i s t y ( 2 ) = 15.e-02,

r b z d o t ’ ,

d t m = l.el0, d t r s t = 1O.Oe9, d t o = 1000.0e-9, d t p = 1.0e-06, d t s l i c = 1000.e-9, i b d y s l i c = 4 , l b l k s l i c = 3 , i j s l i c = 2 ,

i n t t y = ’ e d i t s , l O ’ ,

i n t b o u n d = . f a l s e . , k c o n ( 1 ) = 11, c o n t y p ( 1 ) = ‘log’ p l o t ( 8 ) = ’ n u m v i s ’ , p l o t ( 9 ) = I j o u l h e a t ’ , p l o t ( 1 l ) = ’ s i e ’ , p l o t ( l 2 ) = ’ d i r k e

h i s t n u m = 1, p r o b t y p e ( 1 ) = ’ b z d o t ’ , h i s t x ( 1 ) = 1.05e-2, histnum = 2 , p r o b t y p e ( 2 ) = h i s t x ( 2 ) = 1.05e-2,

p o i n t x (1) = 20.00e-2, p o i n t y ( 1 ) = 2 0 . 0 0 e - 2 ,

p o i n t x ( 2 ) = 2 0 . 0 0 e - 2 , p o i n t y ( 2 ) = 0.00e-2,

P o i n t x ( 3 ) = 1 6 . 0 0 e - 2 , p o i n t y ( 3 ) = 0 .OOe-2,

p o i n t x ( 4 ) = 1 6 . 0 0 e - 2 , p o i n t y ( 4 ) = 2 0 . 0 0 e - 2 ,

p o i n t x ( 5 ) = 1 5 . 0 0 e - 2 , p o i n t y ( 5 ) = 2 0 . 0 0 e - 2 ,

p o i n t x ( 6 ) = 1 5 . 0 0 e - 2 , p o i n t y ( 6 ) = 0 . 0 0 e - 2 ,

p o i n t x ( 7 ) = 1 0 . 0 0 e - 2 , p o i n t y (7) = 0 . 0 0 e - 2 ,

p o i n t x ( 8 ) = 1 0 . 0 0 e - 2 , p o i n t y ( 8 ) = 2 0 . 0 0 e - 2 ,

n b l k = 3, c o r n e r s ( 1 , l ) = 4 , 1 , 2 , 3 , c o r n e r s ( l , 2 ) = 5 , 4 , 3 , 6 , c o r n e r s ( 1 , 3 ) = 8 , 5 , 6 i 7,

$e n d$ i n m e s h

S e n d S e z p h y s

i c e l l s g = 4 , j c e l l s g = 2 0 , tempig = 2 . 5 e - 2 , s i e i g = 2 . 6 7 e 0 5 , gdvlg = 0 . 9 5 , g d v l g = 1 . 0 , donorg = . t r u e . , r a d m o d l g = ” n o n e ” ,

n a m e ( 5 ) = ’ 8 . 8 . 3 ’ , n i g e n = 0 , n i t e r = 3, e q v o l = 2 5 0 0 . ,

e o s m o d l ( 1 ) = ” e x p l o s i v ” , pretend t h e m a t e r i a l i s pbx-9502

s i e i (1) = 0 . 2 6 7 4 e 0 6 , r o i ( 1 ) = 1 . 8 9 4 e 0 3 , a n ( 1 ) = 5 . 5 9 8 4 1 7 , a w ( 1 ) = 1 0 . 9 8 0 0 1 3 , t e m p i ( 1 ) = 0 . 0 2 5 , g m l O ( 1 ) = 1 . 3 7 , csqO(1) = 2 . e 0 6 ,

Send ! Smodtim

tfusi(1) = 5.e06, hfusi(1) = l.eO2, tvap(1) = 6.e06, hvap(1) = l.eO2, tdiss(1) = 7.e06, hdiss(1) = l.eO2, tionize(1) = 8.e06, hion(1) = l.eO2,

eosmodl(2) = “grun”, the material is aluminum tempi(2) = 0.025, roi(2) = 2.7e03, matname (2) = ‘al-new’ , resmodl(2) = ’ tabular’, an(2) = 13., aw(2) = 26.9815, gmlO(2) = 1.136, csqO(2) = 29.e06, tfusi(2) = 5.e06, hfusi(2) = l.eO2, tvap(2) = 6.e06, hvap(2) = l.eO2, tdiss(2) = 7.e06, hdiss(2) = l.eO2, tionize(2) = 8.e06, hion(2) = l.eO2,

gdvlm(2) = O.,

rmingdv(2) = 0.104, delrgdv(2) = 0.001,

gdvlm(3) = O.,

rmingdv(3) = 0.100, delrgdv(3) = 0.001,

eosmodl(3 ) = “idealgas”, matname ( 3 ) = I al-new’ , resmodl(3) = tabular’ I pretend the material is voidium

roi(3) = 1.e-3, an(3) = l., aw(3) = l., tempi(3) = 0.025,

gdvl(3) = 0.95, gdvlb(2,3) = l., gdvlb(4,Z) = l., gridbc ( 2 , 3 ) = f ixedgp’ ,

fill in initial b-theta field binit (1 ) = ’ nocurnt’ , bzi(1) = 0.e-02, rnomfld = 1.e-1, binit ( 2 ) = ‘nocurnt’ , anomres (2) = .false., anomres (3) = .false., bzi(2) = 0.e-02, binit (3) = ‘nocurnt’ , bzi(3) = l.eO,

tmod = 100.e-9,

tflow(4,l) = 2.5e-02, velbc(4,l) = ‘freeslip’, probc(4,l) = ‘wall’,

! !

! Send

! $end ! Sinmesh

Appendix C Listing of MACH2 Modifications File

modvers = ‘x’

if (eosmodl (lblk) .eq. “explosiv”) then

call eosxpl

if (con2on) call eoslreg call eosawan if (eosmodl (lblk) .eq. “idealgas”) then

else if (eosmodl (lblk) .eq. “grun”) then

call eosideal

call eosgrun

else

end if

call eostable

else

*d eos.14,22

*id v9102 *define unicos *d mach2.37

*af , , eostrans *dk eosxpl

cdir$ nolist

cdir$ list

end if

pointer(kp095, qold(O:ip2,0:jp2) common/explcom/xplcritn, xplsie, siecrit data xplcritn/O.l/, xplsie/5.e06/, siecrit/5.348e05/

c--- get new artificial viscosity call numvis (mu, cyl)

c--- it is used in the decision to detonate

include common. h’ include inputcom. h’ include pointer. h’

subroutine eosxpl

c--- explosive model c--- tests whether material is detonated or undetonated material c--- if undetonated, tests whether it should detonate, and if it c--- should, adds the user-determined amount of internal energy c--- to mimic detonation, then converts material to detonated state c--- use gruneisen eos for undetonated and detonated material

end if c------------- use gruneisen eos

c--------------- test for detonation

c------------------ material in this cell is to detonate

do 100 j = 1, jcels

do 100 i = 1, icels

C------------ detonation state is governed by con2 if (con2(i, j ) .lt.0.5)then C--------------- cell contains undetonated explosive if ( j .eq. jcels) then

C------------------ artificial initiator

end if

q(i, j) = 1. qold(i,j) = 2. sie(i, j) = siecrit

lreg(i,j) = nreg(lb1k) lregc = lreg(i, j ) awc (i, j) = awanmlt (lregc) * aw (lregc) anc (i, j ) = awanmlt (lregc) * an(1regc)

if (q(i, j ) .It. qold(i, j ) .and. sie(i, j ) .ge.siecrit)then

sie(i, j) = sie(i, j ) + xplsie con2(i, j) = 1.

ef = sieg + tfusi(1regc) / ef2 = ef + hfusi(1regc)

xyz = one - density(1regc) /ro(i, j) pg = density (lregc) *csqO (lregc) *xyz/

cold curve sieg = 0.5d0 * csq0 (lregc) *

( awc(i, j ) * eosfac * gmll(1regc) )

( awc(i, j) * eosfac * gmlO(1regc) 1

(dv/ (one/density (lregc) -fac*dv) ) **2

eosfac = pm / qe az = one - tsplit / ( tsplit + tiny ) fac = 0.5d0 * ( gmlO(1regc) t 2.d0 ) dv = one / density(1regc) - one / ro (i, j)

ev = ef2 + ( tvap(1regc) - tfusi(1regc) ) / ev2 = ev + hvap(1regc)

bracket disassociation ed = ev2 + ( tdiss(1regc) - tvap(1regc) ) ( awc(i, j ) * eosfac * gml(1regc) )

( (one-fac) + facdensity (lregc) /ro (i, j ) ) **2 * (one - 0.5gm10 (lregc) * x y z )

c------------------ gruneisen equation of state

c------------------ bracket vaporivation

c------------------ bracket fusion

c------------------

C------------------

c------------------ bracket ionization

ei = ed2 + ( tionize(1regc) - tdissilregc) 1 / ei2 = ei + hion(1regc)

eosfac*gmI (Iregc) * (sie (i, j) -ei2) *awc (i, j) / ( az + nfe(i, j ) 1 ( az + nfe(i, j ) )

dtde (i, j) = eosfac*gml (lregc) *awc (i, j) /

p(i, j) = pg+gml(lregc)*ro(i, j)*sie(i, j) csq(i, j ) = (gml (lregc)+one) *gml (lregc) *sie(i, j)

endif if (sie(i, j) .gt.ei.and.sie(i, j) .le.eiZ)then

te(i, j) = tionize(1regc) dtde(i,j) = tiny p(i, j) = pg + gml(lregc)*ro(i, j)*sie(i, j) csq(i, j) = (gml (lregc) +one) *gml (lregc) *sie (i, j)

c------------------ near ionization

if (sie(i,j) .gt. ei2) then

te(i, j) = tionize(1regc) +

( awc(i, j) * eosfac * gml(1regc) )

te(i, j) = tvap(1regc) + az *

endi f

c------------------ near disassociation

c------------------ near vaporization

if (sie(i, j) .gt.ev2.and.sie(i1 j) .le.ed)then

if (sie(i, j) .gt.edZ.and.sie(i, j ) .le.ei)then

eosfacgml (lregc) (sie(i, j)-edZ)*awc(i, j )

te(i,j) = tdiss(1regc) + az * dtde(i, j) = eosfac * gml (lregc) * awc(i, j) p (i, j) = pg + gml (lregc) *ro (i, j) *sie (i, j) csq(i, j) = ( g m l ( l r e g c ) + o n e ) * g m l ( I r e g c ) * s i e ( i , j)

endi f if (sie(i, j) .gt.ed.and.sie(i, j ) .le.ed2)then

te (i, j) = tdiss (lregc) dtde(i, j ) = tiny p(i, j) = pg + gml(1regc) * ro(i, j) * sie(i, j ) csq(i, j) = (gml(lregc)+one)*gml(Iregc)*sie(i, j )

eosfacgml (lregc) (sie (i, j) -ev2) *awc (i, j) dtde(i,j) = eosfac * gml(1regc) * awc(i, j) p(i, j) = pg + gml (lregc) * ro(i, j) * sie(i, j ) csq(i, j) = (gml (lregc) +one) *gml (lregc) *sie (i, j)

endi f if (sie (i, j) .gt .ev.and. sie (i, j) . le.ev2) then

c------------------

s o l i d if ( s i e ( i , j ) . g e . s i e g . a n d . s i e ( i , j ) . l e . e f ) t h e n

c------------------

n e a r f u s i o n if ( s i e ( i , j ) . g t . e f 2 . and. s i e ( i f j ) . l e . e v ) t h e n

e o s f a c * g m l l ( l r e g c ) * ( s i e ( i f j ) - e f 2 ) *awc ( i f j )

d t d e ( i , j ) = e o s f a c * g m l l ( 1 r e g c ) * a w c ( i , j ) p ( i f j ) = pg+gmll ( I r e g c ) * r o ( i f j ) * s i e ( i f j) c s q ( i , j ) = csql ( l r e g c )

e n d i f i f ( s i e ( i , j ) . g t . e f . a n d . s i e ( i , j) . l e . e f Z ) t h e n

t e ( i , j ) = t f u s i ( l r e g c ) d t d e ( i , j ) = t i n y p ( i , j ) = pg + g m l l ( l r e g c ) * r o ( i , j ) * s i e ( i , j ) c s q ( i , j ) = csql ( l r e g c )

r e t u r n end

e n d i f

t e ( i , j ) = e o s f a c * a z *

t e ( i , j ) = t f u s i ( 1 r e g c ) + a z *

t i ( i f j ) = eosfac*gml ( l r e g c ) *

i f ( t s p l i t . e q . 0 ) t h e n

. o r . . o r .

i d e a l g a s ’

) then

e n d i f

1 0 0 c o n t i n u e

c--- s a v e old a r t i f i c i a l v i s c o s i t y do 5 0 j = 1, j c e l s

do 5 0 i = 1, i c e l s

q o l d ( i , j ) = q ( i , j )

5 0 c o n t i n u e

*d eosmat. 1 9 / 2 0

if (eosmodl ( l b l k ) . e q .

eosmodl ( l b l k ) . e q . grun’ % I eosmodl ( l b l k ) .eq. e x p l o s i v ’

gmlO ( I r e g c ) * ( s i e ( i , j ) - s i e g ) * awc ( i f j ) d t d e ( i , j ) = e o s f a c * gmlO(1regc) * awC(i, 1) p ( i , j ) = pg + gmlO ( l r e g c ) * r o ( i , j ) * s i e ( i , j ) c s q ( i , j ) = csq0 ( l r e g c )

( s i e i o n ( i f j ) - s i e g ) *awc (if j )

d t i d e ( i , j ) = e o s f a c * gmll ( l r e g c ) * a w c ( i , j) p i o n ( i , j ) = pg

elseif (eosmodl ( l b l k ) . eq. ’ idealgas’ .or. eosmodl ( l b l k ) .eq. ’ grun’ % .or. eosmodl ( l b l k ) .eq. ’ explosiv’ %

elseif (eosmodl ( l b l k ) .eq. ‘grun’ %

eosmodl ( l b l k ) .eq. ‘explosiv’ then

.or.

elseif (eosmodl ( l b l k ) .eq. ‘grun’ .or. %

eosmodl ( l b l k ) . eq. ’ explosiv’ ) then

*d eossie.34

*d eossie.94

*d eosmat - 4 6’4 7

) then