Acta Astronautica. Vol. 1, pp. 899-908. Pergamon Press 1974. Printed in the U.S.A.
Explosion of a hydrogen ice sphere due to a pulsed fast neutron sourcet J. J. R O S C I S Z E W S K I Air Vehicle Corporation, San Diego, California AND
T. T. Y E H (Received 8 August 1973) Abstract---Using a modified L a x - W e n d r o f f numerical technique, the c o m p u t e r calculations of a hydrogen ice sphere exploding due to an energy deposition by slowed-down fast neutrons radiating from a strong point source are presented. Due to the large energy drop at each collision, the slowing-down process is treated by the approximate method, taking into account only the first three collisions. Calculations include the formation of a shock wave in the air outside the expanding h y d r o g e n sphere. B e c a u s e of the very large neutron energy and high source intensity, the binding and ionization energy of the h y d r o g e n ice is neglected and is treated as a gas with a specific heat ratio of 3' = 3. The neutron source in the numerical example c o r r e s p o n d s to the explosion of the d - t micro-bomb of a diameter of about 1.5 mm. Calculations indicate that in order to recover most of the 14 M e V neutron energy, the radius of the hydrogen ice sphere should be of the order of 20 cm.
Introduction ABOUT 80% of the energy ( - 14 MeV) f r o m the reacting d - t fuel is carried b y the fast neutrons. Therefore, the majority of the energy in a small d - t pellet, ignited by a laser or a relativistic electron beam, will be lost, creating an additional environmental hazard. Also, for any pulsed propulsion of a larger vehicle based on a d - t micro-bomb, it is important to have a larger working mass which could be directed by a thrust-producing device such as a magnetic or geometric nozzle. The most effective absorption of neutron energy occurs due to a collision with particles of equal mass, and therefore the o p t i m u m material is hydrogen. F o r this reason, it is necessary to surround the d - t charge with a low atomic mass solid substance in order to effectively slow down the fast neutron. At a strong neutron source, t e m p e r a t u r e s could be very high, of the order of 1 KeV, and pressures of the order of 107 atm. In such conditions, hydrogen ice was treated approximately as a gas with a specific heat ratio of ~/= 3. The ionization energy terms are small c o m p a r e d to the internal energy. A sinusoidal neutron pulse was assumed with the source located at the center of explosion. To avoid a numerical problem at r = 0, a small spherical cavity was a s s u m e d with zero interface velocity. Because of the large energy loss at each collision, the Fermi-age theory cannot be used. Instead, a p p r o x i m a t e treatt T h i s paper h a s been supported by the Air Force Office of Scientific R e s e a r c h under Contract No. F44620-69-C-0109. 899
900
J.J. ROSCISZEWSKI and T. T. YEH
ment, taking only the first three collisions, was applied. Numerical calculations were conducted using the modified two-step L a x - W e n d r o f f finite difference method [ 1], assuming spherical s y m m e t r y . This method enables a continuous integration including shock waves, which are identified a p o s t e r i o r i as locations with steep density and pressure gradients. The p r o b l e m was treated by assuming a finite hydrogen ice sphere and considering the motion of a shock w a v e in the outside a t m o s p h e r e driven by the expanding hydrogen gas. System of equations C o n s e r v a t i o n of m o m e n t u m , mass and energy can be written in the f o r m W, + Fr = S
(1)
where the subscripts t and r denote partial differentiation with respect to the time and radius, and W, F, and S are the following column vectors:
W=
,o(
u2 V, \ ~-+e+--ot)r vP
)
F= tp
(pu ~ + p ) r ~ t | ~ \ ,._,] \ / u 2 p.u r ~-V,
+
/
0
where u = 3, 2 for spherical and cylindrical s y m m e t r y , respectively; p denotes the gas density, u is the velocity, p is the pressure, e is the internal energy, h = e + p / p , V~ is the ionization energy, a is the degree of ionization, and Q is the net energy added to a gas. This energy is equal to the energy deposited by the collision of the slowing d o w n fast neutrons minus the energy lost by radiation. In order to calculate the energy deposited by neutrons, we shall use an a p p r o x i m a t e simplified treatment: W e shall a s s u m e that the time scale for the h y d r o d y n a m i c p h e n o m e n o n is m u c h larger than the m e a n free collision time for fast neutrons. The mean collision time of 14 M e V neutrons in hydrogen ice is of the order of 4.10 -9 sec and the neutron velocity is of the order of 5.10 7 m/sec. Because the energy drop of the neutron at each collision with the hydrogen a t o m is very large, and it is given b y
only the first f e w collisions are important. The time b e t w e e n the birth of a neutron and its collision is m u c h smaller than the h y d r o d y n a m i c time scale which is of the order of R / a w h e r e R is the radius of a hydrogen ice sphere, a is the speed of sound a = ~ / 3 , k T / M , 3~ is the specific heat ratio, and k is the B o l t z m a n n constant. F o r T = 106 °K, a = 105 m/sec which is m u c h smaller than the neutron velocity. T h e r e f o r e , we shall a s s u m e that the neutron energy is deposited instantaneously,
Explosion of a hydrogen ice sphere due to a pulsed fast neutron source
901
as soon as the neutrons are born. However, the neutron source pulse has a finite time, comparable to the h y d r o d y n a m i c time, and therefore the motion of the medium cannot be neglected during the energy deposition process. The energy source Q(r, t) can be expressed as 0qbi W AE,
Q =
(3)
i=1
where qb, is the neutron flux subjected to the ith collision and q is the number of collisions considered. Denoting N ( t ) the number of neutrons born per unit of time in the center, we shall write the neutron flux subjected to the first collision as
Here a small finite inner radius ro of cavity in the hydrogen ice sphere has been used to avoid the singularity at the origin, and J';~ d r / ~ , denotes distribution function accounting for the variation of the mean free path A,, due to the variation of density of the medium (2~, -- 1/n~,, and ~, ~- ~i(E~) is the scattering collision cross section being a function of the neutron energy E,).
(5)
4IRA, exp \ - Sro A , /
at,
Flux of the neutron subjected to the ith collision is qb,(r,)=
o
~
1-exp
\
A ~.,/J
dr,_,
for//>2
(6)
dr,_,
(7)
and a@,(r,)_ 1 Or, A,
", aqb,_, al, exp o a r , ar,
-
A, /
where dl~ is the element of the path of the neutron subjected to the ith collision, and it is a function of the locations of all previous collisions, l, = l~(r,, r: . . . . . r,). We shall assume that all neutrons after the first collision are located on the cone with the apex angle 0 equal to the mean scattering angle with cos 0 - - 2 / 3 A. For hydrogen A = l, cos 0 -- 2/3 and (see Fig. 1)
d12 =
r2dr2 N//r22 - r, 2 sin 0
(8)
Thus, after the second collision,
0~2(r2) f r2 0 ~ , exp Or-----~ - n~2 Jr Orl o
- o,2
"V/1 - (rl/r2 sin ~)2 ,
dr,
= r2 f ' ~ a@, exp [-(X/r2 2 - r~2 sin 2 i f - r, cos O)/h2] -- A2 Jro Or, ~/r2 2_ rl 2 sin2 ~
(9)
902
J . J . R O S C 1 S Z E W S K I and T. T. Y E H
/,
/
Fig. 1. N e u t r o n collision geometry.
where n ~ t1 ~ n0 inside the integration was assumed for simplicity, and ~
_
1
/~or~
_
1
?/o0" 2
In general, the expressions for l, at i ~> 3 are very complicated. For a very rough approximation, the form similar to Eq. (8) was used for L, and the effect of r, was completely ignored. dl3 =
[
dr3
(10)
Such a simplification is justified by the fact that because of the very large energy drop at each collision, the contribution of neutrons for i/> 3 is small. If the number of density n inside the integration is assumed to be constant and equal to no, we have 0@3(r3) _ r3 ( ' ~ O * z exp [ - ('V/r3_2--r:~ s i n 2 0 ~ r : cos O)[A31 dr~ Or~ A, la~,, Or~ ~,/r, 2 - rz z sin z 0
(11)
This would be the case for a neutron pulse much shorter than the characteristic h y d r o d y n a m i c time. In this paper, because the neutron energy drop at each collision with the hydrogen atom is very large, we shall consider only the first three collisions. The initial and boundary conditions were assumed as follows: u(r)=0 u(ro) = 0 u(ro) = 0
at t = 0 for any time (t >~0) for 0 ~< t ~< t,, where t, is a transition time
where hydrogen temperature T~ reaches a given temperature at the contact sur-
Explosion of a hydrogen ice sphere due to a pulsed fast neutron source
903
face r, b e t w e e n the air and hydrogen ice. This t e m p e r a t u r e c o r r e s p o n d s to the t e m p e r a t u r e at which the hydrogen solid b e c o m e s a gas of a very high density. The calculation, which is based on the two-step, modified L a x - W e n d r o f f method, includes a strong shock wave in the air, an expansion wave in the hydrogen and a contact surface. The calculations, except those near the boundary or near the contact surface, are very straightforward. At r = ro, in addition to u(ro) = O, p(ro) was obtained f r o m the conservation of the total mass of the hydrogen when t ~< 6- For t > t.,, u(ro)= 0, p(ro) was extrapolated f r o m the neighborhood. For t ~< t,, it is assumed that u ( r o ) = 0 , o ( r , ) = po, and there is no shock formed. At t > t,, an isentropic expansion, the contact surface matching conditions and a strong shock relation described below were used (shock tube problem). The isentropic expansion wave (initially a plane wave) gives P~=[1
7y
1 u 3 - u4] 2~'/~H-'
p4
a4
(12)
J
The strong shock relation in the air is: us =
~/
2 p2 T~ + 1 pa
(13)
and the contact surface conditions are p2 = p3 = p,. u2 = u3 = Uc
(14)
Before matching the relations, all the variables at the expansion front, denoted by "4", are obtained f r o m a second- (or higher) order polynomial fitting of the k n o w n data at the neighborhood. After there were enough numbers of mesh space between the expansion front and the contact surface, that is, when t > 6, where t2 is the transition time at I r e ( t 2 ) - r 4 ( 6 ) ] / A r = M (with M = 10 meshes) the expansion region was obtained by the L a x - W e n d r o f f method and the isentropic expansion relation was not used. The strong shock relation Eq. (13), the contact surface condition Eq. (14), and the k n o w n data at the neighborhood of the contact surface were used to iterate for obtaining the contact surface position r,., pressure pc, and velocity uc, and thus all the variables at the contact surface were known. The sample runs were p e r f o r m e d at the following conditions. Initial n u m b e r density of hydrogen ice: h ~ = 4.79 × 1028 m-3 Collision cross sections: crl = 7 x 10 -29 m ~, tr2 = 1.1 × 10-2s m ~, or3 = 3.2 × 10 28 m 2 Step size: Ar = 0.025 1~, Cavity radius: ro = 0.1 )~,
Discussion of results The results of the calculations are presented in Figs. 2 through 8. The energy deposited by slowing-down fast neutrons is presented in Fig. 3. One can see that
904
J.J.
ROSCISZEWSK1
a n d T. T. Y E H
7 1". nH To n,
,~~t~se¢ 9z
=20°K =4.79 x 10ZS,T = :5 =3000K =2.685x102s
E, = 2 2 4 x I0 "'z ,(5levels) O
Po
u~ 5
z
I ~.na,r~zXx ~ ~,\N,,N~ ~
4
~
-I
N =10 sin(~rtlO) u, :: X,/X'=7.45X104 m/sec
N N NN
2
0
I
1
0.5
1.0
I~'F'-'--'"-b
1.5
~ - - -
r
2.0
-
2.5
....
~--
5.0
•
3.5
40
)'l Fig. 2. Total energy distribution due to neutron slowing-down process in hydrogen ice
6.0
5.0
N = 102Ssin(lrt 107) Po = 8 0 . 0 4 9 kg/m 3 U~= ~ / t = 7 . 4 5 x l 0 " m/sec
c~ \
4.0
Qi ,oo--~z$ .0 2.0
/"1/
t
•
0.014~$ec
./
I.C
~ "
/(,.,~:-'~
0
--
0.5
-
q
E
I0
1.5
E
2.0 r/X I
"
"
2.5
~
~
5.0
3.5
Fig, 3. Energy distribution due to first, second, and third collision.
4.0
Explosion of a hydrogen ice sphere due to a pulsed fast neutron source
905
I0 r N = IOZSsin(art I0 r)
10 5 T*K
10+
103
I.,QI/, s - - ~
10 2
0
t
t
]
I
t
I
I
2r/~, 3
4
5
6
Fig. 4. Temperature distribution in hydrogen ice sphere heated by a fast neutron halfsinusoidal pulse [N = 102~sin (10;.n't)] located in the center.
the m a x i m u m energy is deposited near the center of explosion and this energy is rapidly increasing in time. At the radius of r = 3)t,, most of the neutron energy is recovered. The m i n i m u m radius of the hydrogen ice sphere is therefore approximately 20 cm. In Fig. 3 it is s h o w n that the contribution of the third collision is not too large and represents f r o m 2.5% of the total energy (at the m a x i m u m energy) to 15% (at the edge). F o r this reason, an error introduced b y the simplification of the expression representing the contribution of this t e r m and b y neglecting higher order t e r m s should be small. Figure 4 indicates a t e m p e r a t u r e p e a k of a b o u t 107 °K near the center of the sphere, decreasing in time due to expansion cooling. Increase of the neutron source strength b y the order of the magnitude, but with the same total energy, results in a very similar t e m p e r a t u r e distribution as is shown in Fig. 5. Figures 6, 7, and 8 indicate the f o r m a t i o n of two shock waves. One is f o r m e d in the beginning at the center of the explosion and propagates outward. Because of zero velocity in the center (we h a v e assumed zero velocity at r = to), the motion of the shock w a v e results in the simultaneous a p p e a r a n c e of the rarefaction w a v e
906
J.J. ROSCISZEWSKI
and T. T, Y E H
Io8[ N = 1029sin(-~t I08)
107
t • O.042~sec
t.~O. I I
T oIK06 I~,'x~l
,o
t0 ~
10 4 -I
103 0
I0.0
l
I
I
i
I
2
5
4
r/~
v
t
5
6
Fig, 5. T e m p e r a t u r e distribution in h y d r o g e n ice s p h e r e heated by a fast n e u t r o n halfsinusoidal pulse [ N = 1029 sin (108 ~'t)] located in the center,
as in a classical strong explosion problem. H o w e v e r , here energy is released in finite time and is distributed over all the body of the hydrogen. A second system of waves is f o r m e d due to the motion of the interface between the outer shell of exploding hydrogen and the air. This system of waves initially looks like that in a classical shock tube problem. The expansion wave, moving into the hydrogen, results in the cancellation of the first shock wave, leaving the pressure almost constant over the whole volume, occupied initially by hydrogen (Fig. 6). At the beginning, the maximum velocity is reached near the center (Fig. 7). Subsequently, the maximum velocity is dropping and moving to the outer edge of the sphere. Later in time, maximum velocity starts increasing as shown in Fig. 7. The density peak becomes maximum (about twice the solid hydrogen density) at t = 4.0/x sec slightly before the first wave reaches the outer edge of the sphere (Fig. 8). The density distribution between the center of the explosion and the first shock wave for later times (t = 4.0/~sec) is not very much different from the similarity solution for the point explosion (broken line in Fig. 8). The assumption
0
I
2._[._5
4
5
~/m 2)
6
Fig. 6. Pressure distribution at atmospheric explosion of the hydrogen ~phere heated by a fast neutron source.
I0
02
03
plO 4
05
)7)
0
I
2
r
5
4
5
6 ce sphere heated by a fast neutron source.
~ig. 7. Velocity distribution at the atmospheric explosion of the hydrogen
I
1 •
45m
~ec
e-
e~
-i
-1
908
J . J . ROSCISZEWSKI and T. T. YEH
2.4
.....
2.0
4.0 2.0
p 1.6
};~ /q/ t =0.5,~ed
I
I
N = I0 ZSsin(-n-t 107) Po = 80049 kg/rn 3 x I =0.0745
/~~Similarity
solution
r.O
It /
0.8
0
I
2
r/)~
3
4
5
6
Fig. 8. Density distribution at atmospheric explosion of the hydrogen ice sphere heated by a fast neutron source (broken line indicates density distribution corresponding to a similarity solution for point explosion at t = 4.0 p~sec with the same density peak).
of constant density, in calculating of integrals for the second and third collisions, should not be too bad taking into a c c o u n t the fact that, for t = 0.1 /xsec, the neutron source is already turned off. As is evident f r o m Fig. 8, for the pulse time of less than 0.05/~ sec the energy deposition process is essentially uncoupled f r o m the h y d r o d y n a m i c p r o c e s s e s and the a s s u m p t i o n of a constant density even in the expression for the first collision could be made in such a case. F o r t = 10/~ sec the effect of the a t m o s p h e r i c shock and the c o u n t e r - p r e s s u r e is still not too large as shown in Fig. 8. Reference 1. gosciszewski, J. and Gallaher, W., Shock Tube Flow Interaction with an Electromagnetic Field in Proceedings of the Seventh International Shock Tube Symposium, I. I. Glass (ed.), pp. 475-489, University of Toronto Press (1970).