Applications of random-choice method to problems in shock and detonation-wave dynamics

October, 1979

APPLICATIONS OF RANDOM-CHOICE METHOD

TO PROBLEr~S IN SHOCK AND DEJONATIOtI-~/AVE DYNAMICS

by

T. Saito and I. I. G1ass

UTIAS

Report

No. 240 CN ISSN 0082-5255

APPLICATIONS OF RANDO~CHOICE METHOD

TO PROBLEMS IN SHOCK AND DETONATION-WAVE DYNAMICS

by

T. Saito and I. I. G1ass

StibndttedAugus~, 1979

October, 1979 _{UTIAS Report No. 240}

Acknowledge~ents

i

We wish to acknowledge with thanks the valuable discussions witb Prof. K. Takayama, Tohoku l)':p.iversity, Sendai, Japan, while he was on sabbatical leave at this Institute. The ass;i.stance received from Dr. A. J. Chorin and Dr. G. A. Sod in providing uswith their computer programs is very much appreciated.

The financial assistance received from the National Science and Engineering Research Cou!lCil and from the U.S. Air Force under grant No. AF-AFOSR-77-3303 is acknowledged with thanks.

Summary

Although successful numeri cal methods exist for sOlving problems in shock and detonation-wave dynamics, tnere is still a real need of developing new

techniques where the old methods fail to predict important flow properties. For exawple, it has recently been shown in Ref. lthat existing methods fail to predict the interferametrically measured isopycnics in regular and single Mach reflections (let alone complex and doUble-Mach reflections, for which numerical solutions do not even exist). The purpose of the present report is to present eight applications of the Random-Choice Method (RCM) to the solution of problems ~I in shock and detonation-wave dynamics. It is shown that unlike other numerical

methods, the RCM yields sharp-fronted shocks and contact surfaces without re-sort~ng to artificial and perhaps erroneous means of predicting their locations,

"I which depend more on art than science. It is also a very useful method in

showing such fine points as the birth point of the second shock (implosion) wave at the tail of the rarefaction wave in a spherical explosion.

Despi te all these advantages the RCM has yet to be deyeloped to cope wi th problems such as oblique and spherical shock-wave reflections in order to compute the various isolines (pressure, density and velocity) and compare them with available interferometric or other experimental data. For example, isopycnics are much more sensitive indicators of the accuracy of a given numerical method than a comparison of shock shapes (Ref. 1). Undoubtedly, such applications of the RCM will probably take place in the near future, as the need for such . numerical methods now exists.

- - - -- - -

-_

.. - ' - - -.' Tab1e of Contents Acknow1edgements Summary Tab1e of Contents Notation

1. INTR ODUCT ION

2. RANDOM-CHOICE METHOD (RCM) 2.1 General Description of RCM 2.2 Solution of Riemann Prob1em 2.3 Sampling Procedure

2.4 Production of Random Numbers 2.5 Boundary Conditions

2.6 Treatment of Contact Surfaces

3. NUMERICAL RESULTS FOR PLANAR WAVE INTERACTIONS

3.1 Ref1ection of Shock Wave from End Wa11 of a Shock Tube and Subsequent Interactions

3.2 Head-On Collision of a Shock Wave with a Rarefaction Wave 3.3 Head-On Collision of Two Rarefaction Waves

3.4 Contact-Surface Tai10ring

3.5 Shock-Refraction Prob1em at a Stationary Contact Layer 3.6 Shock-Wave Propagation in a Varying-Density Field 4. EXTENSIONS OF RANDOM-CHOICE METHOD

4.1 One-Dimensiona1 Symmetrie Flow 4.1.1 Gen~ral Description 4.1.2 Sod' s Method

4.1.3 Explosion of a Pressurized Helium Sphere 4.2 Reacting-Gas Flow

4.2.1 General Description

4.2.2 Chapwan-Jouguet Detonation

4.2.3 Chorin's Method for Reacting-Gas F10ws 4.2.4 Numerical Results

5. DISCUSSIONS

AND

CONCLUSIONS REFERENCES

APPENDIX A: CONTACT SURFACE MOTION APPENDIX B: PROORAM LISTINGS

TABLES: 1TO 6 AND Al

FIGURES: 1 TO 21 AND Al(a) Al(b) AND A2

iv Page ii iii iv v 1 3 3 4 6 8 9 9 9 10 11 12 12 13 13 14 14 14 15 15 18 18 18 21 24 25 26 Al B1 "

I:lo C D e F i I K m n p p r S t bot T u

u

v V w x Notation speed of sound contact surface

speed of detonation wave total energy per unit volume

function of V, fram Eqs. (1), (2) and (34)

~nteger attached to the mesh point in space inhamogeneous term in Eq. (34)

reaction rate

mamentum per unit volume molecular weight

mass flow (fluid enters a wave fram the left) mass flow (fluid enters a wave fram the right) integer attached to the mesh point in time pressure

sampling point radial distance

state in Riemann problem;

St'

Sr'

S*

time

time increment temperature

particle ve,loci ty shock velocity

solution of Riemann problem, see Eq. (5)

solution vector from Eq. (2)

particle velocity relative to wave front space coordinate

space increl!1ent

z

prQgress parameter -in Eq'. (66)

1, 2, -3 forp1anar. cylindrical,;'spherical :flows, respectively

y _{specific he,a.t ratio ,}

e: _{intern&! ellergy pe;r unit, mass}

e

I

\

density

",

.

'

1. INTRODUCTION

Over the past two decades a number of numeri cal methods have been developed

for the solution of flow problems in nonstationary gasdynamies involving

transition fronts such as shock waves~ _{rarefaction waves and contact surfaces.}

Such fronts occur in shock-tube flows~ _{cylindrical and spherical explosions and}

implosions as weIl as in combustion with deflagration and detonation. In

essence~ _{the set of nonlinea.r hyperbolic partial differential equations of}

motion are replaced by a set of finite-difference equations which are numerically

integrated for a solution. _{Owing to truncation errors an implicit artificial} viscosity is introduced which spreads the contact-surface and shock-wave fronts over several mesh lengths. This is analogous to the action of real viscosity and heat conduction in spreading the shock transition over several mean-free-paths,

and at the contact surface,due to the diffusion of heat and mas~ _{spreading occurs} with the square root of time.

Most finite-difference methods when applied to problems with discontim,llties sllch as shock waves produce oscillations behind them. Von Neumann and Richtmyer _{(Ref. 5)}

introduced an a.rtificial viscosity pressure term into the Lagrangian

form of the gasdynamic equations to get rid of these unwanted oscillations in the

solution. Since then various types of explicit or implicit artificial-viscosity terms have been used depending on the type of equations~ _{Lagrangian or Eulerian}

and the type of finite-difference scheme used. Shock-wave fronts are also smeared by the artificial viscosity term,as weIl as the truncation error~ _{typically over}

several mesh points. Moreover the artificial-viscosity term contains an

arbitrary parameter which must be determined for each particular proble~ _{in order}

to obtain the best results. For a recent review of the subject see Ref.

4.

Unfortunately, the spreading or smearing of such fronts makes it difficult to

know their positions with any precision at a given time. The development ef the Random-Choice" _{Method (RCM) by Glimm, Chorin and Sod (Refs. 2~}

3 and

6),

has

made it possible to overcome these difficulties at the expense of some randomness

in the paths of' these wave fronts. For example, shock waves and contacts~rfaces

are not spread at all and occupy zero zones. Their 10cI;Ltions at any time are not exact but their average positions are. The physical profiles of the head and

