- Email: [email protected]

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Numerical simulation of cloud rise phenomena associated with nuclear bursts Y. Kanarska *, I. Lomov, L. Glenn, T. Antoun Lawrence Livermore National Laboratory, P.O. Box 808, L-231, Livermore, CA 94551, United States

a r t i c l e

i n f o

Article history: Received 11 May 2009 Accepted 11 August 2009 Available online 11 September 2009

a b s t r a c t We present numerical simulations of cloud evolution from nuclear explosions using high-resolution numerical methods. Our numerical approach includes a ﬂuid mechanical model that is a combination of a compressible code (GEODYN) and a low Mach code (LMC). Early stages of nuclear explosions that are characterized by the blust wave propagation are simulated with an explicit code (GEODYN) that solves the compressible Navier–Stokes equations via a high-order Godunov scheme. As soon as the blust wave weakens (10 s) the subsequent cloud rise due to buoyancy forces can be effectively simulated by the LMC code. LMC is an implicit code based on a pressure projection technique, and derived from the compressible Navier–Stokes equations using an asymptotic analysis in Mach number. It analytically eliminates time step restrictions based on sound wave propagation and it is computationally efﬁcient and accurate for simulations of cloud rise dynamics at later stages. We perform a series of cloud rise scenarios ranging from an idealized bubble rise problem to realistic air bursts. We analyze effects of compressible dynamics on cloud evolution at different stages. It is found that in a realistic conﬁguration, interaction of the reﬂected shock wave from the ground with the ﬁreball signiﬁcantly affects cloud evolution, in contrast to the equivalent idealized bubble rise simulations. We validate the code predictions against available experimental data. It is demonstrated that, by providing the initial source from the compressible GEODYN code, late time ﬂow evolution can be successfully simulated with the fast, efﬁcient and accurate LMC code. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The problem of radioactive contamination of the environment has existed for more than 50 years (Glasstone and Dolan, 1977). Nevertheless, the danger of a release of radioactive materials still exists, e.g. from terrorist attacks or nuclear accidents. Taking this into account, decision makers need a tool capable of simulating atmospheric transport/dispersion/deposition of radioactive debris released as a result of a nuclear detonation. Prior to 10 s there exist compressible hydro codes capable of simulating the burst phenomenology (Lomov, 2003). At times greater than 103 s meteorological dispersion models exist (Bradley, 2007) and can accurately predict the debris cloud motion. Currently there are some semi-empirical models (Harvey et al., 1992) that have been employed to predict fallout effects within the 10—103 s time range. However, these models are severely limited because they use sparse data sets and/or data from bursts over desert terrain that may not be applicable to other scenarios. Our goal is to develop a predictive computational tool to accurately specify the evolution of clouds from nuclear explosions for a variety of explosion scenarios. Since the initial phase of the cloud rise is dominated by compressible ﬂow effects it is modeled by an adaptive mesh

* Corresponding author. E-mail addresses: [email protected], [email protected] (Y. Kanarska). 0306-4549/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2009.08.009

reﬁnement Eulerian compressible solver GEODYN (Lomov, 2003). As evolution of the cloud proceeds, the Mach number of the cloud decreases and the time step required by the code to model acoustic waves makes simulation of the late time cloud dynamics with a compressible code quite expensive. Furthermore, Thornber et al. (2008a,b) showed that compressible code is characterized by artiﬁcially high dissipation for low Mach numbers. Therefore, for this time period we simulate the cloud rise using a low Mach number code (LMC) where acoustic waves are removed analytically. The code is based on the incompressible algorithm of Bell et al. (2003). Krispin et al. (1999) simulated gas-particulate cloud rise resulting from nuclear explosions using an incompressible code. Since compressible effects due to cloud expansion with altitude are not included in the incompressible formulation this may lead to unrealistic predictions of cloud evolution. Therefore the incompressible code of Bell et al. (2003) was extended to include expansion/contraction effects due to changes in altitude, as well as handling the large variation in densities found in the post shock explosion environment. This formulation has similarities with the low Mach algorithm of Almgren et al. (2006) proposed for Supenovae modeling. To provide the initial source for the nuclear cloud rise problem we explicitly couple the low Mach code to the compressible GEODYN code. We perform a series of cloud rise scenarios ranging from an idealized bubble rise problem to realistic air bursts. We analyze the effects of compressible dynamics and initial

1476

Y. Kanarska et al. / Annals of Nuclear Energy 36 (2009) 1475–1483

