Iterative numerical solutions and boundary conditions for the parabolized Navier-Stokes equations

Computers& FluidsVol. 13, No. 4, pp. 397.409, 1985 Printedin the U.S.A.

(1045-7930/85 $3.00+ .00 © 1985PergamonPressLtd.

ITERATIVE NUMERICAL SOLUTIONS AND BOUNDARY CONDITIONS FOR THE PARABOLIZED NAVIER-STOKES EQUATIONS M. ISRAELI and A. LIN Computer Science Department, Technicon-lsrael Institute of Technology, Haifa, Israel 32000 (Received 6 July 1984) Abstraet--A new numerical iteration scheme for solving the parabolized Navier-Stokes (PNS) equations is presented. This scheme has all the features and advantages of the successive line over relaxation (SLOR) technique, and thus it can be easily accelerated to get much higher rate of convergence of the global iteration scheme than previously suggested schemes. The choice of appropriate downstream boundary conditions for the PNS and Navier-Stokes equations is discussed in the context of boundary layer simulation. A critical comparison of accuracy and rate of convergence is performed for the flow over a flat plate.

1. I N T R O D U C T I O N

Numerical solutions for the Navier-Stokes (NS) equations, especially for three-dimensional high Reynolds number flows, require a considerable amount of computer time and storage for adequate resolution. Since multidimensional flows are usually solved iteratively, the velocity and pressure fields have to be stored in addition to other (velocity) fields which are required by the iterative procedure. In some cases slow convergence rates were encountered, resulting in very long computations. One approach to cut down the computer requirements is to use approximate equations. In the present paper, we shall consider the parabolized Navier-Stokes (PNS) equations. These equations are derived from the Navier-Stokes equations by neglecting the diffusion effects in the (axial) flow direction (in agreement with the boundary layer approximation theory). Procedures for solving these equations use marching techniques (also through singular points of the boundary layer equations) along the main direction and result in considerable savings both in storage and computational time, if the solution after one marching sweep is acceptable. This would be the case if a well-posed initial-boundary value problem could be formulated. It was shown, however [ 1, 2], that the PNS system has an elliptic nature and therefore the initial value problem is not well posed. However, a well-posed problem for the PNS equations can be formulated by specifying, for the two-dimensional case, two upstream conditions (for the velocities, say) and one downstream boundary condition (for the pressure, say). Difference approximations to the PNS boundary value problem must reflect its elliptic-parabolic nature; the resulting coupled algebraic system must be solved globally (usually by iterative techniques). The multiple marching sweep as a global iteration method has the advantage that the velocity fields are generated during the marching process and only the pressure field has to be stored from sweep to sweep. In cases where the PNS equations are known to be suitable approximation to the NS equations, this approach is attractive. An essential ingredient for this to be true is a proper choice of the downstream boundary conditions. In some cases the I N S equations may be inaccurate, e.g. in regions of nonuniformity, while they may be a fair approximation elsewhere. In such cases the PNS solution can still be used to advantage by incorporation into the numerical scheme for the NS equations on a coarse grid; this is known as the "Booster" method [3]. Numerical experience shows that the rate of convergence of the usual multiple marching sweep deteriorates quickly as the marching stepsize Ax decreases and may eventually diverge. In the present paper we analyze the classical iterative marching procedure by considering its effect on the underlying elliptic equation for the pressure. It turns out that this iterative scheme (for the pressure) is not efficient as compared to other 397

398

M. ISRAELIand A, LIN

known iterative methods for elliptic equations. This analysis is then used as a tool for constructing an effective iterative method for the original system of equations. In Section 2 we present the governing differential and difference equations for the simple Cartesian case, and the basic marching procedure. The global iterative method is described in Section 3. Section 4 describes the analysis of the global iterative methods and the new iterative scheme. Sections 5 and 6 describe the boundary conditions that are reasonable to be imposed on the PNS system, and some results of the present technique are discussed. 2. THE GOVERNING EQUATIONS AND MARCHING PROCEDURE We consider two-dimensional steady incompressible laminar viscous flow along a fiat plate. The governing PNS equations are [1 ]: Ux + vy = 0

(1.1)

L(u) = -Px

(1.2)

L(v) = -py

(1.3)

where the following conservative form for the operator L is used:

L(f)

=

(u f ) +

(v f )

(1.4)

R e O f 2"

Here Re is the Reynolds number based on one of the domaln's dimensions (for example, the plate length). Numerical solutions of eqns (1.1)-(1.3) are obtained by spreading a grid over the computational domain. Let us assume that the grid points are distributed evenly along the x and y coordinates, with the spacing Ax and Ay, respectively. When differencing these equations it should be remembered that the elliptic feature should be reflected [ 1, 4] in the finite difference approximation. In order to be consistent with the boundary layer (parabolic) nature of the flow, the axial gradients of the velocities should be computed using only upstream values, while the elliptic nature is preserved by forward differencing the axial pressure gradient [ 1, 4, 5]. Using the conservative form (1.4), eqns (1.1)-(1.3) are differenced at the point (i, j) as follows: AX

( v i j - v~,j-l ) = 0

ui, j + u;,j-i - ui-l,j - ui-l,j-I + 2 ~ 1 Ax

ui,1 - u,-,, 2 + i - ~ y

[(uv),,j÷,

-

(uv),j_, ]

AX

R,ay z

(Uij+~

,

- 2ui, j + ui,j-i) = Pid - Pi+~.j Ax

(UI~)i,j "3L (Ul))i,j+l -- (Ul))i-l. j -- (Ul))i-I,j- 1 "~- 2 ~

2

"l- - R - ~

(Ui,j+I -

Ui, j -

(1.5)

(1.6)

[(I)i,j+ 1)2 -- (1)/,j)2]