tail of rarefaction waves are perfectly sharp, whereas in other methods they come_{out rounde~.} .

Due to the randomness~ _{the profiles ar~ not smooth but on the}

average very close to the exact values for the rarefaction wave. In case of a shock-tube flow ~ the unifb'_{rm states separated by the contact surface 'are obtained}

exactly without oscillations. Boundary conditions are readily handled by the

RCM. _{The time required for the RCM may be two or threefdld longer than for the}

other methods. However, the number of grid points are fewer for the same resolution.

Consequen~~y, _{the computation time can be much faster.}

It is known that wave-interaction problems can be analyzed by using the method of characteristics. The solution~ _{in this Cl;Lse, is computed with the aid}

of a grid of characteristic lines~ _{which is constructed inthe course of the}

computation. This method is used mainly for a detailed description of the flows. For example~ _{the method of characteristics permits one te determine accurately the}

birth point of secondary shock waves in a flow at the point of intersection of the characteristics of one family. However, if a large number of such shock waves

occur~ _{difficulties would be encountered. Accordingly, the method of characteristics} is usually applied to problems where the number of discentinuities is small.

The RCM·does not resort to any type of finite-difference sCheme to obtain the solutions for wave-interaction prob1ems. A set of states at two adj acent mesh points i and (i + 1), at t ::;:: 0, form a Riemann prob1em (shock-tube prob1em) with the Eulertype of equations of motion, name1y, mass, .mamentum and energy. Fram the method of characteristics it is known that a Riemann prob1em has a se1f-simi1ar solution which consists of four uniform states consisting of the

·two initia1 states, two final states, separated by a contact surface and a nonuniform state, the rarefaction wave. All the thermodynamic and dynamic properties can be calculated exact1y fram an algebraic consideration of the transition re1ations for shock fronts and rarefaction waves (Ref.

34).

Af ter the Riemann problem .. has been solved, one of the five states is Chosen at randam. This Choice makes it possib1e to find a solution subsequent1y, at a time .6.t/2, at a random point P, which is located between the mesh points i and

(i + 1) as shown in Fig. 1. Although P is located in state (2), this solution is assigned to mesh point at (i + 1/2). It may appear that this is in error since the mesh point is really in state (3) in Fig. 1. However, if we realize that the mesh points are the representatives of e1ementa1 regions of the who1e flow field, '(i + 1/2) could have been at point P, i.e., the mesh points are not fixed in physical space. This randamness, on the average, is cancelled out and the correct solution is obtained. Although this sampling procedure gives randonmess in the wave posi tions and the shape of the rarefaction waves, i t is in reality aceeptab1e and the calculated values are exact due to the fact that the method.is free ·fram truncation errors and artificial viscosity.

Many wave-interaction prob1ems were investigated theoretical1y and experi-mental1y ·at .UTIAS. Fini te-difference schemes with artificial viscosi ty were

used mainly for solving the gasdynamic equations in Eulerian or Lagrangian form. In this report, the RCM is first app1ied to same of the practical wave-interaction prob1ems. Although the mathematical justification of the RCM has already been

shown in Refs. 2 and

3,

it is important to check the usefulness and the app1ica-bi1ity of the method to practical prob1ems.

So far, the RCM has been applied to planar, cylindrical and spherical f10ws (Refs.

3

and

6)

as we11 as to combustion f10ws (Ref.

7).

Same very important app1ications of theRCM would be to pseudostationary ob1ique shock-wave reflection prob1ems as we11 as nonstationary cylindrical and spherical

shock-wave reflections. This has not been ful1y accomp1ished to date.

In Chapter 2, the RCM is described for the sake of completeness. Prob1ems invo1ving different.gases with different values of th$ specific heat ratio

7

were not calculated previous1y using the RCM. In this report the program was modified in order to make it possib1e to solve prob1ems invo1ving combinations of gases with different values of

7.

The procedure is described in Appendix A. Chapter

3

deals with several examp1es of simp1e shock-tube f10ws solved by using the RCM. The. results are compared with exact or approximate solutions. Chapter

4

extends Sod':s RCM for cylindrical and spherica1 symmetrie f1ows. Chorin's method for calculating reacting gas f10ws is also brief1y described. The methods are applied to same practical examp1es, which required a modifica-tion of Chorin's program. Discussions and conc1usions can be found in Chapter

5,

and the program 1i'stings are gi ven in Appendix B.

2. RANDOM-CHOICE METHOD (RCM) 2.1 GENERAL DESCRIPTION OF RCM

The equations for an invisc1d non-heat-eondueting one-dimensiona1 flow ean

be written in conservation for.m (Ref. 3):

(1)

where

(2)

and the subscripts indieate differentiation. The total energy per unit volume e, may be written as

1 2

é

=

p€ +

'2

pu

Assw:n1ng the gas is polytropie ,the internal energy per unit mass is gi ven by

~

=

1 p

i - 1

P

".

(-4 )

Letting öt and ~ be time and spatia1 inerements, re spective1y , the solution

of the system of equations, Eq. (1), is to be eva1uated at the mesh points

(i~, nöt) , ~d ( i + 1/2)~, (n + 1/2)öt}. Letting v~ approximate V(i.6x,;

llÖt) "and

v~:î//~

approximate V [(i + 1/2)&, (n + 1/2) öt), it is necessary to find

n+1/2. n n

vi +1/ 2 g~ven _{vi' vi +1 ·}

Consider an initia1 va1ue prob1em for the system of equations given by Eq. (1), with the discontinuous initia1 data,

x ~ (i + 1/2).6x

V(x, nöt)

=

(5)

x

<

(i + 1/2)!$.

then Eq. (1) together with Eq. (5) is called a Riemann prob1em. A method for

solving such a Riemann prob1em wi11 be described in detail in the next seetion.

Let v(x, t) denote the solution of this R;iemann prob1em and let

Sn

be a value

of the random v:ariab1e e equidistributed in

T-

1/2, 1/2j. Then define

v~:îj~

by the solution of the Riemann prob1em at the point [(i + 1/2 + en)!$., (n + 1/2).tit

J

v~:îj~

=

v {( i +

~

+ en ) !$., ( n +

~

) öt" } (6)

name1y, at each time step the solution is first approximated by a piecewise

constant; i.e., Riemann prob1ems are formed with respect· to each pair of every

other mesh points. Time is then advanced exact1y, and the new va1ues on the mesh points in between those used to construct the Riemann prob1ems arechosen by

sampling (Fig. 1). The justification of this method for solving the system of

equations, .. Eq. (1), can be found in Refs. 2 and

3.

2 .2 SOLUTION OF A RIEMANN PROBLEM

The method of solving a Riemann prob1em is now described. Consider the system of equations given by Eq. (1) with the initia1 data,

x

>

0

V(x,

0)

=

( 7)

x<O

The solution at later t.imes wi11 consist of three states; a 1eft. stat~ .S,t, a

right state Sr, a.middle state S* with u

=

u*, p

=

p* and p

=

p* separated from

S,t _{and Sr by waves, Wr ,} W,t, which may be either shock waves or rarefaction waves.

A contact surface (dx/dt

=

u*) separates the g.as initial1y at x

<

0 fram the gas

ihitial1y at x > O. The va1ues of u and p are continuous across the contact

surface while p-and other tbennodynamic quantities, .in general, are not. The

contact surface divides S* into two regions with differing va1ues of densities

P*,t and ~r, but equal constant values of u* and v* (Fig. 2),

Using Godunov's iterative method modified by Chorin and Sod (Refs. 3 and 6)

p,

u, p at the sample point p(e~, 1/2 ~t), - 1/2 ~

e

~ 1/2, are determined as

fo110ws: Define the quantity

M = (p - p

)/(u - u)

r r

*

r

"*

(8)

If the right wave is a shock wave, using the jump condition across the shock wave, we ob tain

where Ur is the velocity of the right shock wave. From the Rankine-Hugoniot conditions one obtains

where M

= (

p) 1/2 rk ( P * ) r Pr r 'f'1 Pr p* - >1 p _r -( + 1 1 .)1/2

rp

('I)

== 2'

'I

+

r

-1 2 2 ( 10) (11)

If the right wave is a rarefact,ion wave, using the isentrapic ],aw pp-

r

= constant

and the constancy of the right Riemann invariant rr

=

2(1P/p)lj2/(r -

1) ~ u, we

find

4

where The f'unction ~~1 . 1}

<

1 1 (12) (13) (14)

is continuous at 1}

=

1, with <1>(1)

=

~l-f1)

=

<1>2(1)

=

r

2 • Simi1ar1y, we de fine . M.t

=

CP.t - p*)