conditions on the cloud evolution. The numerical models and algorithms are described in Section 2. Numerical results are discussed in Section 3. Discussion and Conclusions are presented in Section 4. 2. Numerical model 2.1. The gas dynamics equations The gas dynamics is described by the compressible Eulerian equations that can be written in general form as

@Q þ r F ¼ S: @t

ð1Þ

This equation is written in a vector form and contains ﬁve equations for scalar components of a vector Q. F ¼ ðF 1 ; F 2 ; F 3 Þ are the components of the inviscid ﬂuxes, deﬁned as

Q ¼ ðq; qu; qv ; qw; qEÞT ;

F1

¼ ðqu; qu þ p; quv ; quw; uðp þ qEÞÞT ; 2

Fig. 1.Density 1=c proﬁles for real atmosphere (solid line) and isentropic curves: ðzÞ q0 ¼ C pp00ð0Þ for C ¼ 0:2; 0:25; 0:3; 0:35; 0:4.

F2

¼ ðqv ; quv ; qv 2 þ p; qv w; v ðp þ qEÞÞT ;

F3

T

¼ ðqw; quw; qv w; qw þ p; wðp þ qEÞÞ ; 2

where q is the density, v ¼ ðu; v ; wÞ the velocity vector, E ¼ e þ u2 =2 the total energy, and e is the internal energy. The pressure p could be speciﬁed by an equation of state in a general form as f ðp; V; TÞ ¼ 0, where V is a volume, T is the temperature. The numerical approaches presented here are valid for the general case as well. For simplicity, we consider here the case of an ideal gas equation of state in the form

p ¼ ðc 1Þqe ¼ qRT;

ð2Þ

where c is the adiabatic exponent and T is the temperature. The term S ¼ ðSq ; Su ; Sv ; Sw ; SE Þ represents additional sources and forces that exist in the system, e.g. a gravity term for the momentum equation, and will be discussed later. Substituting the deﬁnition of E in the energy equation (1.5) and deriving kinetic energy from the momentum equation we can get an equation for internal energy as

q

De þ pr u ¼ SE : Dt

ð3Þ

The equation (1) directly represent conservation of mass, momentum and energy and provide a suitable model for a broad range of phenomena, including cloud rise events of interest. The formation of particulate clouds resulting from explosions is a complicated process and it involves both compressible and incompressible ﬂow regimes. The initial phase of the ﬁreball formation is dominated by compressible ﬂow effects and propagation of shock waves. At this stage we solve compressible equation (1) using an explicit higher order Godunov Eulerian hydrocode GEODYN (Lomov, 2003). The numerical integrator of GEODYN is based on the second-order Godunov piecewise parabolic method for gasdynamical simulations of Colella and Woodward (1984) that is extended for multiple dimensions using operator-splitting. The ﬂuxes are obtained by solving a Riemann problem at the cell interface using left and right limited quantities. The integration kernel is embedded within the block-structured AMR framework of Berger and Oliger (1984) in order to maximize resolution in the area of interest. The solution procedure is based on an explicit conservative differencing, with timestep size determined by the CFL condition. The single-ﬂuid algorithm is generalized to treat multiple species by adopting a volume-of-ﬂuid (VOF) approach following Miller and Puckett (1996). The further evolution of the cloud is characterized by weakening of the shock waves. Thus pressure inside the ﬁreball reaches equilibrium with the surrounding atmo-

spheric pressure. As a result, the Mach number decreases and the time step required by the code to model the acoustic waves makes simulation of the late time cloud dynamics with a compressible code quite expensive. Furthermore, in nondimensionless form, the pressure gradient in the momentum equation appears with the factor 1=M2 , where Mach number M ¼ juj=c; u is the local velocity and c is the speed of sound. Obviously, as the Mach number vanishes ðM ! 0Þ the pressure gradient term in the momentum equation (1) has a singularity, which is a well-known problem for weakly compressible ﬂows (Choi and Merkle, 1993; Munz et al., 2003). In this regime, compressibility effects have little inﬂuence on the momentum transfer, since pressure becomes only a weak function of density. To prevent this inaccuracy, the pressure can be decomposed into two contributions (Choi and Merkle, 1993; Bell et al., 2003): pðx; tÞ ¼ p0 ðx; tÞ þ pðx; tÞ, with p=p0 ¼ OðM2 Þ. Here, p0 and p are termed the ‘‘thermodynamic pressure” and the ‘‘hydrodynamic pressure”, respectively. One can obtain the low Mach approximation of the governing equation (1) by following standard procedures of asymptotic analysis and neglecting terms of order OðMÞ and less. With this variable decomposition, only thermodynamic pressure appears in the equations of energy, state and modiﬁed continuity equation. In the momentum equation, the gradient of the thermodynamic pressure vanishes, leaving only the gradient of hydrodynamic pressure. As a result, the singularity is removed, the acoustic waves are eliminated and time step is limited by local velocity juj not juj þ c. Let us consider the transformation of the continuity equation ﬁrst. The continuity equation can be rewritten in the form:

