Numerical studies of a detonation analogue

Cranfield

College of Aeronautics Report No. 8704

March 1987

.^ DELFT

Numerical studies of a detonation analogue

by

L 4 ME 11987

J.F. Clarke, P.L. Roe, L.G. Simmonds

and E.F. Toro

College of Aeronautics Cranfield Institute of Technology Cranfield, Bedford MK43 OAL, England

Cranfield

College of Aeronautics Report No. 8704

March 1987

Numerical studies of a detonation analogue

by

J.F. Clarke, P.L. Roe, L.G. Simmonds

and E.F. Toro

College of Aeronautics

Cranfield Institute of Technology

Cranfield, Bedford MK43 OAL, England

ISBN 0 947767 57 6

£7.50

"The views expressed herein are those of the authors alone and do not necessarily represent those of the Institute."

Contents »

Pa

1. An analogue model for detonation-wave behaviour 1

J.F. Clarke

2. Four numerical methods tested on the analogue 9

E.F. Toro

3. A moving finite element solution of Pickett's Detonation analogue 48

L.G. Simmonds

4. General review of Methods applicable to the detonation problem 78

P.L. Roe

1

-1. An Analogue Model for Detonation-Wave Behaviour A General Conservation Law

Consider a general conservation law, expressed in global form, as follows:

x^ X2

^ ƒ G(x,t)dx = F ( X p t ) - F2(x2,t) + ƒ H(x,t)dx Referring to the sketch, the law states that

'F,^'P/-ot-)

^ f c >

n='^(^. >t)

the rate of gain of a property G in the space between stations x, ^ is due to a net flux F, - F^ of G into the space, plus the contribution from any sources of G that exist there; H is the source strength per unit of the coordinate x.

If F, G and H are all continuous in x, 2 equation (1) can be written in the form X, + h

3G

X, + h ƒ |Y(x,t)dx = F(x^,t) - F(x^ + h,t) + ƒ H(x,t)dx.

^1

^1

where x^ has been written as x, + h; letting h ->- 0 leads to the differential equation

9G ^ 3F u „ However, suppose that there is a discontinuity in G at x = x(t). In

2 -x ( t ) - -x^ : ^ ƒ Gdx + ƒ G dx d t x^ x ( t ) + x ( t ) - Xp

= [ G ] S M / + ƒ ) | | d x 3

X, x ( t ) +

where

[G] E G ( x ( t ) - , t ) - G ( x ( t ) + . t ) , 4 S E ^ x ( t ) . 5

If H is hypothesised to be continuous throughout the region, application of (3) - (5) to (1) in the limit as x, ->• x(t)-, X2 ^ x(t)+ at a fixed time t gives the following result,

[G]S = [F], 6 where [F] is the jump in F analogous to the jump in G defined in (4). Thus

any discontinuity in G must be associated with adiscontinuity in F, and must propagate at a speed S that is directly related to the jumps in G and F, as depicted in (6).

The Detonation Analogue

Pickett's (1985) analogue of the equations that describe the behaviour of real detonation waves can be constructed from the results of the previous

Section. Two "densities" G are involved. The first, written as p, is

subjected to changes in the x,t space as a consequence of the flux p. Since there are no sources of P present in the hypothetical system, p can be

identified with F; with H equal to zero everywhere (2) gives the differential equation

P, . p^ = 0 7 relating p to p. Any discontinuity in p must propagate at a speed S given by

(6), namely

S = [p]/[P]. 8 The second "density" is written as A; it is unconnected with any flux

but is influenced by distributed sources in the x,t field. In this case (2) gives

3

