NASA Contnctor Report 177984 

lCASE REPORT NO. 85-37 

ICASE 

NASA-CR-177984 

19860004488 

STABLE BOUNDARY CONDITIONS AND DIFFERENCE SCHEMES 

FOR NAVIER-STOKES EQUATIONS 

Pravir Dutt 

NASA Contract No. NASl-17070 

August 1985 

INSTITUTE FOR CO~WUTER APPLICATIONS IN SCIENCE AND ENGINEERING 

NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23665 

111111111111111111111111111111111111111111111 

NF010ll 

rr'.!'" ...... , C"'pV L : !,~ ~ k HI" U I 

NOV 2 51985 

.IJ.i;Sl.[Y hCSEARCH CErJTER 
STABLE BOUNDARY CONDITIONS AND DIFFERENCE SCHEMES 

FOR NAVIER-STOKES EQUATIONS 

Pravir Dutt 

Institute for Computer Applications in Science and Engineering 

ABSTRACT 

The Navier-Stokes equations can be viewed as an incompletely elliptic 

perturbation of the Euler equations. By using the entropy function for the 

Euler equations as a measure of 'energy' for the Navier-Stokes equations, we 

are able to obtain nonlinear 'energy' estimates for the mixed initial boundary 

value problem. These estimates are used to derive boundary conditions which 

guarantee L2 boundedness even when the Reynolds number tends to infinity. 

Finally, we propose a new difference scheme for modelling the Navier-Stokes 

equations in multidimensions for which we are able to obtain discrete energy 

estimates exactly analogous to those W8 obtained for the differential 

equation. 

Research was supported in part by the National Aeronautics and Space 

Administration under NASA Contract No. NASl-17070 while the author was in 

residence at ICASE and in part by NSF Grant No. 82-00788, ARO Grant No. DAAG 

29-82-0090, and NASA Grant No. NAGI-506 while in residence at the University 

of California, Los Angeles in connection with his Ph.D. Thesis. 

i 
INTRODUCTION 

For computational problems involving the Navier-Stokes equations, it is 

necessary to limit the domain of computation and introduce artificial boundary 

conditions. Naturally, we would like these boundary conditions to be stable, 

compatible with weak boundary layers, and to remain valid even when the 

Reynolds number tends to infinity. Such a set of boundary conditions were 

proposed by Gustaffson and Sundstrom in [4]. They used energy estimates on 

the linearized Navier-Stokes equations to obtain boundary conditions of 

maximal dissipative type. In this report we define an 'energy' in terms of 

the entropy function for the Euler equations and obtain fully nonlinear 

'energy' estimates from which we are able to extract a family of boundary 

conditions with the above properties. An attractive feature of these boundary 

conditions is that they are easy to implement and can be expressed in terms of 

the physics of the problem. 

The Navier-Stokes equations are an incompletely elliptic perturbation of 

the Euler equations -- which are themselves a hyperbolic system of conservation laws with entropy functions. It was observed by Mock [5] that by 

introducing the gradient of the entropy as a new variable a system of 

hyperbolic conservation laws can be reduced to a symmetric, hyperbolic system 

in terms of this new variable. Further, Harten [5] showed that if the 

dissipative terms in the Navier-Stokes equations are rewritten in terms of 

this new variable then the matrix coefficients of the dissipative terms have 

certain symmetry properties. We are able to show that the augmented mat rix 

formed from these matrix coefficients is, in fact, negative semidefinite. 

This observation is crucial to the energy estimates we obtain for the NavierStokes equations. 
-2-

This leads us to propose a new difference scheme for modelling NavierStokes equations in multidimensions. We are able to obtain discrete 'energy' 

estimates -- which are exact analogs of the 'energy' estimates we obtained for 

the differential equation -- at the semidiscrete level, even for meshes with 

unequal mesh widths. Thus we are able to propose boundary conditions and a 

difference scheme for the Navier-Stokes equations which give a priori 

boundedness of 'energy' for all time. 

This report is organized as follows: In Section 2 we define the NavierStokes equations and obtain the necessary results to derive the 'energy' 

estimates of Section 3. In Section 4 we propose a family of 'stable' boundary 

conditions and relate them to the physics of the problem. In Section 5 we 

propose a new method for differencing the Navier-Stokes equations in 

multidimensions and obtain discrete 'energy' estimates for our difference 

scheme. Finally in Section 6 we obtain stable boundary conditions for the 

difference scheme and conclude by displaying some numerical simulations in 

Section 7. 

2. PRELIMINARIES 

We consider systems of hyperbolic conservation laws of the form: 

(2.1) 

Here q(x,t) 

function of n 

q + 

t 

d 

I 

i=l 

o. 

1S . an n co 1 umn vector 0 f k un nowns, fi(q) is a vector valued 

components, x and 1 d f = (f , ••• ,f ). 
-3-

We can rewrite (2.1) in matrix form: 

(2.2) 

d i 

q + L A (q)qx. = 0 

t i=l 1 

where Ai(q) The system (2.1) is called hyperbolic if the matrix 

(2.3) 

has real eigenvalues and a complete set of eigenvectors for all real wi. 

Following Mock, a scalar function V(q) is an entropy function for (2.1) 

if: 

i) V satisfies 

(2.4) 

is some scalar function called entropy flux in the xi 

direction. 

ii) V is a convex function of q. 

It follows from (2.4) upon multiplying (2.1) by 

solution of (2.1) also satisfies: 

(2.5) 

where F 1 d (F , ••• ,F ). 

d 

Vt + L 

i=l 

o 

that every smooth 
-4-

The Euler Equations of Gas Dynamics 

Description of variables: 

p denotes density, 

u denotes velocity in the x direction, 

v denotes velocity in the y direction, 

w denotes velocity in the z direction, 

m u' m v' and m are the components of momentum in the 

w 

x, y and z 

directions respectively, 

T is the temperature, 

p is the pressure, 

U is the thermodynamic entropy, 

E is the energy, 

R is the universal gas constant, 

Y is the ratio of specific heats. 

Note that we shall use (x, y, z) and (xl' x2' x3) interchangably to denote 

... the spatial vector x. 

We shall also need the following thermodynamic relations: 

