Acta Astronautica Vol. 7, pp. 73-78 Pergamon Press Ltd., 1980. Printed in Great Britain
Constant velocity reactive shock waves for testing numerical methods H. M. S T E R N B E R G Naval Surface Weapons Center, White Oak, Silver Spring, MD 20910, U.S.A. (Received
10 O c t o b e r 1978)
Abstract--Piston driven reactive shock waves with constant velocity fronts are constructed for use in testing numerical methods. Explicit solutions are obtained for underdriven, Chapman-Jouguet and overdriven waves in which the reaction rate is proportional to the fraction of material unreacted. REACTIVE shock w a v e s for which exact expressions for the flow are obtainable will be useful for testing numerical methods, like the one recently developed b y Chorin (1977). The w a v e s in mind consist of a constant velocity shock front moving into a constant state gas at rest, driven b y a piston, with reaction behind the front governed b y (I)
d~f/dt = - K / L
a special case of dl3/dt = - K~O"f (T, p ).
(2)
Here/3 is the fraction of gas unreacted, t is the time, f ( T , p) is a specified function of the t e m p e r a t u r e T and the density p, and K and n are constant. The rate function in eqn (1) is used as an example b y Chorin (1977). More complicated cases, involving similarity solutions for reactive w a v e s with non-constant velocity shock fronts are treated in a previous p a p e r (Sternberg, 1970a). In m o s t of these cases, ordinary differential equations m u s t be solved to get the flow behind the shock front. Construction of the constant shock front velocity w a v e s p r o c e e d s as follows. T a k e the piston path as the particle path f r o m the origin in the distance-time (x, t) plane and denote the variables on this path, ahead of the shock front, and immediately behind the shock front, by the subscripts 1,0 and i, respectively. Let E be the internal energy, p the pressure, c the sound speed, u the particle velocity and D the constant front velocity, with y and q positive constants (y > 1). The entire flow will be k n o w n as soon as the piston path and the values of the flow variables along it are found, since the stipulation that the front m o v e s at constant velocity implies that the flow is steady, i.e. the variables are constant on the lines tl).
(3)
E = P/[P(7 - I)].
(4)
x -
xl = D(t
-
Specify the equation of state of the gas to be
73
74
H.M. Sternberg
T h e s t e a d y flow b e h i n d the f r o n t r e q u i r e s t h a t t h r o u g h o u t the w a v e E - E0 - (1 - / 3 ) q = (p + p0)(l/p0 - l / p ) [ 2 ,
(5)
p ( D - u ) = p o ( D - Uo),
(6)
P -Po
= poD(u
(7)
- uo).
N o w , t a k e P0 = 1, p0 = 1, Uo = 0, a n d w = 1 - l / p , so t h a t , f r o m e q n s (4)-(7), u = w D , p = 1 + w D 2, c 2 = 3 " p ( 1 - w),
( D 2 - 3")w - [(3" + l)12]D2w 2 = (3' - 1)(1 - / 3 ) q .
(8) (9)
F o r t h e c o n d i t i o n s j u s t b e h i n d t h e s h o c k f r o n t , s e t / 3 = 1 in e q n (9). T h e n , f r o m e q n s (8) a n d (9), wi = 2 ( D 2 - 3')/[(3' + I)D2], Pi = 1 + w i D 2, ui = w i D .
The Chapman-Jouguet h a s the p r o p e r t i e s
(10)
(C J) d e t o n a t i o n , f o r w h i c h t h e s u b s c r i p t J will b e u s e d , Dj=u~+cj,/is=0
•
(11)
S u b s t i t u t i o n o f e q n s (11) i n t o e q n s (8) a n d (9), a n d u s i n g p o s i t i v e s q u a r e r o o t s f o r the d e t o n a t i o n c a s e , l e a d s to w~ = ( D fl - 3')/[(Y + 1)D~], p j = ( D fl + 1)/(3" + 1), uj = wjD.r,
(12)
q = ( D 2 - 7)2/[2(72 - I)D~21,
(13)
wj = [ ( ~ / - 1)/3']q[{1 + 23"/((3' 2 - 1)q)} ' / 2 - 1],
(14)
p~ = 1 + ( 3 ' - 1)q[{1 + 27/((3, 2 - 1)q)}lt2 + 1],
(15)
D~ = [(3'2 _ 1)ql2]~t2 + [(3'5 _ l ) q / 2 + 3']in.
(16)
F o r the C J s t a t e , f r o m e q n s (10) a n d (12), t h e c o n d i t i o n s a t t h e f r o n t a r e W u = 2 w j , ptl = 2 p j - 1, Uu = 2 u j .
(17)
W i t h a = w l w i , o" = 2(3" - 1)q/[(3' + 1)D2wi2], a n d e q n (9), ~8 = (or - a + a 2 ) / o ",
(18)
D = [(3, 5 - 1)q/(8cr)] It2 + [(3"5 _ 1)q/(8cr) + 3']1/2.
(19)
and
When
t h e w a v e is
a CJ
one,
~r =
1/4 a n d / 3
= (2or -
1) 1/2.
Waves for testing numerical methods
75
The particle path f r o m the origin (the piston path) is obtained f r o m d x / d t = u = w i a D b y using d x = D w i a [(d[3/da )l(d131dt )] d a .
(20)
Notice that eqns (3)-(20) do not depend on the f o r m of the rate function. When d g / d t = - K / 3 , we have f r o m eqns (2), (18) and (20), Dwi (°'
x I = ' - K --J1
a(1-2a) o-a+a
2 da,
tl = - ( I / K ) ln[(o - a l + a l 2 ) / o ] .
(21) (22)
The integration in eqn (21) depends on the coefficients in the d e n o m i n a t o r of the integrand. L e t s = 4 o - 1. The relation b e t w e e n the sign of s and the shock front velocity is gotten f r o m
4
[(Dfl-Y)2 (D~)')2]
s = 4 o - 1 = (7 + 1)2D2wi 2 L
(23)
D~
where eqns (10) and (13) have b e e n used. Equation (23) can be written s =[I/(D:(D2-y)2)](DDj+'y)(Dj+D)(DDj-~,)(Dj-D).
(24)
The sign of s is thus seen to depend on (DD~ - y)(D~ - D ) . F o r detonation w a v e s D D j - y > 0. This is because, in the p,v plane (v = I/p), the slope of the line containing (P0, v0) and (Pi, v~) must be less than the slope of the tangent to the unreacted Hugoniot at (P0, v0). The possibilities for reactive shock w a v e piston paths for steady flows with d g / d t -- - K g , found with eqns (21) and (22), are listed below. I. T h e C h a p m a n - J o u g u e t
w a v e . D = D~, D D : > y , s = O.
x l = ( D w i l K ) [ 2 ( 1 - a l ) - ln(2al - 1)], ]
tu = - ( 2 1 K ) ln(2al - 1), 112 < al ~< 1. /
(25)
At x t = tl = O, P l = Pu, Ul = uu, p~ = pu. T h e CJ state is a p p r o a c h e d as tl~oo. II. The u n d e r d r i v e n w a v e . D < D:, D D j > y , s > O. xl = (DwiIK)[2(I
- a O - (1/2) lnl(o - al + a t ' ) / o l
+ ~/s {tan-1((2al - D/X/s) - tan-~(l/~/s)}],
(26)
tl = - ( l / K ) ln[(o - al + al2)/o], 1/2 < al ~ 1. The underdriven w a v e can be sustained as a w a v e with a constant velocity front for only a short time. This time will c o r r e s p o n d to a value of / 3 ( > 0) which is a minimum. F r o m eqn (18), this is seen to o c c u r for al = 1/2. Notice (eqn 22) that for
76
H.M. Sternberg
the infinite reaction rate case (K--> oo) the entire underdriven wave shrinks to the origin. The limit line which ends the steady underdriven wave is discussed in Sternberg (1970b). Ill. The overdriven wave. D > Dj, DD~ > % s < O. x~ - - - - ~ [ 2 ( 1 - tx~)-~
x
{
- X / ( - s) (27)
2 a ' - 1 - X / ( - s)l - l n I1 - V ' ( - s)[} In 2al l + x / ( - s ) 1-;,X/(-s)
t, = - ( I / K )
+ X / ( - s)12 < al ~ 1.
ln[(o, - al + al~)/o'],(1
In all of the above cases the flow variables on the piston path are calculated with eqn (8), from
(28)
Ul = w i a l D , P l = (1 - wi¢~1) -l, P l -- 1 + w i a i D 2.
The rate function dl3/dt = - K [ 3 was used in the examples which follow. Note that an additional independent quantity must be specified if the wave is not a CJ wave. Example 1 is a CJ wave with ~/= 1.4, q = 12. From eqns (12) and (16), Dj = 5.075818, pj = 11.151635, u~ = 2, p: = 1.650234, and from eqns ( 1 0 ) , w i = 0.788050, p~ = 21.303270, ui = 4. Figure 1 shows p vs x at t = 0.2 for K = 1, 10 and I00. In an infinite reaction rate case (K --, ~) the piston moves at the CJ particle velocity and the CJ conditions hold between the piston and the front. Example 2 is an underdriven wave with y = 1.4, q = 12 and or = 1/2. Here,
R=|
Pj
'°
i
01,
i
01,
i
,i0
i
x
Fig. 1. S t e a d y r e a c t i v e C J w a v e w i t h q = 2, ~ = 1.4. P r e s s u r e v s d i s t a n c e at t = 0.2, for K = l, 10, 101).
Waves for testing numerical methods
77
D = 3.765872, Ds = 5.075818, wi = 0.751068, Pi = 11.651496 a n d us = 2.828427. I n Fig. 2 p v s x is p l o t t e d f o r K = 1.2 at t = 0.5. T h e t h i r d e x a m p l e is a n o v e r d r i v e n w a v e w i t h ~, = 1.4, q = 1, ~r = 3/16. W i t h t h e a b o v e f o r m u l a s o n e g e t s D = 2.228286, Ds = 2.063951, wi = 0.598367, p~ = 3.971048, u~ = 4/3. I n Fig. 3, p v s x at t = 0.314 is s h o w n f o r K = 1, 10 a n d 100. T h e c o n s t a n t s t a t e b e t w e e n t h e p i s t o n a n d t h e f r o n t f o r t h e infinite r e a c t i o n r a t e c a s e is f o u n d b y s e t t i n g /3 = 0 in e q n (9). T h i s l e a d s to D = 2.228286, p ® = 3.228286, u®= 1, p ® = 1.814143. A final r e m a r k will call a t t e n t i o n to a p o i n t m a d e ( S t e r n b e r g , 1970b) r e g a r d i n g CJ w a v e s w i t h c o n s t a n t v e l o c i t y s h o c k f r o n t s . W h e n t h e r a t e f u n c t i o n d[3/dt = -Kf3nf(T, p) is u s e d w i t h n ~> 1 a n d K finite, all t h e c h a r a c t e r i s t i c s s t a r t at t h e s h o c k
12
10 P $ !
1.0
I
l
/
I
I
1.2
1.4
1.6
1.8
2.0
X
F i g . 2. S t e a d y r e a c t i v e u n d e r d r i v e n w a v e w i t h q = 12, T = 1.4, ~ = 0.5. P r e s s u r e v s d i s t a n c e a t t = 0.5, f o r K = 1.2.
4.0
K=I
3.8
3.6 P
3.4
P~ 32
3.0
I
L
I
o'.,
J
oi,
oi,
x F i g . 3. S t e a d y r e a c t i v e o v e r d r i v e n w a v e w i t h y = 1.4, q = 1, ~ = 3/16. P r e s s m ' e v s d i s t a n c e a t t = 0 . 3 1 4 , f o r K = 1, 10, 100.
78
H.M. Sternberg
front and the steady piston driven wave, with the piston motion given b y eqns (2) and (20), is uniquely determined by the data at the front. Therefore, other steady CJ w a v e s with the same values of y, q, n and K are not possible since they would contradict the uniqueness of this piston driven steady wave. If one attempts to construct a CJ w a v e with a free or fixed b a c k boundary, or with a b a c k boundary consisting of a constant velocity piston moving at less than the CJ particle velocity, and if K is finite and n ~ 1, the resulting w a v e will not be steady and the front will not m o v e with a constant velocity, but will accelerate gradually and a p p r o a c h the CJ detonation velocity as t -~ 0o. It is assumed here that q and K are large enough so that the w a v e does not die out.
References Chorin Alexandre Joel (1977) J. Comp. Phys. 25, 253--272. Sternberg H. M. (1970a) Q. J. Mech. Appl. Math. 23, 77-99. Sternberg H. M. (1970b) Astronautica Acta 15, 359-369.