-\ = r(P,M, , 9

where r is the source term, and its possible explicit dependence upon both

P and X is recognized in (8). Any discontinuities in X must remain fixed in

X, as can be inferred from (6).

Finally, the set (7) - (9) is completed by specifying an "equation of

state" that relates P, A and p algebraically. Specifically, we choose

p = Jp2 + ^QA, 10

where Q is a constant that has a role as an "energy of reaction", with X

playing the part of a chemical concentration or progress variable,

Putting (10) into (7) gives*

P^ + PP^ + ÏQX^ = 0 . 11

If Q vanishes, (11) becomes an inviscid Burgers equation; as a consequence

our system has one of the key attributes of non-linear wave propagation,

albeit that the waves only propagate along paths

^ = P > 0. 12

With X interpreted, loosely, as a chemical progress variable one should

view (9) as an analogue of a species-conservation equation. Consequently x is

behaving as a Lagrangian variable or particle label and not as a spatial

coordinate.

The "reaction" or source term r(p,A) is to be switched on by the passage

of a shock. The shock must obey (8) and hence must follow a path x (t) such that

- d | = S(t) = [p]/[p]. 13

Recall that jumps in A are only permissible at fixed values of x.

(

that

As a consequence [A] at a shock path x (t) must vanish; then (10) shows

[p]

-

K P ' ]

= è {p_ + P^} [P],

14

where p is written for p(x (t)+,t). It follows that

dx^

- ^ - S(t) = i {p_ + P^}. 15

* Equations of this type appear as rational approximations to processes

that occur during weak-wave propagation through combustible gases: cf the

case for which A^ = p described by Clarke (1985).

4

-The source term will be chosen in the present case to mimic some

\/ery

simple high-explosive behaviourjthus

r(p,A) = 2(1 - A)^U(x ft) - X ) ,

₁₆

where U(y) = 1 , y > 0; = 0, y < 0 is the unit step function.

Explicit dependence on p is not a part of our present model, but could

easily be incorporated if required.

Note that, with A equal to zero ahead of the shock, we can assume that

the value of p there, namely p , is likewise zero; p equal to zero will then

follow from (10).

Hugoniot and Rayleigh Line Behaviour

The state equation (10) is the analogue of the Hugoniot relation

'V '•

p = Jp' + iQA . 17

Similarly (8) shows that, in view of the assumptions about A and p ,

.^ : P = Sp 18

is the associated Rayleigh line relation

jC

. The sketch depicts 'i^ : A = 0

and '^ : A = 1 with *i-lines typical of Chapman-Jouguet (CJ) and strong

detonation behaviour.

5

-The subscript -e indicates a final equilibrium value on the "all-burnt" Hugoniot

ly : A = 1. The frozen Hugoniot 7/- A = 0 is the locus of p , p pairs.

Special CJ values are indicated by an appropriate subscript.

The CJ point is found by equating the slope of "J^: A = 1 to the slope of

JLJ when P has the value P^rj' ^^^^

^^'^rc- ^ -

1

(.J = PeCJ = ^CJ = ^P-CJ ' 19

where the last result follows from (15) with p equal to zero. Combining (17) and (18) gives

PeCJ = ^CJ PeCJ = ^PeCJ ^ ^^ ' from which it can be seen that

S c J = Q ^ • 20 Consequently (19) gives

^eCJ = Q^ = ^P-CJ 21 and (18) gives

PeCJ = Q = ^P-CJ • 22 Structure of Steady-State "Detonations".

Fix a coordinate E, to the shock, so that t

C = X - x^ - ƒ S(t)dt = X - X5(t) ,

₂₃

where x is the shock position or time zero. In terms of the variables

E,

,t

equations (9) and (11) become

A^ - SA^ = r = 2(1 - A ) ^ U ( - a , 24

'^ + (p - s)p^ + m^

=

25

If the motion is steady relative to the propagating shock front both A ^ and p. vanish and S must be a constant.

6

-Equation (24) then integrates to give

A = 0, C > 0; 26a

S(l - A ) ^ = (C + S)U(C + S ) , 26b

where the step function U(^ + S) ensures that A remains equal to unity

from the point C = -S at which this value is first reached. The

steady-state version of (25) gives

P -

0,

O

0; 27a

ip2 - PS + èQA = 0, ^ < 0. 27b

(The constant of integration in (27b) is Jp ^ - p S, by evaluating the

left-hand side of (27b) at C = 0- where A = Q; but S = Jp_ and the result

follows.) Alternatively, (27b) gives

P = S + (S^ - Q A ) ^ 27c

since only the strong solution is needed here.

It is useful to reorganise the results in (26) and (27) a little,

as follows. First, (26b) gives

A = 1 - S"'(^ + S)2U(C + S) 28

from which the explicit result

A^ = -2S"'(C + S)U(^ + S) 29

derives. Note that A^^ 0 as C-^ -S at the equilibrium end of the reaction

region. Second, using (28) in (27c), the result for p can be put in the

form

p = S + S ' ' Q ^

{ U

+ S)2 + S^(S^Q"'- 1)}^ 30

when ^ is within the reaction zone, -S < C < 0. Thus

P^ = S " ' Q^ {a + S)' + S ^ ( S ^ Q " ' - l)}~hK + S) 31

and it is clear that p_->-0 as C-*--S, but only for strong waves, for which

S > Spj = Q (see (20)). For CJ waves, which have S = Q , p^- is constant

and equal to 1, so that p has a discontinuity at ^ = -S.

7

-Comments on Conditions Downstream of Detonations

The original pair of differential equations have two sets of

characteristics in the x, t plane, namely

dt

dx

dt

'8p'

32a

32b

where C and P represent the traces of wavelets (propagating along x-positive

only) and particle paths on the x, t plane. The existence of C and P allows

for communication to and fro in the space behind a shock-front; "downstream" on

a particle path P and "upstream" on a wavelet C.

Consider the situation near a steadily propagating strong, or overdriven,

wave, as in the following sketch :

<U«J4.1

^

*- s > cp

- ^ ^ ^

AJt-The whole of the detonation-wave (shock plus reaction zone, whose tail is at

X := 1) is accessible from downstream via the wavelets C. This is because the

value of p at A = 1, namely p , is given by (27c) as follows:

S +(S^ - Q ) ^ S, S > Q^

₃₃

and C has the shape shown in the sketch.

8

-If the wave is a CJ-wave, S is equal to Q^ and C-lines at A = 1 have

the same slope on the x, t diagram as the shock and its associated reaction

region. The situation is then like the one sketched below:

V - i

\«l

_{C:4 5t-#*A^}

_{'. 42*- -d^}

t£A

The shock/reaction-zone combination is now isolated from events that take place

in regions to the left of the A = 1 line. This will be true in all situations

other than those for which a shock is formed in the downstream regions; since

a shock will have sufficient speed to overtake the A = 1; p = Q^ boundary.

The sketch above shows a C line to the left of A = 1 consistent with the

existence of an expansion (Taylor wave) behind the CJ detonation.

References

Pickett, W. (1985) "Detonation in miniature". Chapter 4 in "Mathematics of

Combustion", Ed. J.D. Buckmaster, SIAM, Philadelphia, pp 133-181,

Clarke, J.F. (1985) "Finite amplitude waves in combustible gases". Chapter 5 in

"Mathematics of combustion", Ed. J.D. Buckmaster, SIAM, Philadelphia, pp 183-245,

9

-2. Four Numerical Methods Tested on the Analogue

2,1 Introduction

The ultimate objective in mind is the development of a mathematical model

and a computer code for simulating detonation phenomena in condensed media„ A

basic requirement of the numerical technique for the code is the ability to

accurately represent strongly shocked reactive flows in more than one dimension.

Since the classical finite difference techniques (e.g. MacCormack) are clearly

inadequate for this problem, as illustrated in this report, a very first task

is to survey the current literature in search for possible candidate up-to-date

numerical methods.

A careful process of selection/implementation/development of various

techniques is then required. This report is part of that exercise which can

be significantly speeded up by a suitable choice of a simplified test model.

Such a model is required to possess some of the essential features of the full

physical problem and to pose some of the crucial numerical difficulties likely

to be encountered in practical situations. If as bonuses the simplified model

resembles the Physics and has exact solutions then we are in possession of a

useful object.

The Pickett's detonation analogue described in detail in Chapter 1 fulfils

most of these requirements and is therefore adopted here for the purpose of

testing several numerical techniques.

2o2 Likely Solutions of Pickett's Analogue

Before going into computational details we briefly describe some of the

likely solution profiles of the detonation analogue and point out some of the

potential numerical difficulties. Pigs. 1(a)-(b)show possible solutions for the

variable 'density' against distance at a given time for the case of constant

boundary condition on the left. Fig. l(c)shows the reaction-progress variable for

any of the three profiles. The boundary function pLj(t) can be chosen in such a

way as to generate more complicated flow situations. In all the cases of Figs.

1(a) - (c) the exact solution for the reaction zone can be found. This will

be used to check the numerical solutions.

10

-2.3 Four Test Problems

For the numerical computation we consider four test problems. These

are illustrated in Figs. 2 - 5 .

Problem 1: The initial condition (IC) is the steady exact solution, the

boundary condition (BC) is Pj^ = p and parameters are chosen so

that p = S, the steady shock speed. The profile for p in the

reaction zone happens to be a straight line. (Pig. 2)

The idea here is to compute profiles at various times and superimpose the

exact solution for comparison. Two things are important, namely: (a) the

shock speed and structure, (b) the reaction zone structure, i.e. shape, width.

The ideal numerical method should convect the exact profile undistorted for

no matter how large the times ti, ta, etc. are.

Problem 2: This is in a sense a more realistic test problem. The reaction

is initiated by a flat-topped shock, with parameters chosen so

that the steady solution satisfies p^ = p = S as in problem 1.

(Pig. 3)

The difficulty here is to numerically 'grow' the steady solution (known

exactly) from a flat-topped shock.

Problem 3: This is a slightly more difficult problem. The reaction is

initiated by a flat-topped shock as in Problem 2 but parameters

are chosen so that an unsupported detonation results. Here we

have Pb = P T *^ Pe " ^*

^^^^'

^^

Here one is expected to model well both the reaction zone and the unsteady

'following flow'.

Problem 4: The intial condition is as in problem 2, but with parameters

chosen to give a curved profile for density in the reaction zone.

This is an overdriven detonation with p. > p > S. (Pig. 5)

The new difficulty here is the representation of the curved profile for

density and pressure (from the equation of state).

11

-2.4 Numerical Methods Tested

Here we report on results of three numerical methods, MacCormack,

Godunov and Random Choice (ROM) methods. I have in mind at least three more,

namely. Roe's Method, a new higher order random choice (RCM2) and a hybridised

version of this (RCM3). Work on the second group of methods has already

produced results and will be reported on later. Tn Chapter 3 a moving finite

element method for the analogue is dealt with.

We want to solve numerically a system of hyperbolic conservation laws,

with source terms, of the form

Ut + F(U)^= S(U) (la)

which for the analogue discussed in Chapter 1 one has

(lb)

with p = |(p^ + AQ) (Ic)

The unknowns of the problem are p and X and depend on space x and time t.

Each method is tested on problem 1 (Pig. 2 ) .

2.4.1 MacCormack's finite difference method

This is a 2 step method, second order accurate, well known and easy to

code. The first step of the method is to obtain a provisional solution U*

(predictive step) at a time level t*. This is done by, say, forward flux

differencing system (1) to obtain (source terms ignored)

The second step (corrective step) is to use the predicted values Ut to evaluate

the fluxes and advance the (corrected) solution to time level n+1 preforming,

say, a backward flux differencing and averaging data, i.e.

U =

p X

. F(U) =

Hp2 + XQ)

0

, s =

0

r ( X )

12

-with

U"+^ = U . - ^ I ' p * - F* 1 C3)

1 1 A X I 1 1-11 ^ ''

U. = |(U; + U*) (4)

When applied to Problem 1 the results are as shown in Fig. 6. Note that

the shock is smeared (about 6 zones),overshoots and spurious oscillations develop

and the region near the tail of the rarefaction is slightly distorted. Given

the large number of nodes (50) within the reaction zone the method still

manages to represent most of the rarefaction quite accurately. For the same

reason, wave speeds are also quite acceptable.

Since the shock is smeared, one has to take an arbitrary decision for

reaction initiation. Fig. 6 shows the case of an activation density, DA = 1.0.

Pig. 7 shows the case DA = 5. Improvement for X (the reaction progress variable)

can be observed. Prom these results it is obvious that this method is

in-adequate for the problem at hand.

Before we proceed to describe the other methods we illustrate the concept

of a Riemann problem, which is central to several modern numerical techniques.

2.4.2 The Riemann Problem

The Riemann problem (RP) is the homogeneous problem given by a hyperbolic

system with the special initial condition consisting of 2 uniform states

(Fig. 8) i.e.

D.E.: U^ + P(U)^ = 0 , Xi $ X ^ X2

I.e.:

if X ^

X Q

if X > X^

0

Given the simplicity of the data in the Riemann problem one can usually

solve the problem exactly, although not always in closed form. If general data

is approximated by piece-wise constant functions then one can transform the

general problem into a sequence of Riemann problems. In Pig. 8 one then has

the RP(i, i+1) in a computational set-up. By piecing together these Riemann

solutions after a fashion one can construct the solution to general problems.

13

-A practical example of a Riemann problem is a shock-tube problem. -A

shock tube is assumed, at t = 0, to be filled with 2 gases separated by a

diaphragm, where pressure, density and velocity have constant values pi, p^

and u-] = 0 on the left side and p^, p^, u^ = 0 on the right. (See Pig. 9)

If for instance p > p^ and p > p then the break up of the initial

discon-tinuity generates a wave system as depicted in Pig. 9, with a right travelling

shock, followed by a contact discontinuity and with a left travelling rarefaction.

The mathematical formulation of the shock tube problem leads to the

Euler equations with corresponding initial condition given by the shock-tube

data.

A Riemann problem is a shock tube problem admitting initial velocities

different from zero.

The Riemann problem for the Euler equations has exact solution for the

case of the ideal equation of state and can be found in the current literature

(Chorin, 1976). For the covolume equation of state, the Riemann problem can

still be solved exactly (Toro and Clarke, 1985), but for more complicated

equations of state the mathematical elegance is lost and iterative procedures

must be employed. One can also solve the Riemann problem approximately.

For the problem of interest in this report the Riemann problem is

simple, but the solution is not available in the literature. We solve it here.

2.4.3 Solution of the Riemann Problem for the detonation analogue

The homogeneous version (no source terms) of the governing equations for

the analogue can be written as

"t

* '^(^)h

" °

(5)

The Riemann problem is formed by adding the piece-wise constant data

U =

J

Ui , X .< X^

U , X > x„

(6)

9P

U is given by equation (lb) and A(U) is the Jacobian matrix -^

given by

~P IQ

0 0

14

-whose eigenvalues are

Xi = 0 , X2 = p (8)

The right eigenvectors are

(9)

P L Roe provided a geometric solution to the analogue's Riemann problem

(5) - (6). Here we give the algebraic solution, based on the admissible

(for the present problem) geometric cases.

Case A

Fig. 10a illustrates the Riemann problem in terms of curves of constant

'pressure' p on the p - X plane (the unknowns of the problem). Since we are

constrained by p > 0 and 0 < x ^ 1, only the first quadrant is of interest.

The states left and right (data of the Riemann problem) are given by

points L and R in Fig. 10(a). One can see that

p-, < P^ and p^ > p^ (10)

In addition one knows from the problem that X-, ^ X (always) i.e. the point

R in the Fig. cannot be above L. Point M represents a state that can connect

the given states L and R. The solid line shows the path from R to L through M.

The interpretation is that from R to M p falls smoothly while x is constant

i.e. p falls smoothly. We are in the presence of a rarefaction as shown in

Fig. 10(b).

Pigs. 10(c) - (e) show the behaviour of p, x and p, i.e. the solution of

the Riemann problem.

Since L and M are connected by the condition of constant 'pressure' p,

Pi = H p f + x^Q) = i ( p ^ + x^Q) = p^ (11)

and R and M are connected by the condition of constant A

ei =

1

ip

Q

62 =

1

0

m r

(12)

15

-one concludes that

Pn, = [pf + QC^i -

\)]^

(13)

Thus the state M is completely determined by equations (11) - (13).

Within the rarefaction the solution is continuous and given by

,X, X

(14)

The solution of the Riemann problem for case (A) is thus complete.

Case (B): P^ < P^ , p^ > P^ , x^ > x^

The complete solution is depicted in Fig. 11. The solution consists of

two waves, a 'contact' given by -^ = x^ = 0, as in case (A) and a shock given

by -jr = S with

S =

i

. (Pr - Pm^ ^.

(15)

From the equation of state for p one obtains

S = H P n , + P^) (16)

where

%^ rp^M(x^ - x^)j^ (17)

(X = X and p = p , always). Hence the solution is known everywhere.

Case (C); p.^ > P^ , p^ > p^ , x.| > x^

This case is essentially the same as Case (B) and is represented by the

points L', M and R in Fig. 11.

Summarising, the solution of the Riemann problem consists of 2 waves

(2 real eigenvalues

X],

A2). The left wave H T = '^1 ^ ^ ^^ ^ 'contact'across

which p and A change discontinuously and p is constant. The right wave

rl Y

^ = \2 "'S either a rarefaction (case (B)) or a shock (Cases (B) and (C))„ The

complete solution may be expressed as follows

16

-p„,= [pf * Q(Ai - X^)]^ > p ^ (18)

If P^ > Py. then Pjj, > p^ (shock) (19)

If Pi ^ Pp then p^ ^ p^ (rarefaction) (20)

Now we are in a position of describing two more numerical methods,

2.4.4 Godunov's Method

This 2 step method was first proposed in 1959 (Godunov 1959) and uses

the Riemann problem solution to obtain provisional values for flux evaluation.

Fig. 12 is a physically motivated description of the method for a 3 x 3 system.

Step 1: At time level n one solves all Riemann problems RP(i, i+1).

For each one we otbain values U. f that denote the solution of the Riemann

problem at X = iAX and t = (n + i)AT. This is asigned to the point

(X,t) shown by the squared nodes in Pig. 12. This step generates a secondary

grid with values at the intercel 1 boundaries at half time levels.

Step 2: Fluxes F(U) (eq. (lb)) are now evaluated on these secondary

grid points, at time level n+|, to give

F J : | = FdJ?:!) ^'^^