I(

u.t - u*)

If the 1eft wave is a shock wave,

M.t

=

P.t(u.t - U.t)

=

P*.t(u* - U.t)

where U.t is the velocity of the 1eft shock wave. As on the right,

(15)

(16)

where <I>(~) is defined as in Eqs. (11), (13) and (14).

Eqs. (8) and (15), we obtain E1iminating u* from

Equations (18), '(17) and (10) or (12) are three equations with three unknawns p*, M.t and Mr' These considerations lead to the fo11owing iteratiop procedure: Choose a starting va1ue p'; (or ~, ~), and then compute

M.r

V₊1 , _{M.tV+1, p*v+1}

(v~O) using ,

(19)

(20) (21) (22) 5

,...

Since there is no guarantee that P remains ~O, Eq. (20) is needed. The iteration is stopped when

Then put

M,e

=

M~+l,

M_r

=

M;rl,

p*

=

p~+l

Chorin picked the value of El = 10-6 , E2 = 10-6 in this i teration cycle. To start this iteration cycle either Mr and M,e or P* is needed. From the pqint of view of computation time, the efficiency of the RCM depends on how fast this iteration converges. Chorin obtained better results putting

.as a starting value of the procedure than following the starting scheme suggested by Godunov (Ref.

8).

Godunov mentioned that the iteration may fail to converge in the presence·

of a st rong rarefaction wave. This problem can be overcome by the following

variant of Godunov's procedure. If the iteration has not converged after J iterat-ions, Eq. (20) is replaced by

~1 ~ V

p*

=

_{S max(E1 , p ) + (1 - S)P*} with S = SI =~. In general S will be reset

1

S

=

Sj = 2" Sj-l

after jJ iteration, j

=

1, 2, 3, '," . Chorin noted that the cases jl>2 had never been encountered and the number of i terations required fluctuated between 2 and 10, except at a very few points.

Knowing

M,e,

M

r , p*, Eqs.

(8)

and (15) give

u*

=

(P,e - Pr + M,eU,e + _{MrUr)/(M,e + Mr )} (23)

2.3 SAMPLING PROCEbURE

Having a solution in the middle state S*, the next step is to determine the solution at each mesh point by sampling. There are four basic cases to be sampled.

1)

The sampling point P = (8~, ~

At)

of the contact surface whose inverse plane is (dx/dt)

=

~, and the right i. e. ,

lies to the right slope in the (x,

t)-wave is a shock t)-wave,

As mentioned before 8 is a random number uniformly distributed over the range ·of [-~, ~].

2) The sampling point P lies to the left of the contact surface and the left wave is a rarefaction wave, i.e.,

•

8.6X

<

U*!J.t/2 and p*

<

P

t

3) The sampling point P lies to the left of the contact surface and the left wave is a shock wave, i.e.,

8.6X

<

u*!J.t/2 and P* ~ P

t

4)

The sampling point P lies .to the right of the contact surface and the right wave is a rarefaction wave, i.e.,

-In w~at follow p, u, P denote the sampled solution, and now how they are determined is explained for each case.

Case 1) From Eq.

(9),

the velocity of the right shock wave U_{r can be obtained as}

u

=u +M/p

r r r r ( 24)

If the sampling point P lies to the right of the shock-wave line dx/dt

=

_{Ur ,} region 1 in Fig. 3,

-

-

-p

=

Pr' u

=

_{ur ' P}

=

Pr

If the sampling point P lies to the left of the shock wave, region 2 in Fig. 3,

-

-

-P

=

P*r' u

=

u*, P

=

P* where ~r is obtained from Eq.

(9)

as

P*r

=

_{-Mr/(U* - Ur )} ( 25)

Case 2) The rarefactio~ wave is bounded on the left by the line dX/dt

=

ut-at, and on the right by dx/dt =' ~ - ~, where at and ~ are the sound veloeities in the states

St

and S*, respectively; at is given, as

(26) and a* can be found by using the Riemann invariant rt

=

constant,

If the sampling point P lies to the left of the head of the rarefaction wave, dx/dt

=

Ut -

at, region

5

in Fig. 3,

P

=

Pt' u

=

Ut'

P

=

P_t

If the sampling point P lies to the right of the tail of the rarefaction wave, dx/dt

=

~

-

~, region:;'3 in Fig. 3, set

-

-P

=

P*

t'

u

=

u*' P

=

P*

/

where ~,e -i.s. calculated using the isentropic relation pp-r

=

A (constant), as

(27)

The value of A can be calculated in the state S,e as

.

1/,/

'

A=pp A'

.

,

*

,e' JfI (28)

If the sampling point P lies inside the rarefaction wave, reg~on 4 in Fig.

3,

the slope of the characteristic dx/dt

=

u - a is equated to ~he slope of the line connecting the origin and P" to obtain

.

ü -

a = 28!5x./ t::.t

Also,

,

2a/(r,e, ~ 1) +

u

=

2a,e/(r,e - 1) + u,e

=

rt

SolvingEqs. (29) and (30), obtain

Since a

- 2 ( 28.6x 1,e

~l

)

u = r,e + 1 t::.t + a,e + 2 u.t . = (

_r,ep

-/-)1/2 _{p and}

_:p

₌

_Ar?,e

2

l/ft;t';"l)

p

=

(a

/r,eA)

-

-r,e

l'

= Ap

(29) (31) (32a) (3Zb) Cases 3 and 4 are mirror imag~s of cases (1) and (2) and the same arguments are applied.

2.4 PRODUCT ION OF RANDOM NUMBERS .

The choice of a series of random nunibers 8n (n

=

1, 2, .•• ) det'ermines the behaviour of the solution. If 8 is close ,to -1/2, the values in the left state St, propagate to the right to ( i + 1/2)!5x., (n + 1/2)t::.t), while, if a is close to 172, the values in the right state Sr propagate to the left. Therefore it is important tO 'choose 8n in such a way that they tend to be equally distributed over [-1/2, 1/2} as soon as possible.

Chorin noted that i t is unreasonable to choose a new value of 8 for each mesh point i, at each time level

n,

since there is a finite probability ,that a givenstate S will propagate in both directions. Chorin overcame thisproblem choosing a new ,value of a once for each time level rather than assigning a nèw value of 8 for each i and n. Chorin further improved the method by making the

sequence of random numbers an' reach equidistribution over [-1/2, 1/21 at a faster rate. This is done by combining a random number 6n chosen from the range of [-1/2, 1/2J and a pseudorandom number kn0

Let ml and m2 be two mutually prime integers (ml

< mz) then consider the