(~ -

2 

m

2 

(X - 1) + m!l) (m + 

T u v 

R 2p2 

(E - (m2 

+ m2 + m2l) (y - 1) u v w 

P 2p 

2 2 2 ( (m + m + m l) U u v w log E - 2p - y log p log P - Y log p. 

up to an additive constant. 
-5-

T, P and p will always be restricted to be positive because of obvious 

physical considerations. q will always denote the vector: 

p p 

m pu u 

q m = pv v 

m pw w 

E E 

The Euler equations are of the form: 

where 

pu pv pw 

2 

+ p pvu pu pwu 

f1 f2 2 

+ p puv f3 pv pwv 

2 

puw pvw pw + p 

(E + p)u (E + p)v (E + p)w 

We shall write the Euler equations in operator notation as: 

(2.6) Eq o. 

The Euler equations have a family of strictly convex entropy functions defined 

by 

V(q) -ph(U). 

The preferred entropy function in most physical applications is: 
-6-

V(q) -pU 

(m 2 2 + m + m ) 2) u v w 

-p log (E - 2p + yp log p. 

The entropy flux functions turn out to be: 

(E -

(m2 + 

2 

+ m;J) Fl 

m 

U log u v -m -m + ym log p u u 2p u 

F2 -m U 

v 

It should be noted that the entropy function V(q) is strictly convex but may 

be nonpositive in general. 

Navier-Stokes Equations 

We shall denote the Navier-Stokes equations in operator notation as: 

(2.7) Nq o 

where 

Nq Eq + (-D)q, 

where (-D)q 

We can represent (2.7) in the alternative form: 

3 . 

Nq Eq + 2 (h1 )x. 

i=l 1 
where 

and 

0 

a 

xx 

hi a 

= yx 

a zx 

a 

xx 

where: 

-7-

3 

I 

j=i 

Aij(q)q 

x. 

J 

are as follows: 

u + a v + a 

yz zx 

h3 

e 

xx 

e = 

yy 

a 

zz 

a xy 

e 

yz 

a 

zx 

w - k.£! 

ax 

0 

e 

xz 

e 

yz 

e 

zz 

e u + e 

xz yz 

h2 

v + e zz 

o 

a xy 

ayy 

a zy 

a u + a 

xy yy 

w - k aT 

az 

-2~ --

au A(~ + av + aW) 

ax ax ay az 

-2~ --

av A(aU + av + aw) 

ay ax ay az 

-2~ --

aw A(~ + av + aw) 

az ax ay az 

a - (~+~) yx U ay ax 

e _ (~+ aw) 

zy ~ az ay 

a = _ (aw + au). 

xz ~ ax az 

v+a w-k1! 

zy ay 
-8-

Here the parameters ~,A and k are as follows: 

~ is the shear coefficient of viscosity, 

A is the second coefficient of viscosity, 

k is the coefficient of thermal conductivity. 

Clearly 

(-D)q 

We write the Navier-Stokes equations as 

(2.8) Nq Eq + (-D)q 

where E is the differential operator corresponding to the Euler equations 

and D is the elliptic perturbation (depending on the parameters ~,A and 

k) due to the dissipative terms in the Navier-Stokes equations. 

The entropy function V(q) = -pU for the Euler equations is not positive 

valued in general. Since we wish to interpret the entropy as a measure of 

'energy' of the system, we can normalize V(q) and define a new entropy V(q) 

which has the properties: 

(i) V(q) a Yq 

(ii) V(q) a <=) q = q for some fixed q. 

For the entropy V is not altered by adding to it an ~rbitrary inhomogeneous 

linear function. Hence we define V(q) as follows: 
-9-

(2.9) 

The associated entropy flux functions are given by 

(2.10) 5 av - i i - .L aq. (q)(f (q) - f (q» .• 

J=1 J J 

We choose as our rest state q = (p, m

u

' m

v

' m

w

' E) such that: 

p>O 

u = v = w = 0 

T > o. 

In a later section we shall put further restrictions on p and T. With this 

choice of q we obtain: 

V(q) p[u - U + (y - 1) (u2 + v

2 + w

2) - y + _TT] + (y - l)p 

2RT 

Or more compactly: 

~1 ~ - F (q) = u(V(q) - (y - 1)p) 

v(V(q) - (y - 1)p) 
-10-

We make the change of variables v = V 

q 

terms of v. Here 

a 

b 

v = c 

d 

e 

Substituting these relations in 

(-D)q 

yields 

(-D)q (-D)v 

and rewrite the operator Dq 

where Cij(v) Aij(q)q. Harten [5] observed that the matrix coefficients v 

Cij(v) satisfy the symmetry relation: 

For the energy estimates we shall derive in Section 3 we need to show that 

(2.11) 

3 3 L L (~i)T Cij(v)(~j) ~O IV ~i, ~jEJR? 

i=1 j=1 

in 

Clearly this is equivalent to proving that the augmented matrix C defined by 
-11-

Cll c

12 

cl3 

C = C

21 C22 c23 

c

31 c32 c33 

is negative semidefinite. 

Observe that Cij = (Cji)T C is symmetric. Hence to check that C is 

negative semidefinite it is enough to perform an LU decomposition of C and 

then check that the diagonal elements of the upper triangular matrix U are 

negative. Thus we obtain: 

(2.12) (2.11) holds <=> A + ~ ~ o. 

Typically, for physical fluids A < 0 and ~ ~ o. Further for most 

fluids under nonextreme flow conditions the following relation holds: 

y~ > A + 2~ _> ~. Pr -

Here Pr denotes the Prandtl number. So (2.11) holds under these conditions. 

3. NONLINEAR ENERGY ESTIMATES 

Let n denote a bounded domain in ~ and let an denote the boundary 

of n. Consider the mixed initial boundary value problem in n: 

(3.1) Nq o -... Y x e: n, t > 0 

where N is the Navier-Stokes differential operator with initial condition: 
-12-

(3.2) 

and boundary condition: 

-30. 

(3.3) Bq o Y x e: an, t > 0 

where B is a boundary operator. 

Define: 

(3.4 ) S(t) '" ~ -30. 

= J V(q(x, t»dx. 

n 

We claim that S(t) gives us an estimate of the energy of the system at 

time t. Note that V(q) is a strictly convex, non-negative function of q, 

i.e., 

i) V(q) 2. 0 

ii) V(q) = 0 <=> q q 

iii) Vqq > cI where c > 0, at least in some appropriate physical domain. 

In particular, (iii) implies: 

(3.5) '" - 2 V(q) ~ allq - qll 

for some a > O.Hence we conclude that S(t) has the following properties: 

i) S(t) 2. 0 

-30. 

ii) S(t) o <=> q(x, t) = q Y x e: n 

iii) S(t) ~ k, where k is a constant => 
-13-

~ 

if q(x, t) lies in a domain where an inequality of the form (3.5) holds. 

We wish to study the time evolution of S(t). Recollect that the NavierStokes equations are: 

Nq = Eq + (-D)q 0 

where 

3 3 

(-D)q L L 

i=1 j=1 

Making the change of variables v = V 

q gives: 

3 3 

(-D)q = (-D)v L L 

i=1 j=1 

where the matrix coefficients Cij(v) satisfy the condition: 

(3.6) 

3 3 

L L ~iT cij(u)~j < 0 ¥ ~i, ~j E i> • 

i=1 j=1 

This leads us to the following theorem. 

Theorem 3.1: Consider the mixed initial boundary value problem (3.1)-

(3.3). Let 1 d (~ ' •• "~) denote the outward unit normal to the boundary an 

and let ",1 ",d F , ••• , F denote the entropy flux functions as in (2.10). Then any 

piecewise smooth solution of (3.1)-(3.3) satisfies the energy estimate: 

(3.7) _~ J \ "'F1 r1 + \ \ '" 1 J() dS 

( 

3 0 0 3 3 0 i 0 ) 

L ~ L L Vq l; A q qx

J

o 

do. 

dt an i=1 i=1 j=1 
-14-

Here do is an element of surface area. 

~ 

Proof: We prove it here for the case q(x, t) is 'smooth enough.' We 

have 

(3.8) Nq q + 

t 

Premultiplying (3.8) by V we get q 

o. 

V q + I V (fi(q) )x. + I V ( I A

ij (q)q) o. 

q t i=1 q 1 i=1 q j=1 Xj Xi 

By (2.4), this reduces to 

(3.9) 

3 

Vt + I 

1=1 

From (3.9) we obtain 

dS 

- £ 11 -= dt 

or 

dS 

-dJI -= dt 

By the divergence theorem 

dS 

dt -!nCt 

~i 

F + 

x. 1 

~i 

F + 

x. 1 

~i i F r;; 

o. 

3 

I V ( 3 I Ail (q)q ) J d-:. 

i=1 q j=1 x. X. J 1 

3 

I v T( 3._ I C1J (v)v ) J dx. ~ . 1 X. xi i=1 J= J 

3 3 T i i- ) + I I v r;; c J(V)V

Xj 

do 

i=1 j=1 
-15-

( 

3 3 T i' ) ~ + f L L v C J(V)V

XJ

, dx. 

n i=l j=l xi 

But by (3.6) 

Hence we obtain 

dS < - f I ~ ~Fi ri + ~ ~ v

T 

L. ., L. L. ., 

ri C ij () v v 1 do. 

dt - an i=l i=l j=l Xj 

This finally yields 

(3.7) dS < _ f (~ ~Fi ri + ~ ~ ~V i ij () ) L. ., L. L. r,; A q qx

J

, do. 

dt - an i=l i=l j=l q 

Remark: Estimate (3.7) is a fully nonlinear 'energy' estimate which • 

holds for a broad class of solutions of the mixed initial boundary value 

problem (3.1)-(3.3). In fact (3.7) is simpler to obtain and broader in scope 

than the linearized energy estimates in current use. 
-16-

4. STABLE BOUNDARY CONDITIONS FOR THE NAVIER-STOKES EQUATIONS 

x 

Figure 4.1 

Since the Navier-Stokes equations are rotationally invariant, we can 

consider a moving coordinate frame (x, y, z) where x points in the 

direction of the inward normal and y and z are tangential to an. Of 

course, we reorient u, v, and w so that u is the component of velocity in 

the x direction, v in the y direction, and w in the z direction. Let 

~ denote the outward unit normal to an. Clearly ~ = (-1, 0, 0). 

Then the nonlinear energy estimate (3.7) takes the simpler form: 

(4.1) 

or 
(4.2) 

where 

or, 

(4.3) 

dS < f 

st - an 

-17-

3 

I A

Ij 

j=1 

q • x. 

J 

Substituting the relations obtained in Section 2 we obtain 

(y-l) 

R [0, u v w 

T' i' i' -(I I - T) I 1 

~v hI key - 1) 1 I aT 

= (-T - =) q R T ax 

o 

2 ~ + A (~ + ~ + aw) II ax ax ay az 

(~+ ~) II ay ax 

(~ + aw) II az ax 

21lu ~ + AU(~ + ~ + aw) 

ax ax ay az 

+ v(~ + av) + u(~ + aw) + k aT 

II ay ax II az ax ax' 

(y -=- 1) (21l + A)U ~ + AU(~ + aw) RT ax ay az 

+ v(aU + ~) + w(aw + ~) I. II ay ax II ax az 

So finally we obtain the needed 'energy' estimate: 
-18-

~ < f u(V - (y - 1)p)da + f k(y.=. 1) (T .=. T) aT da 

dt - an an RT T ax 

(4.4) - f (y - 1) U(211 + A) ~ + A(~ + aw)) da 

an RT ax ay az 

- f (y - 1) Vll(~ + aU)da - f (y - 1) Wll(aw + ~)da an RT ax ay an RT ax az 

where 

(4.5) u(V - (y - 1)'P) [

- (y - 1) (u2 + 2 2) T 1 pu u - u + v + w - y + =T 

2RT 

and 

(4.6) v P [if - u + ~!l (u 

2 

+ v 

2 

+ w 

2

) - y + =TT 1 + (y - 1) p. 

2RT 

So far we have left T and p vaguely defined. We are now going to 

specify T and p, but first we have to make a few assumptions about the 

-l. 

solution q(x, t) to the mixed initial boundary value problem (3.1)-(3.3). 

Assumptions: 

-i -i 

i) 3 p. > 03 p(x, t) > p. ¥ x E n, t > 0 m1n - m1n 

(4.7) ~ ~ 

ii) 3 Pmax > Pmin 3p(x, t) < P

max 

¥ x En, t > 0 

iii) -l. -io 

3 Tmin > 03 T(x, t) > T. ¥ x E n, t > O. - m1n 
-19-

Assumptions (i)-(iii) exclude cavitation, freezing of the fluid, and such 

esoteric phenomena as the formation of black holes. Further, it should be 

noted that if (4.7) is valid V(q) is strictly convex and hence, then we are 

able to obtain boundary conditions which ensure that Set) remains bounded, 