ru¼

1 Dq : q Dt

ð4Þ

It can be shown that using equation of state and energy equation (3) the continuity equation (4) can be transformed to form a pressure evolution equation

Table 1 The temperature and density of air equations from the US Standard Atmosphere, 1976 (NOAA) Altitude range (H, km)

Temperature ðT 0 ; KÞ

Density of air ðq0 ; kg=m3 Þ

0–11

288.15 6.5H

11–20

216.65

4:26 6:5 H 1:225 1 288:15 H 11 0:364exp 6:34

Y. Kanarska et al. / Annals of Nuclear Energy 36 (2009) 1475–1483

1477

Table 2 Parameters for the numerical experiments Numerical experiment

Conﬁguration

Dimension

Initial conditions

E1 E2 E3 E4 E5 E6 E7 E8 E9 E10

Compressible Low Mach Incompressible Compressible Low Mach Compressible Low Mach Low Mach Low Mach Compressible + low Mach

2D 2D 3D 3D 3D 2D 2D 3D 3D 3D

R ¼ 400 m; Dq ¼ 0:2; h ¼ 3150 m R ¼ 400 m; Dq ¼ 0:2; h ¼ 3150 m R ¼ 400 m; Dq ¼ 0:2; h ¼ 3150 m R ¼ 400 m; Dq ¼ 0:2; h ¼ 3150 m R ¼ 400 m; Dq ¼ 0:2; h ¼ 3150 m Air burst of 11 kT yield, h ¼ 3050 m (Dixie) Remap from GEODYN (E5) at 15 s Remap from GEODYN (E5) at 15 s, no wind Remap from GEODYN (E5) at 15 s, wind Air burst of 12 kT yield, h ¼ 1200 m (Charleston) remap from GEODYN at 15 s, wind

ru¼

1 Dp

cp Dt

þ

ðc 1ÞSE : p

axisym axisym

axisym axisym

ð5Þ

For the pressure decomposition in the case M 1 the r.h.s of (4) can be approximated as

1 Dp0 ðc 1ÞSE ru¼ þ : p0 c Dt p

ð6Þ

Thus the continuity equation contains only thermodynamic part p0 of the pressure ﬁeld. In the general case p0 can vary in time and space, e.g. Almgren et al. (2006). In the case of a vertically structured atmospheric pressure–density proﬁle p0 ðzÞ and q0 ðzÞ and neglecting thermal conduction, equation (6) can be reduced to

ru¼

w @p0 ; p0 c @z

ð7Þ

or in its equivalent form

r p0 ðzÞ1=c u ¼ 0:

ð8Þ

For an ideal gas, if the background pressure and density are locally isentropic, p0 =qc0 ¼ const equation (8) reduces to well-known anelastic equation (Durran, 1989; Lipps, 1990)

r q0 ðzÞu ¼ 0;

ð9Þ

that is derived under the assumption of small density perturbations. We would like to keep large density variations in the system and prefer to use continuity in the form of (8). In the simulations considered here we assume that p0 ðzÞ is equal to the atmospheric pressure. Fig. 1 shows the density proﬁle in the earth’s atmosphere relative to the isentropic density proﬁle for the same pressure p0 . Despite relatively small deviations between densities, this factor is important for the cloud stabilization that occurs at the altitude where the density inside the cloud becomes equal to atmospheric. Using q0 ðzÞ or ‘‘anelastic approximation” in (9) instead of p0 ðzÞ1=c in (8) assumes that the rate of cloud expansion is equal to the rate of atmosphere expansion with the altitude. In this case the cloud stabilizes at unrealistically high altitudes since only mixing effects lead to cloud deceleration. On the other hand, using p0 ðzÞ1=c in (8) assumes that the cloud expands isentropically and faster than the atmosphere. Therefore, as the cloud rises, the density difference with the atmosphere reduces. This effect, additionally enhanced by mixing, leads to realistic stabilization heights. We will demonstrate it using numerical examples in Section 3. To approximate the atmospheric pressure in (8) we use data for the US Standard Atmosphere, 1976 (NOAA). The temperature variations are approximated with a piece-wise linear function and the density is approximated as exponential and power functions (Table 1). The pressure is derived using the ideal gas low p0 ðzÞ ¼ q0 ðzÞRT 0 ðzÞ. Let’s consider the transformation of the momentum equation. As mentioned earlier, for M 1 the gradient of the thermodynamic pressure vanishes in the momentum equations. If the external forces represents the gravity ﬁeld with the acceleration vector