sequence of integers,

where

ko

is given,

ko

<

lIl2. This will produce a series of pseudorandOlIl numbers. For example, if kc

=

1,

~l

=

3, ~

=

7,

then kl'

k2,

k3' ••• will öe a repetition of a series of integers,

4

.

, 0,

3, 6,

2, 5,

1. If we consider the following

sequence:

k +

(Sn

+ 1/2)

1

e'

= ... n _ _ _ _ _ _ n m... -

'2

, ~ (n

=

1, 2, ..• )

e~ willbe alSQ in the range of

[-1/2, 1/21.

This modif,ied sequenc~ of random , numbers

en

is employeq. in th!;! program. The advantages of using this sequence

en

is described in detail in Ref. 3.

2.5

BOUNDARY CONDITIONS

{ ,

,

Consider a solid-wall boundary at x = b, wi th the fiuid to the left. The bO\mdary conditions are iJlIPosed on the grid point closest to x

=

b, say io~' ~ pseudo right stat~ _{Sr at (io} + 1/2)~ is created by setting

Pi +1/2

=

Pi

-1/2

Q 0 '

In this manner, waves can be reflected at asolid boundary. 2.6 TREATMENT OF CONTACT SURFACES

\

(33)

The Random-Choice Method was not applied to flow problems having combina-tions of gases with ,different specific heat ratios. In this report the program was. modified to handle such problems, thereby enlarging the applicabili ty of the RCM. The position Of contact surfaces are determined during sampling procedure

and consistent with the RCM. The details are given in Appendix

A.

3. NUMERICAL RESULTS FOR PLANAR WAVE OOERACTIONS

The purpose of this section is to show how well the RCM works for nonsta-tionary wave-interaction prob+ems by using several illus'trative examples. Many

ex~les of interactions of shock waves, rarefaction waves, contact surfaces and with solid boundaries are giv!;!n in Refs.

9

and 10.

In order to illustrate the usefulness of the RCM, the fOllowing èxamples

we re calculated: '

1. Reflection of shock wave from end wallof a shock tube and subsequent interactions •

2.

3.

4.

5.

Head-on collision of a shock wave with a rarefaction wave. Head-on collision of two rarefaction waves. _J

Contact~surface tailoring.

Shock-refraction problem at a stationary ',contact l~er.

6. Shock-wave propagation in a vary±ng~density field.

3.1 REFLECTION OF SHOCK WAVE FROM END WALL ,OF A SHOCK TUBE ~D SUBSEQUENT

INTERACTIONS. '

The initial conditions were set in such a way that we could compare the results obtained using the RCM with those obtained by Gurke ,~~ Schwarzkopf

URef. 11), i.e.,

diaphragm pressure ratio ': P41::: 2t:0 test gas - air

driver gas - air

: ' ?'l ::: 1.40

': ?'4 ::: 1.40

Figure 4 was drawn from the data in Ref. 11. Since these results were obtained With the method of characteristics and there is no interactions between shock waves and rarefaction waves, the solutions obtained here are exact. Figure 5

shows our results obtained with the RCM for a mesh number of 720 along the axis of the shock tube. Only the region close to the end wall is dispàäyed. Since this method does not make use of an artificial viscosity, shock waves and contact surfaces remain perfectly sharp. Owing to this outstanding feature, we can

recognize as many as 68 different regions ,without difficulty in this ,particular case. The exact wave traj ectoties in Fig. 4 are superimposed in Fig.

5

as solid lines for ready comparison with the present results. It can be seen that the agreement in wave positions is excellent. At t:::: 2.5 ms, the error' in the position of the contact surface is about 6 mesh sizes, and appears to be the largest error in the regio,n where the exact solutioIi is presented, for comparison. This error corresponds toabout 0.8% of the whole mesh number, i.e., 1.3 cm for a 180-cm long sho'ck tube. The d:iiscrepancy of the wave"positions from the exact solution depends on the mesh si ze and was about 3% when the mesh number was

180. In Table 1, 'values of pressure P, density p, and particle velocity u, obtained by the present method are compared with the exact solution for each region. In '

these tegions, the numerical errors in both pressureand" density are of the order of 10- . , Although they tend to increase as.time goes on, their magnitude~ are small enough to'consider the numerical values as exact. In Table 2 n~erically obtained values or' p, p, u are listed for all the regions which appear on Fig. 5. In Fig.

5

and Table 2 we see that ~he flow patterns get more and more complex and also, as one can expect, the changes in values of p, u, p become smaller. For example, in regions 67 and 68, acrossthem the changes in pressure and density are from 24.690 to 24.694, and from 8.719 to 8.721, respectively.

All the numerical computat~ons were done using an IMB 370. It took 77.78 min to calculate this exampl~ with 720 mesh numbers and 1600 time steps. This fairly long computation time includes the time spent to calculate,the solutions

at mesh points far fram the end wall in which, however, there is no immediate interest. Therefore, if we could simulate wave interactions only in the region close to the end. .. .wall, significant computation time could be saved.

This can be done as follows. First, calculate the values of p, u, p in the regions 2 and 3 for the given initial conditions by usingeither Rankine-Hugoniot relations or this numerical method. Then distribute states 1, 2 and 3

as the initial conditions over the mesh points and by removing the boundary conditions at the other end wallof the shock tube so that no reflected rarefac-tion wave is generated there. This corresponds to a shock tube with an infinitely long high-pressure section. In the same way, one can alw&ys simulate anly those wave interactions that one is interested in, that is, a larger dynamic range for the wave interactions for the same mesh number. It took 3.19 minutes to calculate the same problem with a mesh nunber of 180 and 300 time steps, or from t

=

0 to t = 3.6 ms in this particular case. Many wave interaction patterns are seen in Fig. 5, such as the overtaking of two similarly facing shock waves, the head-on collisions among shock and rarefaction waves, the interaction of a shock or rare-faction wave with a contact surface. Same of these basic interaction problems were simulated in detail and will be discussed subsequently.

3.2 HEAD-ON COLLISION OF A SHOCK WAVE WITH A RAREFACTION WAVE

This interaction problem was investigated by Gould (Ref. 12) both analytically and experimentally. The initial conditions provided for a shock strength PlO' and rarefaction strength P20 of PlO

=

1.96 and P20

=

0.549. This is achieved by

putting the left diaphragm pressure ratio P70

=

4.15 and the right diaphragm pressure ratio P08

=

3.26 with

r

=

1.40.

Figure 6 shows the resulting wave system for this case. As a result of the head-on collision of a forward facing shock wave (fluid particles enter the wave from the right) and a backward facing rarefaction wave (fluid particles enter the wave from the left) , there appears a forward facing transmitted shock wave and a backward facing rarefaction wave. Since the shock wave increases in strength con-tinuously dur ing' the interaction with the rarefaction wave, each fluid particle crossing the shock wave in this period of time will experience a different

entropy jump thereby forming a contact layer or region. Consequently, the region between the transmitted shock wave and the transmitted rarefaction wave consists

of two new unifor.m regions (3), (4) and a contact region. (Note that secondary interactions of characteristics within the contact layer have been neglected.) Pressure p, and particle velocity u, are the same for regions (3) and (4), but density p, temperature T, and entropy S, are different. The exact values of p,

u, p in the regions (3) and (4) can be determined analytically. The wave trajec-tories, however, in the interaction region were determined by Gould in detail, as outlined in Ref. 12. In Fig. 6, the wave tr~ectories were obtained only by averaging the initial and final stage of the interaction without the detailed method of Gould. The results (circled data) obtained by the RCM are superimposed on those obtained by Gould in Fig. 6. In Table 3, the camputed values of p, u, p are compared with those obtained by Gould. As we can see, the calculated values are correct to at least the first four digits and the error is of the order of 10-

4 .

The agreement in the wave trajectories is also very good and the error is about 1% at the mos t .

The computation was carried out with 180 spatial mesh points and it took 0.72 minutes to proceed through 82 time steps, which was long enough to simulate the interaction process appearing in Fig. 6.

In Table 3, as mentioned before, the pressure p and particle velocity u are the same for regions (3) and (4) but the density p is different. The density change across the contact region is as smal 1 as 0.19% of its absolute value.

3.3 HEAD-ON COLLISION OF TWO RAREFACTION WAVES.

'l'his problem was solved by Steketee (Ref. 13). The computer ,simUlation, wal? done byassuming an experiment using a two-diaphragm shock tube. Initially the center part, which -is separated by two diaphragms both from the left part and rig~tpart of the shock tube, is ,filled with a gas at higher pressure than the other two sections so that, after the two diaphragms are ruptured, rarefaction waves propagate.intothis center section and collide. The assumed initial conditions' 'for the numerical caiculat-ions were 1eft section: centre section: right section: Ar (-y

=

1.6667): H2 (7

=

1.4000): Ar (7

=

1.6667): This will give the forward":facing rarefaction the backward-facing rarefaction wave strength

p = p = P = wave

Plo,

0.088118, p = 1. 7457 1.0000, _P = 1.0000 0.024793, p = 0.49119 strength P20' of 0.59498 and of 0.39945.

Af ter the interaction, a uniform state appears between the two transmitted rare-faction waves. Unlike the previous case tpe interaction process is isentropic

and no contact region will be produced. The values of p, u, p in this uniform

region can be determined exactly. The exact values and the results obtained from the RCM are displayed in Table 4. Again the errors in this calculated values are very smalle Figure 7 is the

(x,

t) -diagram for this interaction. The circles are the results obtained using the RCM and lines are those obtained by the method of characteristics (Ref. 13). The agreement is, in gèneral, good especially before the interaction. The agreement of boundaries between the state (3) and transmitted

" rarefaction waves do not look too good. This relatively large error, typically about

3%, includes the error associated with the method of characteristics and we can expect smaller errors for the RCM itself. On the region of penetration, the characteristic lines are no longer straight, and their trajeciDries had to be determined graphically dividing the rarefaction waves into smal 1 segments. In ·

Fig. 7, trajec~ies of characteristic lines were obtained by dividing both rare-faction waves into foursegments. As far as the wavè trajectorles are concerned, since it is not 'pr&ctical to divide the rarefaction fans into too many segments, the method of characteristics is also approximate in this problem.

I t took 1.72 min to compute 180 time steps wi th 180 spati.al mesh numbers. 3.4 CONTACT-SURFACE TAlLORING

In general, when a shock wave reflects at the channel end wall and collides with the oncoming contact surface, the shock wave is transmitted and either a shock wave

or a rarefaction wave is reflected. However, by choosing the initial conditions properly, thElf 'proride a contact surface such that the reflected wave is a Mach wave. As a

result, the gas from the end wall to the contact surface remains in a uniform or tailored state. Figure 8 shows such a condition where the wave between states (5) and (7) is a Mach wave. Since state (5) is stationary (u5

=

0), states (7) and (8) are also stationary (u7 = u8 = 0) and the contact surface is brought to a complete stop. When the flow is tailored, one can get a relatively long test time in the stationary region, which is very useful for aerodynamic and chemical-kinetic studies.

Tailoring in this case was produced by using two different gases with

different specific heat ratios. This problem requires a high degree of accuracy for

a solution and becomes very difficult for those methods with a finite difference

scheme in which contact surfaces cannot be weIl defined. Given the initial test

gas conditions, the diaphragm'pressure ratio P4l and the composition of the driver

gas can be determined by assuming the incident shock strength P2l (Ref. 14). Using

argon as a test gas (rl = 1.667, MI = 39.94), the diaphragm pressure ratio P4l

and the molecular weight of the driver gas (assuming a diatomic gas with

74

=

1.400)

M4 , were calculated for P2l

=

12.30. The calculated results give values for

P4l

=

58.21 and M4

=

7.251 (P4l

=

10.57). This driver gas can be made by mixing

79.87% of hydrogen and 20.13% of nitrogen, for example. The calculated values

and the results obtained by the RCM are compared in Table 5 and in Fig. 8.

It is seen the calculated values of P, u, P in each region are very accurate. The wave

trajec~es are in good agreement until the moment when the incident shock wave

reflects from the end wall, where it leaves with slight delay, causing a spatial

error of about 2%. Although it fluctuates back and forth one mesh size due to the

randomness of the sampling procedure, the contact surface stopped after its inter-action with the reflected shock wave, as expected.

It took 2.38 min to simulate 250 time steps with 180 mesh numbers. 3.5 SHOCK REFRACTION PROBLEM AT A STATIONARY CONTACT LAYER.

The interaction of a shock wave and a stationary contact surface was

investigated analytically by Bitondo, Glass and Patterson (Ref. 15) and experimentally by Ford and Glass (Ref. 16). Af ter normal reflection of a shock wave at a stationary contactsurface, there are two possible cases; a transmitted shock wave, contact surface and reflected rarefaction wave, or a transmitted shock wave, contact surface and reflected shock wave, depending on the incident shock strength and the initial internal energy ratio (acoustic impedance) across the contact surface. Putting a helium layer in air (Air" He "Air), both cases occur at the two contact surfaces.

In the simulation, the incident shock strength was 3.700. In this case, at the

left contact surface

Ct

(Air

I'

He) a reflected rarefaction wave appears and at the

right contact surface Cr (He' Air) a reflected shock wave results. The solution

for p, u, p can be determined exactly and they are compared with the computed results

in Table 6. The interactions are shown in Fig. 9 with the results computed by RCM.

In this particular case, the shock strength was attenuated after passing through the helium layer, from a pressure ratio of 3.700 to 2.971. For further details see Ref. 16, where it is shown that other interactions soon amplify the shock wave

to nearly its original strength. The order or errors both in the calculated values

in the flow variables and in the wave positions are the same as the previous examples. It took 0.70 min to calculate 140 time steps of this problem with 180 mesh numbers. 3.6 SHOCK-WAVE PROPAGATION IN A VARYING-DENSITY FIELD.

When a shock wave propagates in a gas which is stationary but has a certain density distribution in the direction of the wave propagation, the shock wave changes

its propagation speed and strength. For example, by cooling a vertical-channel

end-wall the test gas remains at constant pressure but has a nonuniform density.

Simulation was done with the following initial conditons. The driver gas was

air

(r4

=

1.400) and test gas was helium (rl

=

1.667). The density distribution is

uniform in the driver section. In the test section, it increases linearly such that P

at the end wall is fourfold PI at the diaphragm position. This density change

can be obtained by cooling the test gas from 2980 _{K at the diaphragm to 74.5}0 _K

The results are displayed in Fig. 10. In this particular case, the shock strength increases fram P2l

=

2.24, at the moment when the diaphragm was ruptured, to P2l

=

3.17, when the shock wave has propagated to the end waLL af the shock tube. Corre-spandingly, the l.ocal shock Mach number is also increased. However, since the .local speed of sound is reduced due to the density increase (temperature decrease);

a 0:

$lP,

the propagati.on speed .of the shock wave is reduced. The partic1e

veloci-·ties are alse decreased. The wave diagram is shewn in Fig. 11, where the exact selutien fer the case when the density p is initially uniferm, is superimposed for comparison. As expected, the shock wave and centact surface decelerate. It took 0.503 min te calculate 200 time steps with 100 mesh numbers. It is werthwhile mentiening that this example is difficult te handle using ether numerical methods. Technically, the difficulty arises from trying to keep the region ahead .of the shock wave statienary until its arrival in a finite difference scheme. The method of characteristics also has difficulties to treat this preblem since all the character-istics inside the nonuniform region are not straight.

4 . EXTENSIONS OF RANDOM CHOICE METHOD 4.1 ONE-DIMENSIONAL SYMMETRIC FLOW 4.1.1 GENERAL DESCRIPTION

Sed developed a method .of calculating ene-dimensienal symmetric flows such as cylindrical and spherical sheck waves by combining the RCM with .operator splitting. He gives seme numeri cal results on converging shock wave problems in Ref. 6. Here, Sod's methed will be applied te a spherical-explesion problem and the results will be cempared with seme results obtained by Brode (Ref. 26) using the artificial-viscesity method. The equatiens fer an inviscid, non-heat-cenducting,

radially-symmetric flew can be written in vector form (Ref. 6): V t + F(V)t

=

-I(V) where

(

m/r

)

I(V)

=

(a-l) m2

I

pr m(etp)/pr

where

a

= 2., for cylindric.al symmetry,

a

= 3, for spherical symmetry and

a

=

1, cerrespends te planar cases described in Chapter 2.

Exact selutiens exist for peint-source explosiens in the limiting case when the primary shock wave is infinitely strong (er implosiens) where the pressure in front .of the sheck wave is negligible, by camparison to that behind it, or in the aceustic limit, when the sheck is very weak (i.e., entropy changes are negli-gible) (Refs. 17 te 25). Although such solutiens are useful whenever they are applicable, they cannot be applied te real preblems of finite sources when the shock waves are neither streng na' weak. These solutions do net predict the

impertant features .of the flew that exists in an actual case such as the secend sheck wave and wave interactiens in the vicinity .of the centact surface. To ebtain the se feature s, Eq. (34) must be sel ved • In gener al, there are two maj or difficulties in solving Eq. (34):

1) A singularity exists at r

=

O.

2) The momentUili equation (these~ond component of Eq! (34)

j

caxlnot be expressed in co~servatiqn form.

Consequently, Eq.~J34) ha~ to be solved numerically. Brod,e (Refs."26 - 28.) was

able tG numerically integ:ra1;;e the nonlinear (inhomogeneous ) partial differential equations, ,Eq! . (3~) ,_ :for ~ I,l},lIllber of illlPort~t flow problems invo1.ving sphericü',·e'