Finally, the solution at time level n+1 is obtained via the leap-frog

type scheme

,,n+l _ ..n _ ATTpn+l _ ph+i]

,^^.

"i " ^i AX |_S'+i S"-y

^'^'^'

Figures 14 and 15 show computed results of Problem 1 using this method

using 50 nodes again for the reaction zone. No overshoots are produced by

this scheme and smooth parts of the flow are quite accurately represented, but

discontinuities are badly smeared. Shocks are smeared, in this problem, to

about 5 zones and the discontinuity in derivative at the tail of rarefaction

is also smeared into several zones, and gets worse as time increases. Being

a conservative scheme shock positions (averaged) are correctly computed,

17

-2.4.5 The Random Choice Method (RCM)

This method was first successfully used in practical computations by

Chorin (1976). A readable account can be found in Toro and Clarke (1985).

The method can be implemented in a single or double grid. Here we

briefly describe its implementation on a double grid„ There are 2 steps

(see Pig. 1 3).

Step 1: Solve all the Riemann problems. RP(i, i+1) as done for the

Godunov's Method. Now, in order to extract a value u " ^ from the RP(i, i+1),

say, we point out that at time t = (n + | ) A T and (i - J ) A X -f X ^ (i + J'AX

there are infinitely many solutions, depending on the X-value. The RCM picks

a solution at random in the given range for X. The procedure relies on some

statistical properties of a sequence of random numbers {e|^} in the interval

[0,1]. For details see Toro and Clarke, 1985.

Having obtained a value U*?^? from the solution of the RP(i, i+1) one

assigns it to points (X,t) of a secondary (staggered) grid represented by

the squared nodes in Pig. 13.

Step 2: Now the data is given by the squared nodes (time level n + J ) .

To complete the solution to time level n + 1 one carries out the same procedure

n j_ 1

of step 1. Note that the solution U. is given at the same X-position as

n

for U-, i.e. on the original grid formed by the circled nodes.

This method has some remarkable properties. Pig. 16 shows the computed

solution of a shock tube problem (Euler equations). For details, see Ref. 2.

Note that all discontinuities (shocks, contacts and discontinuities are

derivative) are perfectly sharp. Some randomness inside the rarefaction can be

observed.

2.4.6 Source Terms via Operator Splitting

The complete problem is represented by the non-homogeneous system (la)

i .e.

U^ + [F(U)]^ - S(U)

with the source term given by (lb), i.e.

0

S(U)

18

-The operator splitting technique proceeds in two steps

Step 1 Solve the homogeneous part of the problem with data at time level n i.e.

to get a provisional (predicted) solution Ö at time level n + 1 in a time U^ + [F(U)]^ ^ 0

a provisi step of size AT.

The solution method to obtaine O""*" may be based on the solution of the (homogeneous) Riemann problem as presented in Section 2.4.3 e.g. the Random Choice Method.

Step 2 Solve the problem U^ = S(Ö""'^)

which is an initial value problem for a system of ordinary differential equations. A sophisticated machinery may be put to use to solve this problem

step 1),

problem with data Ü""*" (obtained from solving the homogeneous problem of The simplest solution to this initial value problem is via the Euler (1st order) method, i.e.

U"^l = Ü"^^ + AT S(Ö"+1)

The complete soltuion (corrected) to the inhomogeneous problem is U"^^

For complicated rate equations one may have to use higher order solvers in this step, e.g. Runge-Kutta or even stiff solvers.

For the results reported here we simply use the 1st order Euler method.

Some computed results are shown in what follows.

For Problem 1 of the analogue. Figs. 17 and 18 show RCM results. Apart from some randomness in the shock positions (small) and randomness in the reaction zone (large for the 'pressure') these results are quite satisfactory. Some

19

-2.5 More computed results

Here we show more computed results.

Fig. 19 are results for Problem 2 (see Fig. 3 for definition) using

Godunov's Method. Here reaction is initiated by a flat-topped shock. At

t = 8 the steady (expected) state solution has already been established in the

computations. Again, given the large number of nodes used (50 for reaction

zone) the smearing of the shock due to Godunov's Method is not clearly seen

in the plots. However the smearing of the discontinuity in derivative is

quite apparent. Here one expects Pu = p which is not attained by the method.

Pig. 20 shows RCM results for the same problem. The steady state is

established at about t = 8 as for the previous figure. Shocks are absolutely

sharp as expected from RCM. Also the discontinuity in derivative is very

accurately represented i.e. p, = p . The questionable feature of the results

is the randomness within the rarefaction. This is typical of RCM.

Fig. 21. This is a continuation in time for the RCM computations

of Fig. 20. The steady state solution is nicely preserved.

Fig. 22. This is Godnov's Method on Problem 3 (see Fig. 4 for definition).

Reaction is initiated by a flat-topped shock and from the theory (see Chapter 1)

we expect a solution to an unsupported detonation, a steady reaction zone

followed by an unsteady "following flow" with p, = p.|. < p ^ S.

Fig. 23. This is an RCM solution to the same problem. The difference

between the Godunov and the RCM solutions is that RCM gives 2 discontinuities

in derivative at pj and p which from the theory is correct. Godunov's

Method has given a poor representation of these features.

Pig. 24. This is a Godunov solution to Problem 4 (see Fig. 5 for

defi-nition). The expected solution for large time is one overdriven steady

detonation with a curved profile for the reaction zone.

The numerical experiments on Pickett's analogue have been useful in

providing us with some ideas about the clear limitations of MacCormack (or

Lax-Wendroff), and, to a certain extent, Godunov's Method. RCM performs well if

one is prepared to accept the randomness in smooth parts of the flow.

2.6 Conclusions and further developments

Pickett's detonation analogue has proved very useful to the exercise of

selecting/testing possible candidate numerical methods for the full physical

problem. Of the three techniques tested RCM gives the more promising results.

For the near future we intend to test 3 more methods, namely. Roe's

flux difference splitting method, and two new techniques based on RCM which we

are now developing.

We feel that RCM has the potential for the full problem but research

needs to be carried out into two aspects: the problem of randomness in smooth

parts of the flow and the 2-dimensional extension of RCM. We are already

doing something about the first aspect with quite promising results so far.

Some of these new developments are included in the Appendix to this Chapter,

References

1. Chorin, A. 1976

Random Choice Solution of Hyperbolic Systems.

J. Comp. Phys. 22, pp 517-536.

2. Toro, E.F. and Clarke, J.F. 1985.

Application of the Random Choice Method to Computing Problems of Sol

id-Propel lant Combustion in a Closed Vessel.

College of Aeronautics report NPP 85/16, Cranfield Institute of Technology.

3. Godunov, S.K. 1959

Mat. Sb. 47 (1959), 2 7 1 ; also as USJPRS t r a n s l a t i o n 7226 (1960).

4. Toro, E.F. 1986

A New Numerical Technique for Quasi-Linear Hyperbolic Systems of

Conservation Laws.

College of Aeronautics Report CoA 86/26, December 1986,

Cranfield Institute of Technology, U.K.

21 -CHAPTER 2: Appendix

A New Technique for Pickett's Analogue

Here we present some preliminary results using a new technique based on the random choice method. The details will be reported on elsewhere (Ref. 4, Chapter 2 ) .

The technique can be viewed as a 'random' generalisation of Godunov's method discussed in Chapter 2, Section 2.4.4.

If in equation (21) the fluxes p""*^?, at the intercel 1 boundaries are evaluated at 'randomly' sampled solutions U. f then the new higher order

1 + 2

method is obtained. It may also be interpreted as random version of the Richtmyer two-step version of the Lax-Wendroff method, where the

Lax-Friedrichs intermediate step has been replaced by a Random Choice procedure. Fig. 25 shows results using the new Higher Order Method (HORCM) on

problem 1. Only 10 grid points within the reaction zone have been used. The computed solution is accurate inside the reaction zone, despite of the rather coarse mesh. As all higher order methods, the new technique runs into

difficulties where shocks are present. Here the response of the method are overshoots, the presence of which has a random character. Fig. 26 shows results for the same problem with more profiles on it.

It is well known that a higher order method may be modified to respond adequately to the presence of shocks and other discontinuities. We feel that the new technique discussed here has all the elements to provide a suitable modification. A possible answer is to advance via the Random Choice Method when shocks are involved and use the New Method elsewhere in the flow. The greatest problem then becomes that of finding sound and sufficiently general switching criteria. We are presently working on this problem for

the case of the Euler equations. Preliminary results based on the concept of a non-dimensional length are very encouraging.

For the present problem given by equations (la-c) it is not difficult to take decisions as to which algorithm to use. For expansions use HORCM,

for compressions use RCM. Fig. 27 shows results using the hybridised version of the new method. Compare Figs. 25 and 27. Overshoots have disappeared and shocks are absolutely sharp.

22

-For comparison we have solved the same problem under same conditions using RCM alone (Fig. 2 8 ) . The new method preserves the good feature of RCM (absolutely sharp shocks) and improves the representation of the smooth parts of the flow. Notice that the RCM results, specially for the pressure p, have a considerable amount of randomness within the reaction zone.

Fig. 29 shows results using the hibridised version of the new method for times up to about 10 units. The small dip near the tail of the rarefaction

(end of reaction zone) is a typical higher-order method feature. We are working on ways of dealing with this type of difficulty too.

Finally, Fig. 30 shows results obtained using 20 grid points inside the reaction zone.

We are presently devoting some effort to the problem of understanding some of the properties of the new method on basic problems. Also we are working with the Euler equations, but still on rather empirical (numerical)

reaction zone

unreacted material

x=0

X = . X i X r X

-1 ( a ) Overdriven D e t o n a t i o n , / J j = ^ = S

P P

— X

x = 0

1 ( b ) Overdriven Detonation , P^^-P^^^

P P.

^ - X

1 "

x = 0

1 ( c ) Unsupported D e t o n a t i o n , P^ = /^< ^=S

X

1 ^

reaction t a k i n g place.

0

zero reaction

x = 0

1 ( d ) Reaction Progress Variable

t = 0

t = t,

t = t.

x = 0

Fig. 2 Problem 1

Supported detonation , IC is steady

exact solution, /J, = /J = S.

P.

x = 0

f

9-Fig.3 Problem 2.

Supported detonation , IC is

f lat -topped shock , /J, = /J = 5

p^

p^y

t= ts

s \

Problem 3.

Unsupported d e t o n a t i o n , IC is

flat - topped shock, ^ = /». > /J = S

P

I

t = t,

x = 0

Fig.5

4 j

-Problem 4.

Overdriven detonation , IC is

f l a t - t o p p e d shock, P^^ P, .

'h

\d. •- 9 . 1 . 2 . 0 . e 1 X10 ia 9 2r o CJ e . < Ui Of 4 2 . 3 É - 1 - - ^ ) XI 7 G Ui c ir 3. 2 1 . 3 ( a' 3 PI GURE i

J

/ - ^

/

^

^ ^ / ^ ^ ^ l 2 4 6 8 10 12 D I S T ^ N C E

\ X

^*^*^' ^ k \

\ \

\ \

2 * S 8 10 12 DISTANCE ^^

f

^ ^ ^

k

2 4 6 9 10 12 DrS"rANCE DETONATION ANALOGUE

; nACCORnACK-COnPUTED (DOTTED) AND E<\C^ Q^ 2'5.03 S 5 . 3 3 Cf^L NUnBFR3. "53 DA -^ . '

1

l

1> 14 )

fl

14 < ^ULL ^ -- 1 .3 IG IG 16 SOLU^ ^ONS.

2. 10 12 DISTANCE 16 0 12 DISTANCE 16 Ui et: -D tri UJ a. X10l 7, 5. 4 3.

24

1 0 10 12 DISTANCE 14 16

DETONATION ANALOGUE

"^ICURE

7 :

HACCOR^ACK-COnPUTED (DOTTED) AND E<^CT f FULL) SOLUTIONS