Ui-l,j+l "k U i - l , j - l )

-~y (p~,j- p~.j+~). (1.7)

= 2 Ax

Here a forward difference is used for the Px term in the u momentum equation (1.2) [1], where pj+~.j is treated as k n o w n . This scheme results in an unconditionally stable ("departure-free") marching. During the marching the nonlinear algebraic system is solved at every station i by the full Newton method [ 1]. Obviously, the results after one marching sweep are strongly dependent on the initial pressure field. The scheme (1.5)-0.7) is of the first order in x. It was previously thought that a stable marching scheme must be of the first order in the marching direction because it requires backward differencing for the velocities to use the upstream information, but forward differencing for the streamwise pressure gradient. It turns out that this effect can be

Iterative numericalsolutionsand boundaryconditions

399

achieved by a judicious choice of the placement of the variables to be solved at each 1 station. The choice can be explained most easily by taking v = 0 and ~- = 0 in eqn (1.4) for u, yielding Ux2

= --Px.

(1.8)

The first-order difference scheme (1.6) then becomes (1.9)

Ui, j 2 - - U i - l , j 2 = P i , j - - P i + l , j .

The unknowns are u~,j and p~,j. The first suggested second-order method [6, 7] is Ui+l, f l -

u~,f l = P i j - Pi+l,j

(1.1o)

with the unknowns U;+l,j and P,d. The scheme is centered about i + ½ and is second order. Since the purpose of this work is to analyze the global marching procedure, we will use the classical first-order scheme defined by eqns (1.5)-(1.7). 3. GLOBAL ITERATIVE METHODS Numerical experiments with the present computer code show that the solution after one marching sweep is not close to the final solution of the PNS equations, when the initial pressure field is constructed using the boundary layer assumption py = 0. Since the Px term is forward differenced, some global iteration over the entire solution domain should be performed in order to converge the explicit contribution to this pressure term. The simplest global iterative technique to solve the equations is by multiple marching sweeps with the primitive equations (1.5)-(1.7), where only the pressure field is kept from iteration to iteration [l]. Numerical experiments also show that for certain nets this procedure diverges. The divergence occurs also for the linearized version of eqns (1.5)(1.7). Figure 1 presents the residium of the pressure field as a function of the global iteration's sweep number for a 21 × 11 field. A jump is encountered every l0 iterations (probably related to the arrival of the boundary pressure pulse, traveling at the numerical scheme speed) leading to ultimate divergence. However, convergence was reported with a different mesh and boundary conditions, and also when combining the above procedure with a multigrid technique [8]. It was thought that the replacement of one of the momentum equations by the Poisson equation for the pressure will improve the convergence properties. Indeed, the two resulting schemes converged at a satisfactory rate, but the solution did not satisfy the replaced momentum equation. A successful implementation of the marching technique is derived in the next section. A short and reduced version of the analysis was presented [5]. 4. THE ITERATIVE SCHEME FOR THE PRIMITIVE SYSTEM To motivate our approach, we start with the PNS equations (1.1)--(1.3) and linearize them about a constant state. The linearized system remains (l.1)-(l.3), but eqn (1.4) is changed to