i.e., 

S( t) < k Y t > 0 where k is a constant; then by the discussion 

following (3.5) 

=> (~ IIq(i, t) - 'ill/. di ylz ~ C (another constant) Y t > O. 

Hence we obtain L2 boundedness of the solution for all time. 

Choose 0 < P < P

min such that 

(4.8) v - (y - 1)p > 0 

and 0 < T < T. such that m~n 

(4.9) u - u < 0 

Y q staisfying (4.7). 

It is easy to see that this can always be done. For by (4.6) 

V(q) ~ p[u - u - y + :] + (y - 1)p 

or 

V(q) ~ p [(y - 1) log (~) T - Y +-- log (;)] + (y -

P T 

l)p. 
Clearly 

Hence we get 

T 

T 

-20-

log (~) ~ 1 V T > O. 

T 

Thus if p > p. > 0 we can choose p < p. such that - ml.n ml.n 

V(q) - (y - 1)p > O. 

Hence (4.8) holds. 

Now 

u - u (y - 1) log (~) - log (~). ., p T 

If (4.7) is valid then 

Pmax T . 

U - U < (y - 1) log (--) - log ( ~n). p T 

We can choose 0 < T < T

min such that 

pT. 

(y - 1) log ( m~x) - log ( m~n) < O. 

p T 

Hence (4.9) holds. 

The Navier-Stokes equations are: 

(4.10) Nq Eq + (-D)q 0 
-21-

where E is the hyperbolic differential operator corresponding to the Euler 

equations and D is the incompletely elliptic perturbation, depending on the 

parameters, ll, A, and k, due to the dissipative terms in the Navier-Stokes 

equations. I t is assumed that ll, A, and k are proportional to a small 

parameter E > O. Thus the Navier-Stokes equations may be viewed as an 

incompletely elliptic perturbation of a system of hyperbolic equations. 

The question we are concerned with is: for which boundary conditions are 

the solutions of the perturbed problem (i.e., E > 0) well defined in some 

time interval 0 ~ t ~ TO and are bounded in an appropriate norm as E + O. 

Michelson [8] suggests boundary conditions of the type given below for 

the singular perturbation problem: 

(4.11a) o 

where 

and a = (aI' a 2 , a3 , a 4 ) is a multi-index in some finite set A c ~. 

It is well known that a singular perturbation problem of the type we are 

considering exhibits a boundary layer phenomenon. We wish to choose boundary 

conditions in such a way that a strong boundary layer does not develop. For 

the boundary conditions to be compatible with a weak boundary layer, we need a 

condition of the sort given by: 

(4.11b) 0) o 
-22-

gives a well-posed boundary value problem for the hyperbolic part of the 

Navier-Stokes equation: 

Eq o 

(see [8]). 

Further, from Strikwerda's work on initial boundary value problems for 

incompletely parabolic systems, we know that the number of boundary conditions 

which should be imposed for the Navier-Stokes equations to obtain a well-posed 

problem are: 

5 boundary conditions for inflow boundary and 

4 boundary conditions for outflow boundary. 

At the same time the unperturbed hyperbolic system requires: 

5 boundary conditions for supersonic inflow, 

4 boundary conditions for subsonic inflow, 

1 boundary condition for subsonic outflow, and 

none for supersonic outflow. 

The boundary conditions we are going to impose will be of the form: 

(4.12) g 

where R is a matrix of rank at most 4. 

To get a set of boundary conditions that also works for £ = 0 S must 

be chosen in such a way that Sq = g is a proper set of boundary conditions 

for the unperturbed hyperbolic problem. If Sq = g gives too many boundary 

conditions for the unperturbed hyperbolic problem, the solution will contain 

boundary layer components of the form exp(-x/£) for £ + 0; (see [4]). 
-23-

We now proceed to suggest a set of maximal dissipative' boundary 

conditions for the Navier-Stokes equations which are compatible with weak 

boundary layers. 

Outflow Boundary Conditions 

We need to specify 4 boundary conditions for the perturbed problem. For 

supersonic outflow (-u > c > 0) we need not specify any boundary conditions 

for the unperturbed hyperbolic problem, while for subsonic outflow (c > -u > 

0) we need to specify only one. From (4.4) we have 

dS < J u(v - (y - 1)"P)dcr + J k(y - 1) (T; T) ~~ dcr 

dt - an an RT 

By (4.8) 

- J (y - 1) U«21l + A) ~ + A(av + aW))dcr 

an RT ax ay az 

J (y - 1) (av au) J (y - 1) (aw aU) - vll ax + ay dcr - wll ax + az dcr. 

an RT RT 

J u(V - (y - 1)p)dcr ~ o. 

an 

The boundary conditions we suggest which give a decay of energy are: 
-24-

key - 1) aT 

(Xl T > 0 R ax g1' (Xl 

(211 + A) ~+ ax 

A-+- (av aw) -(X u 

ay az 2 g2 ' (X2 ~ 0 

(4.13) 

( aV au) 11-+- -(X v= 

ax ay 3 

( aw aU) 11 - + - - (X. W 

ax az £+ 

It should be noted that if we put (Xi = 0 we must also put gi = 0 for 

i 1, ••• ,4. Since we want (4.13) to be compatible with the preceding 

discussion, we end up with the following set of maximal dissipative' boundary 

conditions: 

Supersonic outflow 

key - 1) aT 

-- 0 R ax 

(211 + >.) ~+ ax >.(~ + aW) 

ay az 

0 

(4.14) 

11 (~+~) ay ax 0 

(aw + au) 11 ax az 0 
Subsonic Outflow 

(4.15) 

-25-

k( y - 1) !!. = 0 

R ax 

au (av aw) (2j.J + x) - + x - + - - a u = g ax ay az 2 2 

( av + ~) = 0 J.l ax ay 

Inflow Boundary Conditions 

By (4.9) 

Subsonic Inflow 

f pu[U - U]da ~ 0 

an 

Since c) u we have 2 

u < yRT. Hence 

f (y - 1) (pu)u2 da ~ f 

an 2RT an 

(y - 1) (pu)yRTda 

2RT 

=) dS < _ f (y - 1) U(2J.l + x) ~ + x(~ + aW))da dt - an RT ax ay az 

+ f (y - 1) v ( (pu) v _ j.J(av + ~))da an RT 2 ax ay 

+ f (y - 1) w(PU)W _ j.J(aW + aU))da 