explosions. payne, (Ref. 29) solved the" caSe of_ a cylindrical implosion\

Pl"

using "a "

Similar methoq.. Lapidus (Ref. 30) computed a cylindrical-implosion pról>lem in Ca,rtesian

co-ordinate ,in two, space dimensions . Althougl.l tl.lese methods have bèen· used

sucqess:fUlly to exp1ain physical phenom~na,. t~ey have inherent disadvantages in

that the use of artif~cial viscosii;;y smearsout the wave fronts over.several mesh

numbers and ta f~nd their exact 10cation becomes more of an art than a ~science.

*

By using Sod's ~ethod such difficulties inc1uding-1)aod 2) areeliminated:

comp1ete1y.

4.1.2 SOD'S METHOD.

In ·Sod's method, the first step is to remove the inhomogeneous term -!(:V)

from Eq. (34), usi,ng the methodknown as operator sp1itt.ing. Thus the system

(35)

is solved. In this system, the momentum equation cao be expressed in con-servation

form and the 'RCM, described in Cha;p:ter 2, cao be used to obtain solutions. Once

the ,system of equations" Eq. (35), has been solved, the system of ordinary

differential eq~ations given by

Vt=-I(V) (36)

is solved in turn.

In ,the program, Eq. (36) i,s solved using a Cauchy-Euler scheme by utilizing