g ¼ ð0; 0; gÞ, and the reference state values of pressure p0 and density q0 that satisfy hydrostatic equilibrium, i.e dp0 =dz ¼ q0 g, we can derive analytically

@ qu þ r quu ¼ rp ðq q0 Þg; @t

ð10Þ

We do not consider here the transformation of energy equation separately. It is modiﬁed using p0 instead of p in (1.5). If one neglects thermal conduction and there are no thermal sources the density evolution equation (1.1), momentum equation (10) and continuity equation (8) completely describe the ﬂow hydrodynamics. The system can be efﬁciently solved by implicit projection methods, e.g. Chorin (1968), Bell et al. (2003), Almgren et al. (2006). We use the incompressible solver that includes the Boxlib library (Rendleman et al., 2000). The detailed description of the numerical algorithm implementation can be found in Bell et al. (2003), Almgren et al. (2006). An unsplit second-order Godunov procedure of Colella (1990) along wit the MAC projection method (Bell et al., 2003) and MUSCL advection scheme (Colella, 1985) is used. We added several modiﬁcations in the code to simulate nuclear cloud rise events that include extension of the model algorithm for low Mach numbers taking into account the modiﬁed continuity equation (8) and using the surrounding atmosphere pressure proﬁle p0 ðzÞ. As a result the model incorporates expansion/contraction due to changes in altitude and local heating and cooling as well as handling the large variation in densities found in the post shock explosion environment. Additionally, the model was explicitly coupled to the compressible GEODYN code to provide initial source data. The model is capable of simulating fully three dimensional cloud rise scenarios with high resolution ﬁnite difference methods and adaptive mesh reﬁnement in a massively parallel computing environment. 3. Numerical results 3.1. Low density bubble At the ﬁrst stage, we test our numerical approach by simulating low density bubble evolution in different numerical conﬁgurations: compressible, low Mach approximation, and incompressible, for equivalent initial conditions. A spherical bubble of radius R ¼ 400 m is embedded in an air environment at a height h ¼ 3150 m above the Mean Sea Level (MSL). We used the analytically speciﬁed atmospheric pressure p0 ðzÞ and density q0 ðzÞ proﬁles for the US Standard Atmosphere (see Table 1). The density of the bubble is speciﬁed to be one ﬁfth the density of the surrounding atmosphere q0 ðzÞ at the altitude z ¼ 3150 m. The bubble is initially at rest. Parameters for these runs are given in Table 2 (experiments E1–E5). Fig. 2 shows the result of an axisymmetric calculation for a bubble evolution at t ¼ 75; 250 and 500 s in compressible (E1) and low Mach (E2) conﬁgurations. The resolution in both experiments was D ¼ L=N ¼ 9:77 m, where L ¼ 5 km and

1478

Y. Kanarska et al. / Annals of Nuclear Energy 36 (2009) 1475–1483

Fig. 2. Left panel shows vorticity and density ﬁeld evolution for the low Mach experiment E2 at time 75 s, 200 s and 500 s. Right panel shows the corresponding variables in the compressible conﬁguration (experiment E1). Parameters of numerical experiments E1, E2 are presented in Table 2. The cloud rises faster in the compressible code. Difference in simulations can be explained by the fact that the compressible code has higher numerical dissipation in the low Mach limit Thornber et al. (2008a).

N ¼ 512. Both conﬁgurations describe the fact that the buoyancy of the cloud causes the cloud both to rise and to potentially be disrupted by the vortical motions generated by its own rising or by Rayleigh-Taylor and Kelvin-Helmholtz instabilities. At the initial stage the bubble deforms and evolves into a toroidal vortex ring (Fig. 2a). As two large rolls extend out, their convergence under the bubble results in an upwelling of dense air into the bubble, forming a slight high-density tail in the wake. This high-density air is convected to higher altitudes at later times. More dense air

inside the bubble is clearly seen in Fig. 2. Because the density of the expanding cloud does not drop with altitude as fast as the density of the atmosphere, the bubble eventually loses buoyancy and, after a few minutes, stops rising at the so-called stabilization altitude (Fig. 2c). The cloud height HðtÞ in Fig. 3 was computed as an evolution of the center of mass of the cloud:

R 0 q ðx; tÞzdV HðtÞ ¼ R 0 ; q ðx; tÞdV

