A flux splitting method for the numerical simulation of one-dimensional compressible flow

Mathematics and Computers North-Holland

in Simulation

309

XXV (1983) 309-320

A FLUX SPLITTING METHOD FOR THE NUMERICAL ONE-DIMENSIONAL COMPRESSIBLE FLOW

SIMULATION

OF

S.R. MULPURU Systems Analysis Branch, Canada ROE IL0

Whiteshell

Nuclear Research

Estuhlishment,

Atomtc

Energy of Cantrda Limited,

Pinuw,u, Munitohu,

Using the technique of flux vector splitting, it is shown that one-dimensional, inviscid, compressible-flow equations possess a split conservation form. Some attractive features of this form for the design of finite-difference solution schemes are discussed. Based on the split form, two solution themes are designed. One is a first-order accurate ‘upwind’ scheme and the other is similar to the Lax-Wendroff scheme. A hybrid scheme, based on a nonlinear weighting of these two schemes, is demonstrated to yield results superior to either of the two in the solution of an ideal shock-tube problem.

1. Introduction The linear,

first-order,

scalar equation

Y+c Lo at

’ ax

where c is a constant, describes passive convection of an initialf(x) with the velocity c. If the initialf( X) is monotone, then preservation of monotonicity is a qualitative aspect of the solution. Godunov [7] proved that there are no linear, second or higher order finite-difference solution schemes that always satisfy this qualitative aspect (i.e. they generate oscillations), and the most accurate first-order scheme that does satisfy it is the one based on ‘Forward-Time’, ‘Upwind-Space’ (FTUS) approximations of the derivatives. Since the one-dimensional, inviscid, compressible-flow equation set can be cast in a form similar to Eq. (l), it is of interest to extend the FTUS scheme to solve it. It is well known [5] that, for the isentropic condition, the one-dimensional, inviscid, compressible-flow equation set can be decoupled into three individual equations, each of the same form as the linear scalar equation (1). The dependent variables of these equations are the two Riemann invariants and the entropy. The coefficients of the spatial derivatives are the eigenvalues of the hyperbolic equation set u + a and u, where u is the flow velocity and a is the local sound speed. Then the extension of the FTUS scheme to the calculation of the decoupled equations is straightforward. For the non-isentropic case, decoupling of a different type is possible. Using flux vector splitting [ 191, we show in this paper that the one-dimensional equation set possesses a split conservation form and can be decoupled for the purpose of calculation into three individual vector continuity equations, each of the form

(2) The three vectors are eigenvectors corresponding to the eigenvalues u & a and u of the hyperbolic set. Exploiting the split conservation form, we then design both a version of the first-order accurate FTUS scheme and a higher-order scheme similar to the Lax-Wendroff scheme. The FTUS scheme does not generate any spurious oscillations, but its accuracy is poor. The latter can 0378-4754/83/$3.00

0 1983, Elsevier Science Publishers

B.V. (North-Holland)

S.R. Mulpuru / Flux splitting method

310

pose a severe problem when computation cost dictates the use of coarse grids. Higher-order schemes, on the other hand, yield more accurate results in smooth regions but generate spurious oscillations near sharp gradients. These oscillations not only damage the accuracy of the results but also may trigger convergence [9] to physically false solutions. This problem can be avoided by using a weighted sum of a lower-and a higher-order estimate, favouring the higher-order one, but subject to the restriction that the solution be free of spurious extrema. Zalesak [21] designed such a weighting procedure and showed that it was a generalization of the flux-correction procedure developed by Boris, Book and Haine [ 1,2,3]. This weighting procedure is applied to the two schemes presented in this paper. The result is an efficient, monotonic, explicit hybrid scheme with essentially second-order accuracy. The scheme is compared with the well known second-order Lax-Wendroff-Richtmyer two-step scheme [16] for the solution of the one-dimensional shock-tube problem of Sod [ 181. The results demonstrate the superiority of the hybrid scheme.

2. Equations One-dimensional,

inviscid,

compressible-flow

equations

in Cartesian

co-ordinates

and conservation

form

are

where x and t denote

q=

position

[!I;

F=

and time and the vectors are

i,~;#c;~P)~.

Here, p, G, and E are the mass, momentum by the perfect gas equation of state, P=(y-

and energy densities,

respectively,

1)pe

P is given

energy, y is the ratio of specific heats and

pe = E - G2/(2p). Eq. (3) can be written _

and the pressure

(5)

where e is the specific internal

s+

(4)

(6)

in a quasi-linear -

[A]$=0

form: (7)

where the Jacobian matrix [A] is aF/&& The matrix [A] has distinct eigenvalues, u and u _t a. Here, u ( = G/p) is the flow velocity and a (= [ ~l’/p]“~) is the sound speed and these can have only real values for physically relevant situations. [A] can be reduced to the diagonal form by a similarity transformation

[61:

[w’L4Pl=

[Al

(8)

where [A] is a diagonal matrix whose elements are the eigenvalues (u, u + a, u - a) of [A] and the columns of [B] are eigenvectors corresponding to the eigenvalues. These properties satisfy a definition of the hyperbolicity of the system [lo].

S. R. Mulpuru

311

/ Flux .splitting method

2. I, Flux vector splitting Steger and Warming [19] discovered a remarkable here. For equations of state of the form

property

of the flux vector i? It is briefly

reviewed

(9)

P=pf(e)

of which the perfect gas equation (5) is one, the nonlinear flux vector F(+) is a homogeneous function of -degree one; that is F( &) = cyF( +) for any value (Y,and the application of Euler’s theorem on homogeneous functions [ 1 l] to F yields F=

(10)

[A]&

The flux vectors in the second and third spatial dimensions also satisfy relations analogous Further, (10) implies that the F can be split into component vectors, each having an eigenvalue scalar multiple [19]. This can readily be seen by substituting a variant of (8) into (10):

I 1

F= [B][n][B]-‘c$ 1

=h,[B]

[B]-‘&-tX,[B]

0

0

1 1 0

1

0

Pr’cp+~,Pl

to Eq. (10). of [A] as its

I 1 0

0

1

DT’6 tll)

= L,.f, I= 1 where

[B

-‘<-,

I= 1,2,

(12)

3,

(13) and A, = u,

Xz=(U+a),

X,=(24-a),

specifically,

2.2. Splitting

into continuity

equations

We note that the component

vectors f, (Eq. (12)) are eigenvectors

corresponding

to eigenvalues

A, of [A]

S. R. Mulpuru / Flux splitting method

312

because

they satisfy the relation [A]J,=X,J,,

I= 1,2,

One can verify (14) by substituting

[A ] is a diagonal

3. (12) and a variant

matrix of the eigenvalues

Using the definition

(14) of (8) into its left-hand

X, of [A]. Multiplying

side, writing

the two diagonal

matrices

thereby

in (15) one gets

of 6(Eq. (13)), one can write (16) as

(17) =

&.!I,

I= 1,2,

This completes the verification of (14). Summation of (14) for all the eigenvectors

(18)

3.

yields

Also, from (11) and (10) we get

Comparing

(19) and (20), we find that the vector 6 is composed jL

Substituting as

of the eigenvectors

f,

;A. I= I

(21)

(11) and (21) into (3) we can write the one-dimensional

equation

set in split conservation

form

(22) The form of (22) contrasts

[QIE+ at

with the characteristic

[A][QIE=O ax

form (derived

from (7) and (8)): (23)

313

S. R. Mulpuru / Flux splitting method

where (24) The characteristic form (Eq. (23)) was the basis [4,8,14,15] for non-conservative finite-difference schemes. The conservation form, however, guarantees that the computed discontinuities propagate with the correct speed, and it has been preferred in shock calculations [ 131. Further splitting the equation set (22) into individual, vector continuity equations $+&,f,)=O,

l=

1,2, 3

reveals some interesting aspects that have a bearing on finite-difference solution schemes. If each vector f, is advanced in a finite-difference grid over a time step At, at a grid point operation (jJ’“‘=

(L,),(j$,

i, by an

I= 1, 2, 3

where (L,), is a finite-difference operator operations for all the vectors j, gives

for the calculation

of vector f,, then

the addition

of these

(27) I

I

for advancing 6 at i, by virtue of (21). This implies that if all the individual operations L, are stable and possess a certain order of accuracy, so does the operation for the calculation of $. The same result holds if all the individual operations are conservative. Thus, a solution scheme for the calculation of 6 can be built up linearly from the individual solution schemes of its component vectorsf,. The solution schemes for the component vectors need not all be the same. For example, for systems of stiff equations, such as those for calculating laminar deflagrations where the eigenvalues differ greatly, implicit schemes can be used to ease the stringent stability requirements for the solution of those equations associated with the eigenvalues of larger magnitude and explicit schemes for the rest.

3. Solution scheme: Scalar continuity equation It is convenient

to consider

$+jgA(x. t)f] and extend the solution uniform spatial interval calculation of (28) is

first a scalar continuity

=o

equation:

(28)

schemes of this equation to the full equation set later. On a finite-difference grid of AX and time interval At, a version of the first-order accurate FTUS scheme for the

-.ty=g,-l/2_&+I/2 f, 17+’ where if a,+ 1,2 2

0,

if a,+ l,2 < 0,

(30)

314

S.R. Mulpuru / Flux splitting method a ,+ l/2 = ( 0,” + fJ,“+J/2,

(31)

a,” = A; At/Ax

(32)

and f,” denotes f(xo + ihx, t, + n At). Equation (29) belongs to a class of one-step, explicit conservative central-difference operators and can be written in a general form as

schemes

that

use three-point

f, n+1-fin=(C-D),-1,2-(C-D)r+,,2 where C and D represent

(33)

the fluxes that are due to convection

and numerical

diffusion.

They are defined

by CI-c l/2 =tlbf>:+l+ Di+,,2

bf>:J?

= t sign{o,+,,2}b,+,,2{(af):+,

(34) -(d):)

(35)

where 0;+1/2 = w+,

+ 4/2

(36)

and 4, I/Z is a dimensionless

positive quantity whose role will become clear shortly. If X in (28) is considered to be a constant, then, according to the well-known Fourier analysis [16], the linear stability conditions for the schemes defined by (33) to (36) in which a,, ,,2 is a constant, are b,+ 1,21 G 4+,/z

G l/la,+

(37)

i/21

and

la,“1< 1

(38)

and the monotonicity

conditions

[7] are

=Gl/lo,+ I,21

1 $bi+l,z

(39)

and ]a”( < 1.

(40)

The upper limit of b,, ,,2 defines Lax’s scheme [ 121. The smallest value of b,, ,,2 to maintain

monotonicity

bI+ l/2 = 1,

defines

the first-order b r-cl/2 =

(41) accurate

Ior+1,2L

FTUS

scheme, whilce the smallest

value to maintain

linear stability, (42)

defines the second-order accurate Lax-Wendroff scheme [13]. We shall consider the last two schemes further. Although it is not clear if the FTUS scheme preserves monotonicity subject to condition (40) when X is not constant, the results (presented later) show that it does. Also, Eqs. (33) to (36) along with the definition (42) define the Lax-Wendroff scheme only if h is invariant over the time step At. This means that the scheme has second-order accuracy in space but not in time. We label it for future reference as a scheme similar to the Lax-Wendroff scheme. Equations (33) to (36) with the embedded parameter b,, 1,2, provide an efficient setup for the nonlinear weighting procedure described later.

315

S. R. Mulpuru / Flux splitting method

3.1. Solution scheme: Full equation set The class of one-step conservative schemes defined by (33) to (36) for calculating the scalar equation (28) can readily be extended to the one-dimensional equation set (3) because of its split form (Eq. (22)). The schemes are

,$,[m:+’- (6 ,:‘I =

(@:+I

-~)=(c-~),~,,*-(c-~),+1,2

(43)

where

o,+ 1,2

=

sign(ul),+,,2(bt),+,,2[(utf,):_, - (utL):llr

tC

[(d+,

(d+1,2=

(45) (46)

+bXM.

The scheme for calculating 6, defined by (43) to (46), results from the mere addition of the schemes for the calculation of its component vectors 1,. The stability condition for calculating 6, then, is the severest of the conditions for calculating the individual vectors 7,. If all the eigenvalues X, (and thereby a,) are considered to be constant, the linear stability conditions follow from the corresponding conditions (37) and (38) for the scalar equation (28), and are given by Ku,) r+ I,21

G

,+1/2<

(b/l

l/l~,+,,*I3

I=

1,

2,

3

(47)

and ](u,):] < 1,

I= 1,2,

3.

(48)

Subject to inequality (47), the larger the value of (6,),+ ,,*, the more diffusive here the lower limit, @,),+1/z and a constant

Iw+I,*I~

=

intermediate

I= 1, 293,

(49)

value

= 1,

@Jr+,,2

the scheme is. We considered

I= 1,2,3.

(50)

Equations (43) to (46) along with (49) define a scheme similar to the second-order accurate Lax-Wendroff scheme, and with (50) they define the first-order accurate FTUS scheme. A hybrid-scheme, resulting from a nonlinear weighting of these two schemes, is described later. 3.2. Comparison

with other schemes

It is interesting to note that the schemes of Rusanov [17] and Van Leer [20] resemble presented in this work. Their schemes are defined in the present notation by (43) where

Equations

1 At

C ,+I/2

=7

D !+I/2

=1 ‘K,+

the two schemes

-

-_(e;;, Ax

+k;q,

I,2

(+?+I

-8:).

(5 1) and (52) should be compared

(51)

(52) with (44) and (45). By virtue of ( 11) and (32) (5 1) and (44) are

316

S. R. Mulpuru / Flux splitting method

PRKYJFE---

z?

m0.000

0.125

0.250

0.375

0.500

0.625

0.750

0.875

000

POSITION Fig. 1. First-order accurate upwind the symbols the computer solution.

identical.

Equation

(FTUS)

scheme (Eqs. (43) to (46) and (50)). The solid lines represent

(52) can be written,

after substituting

the analytical

solution

and

(21), as

3

D ,+ I/2 =

tK,+,,2

c I=

For Rusanov’s K ,+

t(h):+,

-

(53)

(r,>:>.

1

scheme, l/2

=

Maximum[l(a,),+,,2I,

1(~2),+1~21~

1(~~),+1,21].

(54)

S.R.

--I

Mulpuru

/ Flux

splittmg

method

-r--7-

•I+ +

NALYrIc4L HISHER-ORDER SCHER

0.125

000

0.250

Fig. 2. The scheme akin to the second-order

0.375

accurate

0.500

Lax-Wendroff

0.6‘25

0

4 50

0.875

1.000

scheme (Eqs. (43) to (46) and (49)).

For Van Leer’s scheme,

K ,+

l/2

=

Maximum[

(0:)

,+

,,?,

w,+l/z~

(+%,*].

The schemes of Rusanov and Van Leer treat the calculation of the vectors f, , i,, f, as though they are all associated with the same eigenvalue having the largest magnitude, evidently in the interest of conservatively guaranteeing linear stability. This conservatism results in inaccuracies in the calculation of the vectors not associated with the eigenvalue of the largest magnitude. Since the schemes in this paper use the appropriate eigenvalue for each vector, they are more accurate while still maintaining linear stability.

S.R. Mulpuru / Flux splrttmg method

T104

T -

AMLYrIcAL

0 + *

LAX-WENDROFF TWO-STEP SCHEE

I

I

I

m

PKSSURE-

i

?a t+

3”

! SxK

/

0.1aao

0.;25

Fig. 3. Second-order

accurate

0.253

I 0.625

I

0.375

0.500

Lax-Wendroff-Richtmyer

[16]

0.750

0.875

1 .0@0

scheme (with no artifical viscous terms added).

4. Results

The test problem chosen defined in the problem are P=

l,p=

l,u=O,

P=O.l,p=O.l25,u=O,

is the ideal

shock-tube

problem

of Sod [ 181. The initial

(t = 0) conditions

0 < x < 0.5, 0.5
1.

The space-step (Ax) size used is l/100. The time step (At) is calculated at each time level subject to the condition that the maximum Courant number is 0.75, and the calculation is terminated when the cumulative addition of the time steps is close to 0.15. The numerical solutions at the terminated time are compared with the analytical solutions at the same time in Fig. 1 to 4. The solutions obtained with the first-order accurate FTUS scheme (Eqs. (43) to (46) and (50)) are shown

319

S. R. Mulpuru / Flux splrtting method

r--

--

1

~ N4LYrIcAL 0 + + HYBRID SCHEPE

PRESSURE --

0.375

0.500

0.625

0.750

0.875

1 .QMM

rnSITIoN Fig. 4. The hybrid

scheme resulting

from the nonlinear

weighting

of the two schemes presented

in Figs.

1 and 2.

in Fig. 1. There are no spurious oscillations, but the resolution of the contact discontinuity is poor, and the constant state between the expansion wave and the contact discontinuity is not realized. The scheme (defined by Eqs. (43) to (46) and (49)) similar to the second-order accurate Lax-Wendroff scheme improves the resolution (see Fig. 2) but produces post-shock oscillations. For comparison, the solutions obtained with the well known second-order Lax-Wendroff-Richtmyer two-step scheme [16] (with no artificial viscous terms added) are presented in Fig. 3. The scheme presented in this paper (Fig. 2) yields less oscillatory solutions and produces no spike at the tail of the expansion wave. These solutions also compare well with those of several other methods tested by Sod [18]. The CPU time on a DEC-KI-IO computer for each of the two schemes is 16 seconds. Using the weighting procedure described below, the CPU time is 28 seconds. Zalesak [21] designed a weighting procedure that weights the solutions of a higher-order scheme as heavily as possible with those of a lower-order scheme, subject to the restriction that the final weighted

320

S. R. Mulpuru

/ Flux splitting method

results be free of non-physical extrema. He interpreted the procedure to be a generalization of the flux-correction procedure developed by Boris, Book and Hain [ 1,2,3]. For the weighting procedure, both schemes must be in conservation form and stable, and the lower-order schemes must be free of non-physical extrema. The two schemes (Eqs. (43) to (46), (49) and (50)) presented in this work satisfy these conditions. Moreover, they provide an efficient form for the weighting procedure of Zalesak because they differ only in a mulitplication factor (see (49) and (50)). The weighting procedure for two given schemes is clearly presented in [21], and it is applied without change to the two schemes here. The superiority of the solutions of the hybrid scheme, resulting from the weighting of the two schemes, is demonstrated in Fig. 4. 5. Conclusions The one-dimensional, inviscid, compressible-flow equations possess a split conservation form. This has been shown by using the flux vector splitting. Based on this split form, two solution schemes that belong to a class of one-step, explicit, conservative schemes have been designed. One is the first-order accurate ‘upwind’ scheme and the other is a scheme similar to the second-order Lax-Wendroff scheme. These two schemes are more accurate versions of those of Rusanov and Van Leer. Solutions for an ideal shock-tube problem demonstrate the well-known features that the first-order scheme is very diffusive, while the higher-order scheme generates post-shock oscillations but yields more accurate results in the smooth regions. The two schemes are well suited to the nonlinear weighting procedure developed by Zalesak. The hybrid scheme that results from the weighting of the two schemes essentially eliminates post-shock oscillations while maintaining the higher-order accuracy in the smooth regions. The hybrid scheme can be extended to higher spatial dimensions through time-splitting. Acknowledgment I would like to thank Dr. He&i Tamm

for bringing

the work of Zalesak

to my attention.

References [I] D.L. Book, J.P. Boris and K. Hain, .I.

Comput. Phys. 18 (1975)

[2] J.P. Boris and D.L. Book, J. Comput. Phys. II (1973) [3] J.P. Boris and D.L. Book, J. Comput.

248-283.

38-69.

Whys. 20 (1976) 397-431. [4] M.B. Carver, J. Comput. Phys. 35 (1980) 57-76. [5] R. Courant and K.O. Friedrichs, Supersonic Flow and Shock Waves (Interescience, New York, 1967) Section 37. [6] H. Eves, Elementary Matrix Theory (Allyn and Bacon, Boston, 1966) Section 5.1. [7] S.K. Godunov, Mut. Sb. 47 (1959) 271-306; also Cornell Aeronautical Lab. Transl. [8] W.T. Hancox and S. Banerjee, Nucl. Sci. Eng. 64 (1977) 106-123. [9] A. Harten, J.M. Hyman and P.D. Lax, Comm. Pure Appl. Math 29 (1976) 297-322. [lo] A. Jeffrey, Qunsilineur Hyperbolic Systems and Wuoes Research Notes in Mathematics (Pitman, London, 1976) Section 2.1. [I I] G.A. Korn and T.M. Korn, Muthemuiicul Handbook For Scientists and Engineers (McGraw-Hill, New York, 1968) Section 4.5.5. [12] P. Lax, Comm. Pure. Appl. Math. 7 (1954) 159-193. [13] P.D. Lax and B. Wendroff, Comm. Pure. Appl. Math. 13 (1960) 217-237. Eds., Advances in Computer Methods for Partial Differential [14] B.H. McDonald, in: R. Vichnevetsky and R.S. Stepleman, Equations-III (IMACS, 1979). [15] G. Moretti, Comput. Fluids 7 (1979) 191-205. [16] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial-V&e Problems (Wiley, New York, 1967) Section 12.7 and Chapter 4. [17] V.V. Rusanov, Z. Vych. i Mat. Fiz. I (1961) 267-279; also NRC of Canada Technical Translation 1027 (1962). [18] G.A. Sbd, J. Comput. Phys. 27 (1978) 1-31. [ 191 J.L. Steger and R.F. Warming, NASA-TM-78605, Ames Research Center, Moffett Field, CA (1979). [20] B. Van Leer, J. Comput. Phys. 3 (1969) 473-485. [21] S.T. Zalesak, J. Comput. Phys. 31 (1979) 335-362.