the ,solutions of Eq. (35),

V,

,'

to determine the inhomogeneous term -I in Eq. (36),

namely, the 'system Eq. (36) is approximated by

or ₍₃₇₎

Since theso1utions for the' system Eq. (35) are Gn1y obtained at fntermediate'

points and the scheme of Eq. (37) does 'not require 'valuesat f = 0, the singularity

at the 'uis is eliminated. Theboun'daI.'Y' conditions are proper1y imposed on Eq.'

(35), as ;de~cribed in Sec~ 2 - 5. There is no need'to put boundary conditions on

Eq. (36), since i t is solved only at .intermedi~tepoints.

4.1. 3 EXPLOSION OF A PRESSURIZED HELIUM SI>.HERE.

A set of experiment al results for the explosion of pressurizedg1ass spheres

filled with air ,or heliUI!l into air may be.,found in Refs. 2p, 31 and 32. As an example of. Sod' s methed, one of t.he exp~:rimenta~ conditions,: for a helium explosiqn was used as initial conditions .,fÇlr the num~rical ~a.,lys.is ~

. .

diaphra~ pressure ratio: P41

=

18.25

diaphragm density ratio: P4l

=

2.523

initial particle.velocitie$: uI = u4= 0

(3B)

At the' insta.I}.tof.rupture"the planar-wave'conditions appIy, giving ,a·shoelt. stI,'ength ])21

=

6.497. The shock. wave then ,.c1ecelerates until ft becomes'-a s<:J1;llld' .

wave. The contact ' surface dece1erates .and i t.s motion becomes oscillatorr. The

rarefaction-wave head moves a~ constant speed á.nd reflects at the origin_:-.where it~ ·:

mot ion becomes complex. The tail of the ·rarefaction wave accelerates and an

impIQsion '(s~cond shock) wave originates. on it.

EquatiQn (34) may be reexpress~d in the fo~lowing form by using the method of

characteristics (Ref. 25).

Along a righ~ running characteristics or P-wave,

d.x

₌

u + a

dt

(40)

and a+ong a Ie ft, running characteristic ·or . Q-wave,

(41)

5Q +!:... ~ + (ex _ 1) au

=

0

5t C 5t l'

P

-( 42)

In tqe ,present explosion case, a Q-rarefaction wave results .' The values of Q

in the ,plane' case for.state (4) and (3) 'are 2~/(r4 -1) and 2a3/(r4 - 1) - u3,

respect i vely. In the limit, for a . complete, rarefaction wav~ when a3 .-+ 0,

u3 --+ 2a4!C,,4 - 1) and Q --+ -2a4/( 74 - 1). Consequent1y, Q decreases . fr om-the hea4