ð11Þ

Y. Kanarska et al. / Annals of Nuclear Energy 36 (2009) 1475–1483

Fig. 3. Cloud center evolution in GEODYN, LMC and incompressible conﬁgurations.

where V; q0 ¼ ðqðx; tÞ q0 ðzÞÞ and z are volume, density perturbations and altitude correspondingly. While the general features of cloud evolution exist in both conﬁgurations E1 and E2, there are some important discrepancies. Fig. 2 shows that the cloud density is lower in the low Mach experiment E2 than in the compressible experiment E1. Additionally, the vorticity interface in E2 looks more dissipative than in E1. As a result, the cloud rise speed and the cloud evolution is different in E1 and E2 (Fig. 3). These discrepancies can be explained by the differences in the numerical algorithms in the compressible and low Mach approaches. Venkateswaran and Merkle (1998) show that the dissipation term in the momentum equations has an order of Oð1=MÞ for the compressible equations in the low-Mach limit. Additionally, in recent papers of Thornber et al. (2008a,b), it was demonstrated that the leading order kinetic energy dissipation rate in a ﬁnite volume Godunov scheme increases as 1=M as M decreases. In the particular case of the MUSCL advection scheme with van Leer limiter it was shown that the leading order kinetic energy dissipation rate VL can be written as

VL ¼

Dx2 Dx3 a ð3u2xx þ ð2C 3Þux uxxx Þ; uux uxx þ 24 12

ð12Þ

where u and a are the velocity normal to the cell interface and speed of sound, respectively, C is the Courant–Friedrich–Levy (CFL) number and Dx is the cell length. As it is observed in Thornber et al. (2008a), the dissipation is proportional to the speed of sound and the magnitude of the velocity derivatives squared at leading order.

1479

Thus, any low Mach number features are heavily dampened by the compressible numerical scheme. It was demonstrated that the increase in dissipation at low Mach numbers is not physical, but is a property of the discrete system. This may explain the difference between compressible E1 and low Mach simulations E2. Since numerical momentum dissipation is lower in E2 the small-scale structures and eddies are formed at the cloud interface at the initial stage, whereas the interface is smoother in E1 for the same resolution (Fig. 2). Those structures enhance the entrainment of ambient air in the cloud. As a result, the cloud’s density is slightly lower in E2 than in E1 at later times of the cloud evolution and this leads to the discrepancies in the cloud height evolution in E1 and E2 conﬁgurations (Fig. 3). Fig. 3 shows some oscillations of the cloud near stabilization altitude. After reaching the point of neutral buoyancy the cloud continues to move by inertia. Since it is now heavier than the atmosphere the ﬂuid in the center of the cloud begins to fall back toward ground. This results in the convection of low-density air in the upper atmosphere to lower altitudes. Then the cloud’s momentum causes it to fall past the neutral buoyancy point now into a more dense atmosphere, where it again starts to rise and expand. These oscillations, known as atmospheric gravity waves are seen in Fig. 3 at later times. A parcel of ﬂuid vertically displaced from its equilibrium position in a region of stable stratiﬁcation will oscillate at the Brunt–Vasalla frequency deﬁned as

N2 ¼

g @q q @z

ð13Þ

We estimated the period of oscillations at z ¼ 8 and z ¼ 10 km and derived values of 185 and 180 s, respectively. These values are close to the 190 s period of oscillations observed in the simulations in Fig. 3. We also performed similar simulations in the incompressible conﬁguration (experiment E3 in Table 2). Since the compressible effects due to the atmospheric background stratiﬁcation are not included in this conﬁguration the cloud does not expand with the altitude and the incompressible code predicts unrealistically low stabilization altitudes for high altitude bursts Fig. 3. We performed several additional runs with higher resolution to clarify how resolution affects the results. Snapshots of the vorticity ﬁeld at time t ¼ 140 s for two different resolutions D ¼ 2:44 m and D ¼ 9:76 m are shown in Fig. 4. We see that the cloud position remains unaltered as the spatial resolution changes, and the largescale features of the ﬂow, including deformation of the interface and vortex ring development, are unaltered. However, with increasing resolution more small-scale features appear at the cloud interface in the low Mach code (Fig. 4a). These instabilities

Fig. 4. Vorticity and density distributions in Low Mach (E2) and GEODYN (E1) simulations for two different resolutions D ¼ 9:76m and D ¼ 2:44m at time 140 s.

1480

Y. Kanarska et al. / Annals of Nuclear Energy 36 (2009) 1475–1483