- s T ^

^

^ , •

-O

i*l x,

Fig 8 Data for a Riemann problem

exemplified by density.

rarefaction

contact surface

shock

t = t.>0

Pi . P[.^[-0 \ P r . / ° r . U r = 0

_{t r O}

Fig.9 Example of a Riemann problem :

A shock-tube problem (ui = Ur = 0).

A = l

A=Ai

A=A,= A„,

( b )

\ p _{" l}

Prr. J

A l 1 ^(^^r

(c)

(d)

( « )

Fig. 10. Solution of the Riemann problem for case (A).

i

/

/

, A L

R

L'

0

( a )

= S

A , I -09 Pr _{Am = Ap} m Pr

( c )

3Co

( d )

3Co

( e )

n t l

n*l

AT n+S

n

Fig.12 Illustration of Godunov's method.

Fluxes evaluated on secondary

grid points at intercell boundaries.

solution sampled here

n + 1

n J/2

solution

sampled here

n

i-1

i * l

Fig. 13 Illustrations of Random Choice Method

Solution of Riemann problem sampled

in random fashion and values

1 0 . 8_ >-1— i - i _ w 6 . z UJ Q • * -2 . 0 . £

--'-"'>^V^i

^^^.-'-'''''^.^-'YO-q'''^

-'"'^ ^^^^^^^ ^^-^^^ ^ ^ ' ^ ^„„^•'•'''''^ ^ ^ ^ . ^ ' ' ' ' ' ' ^ ^ • " ^ ' ^ J ' " " " ^ ^ .^f^"''''''''^ ^^^^^ ^.^-'^^^ ^-^^-"'^

1

1

1

1 1 1

I i iL 1 1 1 • 1 1 1 I ) 1 2 3 4 5 G 7 8 DISTANCE X 1 0 - 1 10^ 8 . 2 O t - _ O G. < UJ _or ^_ 2 . 0 ^ " ^ ^ ~ ^ ~ ~ " ^ \ ^ " ^ ^ \ ^

X X

\ \ ^ _{^ \ \} \ \ \

V \

_{X X}

\ A. " ^ \ -X \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ " ^ \ 1 1 1 ï 1 T . i 1 0 1 2 3 ' ! 5 G 7 8 DISTANCE X 1 0 l 5_ i . lil oc =3 ^ V) 3 . UJ 2 . 1 . 0 .

e

- ' " " " " " ^ ^ - " ^ ' ' " " ^ ^ . • ^ ' ^ ' ^ ^ ^ ^ - - ^ -''"" ^-"^^^"^ ^^^''^^^ ^ . / - ^ ^^•^ ^^r^ L , , X ^ It - ' ' ' ' ' ^ ..,1''^^'^ .^f'"'''''^ 1 •^*' ..^"""''''''^ .^-""''''''^ I ^^•^ ^^.„^--^ ^^„.'''^ 1

^^^ ^--"^^''''^ ^.^""'''''''^ ' 1

^ ^ ^ " ^ ^ ^ ^ ^ ^ - " " ' ' ^ ^^^^.^f^'"'^^^ [ I 1 I \ i 1 1 1 1 I 1 I I ) 1 2 3 4 5 G 7 8 DISTANCE DETONATION ANALOGUE • 1 1 9 10 T - •• I 9 10 9 10

FIGURE IV :GODUNOV-COriPUTED SOLUTIONS FOR I N I T I A L VALUES AS

1 0 , 8. o G, 7 UJ Q •*0 35 •*0 ••S DISTANCE 55 G3 DISTANCE "G0 DETONATION ANALOGUE

FIGURE <S :GODUNOV-COnPUTED SOLUTIONS FOR I N I T I A L VALUES A S " GIVEN BY DOTTED PROFILES. T F = 1 3 . 2 0 4 , Q= 2 5 . 3 3 , CFL NUf1BER^3.53

COMPUTED SOLUTION BY THE RANDOM CHOICE METHOD WITH VAN DER CORPUT SAMPLING.

•—c z CM 0.0 0'.2 0 . 4 O.C X A X I S O.t 1.0 _{0.4 0.6 O.B} HpSM X A X I S o J-p» t n UJ CL CL ._ O 0.0 0.2 I I I 0.4 0.6 X A X I S —I r-0.8 1.0 u i z cc _ 1 • 1 I 1 1 1 1 I I 0.0 0.2 0.4 0.6 0.8 1.0 X A X I S

Fc^.

\ b

SHOCK-TUBE TEST PROBLEM! COMPUTED (SYMBOL) AND EXACT ( L I N E ) SOLUTIONS.

10 12 DISTANCE IG X10 10. w w UJ 8. Qi. IS) O Qi tL Z

H G.

I -u < UJ 10 DISTANCE ie

' DETONATION ANALOGUE

FIGURE l>: COMPUTED SOLUTIONS BY THE RANDOM CHOICE METHOD.

Q= 25.00 S= 5.00 CFL NUMBER=0.500

8.

Ë ^

z UI o ^. 0 DISTANCE IG X10' UJ V) UJ T 2 . 10 12 DISTANCE 16 DETONATION ANALOGUE

FIGURE I? : COMPUTED SOLUTIONS BY THE RANDOM CHOICE METHOD Q= 2 5 . 0 0 S= 5 . 0 0 CFL NUMBER=0.500

35 40 DISTANCE 45 - r 50 DISTANCE 60 DISTANCE DETONATION ANALOGUE

FIGURE 20 :RCM SOLUTIONS EVERY 100 TIME STEPS FOR I N I T I A L VALUES AS GIVEN BY DOTTED PROFILES. Q= 2 5 . 0 0 S= 5 , 0 0 CFL NUMBER=0.50

ie. (/) G. z UJ Q 4. 2. 0 55 60 65 70 —r— 75 80 85 90 95 DISTANCE 100 105 110 115 90 95 100 DISTANCE 110 115 X101 6 , —r— 75 80 55 60 65 70 85 90 95 DISTANCE 100 105 110 115 DETONATION ANALOGUE

FIGURE 21 :RCM SOLUTIONS EVERY 100 TIME STEPS FOR I N I T I A L VALUES AS GIVEN BY DOTTED PROFILES. Q= 2 5 . 0 0 S= 5 . 0 0 CFL NUMBER=0.50

10, 8_ >t -co G z UJ Q DISTANCE X10' X10"l 10. 8. z o o G. < UJ 4.. DISTANCE X10l 5, UJ O:! v> 3. UJ Q. DISTANCE X10'

DETONATION ANALOGUE

FIGURE 22 .-GODUNOV-COMPUTED SOLUTIONS FOR INITIAL VALUES AS

12, 10 I- 8 <r z Q G-24 0 —r-8 - I — 14 10 12 DISTANCE X10' )G X10 10. z o »—. I— c G. < • UJ Of 4 . 0 L 8 10 12 DISTANCE I * IG X10I X10l 6 , 5. UJ or :D tp to UJ Of Q. 24 1 0 2 10 14 I DISTANCE X10' 16 DETONATION ANALOGUE

FIGURE 23 :RCM SOLUTIONS EVERY 200 TIME STEPS FOR I N I T I A L VALUES AS GIVEN BY DOTTED PROFILES, Q= 2 5 , 0 0 S= 5 , 0 0 CFL NUMBER=0.50

\ i 10. 9 . Cl' a C * , 2-0.

. A

/

r

1 1

1

1

1

1

/

A

/ / ^ /

A A

/ /

A

/

A

/

A

/

A

/ 10 15 20 25 30 35 40 45 DISTANCE 50 55 60 X10'' 10. z o .-. u C-UJ 8. G. 4 , 2. 0 . 1

1

J 10 15 20 25 30 35 40 *5 DISTANCE 50 G0 X10' 6-, 5. ^ 4 cr IT UJ ^ 3 10 15 20 25 30 35 40 45 DISTANCE 50 55 60

DETONATION ANALOGUE

FIGURE XH .'GODUNOV-COMPUTED SOLUTIONS FOR INITIAL VALUES AS

—T-2 8 10 12 DISTANCE 10 12 DISTANCE X10f 12 10 12 DISTANCE -e—I 16 16 T — I e •€> 1 16

DETONATION ANALOGUE, PROBLEM 1

FIGURE 25 : HORCM COMPUTED SOLUTIONS CSYMBOL) AND EXACT SOLUTION

GIVEN BY FULL-LINE PROFILES. TF= 1 . 9 2 1 , Q= 2 5 . 0 0 , CFL NÜMBER=1.00

ie DISTANCE 16 n 16 10 12 DISTANCE X10l 6., UJ tn V) UJ o: o. —1 16 DISTANCE

DETONATION ANALOGUE, PROBLEM 1

FIGURE 27 :HYBRID METHOD SOLUTIONS (SYMBOL) AND EXACT SOLUTION

GIVEN fcn FULL-LINE PROFILES. TF= 2,02G, Q= 25.00, CFL NUMBER=1.00

4. • * — T - 1 — » -12 2 10 DISTANCE 14 16 10 12 DISTANCE 16 1. —r-2 •e—r-10 12 DISTANCE 14 16

DETONATION ANALOGUE, PROBLEM 1

FIGURE 28 :RCM-COMPUTED SOLUTIONS (SYMBOL) AND EXACT SOLUTION GIVEN BY FULL-LINE PROFILES. TF= 2 . 0 1 1 , Q= 2 5 , 0 0 , CFL NUMBER=1,00

10 15 20 25 30 35 40 45 DISTANCE

50 55 60

35 40 45 DISTANCE

DETONATION ANALOGUE, PROBLEM 1

FIGURE 29 : H Y B R I D METHOD SOLUTIONS (SYMBOL) AND EXACT SOLUTION GIVEN BY FULL-LINE PROFILES. T F = 1 0 . 0 3 7 , Q= 2 5 , 0 0 , CFL NUMBER=1.00

10 12 DISTANCE 16 10 12 DISTANCE 16 10 12 DISTANCE 16

DETONATION ANALOGUE, PROBLEM 1

FIGURE 30 :HYBRID METHOD SOLUTIONS (SYMBOL) AND EXACT SOLUTION

GIVEN BY FULL-LINE PROFILES. TF= 2,023, 0= 25.00, CFL NUMBER=1,00

- 48 - : ^

3. A Moving Finite Element (MFE) Solution of Ficketts Detonation Analogue 3.1 Introduction

The present analogue which is required to possess some of the important features of the physical problem is used to test the MFE method. The MFE method introduced by Miller and Miller (1981) and Miller (1981) has been used with considerable success for the solution of scalar hyperbolic equations in one space dimensions.

The MFE method is based on the measure of the error in the solution; this being the optimal feature of the method suggests that good approximations are obtainable on coarse grids.

Recent work on the MFE method by Wathen and Baines (1985) Baines (1985), Edwards (1985) have shown different forms of the MFE method with certain

advantages in each case. In particular we have the local MFE method and the Mobile Element Method (MEM) which is derived from the local MFE method. We give a brief description of the Standard MFE, the local MFE and the MEM which has been implemented.

In Chapter two, four test problems were used to test various finite difference techniques. Here we investigate the difficulties which maybe encountered when applying the MEM to the steady-state problem (problem (1),

for which there exist exact solutions, see Chapter one) and the unsteady problems (problems 2 and 4 ) .

3.2 The Standard MFE Method For an equation of the form

^ ^ f ( u ) , = 0 1 we seek a piecewise linear function v of the form

N where

3-! = 3.(t), j = 1(1)N, are nodal amplitudes

a. = a.(x,s), j = 1(1)N, are linear basis functions of local compact support and s = s(t) is a time-dependent vector of nodal positions s. , j = 1(1)N.

49

-The piecewise linear function and the a. basis function can be represented as in figures (3.2.1) and (3.2.2)

J-1

^j ^ X

> 1

••x Figure 3 . 2 . 1 . F i g . 3.2.2 P a r t i a l d i f f e r e n t i a t i o n of (2) w i t h respect to time y i e l d s da. 9a. •J By the chain rule

3a. 9a.

—J-

= y

s. — i

9t ^. "i 8s.

substituting (4) into (3) gives

Interchanging the order of summation gives 9a-_i

al = r ^ft * ïj

'i \

^• aiT

3v

"^ i= y^fj^^j'

50

-The function g. may be regarded as a second type of basis function. In fact

J

both have the same support and indeed for linear elements it can be shown by direct differentiation (see also Lynch (1982)).

In the case of linear elements equation (9) can be written as

3. = -ma.

10

where m represents the local slope of the solution. We note that the 3- is a discontinuous linear basis function which can be represented as in

figure 3.2.3

B.(x,s)

^ X

Figure 3.2.3

In the approximation space spanned by the basis functions a. and 3- we

J J

define the Ip norm of the residual as V, + f,(v)

The Lp residual which is a measure of the error in the approximation space is minimized with respect to the variables a. and s.. This leads to the

J J

extended system of Galerkin equations

j = 1(1)N. <v^ + fx(v), a.> = 0

<v, . f ^ ( v ) , 3 j > = 0

11

12

where <.,.> is the L^ inner product. Substituting for v^ from (7) we obtain the set of MFE equations

51

-14

15

where

y - i32'S^;a2,S2;... ;a|^,S|^}

A(y) is the MFE matrix, which is square and symmetric and consisting of

inner products of the basis function a and 3 in 2 x 2 blocks

<a.,a.><a.,

&.>

'h

a-><3,.. S.>

and g(y) is a vector arising from the terms involving f(v) in (12).

The MFE set of equations can then be solved for the parameters

a., s-, j = l(l)N) using the preconditioned conjugate gradient method,