an RT 2 ax az 

+ f .!. [ Cf - T) k(1 - 1) aT + (1

2 - 1 + 2) ( )T] d 

an T T R ax 2 pu a. 
-26-

The boundary conditions we could specify are: 

(2)..1 + A) ~+ ax A(~ + aw) ay az 0 

(4.16) (av aU) )..I-+--av ax ay 3 g3 

(aw aU) )..I-+--aw ax az 4 g4 

key - 1) aT 

as T R ax gs' 

Here we must impose the conditions 

2 

> (y - y + 2) as 2 gl + £ for some suitable £ > 0 

depending on T. 
-27-

Supersonic Inflow 

dS < f pu[u - U]dcr + f (X - 1) u[~ u -I (2lJ + A) ~ + A(av + aw)ljdcr dt an an RT 2 ax ay az 

+ f (X - 1) v[<PU) v - Jl(av + ~)dcrj+ f <X - 1) w[<PU) w - (aw + ~)jdO' an RT 2 ax ay an RT 2 Jl ax az 

+ f l..[<T - T) k(X - 1) aT + (PU)TjdO. 

an T T R ax 

The boundary conditions we could specify are: 

<4.17) (av + _aU) Jl ax ay - a3 v = g3 

( aw aU) lJ-+- -a w=g ax az 4 4 

Here the following conditions must hold: 

and as > gl + E for some suitable E > O. 
-28-

Discussion of the Boundary Conditions 

The outflow boundary conditions we have proposed are of the following 

type: 

(4.18) aT 

ax o. 

(4.18) says that the computational boundary corresponds to an insulated wall 

and there is no conduction of heat across it. 

(4.19) o 

(4.19) asserts that there is no shearing of the fluid either in the normal or 

tangential direction against the computational boundary. 

For an inflow boundary one of the boundary conditions we have specified 

is: 

(4.20) pu g. 

In other words, we specify the momentum influx across the computational 

boundary. 

For the temperature component of the state vector, we have the boundary 

condition: 

(4.21) k(y - 1) aT 

R ax as T 
-29-

(4.21) has an interesting physical interpretation. For injection through a 

porous wall into the main stream (sometimes called transpiration) the injected 

fluid may be, say, a coolant at a temperature considerably different from the 

wall, and one needs to consider an energy balance at the wall. A good 

approximation for coolant injection is to use the boundary condition proposed 

by Roberts [in Truit (1960, chapter 11)]. 

(4.22) Injection: k 2.! '" u c (T - T ) ax Pw w p w coolant 

where is the momentum flow of coolant per unit area through the wall, 

Tw is the temperature of the wall, Tcoolant 

;;Iw is the temperture gradient 

information refer to [12, chapter 1.4].) 

coolant and 

is the temperature of the 

at the wall. (For more 

By choosing as and gs appropriately (4.21) can be geared to satisfy 

(4.22). The physical effect of such a boundary condition would be to make the 

temperature of the fluid within the computational boundary stabilize at 

Tcoolant over time. When we discretize the boundary conditions, however, we 

can choose gs so that the discretized version of (4.21) is of extrapolation 

type. This has the effect of allowing the system to evolve as if there were 

no boundary. 

For the velocity component of the state vector, the boundary conditions 

we have specified are: 
-30-

(2~ + A) ~ + A(~ + aw) - a 2 u 

ax ay az 

(4.23) ( ~-+- aV aU) -a v 

ax ay 3 

( ~-+- aW aU) -a w 

ax az 4 

For one-dimensional fluid flow the boundary conditions (4.23) would take 

the simple form: 

(4.24) (2~ + A) au 

ax 

For perfectly diffuse reflection of a gas or fluid against a solid wall 

we have the boundary condition (after Maxwell in 1879): 

(4.25) u 

w 

'" R.(~) ax w 

where Uw is the velocity of the fluid at the wall, R. is the mean free path 

of the gas or fluid at the wall and ~) is the velocity gradient at the 

ax w 

wall (see [12]). 

By choosing and appropriately we can choose (4.24) to satisfy 

(4.25). Since we want the boundary conditions to be nonreflecting, however, 

this would be a bad choice indeed. If we let the computational boundary move 

with the same velocity as that of the fluid at the boundary interface, we 

should not have any diffuse reflection. And, in fact, we can discretize 

(4.24) to achieve this exactly. The discrete boundary conditions so obtained 

are of extrapolation type and the effect of the computational boundary is 

minimized. 
-31-

A last remark is that the boundary conditions proposed are of maximal 

dissipative type and such boundary conditions are intrinsically radiative. 

5. STABLE SEMI-DISCRETE DIFFERENCE SCHEMES FOR THE DIFFERENTIAL EQUATION 

The stable difference scheme we are going to propose is valid for the 

Navier-Stokes equations in multidimensions. However, to avoid getting drowned 

in notation we restrict the discussion that follows to Navier-Stokes equations 

in one space dimension. 

Consider the incompletely parabolic system obtained from a hyperbolic 

system of conservation laws with entropy: 

(5.l) qt + f(q) + (A(q)q ) x x x 

o. 

Let V(q) be a normalized entropy function as in (2.9) for the 

hyperbolic system of conservation laws: 

(5.2) 

and let v

T 

= V 

q 

q + f(q) = 0 

t x 

We can rewrite (5.1) as: 

(5.3) q + f(q) + (C(v)v ) t x x x 

where we assume C(v) ~ O. 

o 
-32-

Divide the x axis into a mesh of points 

and ~x is the mesh width. 

q.(t) = q(x., t). J J 

Let Consider 

{x. } where 1. • IN 

1.e: 

the semidiscrete 

i~x 

difference 

approximation to the hyperbolic system of conservation laws (5.2): 

(5.4) 

where rand s are integers > 0 and h(qj-r,·.·,qj+s) is the numerical 

flux function corresponding to the flux function f(q), i.e., h(q. , ••• ,q. ) J-r J+s 

satisfies the consistency condition: 

(5.5) h(q, ••• ,q) = f(q). 

Let the order of accuracy of the semidiscrete difference approximation be a 

where we can take a > 1. 

Further, suppose (5.4) satisfies a semidiscrete entropy inequality: 

(5.6) 

where H(q ••• q ) is the numerical entropy flux function corresponding j-r' , j+s 

to the entropy flux function F(q), i.e., H(q. , ••• ,q.+) J-r J s satisfies the 

consistency relation: 

(5.7) H(q, ••• ,q) F(q). 
-33-

We propose the following semidiscrete approximation for the Navier-Stokes 

equation (5.1): 

(5.8) 

where 

It is easy to see that the order of accuracy B of (5.8) is given by 

B = min(a, 2). Hence if (5.4) is a second-order accurate approximation to 

(5.2) then (5.8) is a second-order accurate approximation to (5.1). In any 

case, (5.7) is at least a first-order accurate approximation to (5.1). 

Remark: The only second-order accurate semidiscrete difference scheme 

for a hyperbolic system of conservation laws which satisfies an entropy 

condition of the type (5.6) of which we are aware of is one proposed by Osher 

in [10]. 

Lemma 5.1: The semidiscrete difference approximation (5.8) to the 

Navier-Stokes equations (5.1) satisfies the entropy inequality: 

(5.9) 

Fix ~t > O. Let 

q.(tot) * 

J 

~t qJ' - AX ~ - (h(q. , ••• ,q.+ )). L> J-r J s 

Let E q (t) denote the solution to the semidiscrete difference scheme 
(5.10) 

with q~(O) = q .• 

J J 

entropy inequality: 

-34-