Fig. 5. Deviation of pressure from atmospheric, p p0 ðzÞ, in compressible GEODYN simulations at 5, 15 and 85 s is shown by color scale. Deviation from the low Mach constraint in the continuity equation D ¼ wru @p0 for compressible GEODYN simulations is shown by contour lines: D ¼ 0:1 is represented by the black bold line and D ¼ 0:01 p0 c @z

by the red dashed line. (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 6. Vorticity ﬁeld in 2D axisymmetric simulations (left panel) and 3D simulations with wind (right panel) of the Dixie event.

are less pronounced and the cloud interface is smoother in simulations with the compressible code at the same resolutions (Fig. 4b). As mentioned earlier, these differences are explained by increased numerical dissipation for the explicit Godunov scheme in the low Mach limit. In these simulations, there is no additional explicit dissipation mechanism such as viscosity or material diffusion. We assume an implicit subgrid parametrization of numerical dissipation on small scales following ideas of Boris et al. (1992). The bubble is susceptible to growth of both Rayleigh-Taylor and Kelvin-Helmholtz instabilities. While RayleighTaylor instabilities disrupt the bubble on small scales near the center of the bubble, the disruption of the entire bubble is caused by

induced vortical motions of the ﬂuid. The shear along the side of the bubble also induces the Kelvin-Helmholtz instability that generates secondary instabilities as the evolution progresses. A particular feature of the cloud rise is the jet formation in the middle of the cloud. A similar jet was found in numerical and laboratory studies of the unstable growth of an argon bubble subjected to a planar shock wave in Ranjan et al. (2005). However the jet itself is rarely observed in real explosive clouds. This can be explained by that fact that the jet is clearly seen in the vorticity ﬁeld, however its material concentration is very low. Additionally, as shown below, the effect of wind is to destroy the cloud symmetry and the jet formation (Fig. 6).

Y. Kanarska et al. / Annals of Nuclear Energy 36 (2009) 1475–1483

Fig. 7. Cloud center evolution in numerical simulations of the Dixie event in numerical experiments E6–E9. Measured values for the top and bottom of the stabilized cloud (Hawthorne, 1979) are shown by squared points.

1481

Fig. 8. The top and bottom of the stabilized cloud for the Dixie event from different sources: Measurements (Hawthorne, 1979), 3D ‘‘realistic” scenario (E8), 2D ‘‘realistic” scenario (E6), 3D ‘‘idealized” scenario (E4).

Fig. 9. Comparison of the cloud evolution in time: simulations (right panel) and data (left panel). The gas phase is represented by the light gray color. The dust lofted from the surface is shown by the dark color.

1482

Y. Kanarska et al. / Annals of Nuclear Energy 36 (2009) 1475–1483

the cloud top break up quickly in 3D simulations (Fig. 6). Additionally, we performed 3D simulations with wind using the wind proﬁle from Hawthorne (1979). The effect of wind is to destroy cloud symmetry and jet formation (Fig. 5). Due to additional mixing and shear instability the cloud height in 3D simulations with wind is slightly lower than in simulations without wind. In both cases the simulations are close to measured data from Hawthorne (1979) (Fig. 7). The comparisons for the top and bottom of the stabilized cloud in different numerical conﬁgurations and available data are summarized in Fig. 8. 3.3. Low air burst

Fig. 10. The top and bottom of the stabilized cloud for Charleston: 3D ‘‘realistic” scenario (E10), 3D ‘‘idealized” scenario (E9)

3.2. High air burst We performed simulations of the Dixie event, a high altitude (6022 ft) air burst of 11 kT yield. The initial stage of the ﬁreball formation was simulated in the GEODYN code using a 2D axisymmetric conﬁguration (E6). A high energy pill of equivalent energy (11 kT) was used for the initial condition in the GEODYN code. Fig. 6 shows the shock wave formation during the ﬁreball expansion and reﬂection of this shock wave from the ground. After about 15 s the pressure inside the ﬁreball becomes equal to atmospheric. The velocity at this time reaches the value of 100 m/s and the corresponding Mach number is about 0.3. Thus, ﬂow conditions are suitable for the low Mach approximation. We map the velocity, density and particulate tracer ﬁelds from GEODYN at 15 s into the LMC code as the initial conditions (experiment E7). The later stage of the cloud rise was simulated in the LMC code up to 103 s. We performed additional runs where we mapped solution from GEODYN at earlier and later times. At mapping times smaller than 15 s compressible effects become important and cause discrepancies in the ﬁnal ﬂow ﬁeld at later times. Small difference in results is found if we use mapping times larger than 15 s. Thus 15 s was chosen as an appropriate time for mapping from the compressible to the low Mach code. Comparison between 2D GEODYN and LMC simulations (E6 and E7 numerical experiments) is presented in Fig. 7 and shows that the cloud rise speeds are close in both conﬁgurations. Small differences between compressible and low mach simulations can be explained by the dominant effect of the initial vorticity deposition due to the interaction of the reﬂected shock wave from the ground with a ﬁreball. The reﬂected shock wave is seen in Fig. 5 at time 5 s. During this interaction, vorticity is deposited due to the misalignment between the pressure gradient Dp across the shock wave and density gradients Dq across the bubble interface ðDqxDpÞ which later leads to the formation of the primary vortex ring. This mechanism was not included in the idealized simulations and explains difference between the idealized and realistic air bubble simulations (Fig. 7). Despite the fact that the buoyancy difference in ‘‘idealized” (E5) and ‘‘realistic” simulations (E8) is the same, the cloud heights are different. In contrast, the mapped solution from GEODYN to LMC includes the ﬂow ﬁeld after the interaction of shock wave with the ﬁreball, and causes a small difference between compressible and low Mach simulations for experiments E6 and E7 at later times (Fig. 7). As mentioned earlier, the jet at the top of the cloud in our 2D simulations appears too strong and it is rarely observed in real explosive clouds. However the vortex structures at

Finally, we present simulations of the Charleston event (numerical experiment E10). This event belongs to the set of low air bursts that are characterized by the formation of a dust stem. Dusty boundary layers are a common feature of shock interactions with realistic surfaces (Kuhl et al., 1990). The dust appears to be swept out from the surface at the initial stage of the cloud formation due to the close proximity of the cloud to the surface and be involved into the cloud at later times. Data for initial conditions for simulations and wind proﬁles were taken from Hawthorne, 1979. We used the same methodology as described above. The ﬁrst 15 s was simulated in GEODYN and then the solution was mapped into LMC. We assume that the dust layer height is proportional to the shock wave velocity. The initial homogeneous dust distribution with a maximum height of 2 m at the surface was estimated from GEODYN velocity ﬁeld. The dust phase evolution was simulated using transport equation, where dust was advected by the gas ﬂow ﬁeld. This approach is somewhat simpliﬁed since the detail simulation of the dust particles requires multiphase approach, e.g. Lomov (2003). But this is out of scope of the present paper where we focus on the hydrodynamical aspects of the cloud evolution. Results of simulations and comparisons with available data are presented in Fig. 9. The gas phase is illustrated by the light gray color and the dust phase lofted from the surface is shown by the dark color. The shock wave formation during the ﬁreball expansion computed in the compressible code is shown in Fig. 9a at time 5 s. After interaction of the ﬁreball with the surface the dust layer is formed. The dust lofting process and the dust stem formation are clearly seen at times 65 s (Fig. 9b) and 215 s (Fig. 9c). The qualitative agreement with available data is very good. The Charleston burst is closer to the surface than the Dixie burst, so that the effect from the vorticity deposition (due to the interaction of the surface wave reﬂected form the surface with the ﬁreball, described in the previous paragraph) is stronger for the Charleston event. And, indeed, discrepancies between idealized and realistic simulations are more pronounced for the Charleston event (Fig. 10) than for the Dixie event (Fig. 8), that conﬁrming, once again, the importance of the vorticity deposition mechanism for cloud evolution. 4. Conclusions We analyze the cloud rise phenomenon associated with nuclear bursts using high resolution computer codes and available data from shots. Our numerical approach includes a ﬂuid mechanical model that is a combination of an explicit compressible code and an implicit Low Mach code. We analyze effects of compressible dynamics on the cloud evolution. Signiﬁcant discrepancies in the cloud rise speed and interface structure for idealized low density bubble conﬁguration in compressible and low Mach simulations are found. These differences are explained by increased numerical dissipation for the explicit Godunov scheme in the low Mach limit. This behavior is eliminated in the implicit low Mach approach using the pressure decomposition technique. In a realistic conﬁgu-

Y. Kanarska et al. / Annals of Nuclear Energy 36 (2009) 1475–1483

ration, interaction of a reﬂected shock wave from the ground with a ﬁreball creates initial vorticity deposition that affects signiﬁcantly the cloud height prediction. Comparisons between compressible and low Mach simulations show small differences (less than 10%) at later time cloud evolution if this effect is included. The proposed numerical approach is validated against available experimental and empirical data and shows good agreement. Acknowledgments This work performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. References Almgren, A.S., Bell, J.B., Rendleman, C.A., Zingale, M., 2006. Low Mach number modeling of type Ia supernovae. I. Hydrodynamics. Astrophys. J. 637, 922–936. Bell, J.B., Day, M.S., Rendleman, C.A., Woosley, S.E., Zingale, M., 2003. Adaptive low Mach number modeling of nuclear ﬂame microphysics. J. Comp. Phys. 195, 677– 694. Berger, M., Oliger, J., 1984. Adaptive mesh reﬁnement for hyperbolic partial differential equations. J. Comp. Phys. 53, 484–512. Boris, J.P., Grinstein, F.F., Oran, E.S., Kolbe, R.L., 1992. New insights into large eddy simulations. Fluid Dynamics Res. 10, 199–228. Bradley, M., 2007. NARAC: an emergency response resource for predicting the atmospheric dispersion and assessing the consequences of airborne radionuclides. J. Environ. Radioact. 96, 116–121. Choi, Y.H., Merkle, C.L., 1993. The application of preconditioning in viscous ﬂows. J. Comp. Phys. 105, 207–223. Chorin, A.J., 1968. Numerical solution of the Navier–Stokes equations. Math. Comput. 22, 745–762. Colella, P., 1985. A direct Eulerian MUSCL scheme for gas dynamics. SIAM J. Comput. 6, 104–117. Colella, P., 1990. Multidimensional upwind method for hyperbolic conservation laws. J. Comp. Phys. 87, 171–200.

1483

Colella, P., Woodward, P., 1984. The piecewise parabolic method (PPM) for gasdynamical simulations. J. Comp. Phys. 55, 174–201. Durran, D., 1989. Improving the anelastic approximation. J. Atmos. Sci. 46, 1453– 1461. Glasstone, S., Dolan, P.J., 1977. The Effects of Nuclear Weapons. third ed.. US DoD, US DoE. Harvey, T., Serduke, F., Edwards, L., Peters, L., 1992. KDFOC3: A Nuclear Fallout Assessment Capability. UCRL-TM-222788. Hawthorne, H.A., 1979. Compilation of Local Fallout Data from Test Detonations 1945–1962 Extracted from DASA 1251, vol. I. Continental US Tests. Krispin, J., Potts, M., Brown, B., Ferguson, R., Collins, J., 1999. High-order Godunov schemes for multiphase gas-particulate ﬂowﬁeld modeling and simulation. DTRA-TR-99-39, pp. 1–41. Kuhl, A.L., Chien, K.-Y., Ferguson, R.E., Collins, J.P., Glaz, H.M., Colella, P., 1990. Simulation of a turbulent, dusty boundary layer behind a shock. In: AIP Conf. Proc., vol. 208. pp. 762–769. Lipps, F.B., 1990. On the anelastic approximation for deep convection. J. Atmos. Sci. 47, 1794–1798. Lomov, I.N., 2003. Simulation of dense and dilute multiphase compressible ﬂows with Eulerian–Lagrangian approach. In: Proc. of 6th Int. Conf. on Multiphase ﬂows. Miller, G.H., Puckett, E.G., 1996. A higher-order Godunov method for multiple condensed phases. J. Comp. Phys. 128, 134–164. Munz, C.D., Roller, S., Klein, R., Geratz, K.J., 2003. The extension of incompressible ﬂow solvers to the weakly compressible regime. Comput Fluids 32, 173–196. Ranjan, D., Niederhaus, J.H.J., Oakley, J.G., Anderson, M.H., Greenough, J.A., Bonazza, R., 2005. Shock-induced instabilities on a spherical bubble. In: Proc. 25th Int. Symp. on Shock Waves. Rendleman, C., Beckner, V., Lijewski, M., 2000. Parallelization of an adaptive mesh reﬁnement method for incompressible ﬂuid ﬂows. In: International Parallel Distributed processing Symposium. Thornber, B., Drikakis, D., Williams, R.J.R., Youngs, D., 2008a. On entropy generation and dissipation of kinetic energy in high-resolution shock-capturing schemes. J. Comp. Phys. 227, 4853–4872. Thornber, B., Mosedale, A., Drikakis, D., Youngs, D., Williams, R.J.R., 2008b. An improved reconstruction method for compressible ﬂows with low MAch number features. J. Comp. Phys. 227, 4873–4894. Venkateswaran, S., Merkle, C.L., 1998. Evaluation of artiﬁcial dissipation models and their relationship to the accuracy of Euler and Navier–Stokes computations. In: Lecture Notes in Physics. 16th Int. Conf. on Numer. Methods in Comp. Physics, vol. 4(12). pp. 427–432.