[Q

=

2a4!( 74 - 1)

J

to the ·tai1 of the raref'action wave~ Since the ,Q-r.arefa<;:tion wave is isentI,'opic, Eq .. (42) gives fR/5t_

=

-(ex - I)au/;-,that is., as one prQceeds '

along a . Q characteristic Q:decrease. This can. only happen i;f. tpe charact~ristic.

accelerates Jn the· direct ion of the tail of the ,rarefac~ion w~ve; In other wqrds,

~s time goes on, the Q-rarefaction wa'1(e continuesto get' stron~er untf1 i t is

engU+fed .by·the second shock wave.

The pressure profiles obtained by Sod' s method are dispJ,ayed in-Fig,~ ,12.

The numerical res,ults are drawn from theoutput of the computer without sinoothing (except in,Fig. 15, the explosion-wave ,diagram). The results exhibit the general characteristics of a Q-rarefaction wave for an explosion problem very weIl.

Ini tially the-planar rarefaction strengthp34 T 0.356. As time goes on~~ the

absolute values of pressures at th~ tail of therarefaction wave continue to decrease until,the second shock Wave i~ generated around the time number 10

which corresponds -to, 26.5!-ls aft er the glass sphere ruptures. At-time nUI11ber 10,

the rarefactio~ strength is increased to 0.142. The density profiles are shown in Fig. 13. They are 'lui te similar to the pressure profiles except for the discontinuity

across the contact surface. Again, the lowered density ratios or increased rarefaction-wavestrengths are observed. The formation of a second .shock wave and the, deceleration of the primary shock wave can be explained by considering the path 9f the contactsurface as similar to that of a piston. The contact s~face

is decelerated as a result of the spherical n~ture of the flow and sends outahea4 of itself rarefactionpulses (P-characteristics) that overta.ke and deca.y the ,primary shock wave. Howeve~, behind it, compression pulses (Q-cha,racteristics) are sent out which overta.ke to form a second shock (implosion) wave along the tail of the rarefaction wave. (It is worth noting that in the planar case such pulses run parallel to the tail of the wave, but in cylindrical and spherical flows, they

collide with the tail of ,the rarefaction wave.) This 'explains why the path of,second shock is connected with the tail of the rarefaction wave. Although the 'second

shock wave is a backward facing wave, it :isinitially very wea.k and propagates outward at first owing to the high positive particle velocity. However, as it gains strength, it overcomesthis counter flow and finally implodes on the origin with 'unlimited strength idealJ,y. In Figs. 12-a and 13, wecan see thatthe strength of the 'primary shock, Wave decays as it proceeds and how the second shock wave is generated. In the forego~ng figures at ab out 26.5 !-ls after rupture (time number 10), noticeable discontinuities were found that grew into a second shock wave, and were defined as the birth point of the second shock wave. The

particle-velocity profiles are shown in Fig, ,14. In the early stage of the explosion, all partieles move putward (positive velocity). As the contact 'surface decelerates, compression, Waves (Q-:cha,racteristics) and rarefaction waves (P-characteristics) decelerate the partiele velocity creating a.large negative velocity range. In this particular case, i t;l;'eaches' almost -0. 8al (at time number 46, 120 J..l s after rupture). This negati ve particle ve,loci t,y is reversed by the reflected second shock up to small positive values (time number 56). In Figs. 12, 13 and i4, the profiles of p, p and u have complex struqtures around the origin. This is

mainly-due to the 'randomness,of theRCM in the rarefaction wave. Unlike similarly-facing shock waves in the 'planar case, the second shock wave does not overta.ke the primary shock wave because of the lower sound speed and the lower or even opposing particle velocity behipd the p?imary shock wave; Therefore, the primary shock wave and the reflected second shock wave slow down'continuously to become Mach waves. (It is worth noting that thisis the principle behind travelling wave '

sonic-boom,simulators; ,Ref. 33.)

The Wave diagram is shown in Fig. 15. Isobars, isopycnics and isotachs ( constant velocity lines) are ,shown in Figs. 16, 17 and 18, respecti vely. In the figures, the decreases in pressure and density behind the incident 'and reflected rarefaction waves 'ar~ 'lui te apparent.. The pressures and' densi ties' in front of the imploding shock wave ne~ the origin are very small. Tbis imploding shock wav~

also induces large negati ve velocities as it approaches the -origin and a line of zero particle velocity is ,Been. I t is also apparent that the second shock wave is

produced, along the taU. of thera,refaqtion wavE;, prepagates outward at first, 'then

starts. t0 imp10de and increases in stre~gthto unlimited values 'as it hits the.i· origin .

and refJ,ects. .The origin becomes ·a singulari ty antl the c0ntinuum eciu~tions

break-down.. TransportprQpert.ies· would keep. the· thermod;Ynamic' propert_:i:es finite at the

implósl,0n f0CUS ..

REACTING~GAS FLOW

t .$

The RÇ}1'was applied ;to reactiI!g gas f10W by Ch0rir,r (Ref. 7) and wa:s' sho~ to be'

càpable of _ handlj,.ng t±jûe-depende~t . detoJ;lation 'and deflagration waves. wiihfinite

-and in fini te reacti0n rate. Chorin ,emp4asized an important ' advantage of usiI;lg

his method since" the interaction ef _ the- flow ~d the, chemical reaction can b~

taken into account 'whén the Riemann:problem i~ solved, even when .the time scales of

the 'chemistry .. and the fluid flow are very different; As aresult" the ',basic

conSierv~tio~ laws are' satisfied at the -end of e~ch time step. 11' ,the chemical

reactiens 'and the -gas.f10w were to be taken ,into account in separate fracti0nál

steps, the ,basi~ conse.rvation laY[s may be violated at the end of each hydrodyr).amic

st'ep, t~us either ,inducing unWanted ·oscilla:tions and waves,-or requiring ti-me ~.steps

e.maJ,l enqugh.fol;' all chaI).g~fl! to 'be very. gradual, usually a costly remedy. Chorin' s .

main 'oQject seemed to be placed on develo-pi~g the technique :andto: show ·the ,us'efullness

of. the Jnethqd. Th\l-s, inhis program, _ the:r,e . iS',M unr~asonable aasumpti0n.; io

=

1.1

=

7,

i. e." the. specific he'a~ ratios do not .change 'ac~oss . the ree.-ction front., Althoug~, as hest~tes',in his paper, the case

ra"

71 is 'more,diff'icult only beeause of"adait':t~n,al

algebra. We -feIt -tha:t this assumptionmust be 'impr0ved. Here, ,as a firststep, the-,

programwasmodi~_ie,d to have changes in the, spec:ific heat ratio ''So-that

10

f

~

at the detonation-wave front. The,RCM was then applied to -a ,detonation wave in a

2H2.+ 02 gas mixture.

4.2.2 CHAPMAN-JOqqUET DET?NATION

-A Chapman-Jouguet detenatien is describeq briefly in this section for the"

convenien,ce of explaining Cherib.' s' method of calculating reactiI).g-gas ' flow. The

equations to be solved are· (Ref~ 7 and 34)

where

1 2

e = p€ + '2 pu

Ei is the 'intern~l en~rgy per unit mass ~d expressed as,

18

(43)

(44)

( 46)

- - - ---

-€.

=

1 g

~ .., - 1 p ( 48)

q is the energy of· formatien which can be released thrQugh. chemical reacUén,... Ip,

this _present se_ction, it is assumed tha:t part of g is relea·sedin~tan~:aneQus-ly:/

in an infipitely· thin reEloction zone; Deneting .the unburned -and-the: burned gas 'with

o

and 1, reapectively, we haye

1 Po -€o

_

.

- + go

"'0 -

1 p 0 1 P1 + gl €1

=

-"'1 - 1 ~ ( 50)

AssU4lling that the-unburned gEloS .is . on the. right, and,. le~ting D be the velocity of the reacting zone,Eqs. (43) and (44) can be·written in ct':meervation form

_whe,re

w

=

u - D,

o 0

From these. rela~ions ene relii.d.ily deduces.

=

-M

if-

=

-(p - p ) /(,. - ,.)

o 1 0

where ,. is· the specificvelume., 1/ p. Conservation of energy is expressed by

1

€1 - €o -

'2