provided that the 2 x 2 blocks are non-singular. There are in fact two types

of singularities associated with the MFE matrix. Miller and Miller (1981) and

Miller (1981) have introduced penalty functions, which they use to force some

conditions in the parameters a-, s. (j = 1(1 )N) to prevent such singularities.

In this report we deal with the singularities in a different way, which is

discussed in section 3.4.

Finally, once the parameters a. and s. are obtained we solve for a. and

J J J

s. using a time-stepping integration method.

J

3.3 The Local MFE Method

Consider the equation (7) which is used to approximate u.. Using

equation (9) we can rewrite equation (7) as

9v v^.; ^v- > IC

^ = ^ ( ^ j - 9 7 ^ j h 1Ö

Now the local element form of (16) can be expressed as

2

=

I I.

2 ^i^ki

k " i = l ^ ^ k

where

\i/^\ - \ \ - ( ^ \ ^ k i ' ''ki ^^^ '^''' -^yp^

element basis functions shown in figures 3.3.1(a) and (b) and w.,, w.^ are

weighting velocities for each element k, where the element numbering is

defined as k = j-i, (k

=i,...,

N +J)

52

-\l^\'-^

1--^x

Figure 3.3.1 (a) 0|^2(^'§) * • x Figure 3.3.1 (b)

Within each element, k, the local MFE method proceeds to minimize the Lj, residual

\

- fx(^)

k^'\i \i^'\''.^\Ai^

with respect to the W. .. This lead to the elementwise Galerkin equations

^li^ki/A^k <^ki'^i> = -<'^ki-Vl

^ = 1 Wk2/^^k <^k2'^ki> = -<'^k2-VJ

•k = J(l)

H-i

53

-where it is clear that we have a 2 x 2 matrix system of equations,i.e.

^k ^ ^ ^k , k -|(1) N - K

19

where C,. =

\l^'^\' \2/^\^

< ^ 1 ^ l > < ^ k l ^ 2 > " < ^ 2 ^ 1 > < ^ 2 ^k2> , k , k 6 k 2 1 1 2 r(l) N-i , |(1) N-l , , k = |(1) N-^ ,

20

21

and

h -

^<^l'-f(^)x><V-^(^)x>^

k = f ( l ) N-è .

^^k ~ ^k2 " ^kl 22 23 The inner products in (22) can either be integrated exactly or integrated approximately, this will depend on f(v) , Once b. is known, we then

A ~ K

proceed to solve (19) for w. which give information about the element motions. The information is then transferred to the components of the nodal velocity a. s. via (17) which give two equations for each node i.e.

1 - A./As.

1 -

' j ^ j . !

'3 ^•-i,2/^^j

which is usually written as M .y. = w . where 1 -A./As. J J 1 -A. T / A S . T , y = a. J s .

24

25

26

54 -^j =

'y

w. _J-^,2 /As,

VLl/ASj^^

^j - ^j-1' ^^j s. - S

_{j - 1}

27

The 2 x 2 system can be solved directly for a. s. provided

J J

M. ^ 0. The solutions are of the form,

. "^.i.l^'^^1.l"".i-i.2^''^i) - "^./'^^i'(".i.Ll^'^1'

C,^ ?^ 0 and

28

^'^J-L2/^^j -

V L I ^ V ^

j - (A.^i/ASj^^ - Aj/ASj) s- =

29

Finally a- and s. are obtained by a time-stepping integration method. Note that while both the global and local methods give the same solution, the immediate advantage of the local MFE method is that it only requires when possible, the inversion of 2 x 2 matrices.

3.4 Singularities

For linear elements in any number of dimensions the global MFE matrix A can be expressed as M^CM

₃₀

where

\ °

0 \ .

, c =

f\ °1

0 \

and M. and Cj^ are defined in the previous section. Both C and M are diagonal block matrices (apart from boundary effects).

If M or C is singular then A is singular. The singularity of M occurs

if

A./AS. - A.^^/As .^^ . 0

₃₁

This type of singularity corresponds to two adjacent slopes in the solution being equal. This is referred to as "parallelism" or "collinearity" which can be remedied (see Baines 1985).

55 -The singularity of C occurs if

^k2 - ^kl ^ °

32

This is Indicative of a formation of a shock due to node overtaking. This singularity is remedied by replacing two of the equations in (24) with the shock condition (i.e. jump condition or internal boundary condition).

^k2 " ^kl

^k2 " ^ 1

33

so that the new configuration (see figure 3.4.1(b)) moves as a shock

^kl ^ 2

'kl

'k2

L : R i ^j-1 'kl ^j "^k2 ^j+1 (a) (b) Figure 3.4.1 Formation of a shock

Replacing (24) by (33) we find that since Sj^^ is known, only one equation is now needed in the element L (see figure 3.4.1(b))

\ l =^^-Ul

' ^-l,L^kl

34

and similarly for the element to the right of the shock or collapsed element

\ 2 = ' ^ j + U R ^ ^ + l . R 'k2

35

When a second node attempts to overtake a pair of shock nodes the same procedure can be used but this time we usually delete a node.

For each time level t" we compute the smallest time requiredfor s,, to overtake s. ^j. This time step is then compared with an initial fixed time step, DT , say. Node overtaking is then prevented by taking the smallest of the two.

56 -3.5 The Mobile Element Method

First we rewrite the scalar equation (1) using the convective derivative Du

Dt - ^t ^ ^ ^

36

so that Du

Dt - ^^ ^ ^x = °

37

It is this form of the equation that a solution is sought and one which we shall make reference to in this section.

Now instead of solving the 2 x 2 system (25) for a. and s. (j = 1(1 )N' we rewrite the system in the form

(AS. (As.,

j+1

a. J (Wj_^^2.A.é)(As.)

(^j+Ll

'

Vl^^^'V^'

38

The square root factors are introduced in order to maintain a conservative scheme. This is an overdetermined system for a.. To solve this system in terms of s as a parameter, we may minimize the norm.

2 2 D^La. - D"^w with respect to a .. Where D = As. 0

''i

"j-i,2 *

"f

w • + A.. j+i,l j+1,

39

This gives

^'^j-l,2^^j+Ll^

. I V i ^ V

^J ^ (ASj + As.,^)

'

TAiT^Ai-TT ^j

40

Now all that remains is to choose the s.. For the equation (40) using approximate integration on a uniform grid, we have

-9f W . 1 o -t- W . 1 1 J - L 2 j-i-J,l As, + As

j+1

;ix

+ 0(As'

41

57 -also A . , -I-J+1 As . 1 - As. J+1 J 9u 9x O(As^)

42

Equation (40) therefore represents a discrete analogue of equation (37). A natural choice for the s is obtained if we consider characteristic

methods. We start off by requiring that KT; = 0 which would mean in the corresponding discrete version (40) that a. = 0, so a choice which is in close correspondence to characteristic description is

(w

_{+ w.}

j-J,2 " "J+i.l

TÏ-TTÏT]—

)

li/i^'X*^^)

and

a . J 43 44 We note that the same result can be achieved by minimizing

^1

which is trivial in the scalar cases. In the case of systems we cannot expect to obtain a "magic" formula for s. which would make a. zero. So instead we have to

. 2 •^ J

minimize a. with respect to s. for an approximate characteristic description. J 2 J

Further we note that this method will only fail if As. = 0. This shock type of singularity is remedied in the way described in section 3.4. The MEM is therefore able to cope with linear data more naturally than the MFE methods. 3.6 MEM for Systems of Equations

We consider systems of hyperbolic equations and seek an approximate

solution using a single grid. For simplicity we give a summary of the analysis for a two component system (see Edwards (1985)). Using (1) we may write the two component system as

u ^ " . f(u(l>. u < ^ ' ) , = 0

As before, an elementwise minimization over w^• , w^• gives

45

c(i) ^(1) . b(l) ^k -k -k c(2) ^(2) _ t,(2) ^k -k " -k 46 47

58 -For a single grid

D^M.y, = D^w

J-J

-J

48

which is an over-determined system for y .

J where As. J 0 0 0 0 As 0 0

j+1

0 0 0 0 0 0 ' ^ j + 1

49

1 0 1 0 0 1 0 1 - A ( ^ ) / A S . J J

-A(1}/AS. ,

J+1 J+1 - A ( . 2 ) / A S . J J

-A(^J/AS. T

J+l J+1

^r'

f'

s . J ''j-è,2 j+i,l

w

_j-i,2 (2) :,(2) VI _i-:-kl

Again the square root factors are introduced to give a conservative scheme. We regroup terms involving s in (48) on the right hand side of (48). We now have a different over-determined system to solve, namely.

D ^ a = D'^w

where D is defined in (49) and

50

a =

J

1

1

0

0

0

0

1

1

"1?1.2 ^ f'^

53

-The usual minimization procedure of the Lp norm of (50) with respect to a gives the normal equations

L^DLa = L^w

₅₁

where a can be written explicitly in terms of Si,i.e.

52

.._^?l^Li),i!l!Mi)|3

j - -(As + As 71 (As^ + As ) j

a(2)

53

•fl) Now it is clear that there exist no choice of s^which will give both a^ '

• (2) "^ and a\ ' zero. However we can minimize these two parameters by minimizing

2

(P)

,(P)

_I

_{(a^M' . (P = 1.2)}

₅₄

with respect to s-, which results in

- I '

s . = J

P=l

<'-i:i"^n.2^*!:L)

i.j^y^

2 P=l (P = 1,2)

55

Generalising (41) and (42) we see that equations (52) and (53) are a discrete analogue of the differential equation system

,(P) 3f(P) . 3,(P) Du'

Dt

- Ü ^ ^ S ,(P = 1,2)

56

and also the iiininizrtioii orocociuro is the discrete analogue of minimizing Du Dt

(P) 2

_9f 9x

(P)

+ s 9u _9x

(P)

, (P = 1,2)

57

with respect to s.