C{f)'~x(Of)+~(~'f)y i,

Oy 2

(4.1)

where 0 and l? are constant reference velocities. The next step is to discretize the equations only in the x direction to obtain: Um -

Urn_ 1

Ax

+ (vy). = 0 em÷l

D(Um) = D(Vm) - - ( e y ) .

B

Ax

P~

(4.2)

(4.3) (4.4)

400

M. ISRAELIand A. LIN 2.0-

15.

a,_ 1,0<3

-

0.5-

OD-

o

-16

z6

sb

X/I..

Fig. 1. Pressure residue vs the global iteration number. where Urn, I'm and Pm are functions o f y. Here D(fm) is the semi-d_iscretized form o f / : ( f ) at the marching station m. T h e semidiscretized system should be diseretized also in the y direction before solution is attempted, but since the speeitic form of this disc~etization is not important for the following argument, we postpone this step for the sake of transparency. The marching iterative procedure assumes that u,k(y), F,k(y) are known as well as Pmk-' for m = 2, 3, 4, • • • M, where k is the current iteration index. Therefore, the marching scheme for m > 2 is: Vmk -

Vm_ k

+ (V~,),,, = 0

hx

D(Um k)

=

- -

em+l k-I Ax

- -

em k

D(Vmb -(P,%. =

(4.5) (4.6) (4.7)

We now apply D to eqn (4.5) and differentiate eqn (4.7) with respect to y. Elimination o f the V terms between eqns (4.5) and (4.7) gives:

D(UI k - Urn-?) = a x ( e . b , . .

(4.8)

By substitution of eqn (4.6) into eqn (4.8) we obtain: Pm+l k-I - P m k - P m k-l -~ Pm-i k "~- A x 2 ( p ~ k ) m = O.

(4.9)

It follows that the marching scheme for the primitive system (4.2)-(4.4) can be viewed

Iterafive numerical solutions and boundary conditions

401

as a line iterative scheme for the semidiscretized Laplace equation; indeed, upon convergence eqn (4.9) will become: P,,,+I - 2Pro + P~-l + (Pry),, = O. Ax 2

(4.10)

In order to find out the rate of convergence of eqn (4.9) to the final state (4.10) we Fourier transform eqn (4.9) in y assuming appropriate boundary conditions in that direction: e m k = Zme lny,

]amk-I = Zme lny

(4.1 1)

where 12 = - 1 and n is the Fourier wave number. After substituting these definitions into eqn (4.9) we obtain: z..+, -

-

z.

+ 2.,_, -

ax2n

2.. = o.

(4.12)

Transforming in the x direction we define again (4.13)

Z = A e t~'x

where 0 is the wave number in the x direction. By substitution of this definition into eqn (4.12) we obtain:

_1: e_5 I [=[l+ x2n 2-

e-'° I

(4.14)

This means that all the long waves (with small Ax2n2) in the cross-flow direction are only weakly damped irrespective of their structure in the marching direction. In particular, the n = 0 modes which exist for derivative boundary conditions in the cross-flow direction are not affected at all by the relaxation, i.e.

On the other hand, the well-known SLOR scheme for the Poisson equation gives (after the same Fourier transformations): Z,,,+, - (2 + Ax2n2)Z,,, + Z,,,_, = 0

(4.15)

and [q _ e_lo [

q = 2 + Ax2n 2 ~ 2

(4.16)

and also = q2 + 1 - 2q cos 0 "

(4.17)

This quantity is less then 1 for all acceptable q's and cos 0 < 1. Most waves are strongly damped and only the longest waves in both directions are weakly damped by the iteration. This behaviour was used to accelerate the convergence as is done by the SLOR technique, Chebyshev acceleration or multi-grid methods. The question is how to generate an equivalent relaxation scheme for the primitive system in the marching form. This means that we may add terms which can be evaluated during the marching process but should vanish upon convergence. A rational approach to the construction of the relaxation scheme is to retrace

402

M. ISRAELI a n d A. LIN

backwards the steps of the derivation of the discrete Laplace equation from the discrete primitive equations. We start from the SLOR equation (4.15): Z,,,+l - 2Zr,, + Z,,,-t = Ax2n2Zra

which we inverse Fourier transform with respect to y to obtain: Prn+l k - I _ 2pro k + prn_l k =

_Ax2(pyyk),..

Now we substitute eqn (4.8) for the right-hand side of the last equation to obtain 2pro k + P m _ l k = - A x D ( U m

Pm+l k-I -

k -- Um_lk),

which can be written as _ ( p r n + l k - , _ pmlC) ..~ (prnk-I _ pink) ..~ ( p m k - I _ prn_l k) --_ A x D ( V r n _ U r n - l ) .

Adding the equations from m = 2 and using the linear form of D we obtain m

_(prn+lk-~ _ prnk) + E (Pi k-~ - p k) + ( p ? - ,

_ p k) = A x D ( U m - UI)

iffi2

but from eqn (4.3): p2k-1 __ p1 k = - A x D ( Ulk).

(4.18)

Therefore, we obtain for m > 2 D(Um) =

Pm+~k-l -- Pink + -

Ax

1 ~ ( p k - I _ pik).

(4.19)

~ x ~-2

Equation (4.19) contains all the modifications required in order to convert the iteration scheme of eqns (4.5)-(4.7) into a scheme equivalent to the SLOR scheme for one Laplace equation with fl = I. We see that in this approach only the x momentum equation is modified. The new added term can be generated easily during the marching process and is inexpensive in storage (one extra line vector) and computation (one substruction per grid point). In what follows we will rederive eqn (4.19) in a more general way and introduce the over-relaxation parameter fl > 1. In practice we will use difference approximations and boundary conditions also in the y direction and the resulting scheme may not be amenable to the discrete analogue of the Fourier transform. It is therefore worthwhile to generalize the previous approach by using the matrix finite difference formulation. Let the vectors Urn, Vm and Pm contain the N values of the corresponding variables on the mth line (x = constant) of the marching sweep (including the specified boundary values). The u-momentum equation (1.5) can be written in the form: em+~ -

Prn = - - A x D ( U r n ) .

(4.20)

On the other hand, elimination of I'm between the continuity and the v-momentum equations will result in: FPm = Rm - R,,,- t

(4.21)

where F = Ax21(O2/Oy2). Substructing successively the u-momentum equation (4.20) and using eqn (4.21) gives: Prn+t - (21 + F)Pm + Prn-t = 0,

m = 3, 4, • • •

(4.22)

Iterative numerical solutions and boundary conditions

403

which is Laplace's equation. The first equation o f (4.20) can be used as a derivative condition at the left (inlet) boundary, namely: /'3 - / ' 2 = R2.

(4.23)

We now apply the SLOR scheme to the last two equations (ignoring temporarily the downstream boundary condition) to get the downstream marching form: -P~

+ P3 k-~ = R2

Pr,,-t k - (21 + F ) P * + Pm+l k-I = 0;

(4.24) m = 3, 4, - • •

(4.25)

where Pmk = t i P * + (1 - f t ) P m k - l , ft is the over-relaxation factor and the superscript denotes the iteration sweep number. In order to recover the primitive variable formulation, we relate the velocity field in eqn (4.2 l) to the starred pressure field, i.e. F P * = Rm - R,,_~.

(4.26)

Substitution in eqn (4.25) gives: Pr,,-i k -- 2 P * + e m + l k - I ~-- R m - R , n - l ,

m = 3, 4, . . •

(4.27)

Successive summations o f eqns (4.24) and (4.25) gives: P,n+l k-l -- P * = Rm + S,n;

m = 2, 3, 4, • • •

(4.28)

which is the primitive variable marching form o f the u-momentum equation. The source term Sm in eqn (4.28) satisfies: Srn = Srn-I + ( P * - P,n k - l ) + ( P * - t - P,n-lk),

m = 3, 4, • • •

(4.29)

with $2 = 0. It can be seen that S,n vanishes upon convergence. The computational form o f eqn (4.28) for m = 3, 4, • • • is: -2P*

= R,,, +

(4.30)

Sm = S,n-1 -- P,n+l k-1 + 2P'm-1 - P,,,-tk;

Thus, the theory o f over-relaxation can be applied exactly to the constant coefficient case o f system (1.5)-(1.7). For the nonlinear ease this theory can serve as a guide to the choice o f ft. Alternatively, one can choose fl = 1 and apply the multigrid procedure [7, 9]. 5. BOUNDARY CONDITIONS The problems considered here simulate boundary layer like flows, where the inflow conditions (at x / L = 0) are assumed to be known. The edge conditions should be determined by viscons-inviscid interaction with the potential flow. However, we assume here that the edge conditions can be determined from local relations and concentrate on the end (downstream) conditions. These conditions should allow the flow far downstream to reattach and behave like a boundary layer, and to have a small upstream influence. Such conditions should apply to the PNS and NS systems. We start with the I N S system, where only a downstream boundary condition for the pressure is needed. Since py is small in boundary layers, a possible condition would be: p ( L , y) = p ( L , oo)

(5.1)

which is equivalent to the omission o f the normal pressure gradient in the normal

404

M. ISRAELI and A. LIN

m o m e n t u m equation (1.7). But this is not a trivial assumption since the p~, term is of the same order o f magnitude as the other terms in eqn (1.7). Thus this boundaD' condition should be implemented very carefully. In the last marching station rn where the end conditions are applied, we force Pm,j = Pro,n, 1 -< j < n. But when solving the other flow variables on this station, we get the velocities from the continuity and the normal m o m e n t u m equations, where the axial m o m e n t u m equation serves to solve the axial pressure gradient at the last station (which is really not needed). It can be shown for the boundary condition U = V = 0 at y = 0 the following holds: O~V

0y"

-0

at

y=0,

v= 1,2,3, -..

Therefore at the last station we have: P= P~,

V = O,

U = u(y)

which is a very "tough" boundary condition. A more suitable boundary condition can be derived by writing the px term in eqn (1.6) in the form: p x ( L , y) = px(L, oo) -

p~ty.

Since the integral term is small compared to other terms in eqn (1.6), we obtain for x = L the end condition: L ( u ) = - p x ( X = L , y = or).

(5.2)

Here the velocities are solved from the continuity and the axial m o m e n t u m equations for the given Px, and the pressure is defined by the normal m o m e n t u m equation. In many cases, especially for compressible flows over complex configurations, the downstream conditions are known from some special measurements, like those of the streamlines. Here we will consider a very c o m m o n boundary condition of this type; the directions o f the streamlines at the outlet are known, i.e.: V U" The solution procedure in this case is as follows: some values for the Px term are assumed at the outlet station, and by using the equation Px = v U ~ + U2ar

(5.3)

these values are corrected. It should be noted that some o f the equations have to be applied at the half axial interval at the downstream end station. The boundary conditions for the NS equations that will best approximate the present situation and result in a well-posed problem are not obvious. Motivated by boundary layer scaling, we choose eqn (5.1) or (5.2) as one of the downstream conditions. The choice of the second boundary condition is motivated by vorticity transport considerations; the vorticity is advected downstream by the main flow while the back diffusion is small. It follows that specifying an incorrect downstream vorticity will modify the flow only locally, in a thin boundary layer near the exit. In such a boundary layer the scale will be o f the order o f the inverse o f the Reynolds number, and the main contribution to the vorticity will come from the x derivatives. We expect that specification of v ~ = 0, u ~ = 0 at the exit will minimize the local generation of the vorticity. This amounts to use the PNS equations at x = L together with eqn (5.1) or (5.2).

Iterative numericalsolutionsand boundaryconditions

405

Returning to the initial conditions, it was first proposed that the u and v velocities will be given at the first station. The requirement on v is not physical and may unnecessarily affect the results. Some other initial conditions were also considered: 1. To impose the u velocity as well as the angle of the flow in the first station ~, defined as e = v/u. 2. To impose the u velocity at the inlet as well as the Bernoulli condition which is imposed iteratively. 3. To impose vx -- 0 at the inlet. The last possibility is illustrated in Figure 6. 6. RESULTS First we will establish the need for the pressure field iteration procedure. For a 41 X 41 grid, Fig. 2 shows that the values of the maximum v velocity after the first sweep are three to five times smaller than those of the converged solution. This figure depicts also the effect of the downstream boundary condition. The maximum v velocities for the two boundary conditions (5.1) and (5.2) agree initially; however, the difference grows to about 27% at the plate's end. The maximum difference in the wall shear is about 15%. The edge is assumed to be in the potential region where two conditions out of the following three can be chosen: (1) specified streamwise velocity; (2) static pressure and (3) stagnation pressure. It turns out that the difference between the various choices for the edge condition is small (since the magnitude of the square of the edge normal velocity is negligible). The location of the maximum normal velocity Vm~, can be considered to be the edge of the viscous region or the contact point between the lower and the middle deck in the

020.

015,

01(>

/

........... .;.....

---"

/-,\ 0.05.

\

OOO 0.00

025

Legend: Curve

0.50 X/L

075

tOO

OQD

BaSlUS

~-&~,- eqn (51)

o oo

eqn (5.Z

-

Fig. 2. Converged/)IX profiles.

~

One

Sweep

406

M. ISRAELIand A. LIN 0.8-

0.7"

0.6-

05-

o4J o.3=

0.2-

o.1

00 t

o.oo . .

. . . . . .

o.25' .

o.75'.

. . . . . . . . o ~. o . . . . . . .

. . . . . .

i6o

X/L Legend: Curve

D o O eqn(51)

-,~A~-eqn(5.2)

Fig. 3(a). ConvctBed Pm,~ profiles.

09

08

Q7

",

06-

',,

05",%

04-

0.3-

0.2.

1

0.1.

~

o,, B'"

0.O 000

0.02

0.04

Legend: Curve

006 ~

Q08 P eqn(51 )

Fig. 3(b). Convergedp profilesat x / L

0.10 -~-~ =

0.25.

012 eqn(5.2)

014

407

Itemtive numerical solutions and boundary conditions 1.0

T

""-,

•

,

-

// ...... ~...... :..,,.,

;

0.9

, \\

",,

!

" 0.8.-

,

I ,

, ,

i',./;" \/ 1.0

1.1

1.2

1.3

1.4

1.5

t.6

1.7

t8

1.9

Over reioxotion Porometer Legend Curve

--L 2 t e l l eqn{5.2) -~-NL 21sll eqn(5.2) - ¢ ~ : - N L / L 2t=3 eqn(5.2)

-¢3-~-~- L 21e3 i ! q NL 21e3

eqn(5.1 ) eqn(5.1)

Fig. 4. Rate of convergence for coarse and fine grids; R = 100, L = linear and N L = nonlinear.

0.13 -,

,~'

/

0.42-

/

o~.

/"

,,C7//'~Tj

o. ll

0.09

..............

0 0 5 -I O.0

0.1

0.2

Legend: Curve

0.3

0.4

0.5 X/L

= = eqn t ~ . ~]) ~ - e.~,~-e-

0.6

0.7

0.8

0.9

-~~ e , , - e qeqn n (5.2) ( 5 . 1 ) PPNS NS

Fig. 5. Comparison of u,=,~ for the PNS and NS equations.

1.O

408

M. ISRAEU and A. LIN

triple deck theory. Thus, h is reasonable to compare the Vm,x values of the PNS equations to those o f the Blasius solution. It can be seen from Fig. 2 that ~m~, of the Blasius solution changes more rapidly with x than that of the PNS equations, Figure 3 compares the pressure overshoot for the boundary conditions (5.1) and (5.2). There is a m a x i m u m in both pressure profiles due to the viscous interaction, and this m a x i m u m is closer to the wall than that of v. We note that the considerable difference in the pressure for the two boundary conditions is felt throughout the whole computational domain. Figure 4 shows the average rate of convergence as function of the over-relaxation parameter 9, for problem (1.1)-(1.3) and its constant coefficient (linear) version. The two sets of curves correspond to the end conditions (5.1) and (5.2). For the linearized version with Dirichlet end conditions the dependence on 9 agrees with the theory o f SLOR for Poisson equation. For the nonlinear case the curve is shifted to the right while keeping a similar shape. There is a strong dependence on Ax and a weak dependence on Ay (compare the curves for the 21 × 3 grid to those of the 21 × 11 grid); thus, low resolution runs can be used to estimate 9. The end condition (5.2) gives results similar to those of the previous case, but the rate of convergence is smaller. Figure 5 can serve as a justification for the use of the PNS equations. It can be seen that even for R = 100, the difference between the NS and the PNS solutions is quite small away from the end stations. The effects o f the mesh can be seen in Fig. 6. It can be seen that the numerical scheme for both the NS and the PNS systems is much more sensitive to the spacings in the normal direction than in the axial direction. This effect is very pronounced at 1he leading edge and diminishes as x --~ oo.

015

*-~"

. . . . .

~-.....

,~ '\

!/

0~2

'*"*~

'\,,

009o.

,/'~"" ,c

o8.

+ ,,

0.07. , ........ 0(3

, O~

....

, ......... 02

, ..... 0.3

,-

,

0.4

.....

0.5

,

-, . . . . . .

0.6

07

i-- 08

l 0.9

-, 1.0

X/L Legend:Curve

~ r~ = •e - , e - - e . -I,4-

NS : 11e2I NS:dte41 PNS,21e21

Vx.O

Ot

~-~~ * - ~

i

NS = 21121 PNS,21e21,eqn(5.2) PNS,41e41,eqn(5.2)

x'O

Vtg. 6. Variation of the maximum velocity for the NS and PNS systems with different meshes and boundary conditions.

Iterative numerical solutions and boundary conditions

409

Acknowledgements--This research was supported by Air Force contract no. 77-3405, and by the Stiftung Volkswagen Werke. Partially supported by the Technicon VPR Fund--Lawrence Deutch Research Fund, Grant No. 121-622 (1982).

REFERENCES 1. S. G. Rubin and A. Lin, Marching with the parabolized Navier Stokes equations. Isr. J. Tech. 18 (1980). 2. M. Israeli, V. Reitman, S. Salomon and M. Wolfshtein, On the marching solution of the elliptic equations in viscous fluid mechanics. Proc. 2nd Int. Conf. Numerical Methods in Laminar and Turbulent Flows (Edited by C. Taylor et al.) Venice (1981). 3. M. Ungarish and M. Israeli, Improvement of numerical schemes by incorporation of approximate solutions applied to rotating compressible flows, Proc. 7th Int. Conf. on Num. Methods in Fluid Dynamics, Stanford, Lecture Notes in Physics 141, Springer-Verlag (1980). 4. A. Lin and S. G. Rubin, Three dimensional supersonic flow over a cone at incidence. AIAA J. 20(11), 1500-1507 (1982). 5. M. Israeli and A. Lin, Numerical solution and boundary conditions for boundary layer like flows. Proc. 8th Int. Conf. on Num. Methods in Fluid Dynamics. Aachen, West Germany, Springer-Verlag, pp. 266-272 (1982). 6. M. Israeli, NASA Lewis seminar, July 1982. 7. M. Israeli and M. Rosenfeld, Marching multigrid solutions to the parabolized Navier-Stokes (and thin layer) equations. 5th GAMM Conf. Num. Fluid Dynamics, Rome (1983). 8. S. G. Rubin, Incompressible NS and PNS solution procedures and computational techniques. Von-Karman Institute Lecture Notes (1982). 9. S. G. Rubin and D. R. Reddy, Analysis of global pressure relaxation for flows with strong interaction and separation. Computers Fluids 11(4), 281-306 (1983).