("0 - 1'1) (po + _{P1 )}

=

0 8uQstituting Eqa. (49)and (50) into Eg. (55), we have

(51)

( 52)

(53)

(54)

( 55)

~here 6. <gl - go· In the. (1"1. Pl)- plane t~e 'lines through ("0, ,Po) tangent ,to the curv~expressed by -Eq. (56) are called the Ray1eigh .line-s.--(Fig,. 19.). Their peints , _{of tangency, 81 and 82 are called the Chapman-Jouguet} · (-crr)' point;a'i,. ) The

uppe~ portion .of the éurve 'cerresponds to' d'eton'atio~'s; - the 'portio~' abov~ Sl

to strong detonations and the portion helow to weak detonations. The ~ower part of the curve corresponds to deflagrations. A portion of the curve is amitted.as it corresponds to physically unacceptable condi tions where

M2

< 0.

The veloci ty and strength of astrong detonation are entirely determined by the state of the

unburned gas in front of the detonation and one quantity behind the detonation just as for shock waves. Let Po, Po and Ua be given, as well as Pl' and assume the unburned gas lies to the right of the detonation, then from Eq. (56)

where and 2 Il_i y_ ,.;;;;; 1 1

=

----=-

_{7- + l '} 1 (57) i

=

0, 1 (58)

The states on the curve located between the CJ point Sl and the line 7

=

_TO correspond to weak detonations. In this region, a CJ-detonation is followed by

a rarefaction wave. '

In what follows an explicit criterion for determining whether a detonation will be astrong detonation or a CJ-detonation is described. It is shown in Ref. 34 that at Sl, the velocity IWl

t

= al where al = (n Pl/PJ)1/2 is the sound speed, i.e., a CJ-detonation moves with respect to the burned gas with a velocity equal to the velocity of sound in the burned gas. This fact is used to determine the density PCJ, velocity UCJ' and pressure PCJ behind a CJ-detonation. Fram Eqs. (51) and (54)

(60)

/

~herefore we get

(61)

Fram Eqs. (56) and (61) we have,

where

b

-Therefore

. .

-Given PCJ, PCJ

=

P1 = TJ.-l CM be obtained from Eq. (61). Since M

= -

P.l wl

and wl = -al' we find

The velocity of the CJ-detonation wave DCJ is found from

P (u - D ) _{o 0 CJ}

=

-M

(64)

Suppose ul' the velocity of the burned gas, is given. If ul

<

UCJ a CJ-detonation appears, followed by a rarefaction wave. If ul

=

UCJ a CJ-detona-tion appears alone, and if ul

>

UCJ astrong detonation takes place.

An outline of Chorin's method for calculating reacting-gas flows is described in the next section.

4

.2

.3

CHORIN I S MFJrHOD FOR REACTING-GAS FLOWS

In this section, Eq. (47) is replaçed by

€

=

1

E

+ Zq

., - 1 P (66)

where Z is a progress parameter for the reaction, and q is the total availab1e binding energy (q ~ 0). Z is assumed to satisfy the rate equation,

where K

=

0 if' T

=

g

< T

P - 0 K = K if' T =

g

>

T o 'P 0

(68)

To is the ignition temperature and

Ko

is the reaction rate. Chorin has solved Eqs.

(43), (44), (45), (46), (66), (67)

and

(68)

using the RCM and has shown that a Riemann type problem can be solved even when def'lagrations and detonations are included along with shock and raref'action waves. In the present case, the initial data f'or a Riemann problem has the f'orm of'

and

Sp} P

=

Pt'

p

=

p

t'

u

=

Ut'

Z

=

Z t) f'or x

~

0

-

l

s

(p = P , P = P , U = U , Z = Z ) f'or x

>

0

r r r r r "

When there is no chemistry (Ka =

0),

Z

=

constant and Eq. (69) reduces to Eq.

(7) and its solution is given in Chapter 2.. In case

Ka

f

0, right and lef't waves may now be CJ or strong-detonation waves as well as shock and raref'action waves. Chorin has incorporated these possibilities into the solution of' the Riemann problem and deveJ.oped a computer program to calculate tbis problem.

The state Sr remains a constant state; ur and Pr are must change at constant volume (and thus can do no work). can be f'ound by integrating Eqs .•

(67)

and

(68),

with z(o)

Zr + _{5Zr , 5Zr}

<

O. The new pressure is wri tten as

f'ixed. rhe energy in Sr The cha.ne;e,_{6Zr in Zr} = Zr and Z( ~/2)

=

(70)

In what f'ollows, the superscript new, is dropped. Similarly, Zj changes to

Zt

+ 5Z,e, and a new Pt is f'ound using the analog of' Eq. (70).

In S* the values' of' Z dif'f'er from the values Zr + EZr,

Zt

+ BZJ' Let

Z*.t

be the value of' Z to the lef't of' the slip line (dx/dt =~) and let Z*r be the

val~e of' Z to the right of the slip line. The dif'f'erence in energy of' formation

across the right wave is N

=

[z*r - (Zr + _{5Zr ) }q, and across the lef't wave it}

is

b.t

=

{Z*t - (Zt

+ fJll) Jq. Iteration Will be carried out on the v:alues

Z*t,

Z*r' 4l,~. In the f'irst iteration, E~. (19) to (22) are iterated with Z*r

=

Zr + 'öZ_{r ,}

Z*t

=

Zt

+ 5Zt, and thus l:::....t =

!::sr

=

o.

When this iteration has . converged, a new pressure p* is given, and new densities ~t, P*r can be f'ound f'ram

Eqs.

(9),

(16) or the isentrapic equation of' state. New temperatures T*~ = p*/~t, T*r

=

p*/ P*x' are evaluated, Egs.

(67)

and

(68)

are sOlved, and new values

Z*.t,

Z*r'

At,

b._{r are f'ound.} I f Ar

> 0 the right wave is either a shock or a raref'action

wave, and if'

b.r

<

0 the rIght wave is either a CJ-detonation f'ollowed by a rare-f'action wave or astrong detonation.

Let 'I.l.JE. be the velocity in S*. Gi ven ~,

!::sc,

we can f'ind the "Vel,oci ties

UcJt' UcJ_r bebind possible CJ-detonations on the right and lef't. If' ~ ~ UcJ r

the right wave is a CJ-detonation fo11owed by rarefaction, and if u*

>

UcJ r the right wave is a st rong detonation. The CJ state is unaffected by S* Tsince it depends only on Sr) and as far·as the Riemann solution is concerned it is a fixed state. If the right wave is a CJ-detonation, Mr is redefined as

where P_CJ is fram Eq. (65) and then

Mr

=

(PCJPCJ) 1/2 _{cJ>2 (p:JPCJ) , P*/PCJ}

~

1 (71)

cJ>2 is defined in Eq.

(13).

If the right wave ~s astrong detonation, from Eq.

(59),

where 2 2 _{(Pr -}

P*)(~l

Pr + P*) (cJ>3)

=

- - - 2

---.;.;~---.;~---(~

-

~12

) Pr +

(~12

l)p* -~o (72)

Simi1ar expressions occur on the 1eft. A second i teration starts with M.t,

Mr

fram the previous iteration, and written out in ful1, appears as follows:

V~O

U

*

v

= (

.

+ MV + MV )/(MV + MV)

P

_{t -}

Pr ~

Ut

_{r ur} ~ r where

if right wave

=

CJ-detonation,

otherwise,

if 1eft wave = CJ-detonation,

=

(P.e' P_{t ,} u.e) otherwise,

v+1 _ 1/2 V+1

M

r - Pr cJ>3(Pr~' Pr' P* ) if right wave = strong detonation,

- ( n

)1/2~(

V+1/) o therwi se - Pr-r ~ P* Pr