dnE I::J.- E E .::...:1.- + - h(q ••• q ) dt I::J.x j-r' , j+s 

Then q. * (I::J.t) J 

o 

Since (5.10) satisfies the 

¥ sequence {I::J.tkl 3 I::J.t k + 0 3 a sequence {Ekl 3 Ek + 0 and the discrete 

entropy inequality: 

(5.11) 

V(q~(I::J.tk)) - V(q.) 

J I::J. J + ~- [H(q. , ••• ,q.+)] < Ek tk uX J-r J s -

holds. 

We now define a discrete difference approximation for the Navier-Stokes 

equation: 

Let qj(t) denote the solution to the semidiscrete difference scheme: 

Then it is obvious that q. (!::J.t) 

J 

entropy function 

q. (I::J.t) + O(I::J.t 2). Define an approximate J 

W(q.(I::J.t)) - '" * T I::J.t = V(q.(I::J.t)) - v . .- I::J. 

J J J uX 
-35-

By (5.11) 

Atk 

W(qJ.(Atk » <V(q.) --A - [H(q ••• q )] J Ax j -r ' , j +s 

or 

(5.12) 

By a Taylor series argument: 

And since q. (At) 

J 

W(q. (lit» 

J 

by letting Atk + 0 in (5.12). 

~ 2 V(q.(lIt» + O(lIt ). J 

Define a discrete version of the energy S(t) by: 

N 

(5.13) S(t) 2 

j=l 

We can now easily obtain a discrete version of the energy estimates of Section 

2. (Henceforth we shall denote h(q ••• q ) j-r' , j+s and 
-36-

Theorem 5.1: Consider the semidiscrete difference approximation: 

(5.8) 

and let 

dq. 

--=.J.. = dt 

S(t) 

N 

I 

j=l 

flxV(q.). J 

Then the following 'energy' estimate holds: 

(5.14) 

Proof: By Lemma 5.1 

Hence 

dS 

-= dt 

dV(q. ) 

J < 

dt 

using summation by parts. 

C < 0 => 

Since 
-37-

Hence 

Remark: The generalization of this theorem to multidimensions is 

obvious. 

We now extend these results to uneven meshes, though restricting (5.4) to 

the Godunov scheme. 

The Godunov Scheme 

Consider a grid of points {x.}. lli and define ~X. J J€ J Define 

x. 1 + x. J- J 

2 and OX. 

J 

~X'l+~X' J- J 

2 

6x. 1 6x. J- J 

I" ~I" "I I I I I I 

I #1 1 

~ qj+1 q. 1 qj J- . . x. 1 x. x.+ 1/ xj + J- ~- 1 /2 J ~ 2 

' ... .' OXj 

Figure 5.1 
-38-

Consider the hyperbolic system of conservation laws: 

(5.15) with Riemann initial data: 

Then the Godunov flux 

(5.16) 

{ 

q. 

q(x, t)lt=o = J 

qj+l 

f~ J 

f(q(x 

for x < ° 

for x > 0 • 

is defined as: 

0, t 

where q(x, t) is the solution to the Riemann problem (5.15). 

The semidiscrete version of the Godunov scheme is: 

(5.17) 

where is defined in (5.16). The semidiscrete Godunov scheme 

satisfies the entropy inequality: 

(5.18) 

where 

and is defined as: 

dV(q.) J 

dt 

is the numerical entropy flux for the Godunov scheme 

0, t 
-39-

where q(x, t) is the solution to the Riemann problem (5.15). 

A natural modification of (5.8) for uneven meshes is: 

(5.19) 

Clearly (5.19) is a first-order accurate semidiscrete difference approximation 

to (5.1). We claim that (5.19) satisfies the entropy inequality: 

(5.20) dV(qj) 6,- G T 6,- 6.+ v. 

d + -1:- (F (q., q.+l») + v. -1:- (C(v.+lf ) J) < o. 

t uXj J J J uXj J 2 6.Xj -

The proof is exactly the same as in Lemma 5.1. 

Define a discrete version of the 'energy' S(t) for the uneven mesh 

{x.} by: J 

N 

(5.21) S(t) 2 

j=l 

ox. V(q.). J J 

Then we obtain the following theorem. 

Theorem 5.2: Consider the semidiscrete difference approximation: 

(5.19) 

and let 

~ = 6.- (f~) dt - oXj J 

N 

S(t) = 2 

j=l 

ox. V(q.). J J 
-40-

Then the following 'energy' estimate holds: 

(5.22) 

Proof: The proof is exactly as in Theorem 5.1. 

Remark: Theorem 5.2, though a simple extension of Theorem 5.1, will prove 

very useful in proposing stable boundary conditions for semidiscrete 

difference approximations to the Navier-Stokes equations. 
-41-

6. STABLE BOUNDARY CONDITIONS FOR THE DIFFERENCE SCHEMES 

",t.x'~ t.x ~ ~ t.x ~ ~ t.x ~t.x'-.! I I 

I I I I I I I I 

I I I I I I I I 

I I I I I I I I 

I I '~ I I I I I 

( 1 

• • > 

Figure 6.1 

Suppose we have the boundary condition at x = 0: 

a aT T 

ax g where a can be arbitrarily small. 

Then a discrete version of this boundary condition would be: 

=) T = o 

aT

1 - gt.xo 

(a + t.xo) • 

So if t.xo » a =) TO + - g and hence could give rise to a steep 'numerical 

boundary layer.' To avoid this we should have t.xO of the same order of 

magnitude as a or smaller. If, however, we choose an even mesh this becomes 

computationally infeasible since a can be arbitrarily small. 

To overcome this difficulty we propose the following: Divide the 

interval [xO' ~+1] into N + 2 points xO' x1 ,···,xN, ~+1. As before let 

t.xJ = xJ+1 - xJ. Choose 
t.x. 

1. 

-42-

t.x for i 1, ••• ,N-1 

where we may choose t.x' < < t.x. (See Figure 6.1.) So the mesh points are 

evenly spaced in the interior, but Xo is close to xl and xN+1 is close 

to xN· 

Define ox = t.x + t.x' Henceforth we shall restrict to be 2 q 

q (p, mu ' E)T. 

The semidiscrete difference approximation for the mixed initial boundary 

value problem is tailored according to (4.19): 

dq1 t.- G + _1_ [ t. + vI t. 

+ Vo] --+- [f (q1,q2)] C(v3/ 2) t.x - C(v 1/2 ) = 0 dt ox ox t.x' 

dqJ t.- G t.-

[C(VJ+ 112) 

t. + VJ 1 --+- [f (QJ,QJ+1)] +- t.X

J = 0 for j=l, ••• , N-1 dt t.x t.x 

dQN t.- G 1 [C(VN+ 112) t. 

+ vN t. + vN- 1 ] --+- [f (qN,qN+l)] +- t.x' C(vN_ 112 ) t.x dt ox ox 

(6.1) with initial condition: 

and boundary conditions: 

where Band B are boundary operators. 

0 
-43-

In general, B and 13 may be an underdetermined set of boundary 

conditions. For example, B(qO' ql' q2) may be of the form: 

and is unspecified. 

If B is a 'legitimate' boundary operator we can connect ql to qo by just 

two waves, say a 2- contact and a 3- shock or rarefaction, such that 

B(qO' ql' q2) = 0 holds. We make this statement more precise by examining 

the boundary Riemann problem. 

Suppose we have a boundary (a wall) given by x = st. We want to give a 

number, m, of nonlinear boundary conditions depending on q at the boundary: 

h(q)/x=st 0 

h: JRl ... JRIl, he:C"" 

(6.2) 

and initial condition: 

q(x, t)/t=O = qr (a constant state) for x > O. 

The underlying differential equation is 

(6.3) o 

where 
-44-

A(q) af 

aq 

We want to construct similarity solutions for this initial boundary 

problem. Clearly, for the problem to be well-posed we have to choose qr so 

that 

where Xi(q) are the eigenvalues of the matrix A(q) and m = n - k. 

We can connect qr to the boundary by n - k waves. The rarefunction 

and shock curves of the jth family determine a wave curve: 

«q E

r ), E . > 0 r' J 

Wj(q E j ) 

= sj(q, E.), 

r' 

E . < 0 

r J J 

where Wi is twice continuously differentiable in both its arguments (see 

[2] ) • 

Given qr and sufficiently small parameters Ek+1,· •• ,En we define 

~n ~k+l 

states q , ••• ,q inductively by: 

~n n 

q W (q , E ) r n 

~j-l Wj q (qj , E . ) for j n, ••• ,k + 2. 

J 

Set: 

qR, _ ~k+l ( - q qr' E: , ••• ,e:). k+l n 
-45-

Then the boundary condition h(q) o becomes m equations in the n - k 

unknowns 

By the implicit function theorem we will be able to solve for in 

terms of in a manner if m = n - k and the Jacobian has rank 

m. Differentiating we get the local solvability condition: 

(6.4) o~ (r~+l, ••• ,tn) has rank m = n - k 

oq I I q=qr 

where in the ith eigenvector of the matrix A(q) corresponding to the 

eigenvalue Ai. The construction is described in Figure (6.2). 

t 

Figure 6.2 

-n+l 

qr - q. 
-46-

(6.4) is also a necessary condition for the linearized initial boundary value 

problem to be well-posed. 

We relate this now to the boundary operator B(qO' q1' qZ) = O. Let 

denote and denote corresponds to 

the m boundary conditions: 

(6.5) h(q) I x=O o 

and the initial condition for the initial boundary Riemann problem (6.2) is: 

q(x, t) I t=O 

So B(qO' q1' q2) is a 'legitimate' operator if (6.4) holds. Assuming that 

it is, we can connect to a state such that 

holds. 

Returning to the difference scheme (6.1) we can now define the Godunov 

flux using this construction. Suppose, at 

and we have the boundary operator B(qO' q1' q2) = O. 

Consider the Riemann initial boundary value problem: 

q(x, t)lt=O for 

and boundary condition: 

h(q) = 0 for x o 

as in (6.5). 
-47-

f>x !:.xl2 

,. ~ ., I I I I I I 

qo dJ~ ql I Xo 1'\ xl ox • Ax+Ax' ~31 I 2 

I 2 I 

'-f ~ 

Figure 6.3 

Then if (6.4) holds we can join ql to a state qo in a small enough 

neighborhood of ql so that B(qO' ql' q2) = o. 

In general, the number of boundary conditions, m, specified by the 

boundary operator B will be more than the number of boundary conditions 

needed for the hyperbolic Euler equations since we are solving the NavierStokes equations, Le., m > n - k. Hence, when we connect ql to a state 

qo lying to the left of the line x = x~2 some of the waves will lie in the 

region x < x 1/2 • 

boundary. 

This corresponds to letting waves be radiated at the 

We know that we can connect ql to a state qo such that B(qO' ql' 

q2) = 0 if qo and ql are sufficiently close. We chose an uneven mesh 

precisely for this reason. We set llx = llx = llx'" o N and for 

i = 1, ••• ,N - 1 where llx'" could be arbitrarily small. In fact, by choosing 

llx'" small enough we can make qo lie as close to ql as we want. This has 

the effect of making the waves in the cell bounded by x = XU and x = Xl 

very weak and hence we conclude that 
-48-

Recollect that the energy estimate we obtained for an uneven mesh was of 

the form: 

(6.6) 

Since as f,x' -+- 0 the estimate (6.6) reduces to: 

(6.7) 

for sufficiently small f,x'. 

One further point that should be noted is that when we take a fully 

discrete version of our semidiscrete scheme the CFL condition we have to 

impose to prevent wave interactions is: 

(6.8) 

where ~(A) is the spectral radius of the matrix A(q). Clearly (6.8) is not 

an unduly restrictive condition. 

For the numerical simulation presented in the next and last section, we 

nondimensionalize the Navier-Stokes equations. 

be a reference length, Pf a reference density, and uf a 

reference velocity. Define: 

Let L 

x * x/L, 

* t = tuflL, 
-49-

* 2 

p = p/Pf uf ' 

* 2 T TR/uf , 

* P P/Pf' 

* 2 E E/Pf ufo 

We also need the following parameters, which occur extensively in fluid 

dynamics: 

Here Pr and Re 

Pr 

Re = 

]..Ic 

~.....:]..I:;..R-!y~ = --E. 

k(y - 1) k 

are abbreviations for the Prandtl number and Reynolds 

number respectively. Henceforth we drop the for notational 

convenience. For our numerical simulation we use Stokes' assumption: 

The nondimensionalized Navier-Stokes equations then take the form: 

P + (pu) 0 

t x 

4 ~2 T 

() [ ( ) ] a eu aU) + y CJ E t + E + p u x - 3 Re a x ax 7(-y----1:-.):-;P:"'"r--::::R-e ax 2 
-50-

Rewriting the dissipative term (-D)q in terms of our new variable v = V , q 

we obtain: 

0 0 0 

(C(v)v ) 0 -4T -4uT a 

x x 3Re 3Re ax 

-4uT 2 yT2 

0 -4u T 

3Re 3Re (y - ORe Pr 

Our nondimensionalized boundary conditions take the form: 

Supersonic Outflow 

Subsonic Outflow 

Subsonic Inflow 

au 

-= 0 

ax 

aT 

-= 0 

ax 

aT 

- = 0 

ax 

au 

-= 0 

ax 

y aT _ a T = g where 

Re Pr ax 3 3 

where 

2 (y - y + 2) 

a3 > 2 gl 

0 

u 

-

T 

-1 

T 
-51-

Supersonic Inflow 

y ~- Re Pr ax a3 T where 

We now propose a discrete version of these boundary conditions which are 

of extrapolation type and hence minimize the effect of the computational 

boundary. 

Supersonic Outflow 

We can connect qo to q1 without any waves at all by putting 

Subsonic Outflow 

The boundary conditions proposed for supersonic outflow would work 

equally well for subsonic outflow. This corresponds to choosing a 2 o in 

the subsonic outflow boundary conditions. If we choose the 

discretized boundary conditions would be: 
-52-

TO T1 

where we could put 

4 (u2 - u ) 1 

g1 - 3Re !::J.x - u2 u2 • 

Then 

(4/3Re) - u

2 !::J.x')u1 - g1 !::J.x' 

Uo (4/3Re) 

Clearly Uo + u1 as !::J.x' + O. We choose by joining 

wave interactions. 

Subsonic Inflow 

!::J.x' 

!::J.x 

to qo by two 

(T2 

- T

1

) 

g3 = -Y- - u3 T2• Re Pr !::J.x 

2 Here u3 > (y - y + 2)g1/2. 

SUEersonic Inflow 

(pu)O = g1 (pu)1 !::J.x' 

= --[(pu)2 - (pu) d !::J.x 

4 (u1 - uO) 4 (u2 - u 1) 

:nre !::J.x> - u2 u1 g2 = JRe !::J.x - u2 u2 

(T1 - TO) (T - T 1) y - u3 T1 Y 2 

Re Pr !::J.x> g3 Re Pr - u3 T2• !::J.x 
-53-

As mentioned earlier all the boundary conditions are such that qo + q1 

as ~x' + 0 and they are a discretized version of the boundary conditions for 

the differential equations, i.e., 

e:R .l9. + Sq ax g. 

It is easy to verify from (6.7) that for ~x' small enough we get bounded 

growth of the discrete version of the energy S( t) in time by choosing a2 

and appropriately. 

7. NUMERICAL RESULTS 

Throughout the simulations we use the following parameter values: 

Pr 0.7 

y = 1.4 

~x 0.1 

~x' 0.000001 

~t 0.01. 

Simulation 1: For the first simulation we use a very large value of the 

Reynolds number 

We run the program with Riemann initial data: 
-54-

PL = 1.0 

u

L = 3.0 

PL 1.0 

and 

PR 0.4734821 

u

R 2.1393370 

PR 0.3333333. 

At the left boundary our boundary conditions correspond to supersonic 

inflow. At the right boundary we have supersonic outflow boundary conditions. 

The solution to the Euler equations with initial data corresponding to 

the above Riemann problem is a 3 shock moving to the right with a speed 

3.774567. 

The numerical simulations bear this out very well. The boundary 

conditions are radiative and allow the shock to pass through. Further for 

long time the solution stabilizes to a constant state. 

ACKNOWLEDGEMENTS 

The author is grateful to Professor S. Osher for his valuable help and 

encouragement and would like to thank Professor Ralston and Professor Bube for 

helpful discussions. 
P 

3. ij-l 

3.2] 

3. 0 ~ 

2.8 

2.6 

2.ij 

2.2 

2.0 

1.8 

1.6 

l.ij 

1.2 

-55-

T=O.O 

- - - = ~ - ~ 

1.,1,.11",11",111,111""111",,111 .. ,. I 

0.0 0.2 O.ij 0.6 0.8 1.0 

X 

Figure 7.1 
-56-

T=O.05 

P 

3. ~~ 3.2 

3.0 

2.8 

2.6 

2.l! 

2.2 

2.0 

1.8 

1.6 

lot! 

1.2 

1.0 -. 

- 0.8 

0.6 

O.L! 

~: ~] ~J~i~i'~"~"~i~iJ~i'~i~"~"~"~J~"~"~"~i~"~Ji~"~i~"~"~iJ~i~"~"~"~i~iJ 

0.0 0.2 D.L! 0.6 0.8 1.0 

x 

Figure 7.2 
P 

3.4: 

3.2 

3.0 

2.8 

2.6 

2.4: 

2.2 

2.0 

1.8 

1.6 

1.l! 

1.2 

-57-

T==O.l 

1. 0-J?---==---=--9---=, 

0.8 

0.6 

O.l! 

0.2 

El 

O.O~~~~~~~~~~~~~~ 

0.0 0.2 0.4: 0.6 0.8 1.0 

X 

Figure 7.3 
P 

3.l! 

3.2 

3.0 

2.8 

2.6 

2.l! 

2.2 

2.0 

1.8 

1.6 

loll: 

1.2 

-58-

T==0.2 

:: ~-j3----EC=---eO!--iOa--t:C:---a3~O:;-'-::l~ 

0.6 

O.l! 

0.2 

O.O~'I~'~"~"~'~"~'I~'~"~"~"~'~'I~"~'~"~"~'~'I~"~'~"~"~'~'I~"~"~'~"~"~I 

0.0 0.2 O.ll: 0.6 0.8 1.0 

x 

Figure 7.4 
P 

3.4: 

3.2 

3.0 

2.8 

2.6 

2.4: 

2.2 

2.0 

1.8 

1.6 

1.l! 

1.2 

1.0 

0.8 

0.6 

O.l! 

0.2 

0.0 

0.0 0.2 

-59-

T=1.5 

0.4: 0.6 0.8 1.0 

X 

Figure 7.5 
RHO 

3.4: 

3.2 

3.0 

2.8 

2.6 

2.L! 

2.2 

2.0 

1.8 

1.6 

1.2 

1. D-G----f;3---6:J 

0.8 

0.5 

D.L! 

0.2 

-60-

O.O~~~~~~~~~~~~~~~ 

0.0 0.2 0.4: 0.5 0.8 1.0 

x 

Figure 7.6 
-61-

T=O.Ob 

RHCl 

3.L!-

3.2 

3.0 

2.8 

2.6 

2.4: 

2.2 

2.0 

1.8 

1.6 

1.4: 

1.2 

1.0 ~ ~ 

~ ~ 

0.8 

0.6 

0.0 0.2 0.4: 0.6 0.8 1.0 

x 

Figure 7.7 
-62-

T== 0.1 

RHCl 

3.4 

3.2 

3.0 

2.8 

2.6 

2.4 

2.2 

2.0 

1.8 

1.6 

1.4 

1.2 

1.0 - - - - - 0.8 

0.6 

O.L! 

0.2 

0.0 . I 

0.0 0.2 D.L! 0.6 0.8 1.0 

X 

Figure 7.8 
-63-

T=0.2 

RHO 

3.4 

3.2 

3.0 

2.8 

2.6 

2.4 

2.2 

2.0 

1.8 

1.6 

1..4 

1.2 

1.0 

0.8 

0.6 

O.l! 

0.2" 

O. O~ .'~I~ii~I~"~"~I"il~ii~I""~iI~I"''~"~I"il~"~I"li~ii~I~"~I~li~i:~I""~ii~I~IiJ 

0.0 0.2 0.4 0.6 0.8 1.0 

x 

Figure 7.9 
-64-

T= 1.5 

RHCI 

3. l!1 

3.2 

3.0 

2.8 

2.6 

2.l! 

2.2 

2.0 

1.8 

1.6 

1.l! 

1.2 

1.0 ..... - - =::J - - -

0.8 

0.6 

O.l! 

0.2 

0.0 I . I' .. 

0.0 0.2 D.l! 0.6 0.6 1.0 

X 

Figure 7.10 
-65-

T==O.O 

U 

3.~j 

3.2 

3.0 

2.8 

2.6 

2.Li 

2.2 

2.0 

1.8 

1.6 

loLi 

1.2 

1.0 

0.8 

0.6 

O.I.! 

0.2 

0.0 

0.0 0.2 O.Li 0.6 0.8 1.0 - x 

Figure 7.11 
U 

3.l! 

3.2 

3.0 

2.8 

2.6 

2.l! 

2.2 

2.0 

1.8 

1.6 

lot! 

1.2 

1.0 

0.8 

0.6 

O.l! 

0.2 

0.0 

0.0 0.2 

-66-

T=O.05 

O.l! 0.6 0.8 1.0 

X 

Figure 7.12 
U 

3.ll: 

-67-

T=O.l 

3.2 

3.0~~~~~~~ 

2.8 

2.6 

2.ll: 

2.2 

2.0 

1.8 

1.5 

1..ll: 

1.2 

1.0 

0.8 

0.6 

O.ll: 

0.2 

o. 0 ~M'"I"'T"T"T"'I'"I''TT'T""......,.-r'''!"'T'T'T'''''''''''r''T'''P''I''"'M''':''''I'''''''''''''~I'''T'''!''T"~...........-r0.0 0.2 0.4: 0.6 0.8 1.0 

X 

Figure 7.13 
-68-

T==0.2 

u 

3. ~j 

3.2 

3.0 -- -

2.8 

2.6 

2.4 

2.2 

2.0 

1.8 

1.6 

1.4 

1.2 

1.0 

O.B 

0.6 

O.l! 

0.2 

0.0 

0.0 0.2 0.4 0.6 0.8 1.0 

X 

Figure 7.14 
-69-

T=1.5 

U 

3.L! 

3.2 

3.0~~~e-~-9~~~~~~~~ 

2.8 

2.6 

2.L! 

2.2 . 

2.0 

1.8 

1.6 

l.L! 

1. .2 

1.0 

0.8 

0.6 

O.L! 

0.2 

O.O~~~~~~~~~~~~~~ 

0.0 0.2 O.L! 0.6 0.8 1.0 

X 

Figure 7.15 
-70-

Simulation 2: We use a low Reynolds number: 

Re = 500. 

Once more we run the program with Riemann initial data: 

1.0 

1.0 

and 

1.625000 

0.3798263 

2.000000. 

At the left boundary we have subsonic inflow boundary conditions and at 

the right boundary the boundary conditions correspond to subsonic outflow. 

The solution to the Euler equations with the above Riemann initial data 

is a 1 shock moving to the left with a speed = 0.6124516. 

The numerical results, once again, have all the desirable properties we 

observed in Simulation 1. 
RHCl 

2.2 

2.0 

1.8 

1.6 

1.4: 

1.2 

1.0 

0.8 

0.6 

0.4: 

-71-

T=O.O 

--

0.0 0.2 0.4: 0.6 0.8 1.0 

X 

Figure 7.16 
RH(j 

2.2 

2.0 

1.8 

1.6 

1.4 

1.2 

0.8 

0.6 

-72-

T=0.2 

0'lJ:]' 

0.2 . 

~i~J~"~J~J~"MJ~jriJMj~j~"MJTJ~"MITJ~jjMjTj~"~jTj~lj~jTj~"~jTj~"~ITj~"~iTi~"~jTrjl 

0.0 0.2 0.4 0.6 0.8 1.0 

x 

Figure 7.17 
RHt:! 

2.2 

2.0 

1.8 

1.6 

1.2 

0.8 

0.6 

o.~ 

-73-

T=0.4 

0.0 0.2 O.~ 0.6 Q.8 1.0 

X 

Figure 7.18 
RHCl 

2.2 

2.0 

loB 

O.B 

O.S 

-74-

T=0.8 

0.0 0.2 O.~ O.S O.B 1.0 

X 

Figure 7.19 
-75-

T=1.5 

RHD 

2.2 

2.0 

1.8 

1.6 -

- - - - :3 - - - - -

1.14 

1.2 

1.0 

0.8 

0.6 

0.14 

0.2 . I 

0.0 0.2 0.14 0.6 0.8 1.0 

X 

Figure 7.20 
U 

2.2 

2.0 

1.8 

1.6 

1.2 

0.8 

0.6 

-76-

T=O.O 

0.0 0.2 o.ij 0.6 0.8 1.0 

X 

Figure 7.21 
U 

2.2 

-77-

T=0.2 

= El 

0.0 0.2 O.ij 0.6 0.8 1.0 

X 

Figure 7.22 
-78-

T=O.4 

U 

2.2 

2.0 

1.8 

1.6 

U] 

1.2 

1.0 --

0.8 

0.6 

:: :t " ." ",'. , ." .I" ". ,. I' .: •• ,,6, . , :" .. ~ 0.0 0.2 O.L! 0.6 0.8 1.0 

x 

Figure 7.23 
U 

2.2 

2.0 

1.8 

1.6 

1.2 

1.0 

0.8 

0.6 

0.4: 

-79-

T=0.8 

= - ~ 

0.0 0.2 0.4: 0.6 0.8 1.0 

X 

Figure 7.24 
U 

2.2 

2.0 

1.8 

1.6 

1.2 

1.0 

0.8 

0.6 

-80-

T=1.5 

o. ~t -- --- -= El 

o. 2 ~ ,:, , , ,:" ,,:, , . , i ' , , ,~ , , ,~, , ,,: "'~'" ,~, , , , i 

0.0 0.2 O.ij 0.6 0.8 1.0 

x 

Figure 7.25 
P 

2.2 

2.0 

1.8 

1.6 

1.2 

0.8 

0.6 

O.l! 

-81-

T=O.O 

0.0 0.2 O.l! 0.6 0.8 1.0 

X 

Figure 7.26 
F 

2.2 

2.0 

1.8 

1.6 

loY: 

1.2 

1.0 

0.8 

0.5 

O.Y: 

0.2 

0.0 0.2 

-82-

T=0.2 

O.ij 0.6 0.8 1.0 

X 

Figure 7.27 
P 

2.2 

2.0 

1.8 

1.6 

1.2 

0.8 

0.6 

O.L! 

-83-

T=O.4 

0.0 0.2 O.L! 0.6 0.8 1.0 

X 

Figure 7.28 
-84-

T=0.8 

P 

2.2 

2.0 

1.8 

1.6 

1.l! 

1.2 

1.0 

0.8 

0.6 

o.~ 

0.2 

0.0 0.2 o.l! 0.6 0.8 1.0 

x 

Figure 7.29 
P 

2.2 

2,0 

1.2 

1.0 

0.8 

0.6 

-85-

T=1.5 

:: :J .... ,,.." ,~, ,,..,-,.,....,.,..,.." '-, ...... ,.,....,., ·..,..,",..,.'.,....,.".,.....I'~· ..... ",..,.,,-.. ',..,.'..... 1'.,...,"..,..,',,..,., ..... , ,,....,., ,""'" ,,..,-,.,....,.".,.....,.~ 

0.0 0.2 O.L! 0.6 0.8 1.0 

x 

Figure 7.30 
-86-

REFERENCES 

[1] S. Abarbanel and D. Gottlieb, Optimal time splitting for two- and threedimensional Navier-Stokes equations with mixed derivatives, J. Compo 

Phys., Vol. 41, No.1 (May 1981), pp. 1-33. 

[2] J. B. Goodman, Initial boundary value problems for hyperbolic systems of 

conservation laws, Center for Pure and Applied Mathematics, University 

of California, Berkeley (September 1982), PAM-99. 

[3] S. K. Godunov, A finite difference method for the numerical computation 

of discontinuous solutions of the equations of flow dynamics, Mat. Sb., 

Vol. 47 (1959), pp. 271-290. 

[4] B. Gustafsson and A. Sundstrom, Incompletely parabolic problems in fluid 

dynamics, SIAM J. Appl. Math., Vol. 35 (1978), pp. 343-357. 

[5] A. Harten, On the symmetric form of systems of conservation laws with 

entropy, ICASE Report No. 81-34, (October 1981). 

[6] A. Harten and P. D. Lax, A random choice finite-difference scheme for 

hyperbolic conservation laws, SIAM J. Numer. Anal., Vol. 18 (1981), pp. 

289-325. 

[7] P. D. Lax, Hyperbolic systems of conservation laws and the mathematical 

theory of shock waves, SIAM Regional Conference Series in Applied Math., 

No. 11 (1978), SIAM, Philadelphia. 
-87-

[8] D. Michelson, Initial-boundary value problems for incomplete singular 

perturbations of hyperbolic systems, Proceedings of AMS-SIAM Summer 

Seminar, B. Enquist, S. Osher, R. C. J. Somerville, eds., La Jolla, 

1983. 

[9] J. Oliger and A. Sundstrom, Theoretical and practical aspects of some 

initial boundary value problems in fluid dynamics, SIAM J. Appl. Math., 

Vol. 35 (1978), pp. 419-446. 

[10] S. Osher and S. Chakravarthy, High resolution schemes and the entropy 

condition, SIAM J. Numer. Anal., Vol. 21 (1984), pp. 955-984. 

[11] J. C. Strikwerda, Initial boundary value problems for incompletely 

parabolic systems, Comm. Pure Appl. Math., XXX (1977), pp. 797-822. 

[12] F. M .White, Viscous Fluid Flow, McGraw-Hill, 1974. 
1. Report No. NASA CR-177984 I 2. Government Accssion No. 3. Recipient's C.Ulog No. 

ICASE Report No. 85-37 

4. TItle and Subtotle 5. Repon o.t. 

STABLE BOUNDARY CONDITIONS AND DIFFERENCE AUl!ust 1985 

SCHEMES FOR NAVIER-STOKES EQUATIONS 6. Performing OrganiZition Codl 

7. Authorh) 8. Performing Organiz.tion Repon No. 

Pravir Dutt 85-37 

10. Work Unit No. 

9. Pe~orminf. Orl/ilniZa¥on Nae;e and Addresa AI' . nst tute or omputer pp ~cat~ons in Science 

and Engineering 11. Contr.ct or Grant No. 

Mail Stop 132C, NASA Langley Research Center NASl-17070 

Hampton, VA 23665 13. Type of Report and Period Covered 

12. Sponsoring A~ncy Name and Address Contractor Report 

National Aeronautics and Space Administration 

Washington, D.C. 20546 14. Sponsoring Agency Code 

505-31-83-01 

15. Sur,ementary Notes 

angley Technical Monitor: Submitted to SIAM J. Numer. Anal. 

J. C. South Jr. Additional Support: NSF Grant No. 82-00733 and 

Final Report ARO Grant No. DAAG 29-82-0090 and NASA Grant 

NAr-l-<;Oh 

16. Abstract 

The Navier-Stokes equations can be viewed as an incompletely elliptic 

perturbation of the Euler equations. By using the entropy function for the 

Euler equations as a measure of , energy , for the Navier-Stokes equations, we 

are able to obtain nonlinear , energy , estimates for the mixed initial boundary 

value proble~. These estimates are used to derive boundary conditions \vhich 

guarantee L boundedness even when the Reynolds number tends to infinity. 

Finally, we propose a new difference scheme for modelling the Navier-Stokes 

equations in mu1tidimensions for which we are able to obtain discrete energy 

estimates exactly analogous to those we obtained for the differential 

equation. 

17. Key Words (Suggested by Author(s)) 18. Distribution Statement 

energy estimates 64 - Numerical Analysis 

gradient of entropy variables 

stable boundary condit ions 

stable difference schemes Unclassified - Unlimited 

19. Security aassif. (of this report) 20. Security Classif. (of thIS Plgel 21. No. of Pages 22. Price 

Unclassified Unclassified 89 A05 

For sale by the National Technicallnformatron Service. Springfield, Virginia 22161 

NASA·Langley, 1985 
End of Document