Pressure Gradient method for solving incompressible Navier-Stokes equations with curvilinear coordinate system

Math1 Comput. Modelling, Vol. 14, pp. 132-136, Printed in Great Britain

Pressure

Gradient

Equations

Method for Solving Incompressible

with Curvilinear

Cha-Hsiang

08957177/90 $3.00 + 0.00 Pergamon Press plc

1990

Coordinate

Navier-Stokes

System

Tan and Michael Pecht

Department of Mechanical Engineering, University College Park, Maryland 20742, U. S. A.

of Maryland

.J. C. Duh Sverdrup Technology, Inc., M.S. 500-217, Cleveland, Ohio 44135, U. S. A.

Abstract.

A nonstaggered

using primitive

variables,

Pressure

Gradient

(PG) method, one of the fmite difference

equations

based on a curvilinear

coordinate

uses pressure

instead of pressure

gradient

contribution

of this method is its adoption

troublesome

calculation

cases.

of the pressure

near the boundary.

terms is relatively buoyancy

The comparisons

nonstaggered

grid which eliminates

is that the discretization

and 3-D unsteady

the

for

error of the convection

flow in the high Reynolds

have been studied, including

the present numerical

to conventional

The most prominent

and any special treatment

of the convecting

schemes

Navier-Stokes

The PG method, in contrast

conditions

advantage

incompressible

as the field variable.

of a modified

to the magnitude

convection, between

system.

boundary

Another

insensitive

A variety of test problems

transient

Lewis Center

is applied to solve the 3-D, unsteady

methods,

velocities

NASA

number

Couette flow, flow past a cylinder,

cavity flow driven by shear and body forces. results and the existing data are shown to be in

good agreement. Kevwords,

Pressure

Gradient

method;

Curvilinear

Coordinates;

Navier-Stokes;

Finite-

differencing.

satisfy the discretired

INTRODUCTION

unrealistic. In solving incompressible primitive

variables,

identified, pressure

Navier-Stokes

two major difficulties

with

programming

the

coordinate

is unavailable;

strategies

spurious

pressure

have led to the development

for iteratively

construction

of various

undesirable

spurious

area include Harlow

computing pressure

Previous

(1972), Maliskar

system.

works in this

equations

and Raithby

(1984),

based on a curvilinear

computing

coordinate

and the accuracy

PG method,

a variety of test problems

flow problems

have favored

gradient

system.

To

of the proposed

are studied, including

transient

cavity

(PG)

buoyancy

driven flow.

MATHEMATICAL

field.

FORMULATION

the by

When a curvilinear

It is well known that the use of a

pressure

and 3-D unsteady

Gradient

Navier-Stokes

It has been a trend

grid usually causes an oscillating

or even an oscillating

Pressure

into two types,

(MAC) grid, which was first developed

and Welch (1965).

nonstaggered

grid.

ago that most of the schemes for

incompressible

use of a staggered Harlow

can be categorized

in high Re number flows.

the feasibility

Couette flow, flow past a cylinder,

grid and nonstaggered

since two decades

(MAC) grid could cause

to be less accurate

demonstrate

convection, staggered

with a curvilinear

it has been pointed out by Tan

In the present paper, a nonstaggered

and Tan (1988). Grid arrangements

of problems

Moreover,

grid in

This merit is identified

method is applied to solve incompressible

the

and Welch (1965), Chorin (1967),

and Spalding

on computation

the computation

modes.

and the

for eliminating

modes.

is very attractive.

but are physically

of nonstaggered

of various

the pressure

grid networks

equations

the simplicity

and Duh (1989) that the staggered

the other is that a proper grid

is needed to eliminate

These difficulties

especially

have been

one is that an explicit equation governing

arrangement

Patankar

equations

governing

However,

pressure

dimensionless

field,

coordinate

governing

can be written as follows:

The results 732

system is introduced,

equations

of incompressible

the flow

733

Proc. 7th Int. Conf on Mathematical and Computer Modelling

f$(JUjg;O

(1)

Fl

, where ci is the boundary-fitted curvilinear coordinate system; J and giJ are the metric coefficients as defined in Thompson and Warsi (1985); Ui is the connavariant velocity component; PO is the pressure gradient term, and f,t,is the

(small), forcing too much (little) mass out from one control volume to the other. Subsequently, the pressure gradient for the next iteration should be reduced (increased). This iteration process is repeated until the mass flux is balanced between two neighboring control volumes. It is seen that Eq. (5) performs this strategy. To accelerate the rate of convergence, the velocity correction equations (Patankar and Spalding, 1972) can be applied directly with the use purturbated pressure gradients.

$-

B

I I

P

P

P

body force. It should be noted that variable $ denotes the covariant velocity component in the present paper. Solution Procedure The Navier-Stokes equations for incompressible flow are governed explicitly by the pressure gradients, not by the pressure itself. If the condition of compatibility for the pressure Vx(VxP)=O

(3)

and the continuity constraint (a) Q stands for u and v v*o=o

(4)

7 -8i-gi -8

are combined, the computation of pressure is no longer needed. Following this concept, we propose an iterative procedure to compute the pressure gradients instead of the pressure itself. This procedure is further facilitated by using a grid that, being staggered originally, is modified into a nonstaggered one shown in Figure 1. As a result, the proposed pressure gradient (FG) method enjoys the simplicity of the nonstaggered grid and, at the same time, avoids the contamination of the spurious pressure modes, the computation of the pressure boundary conditions, and the violation of mass balance. To introduce the pressure gradients into the continuity constraints, the modified form for the artificial compressibility method (Chorin, 1967) can be written as:

,

-81 -z-

-$ (b) $ stands for

(5)

where yi is an artificial time step in terms of metrics. During the calculation, Eq. (5) is iterated to obtain the correct pressure gradient fields and ensure that the continuity equation is divergence-free. This iterative procedure follows the proper physical trend in correcting the net mass flow crossing the interface between two neighboring control volumes. If the net mass flow is positive (negative), it implies the pressure gradient at the interface is too large

8-

u

and

V; 8

P

I

stands for

U, v, Wax

and

aPlaY

‘d is modified into a VI-I nonstaeeered one

FIG. 1. An onelnallvd

Test Problems

Case1: Couette The inner cylinder with radius rl rotates about its axis at a constant angular velocity w. The outer cylinder with radius r2 is stationary. The problem is schematically shown in Figure 2. The x and y coordinates and the velocity

Proc. 7th Int. Conf. on Mathematical and Computer Modeiling

734

FIG. 2. Schematic of Couette flow uroblem

components are nondimensionalized by using rt and rlw, respectively, whereas the pressure is nondimensionalized by using (rlw)*. The Re number is defined in this case as rt*w/v. The analytical solutions to this problem can be found in white (1974) as: u=(l-k2r2)y,and

v=

(1 - k5 r2

-(l-k2r2)x

3 the flow vast a cvlinder nroblem

(6)

(1 - k5 r2

, where r* = x2 + y* and k = rtjr2. Case 2: Flow Past a Cvlinder The nondimensional variables x, y, u, v and p for this problem are defined based on characteristic quantities D, u, and pum, respectively; Re is defined as puJ/p, where D is the diameter of the cylinder, and pm, u, are free stream pressure and velocity. The transformation domain is indicated schematically in Figure 3. The cylinder is surrounded by a circular domain of large but finite diameter @, = 12.5 D). This domain is covered with a non-uniform grid distribution shown in Figure 4. FIG. 4. An O-tvoe mid distribution used for

Case 3: Transient Buovancv Convection in a Sauare Cavity The continuity, Navier-Stokes, and energy equations are formulated using a conventional Boussinesq approximation. The initial conditions are u = v = 0 = 0 at t = 0 everywhere. The boundary conditions are

u+

0

u=v=c1,8=0~at

al

y=o, x=0,

u=v=

ae

&5;=0

u=v=o,e=e,at

Case 4 3-D Unsteadv Cavitv Driven Flow

aty=l,

(7) x=1.@)

the flow nast a cvlinder uroblem

This test problem is the flow driven by an artificial body force which has the following components in the Navier-Stokes equations: fu = 2 x2 y z + 4 x3 y* z2 t* - 2 y z t (2/Re + l),

(9)

f,=-xy2z+x2y3z2t2-2xzt(l-l/Re),

(10)

f, = - x y z* + x2 y* z3 t* - 2 x y t (1 - l/Re).

(11)

735

Proc. 7th Int. Co& on Mathematical and Computer Modelling

The selection of these artificial source terms provides an exact solution to the governing equations, the exact solutions are: u=2xyzt,v=-xyzt,w=-xyzt,

(12)

p,=-2yzt,p,=-2xzt,p,=-2xyt.

(13)

Literature

L sep

SOP

0.56625

52.2

1.9

___---

_____- -

53

2.1

I.0015

0.5275

53.7

2.515

I.05

Kawaguti

e

CD1

cDP 17

I

Taneda (experiment) Kawaguti

For computational stability, we select second-order upwinding scheme and finite volume approach to discretize the convection terms. The convergence criteria am: (maximum relative error of velocity components) < 0.01, and (maximum residual of continuity equation) < 10m5.To check the accuracy of the scheme, we define the mean error for

and

Jaln

I

Present

1

0 53

1.03

result

52

I .95

dependent variable $ as:

(14)

, where N denotes number of grids used. In case 1, Table 1 gives the standard deviations of u and v for different values of Reynolds number. It should be noted that the exact solution for this problem is independent of Re number. From the results presented in Table 1, the numerical results are seen to be in good agreement with the exact solutions. Table 1 The mean errors on Couette flow uroblem

ulot for the flow vast a cvlinder nroblem

Qble 3 Comparisons of Bansient buovancv-driven flow (Steadv-state solution) In case 2, the quantities that are presently used to verify the physical phenomena are the wake length, L,,,*e; the separation angle Gsep;the pressure drag, Cot,; the friction drag, and the total drag, CD Table 2 summarizes these results at Re=40 and compares with those reported in the literature (Kawaguti, 1953, Taneda, 1955, Kawaguti and Jam, 1966, Son and Hanratty, 1969, Jordan, 1970, and Kwak, etc., 1984). The results obtained from the PG method are generally in good agreement with those reported in the literature. The flow in the form of velocity vectors is shown in Figure 5. A standing eddy has formed at the rear of the cylinder. The velocity inside the eddy is very small.

Literature

Gr

u max

V

max

Vahl Davis

1,000

3.649

3.697

Present

1,000

3.483

3.505

Vahl Davis

10,000

16.178

19.617

Present

10,000

15.I77

18.45 1

result

result

In case 3, the computation is performed on a 21 by 21 uniform grids with Re=l, G1=103, 104, Px=O.7, GH=l and

Proc. 7th Int. Conf. on Mathematical and Computer Modelling

736

B&I. The distribution of the vertical velocity component v along y=OS axis for several time steps are plotted in Figure 6. It has been found that the flow recirculates clockwise due to the buoyancy effect and the flow reaches the steady state around t=l . 1. The values of u maxi vmax are given in Table 3. These values compare very well with the results of de Vahl Davis (1983).

” L.B

0

-

T-t.2

FIG. 6. The distribution of v alone v=O.5 for the

In case 4, the numerical results obtained by the PG method is identical to the exact solution given in Eqs. (12-13). Such identity is achieved because the second-order upwind scheme on spatial discmtization and the first-order explicit scheme on transient term are used. Conclusions The proposed PG method has been applied to various fluid flow problems with curvilinear coordinate system. The numerical results have been compared with analytical solutions, experimental data, and/or numerical results available in the literature. The test problem discussed in this paper by no means cover all possible aspects of fluid flows. However, the success in solving these problems enhances our confidence in the PG method.

Chorin, A. J. (1967). A Numerical Method for Solving Incompressible Viscous Flow Problem. J. m, 2, pp. 12-14. Chorin, A. J. (1968). Numerical Solution od the NavierStokes Equations. -ComD., 22, pp. 745-762. Courant, R. (1957). Calculus of Variations and Supplementary Notes and Exercises. revised and amended by J. Moser, New York University. de Vahl Davis, G. (1983). Natural Convection in a Square Cavity: A Comparison Exercise. h J. Num. Meth. in w, 3, pp. 227-248.

Hadow, F. H., and Welch, J. E. (1965). Numerical Calculation of Time-Dependent Viscous Incompressible HOW of Fluid with Free Surface. Phvs. Fluid, 8, pp. 2182-2189. Kawaguti, M. (1953). Numerical Solution of the NavierStokes Equations for the Flow around a Circular Cylinder at Reynolds Number 40. J. Phvs. Sot. &Q, 8, pp. 745-757. Kawaguti, M., and Jain, P (1966). Numerical Study of a Viscous Fluid Flow Past a Circular Cylinder. J. $0~. Japan, 21, pp. 2055-2062. Kwak, D., Chang, J. L. C., Shanks, S. P., and Chakravarthy, S. R. (1984). An Incompressible Navier-Stokes Solver in Three-Dimensional Curvilinear Coordinate System Using Primitive Variables. AIAA 22nd Aerospace Aciences Meeting, Reno. Maliskar, C. R., and Raithby, G. D. (1984). A Method for Computing Three Dimensional Flows Using NonGrthogonal Boundary-Fitted Coordinates. u Numerical Methods Fluids, 4, pp. 519-537. Patankar, S. V., and Spalding, D. B. (1972). A Calculation Procedure for Heat, Mass and Momentum Transfer in Three Dimensional Parabolic Flow. Int. J. Heat Mass Transfer, 15, 1782. Son, J. S., and Hanratty, T. J. (1969). Numerical Solution for the Flow around a Cylinder at Reynolds Number of 40, 200 and 500. J. Fluid Mech, 35, pp. 369-386. Tan, C.-H. (1988). A Numerical Scheme for Solving 3-D, Unsteady, Incompressible

Navier-Stokes Equations.

Ph. D. Dissertation, University of Maryland. Tan, C.-H., and Duh, J. C. (1989). Comparison between Pressure Gradient Method and MAC Method on High Re Calculation. Proceedings 5th Int. Symposium on Numerical Methods in Engineering, Switzerland. Tan&t, S. (1955). Experimental Investigation of the Wake behind Cylinders and Plates at Low Reynolds Numbers. a Sot. Jm, 11, pp. 302-307. Thompson, J. F., Warsi, 2. U. A., and Wayne Martin, C. (1985). Numerical Grid Generations, North-Holland. White, F. M. (1974). Viscous Fluid Flow, pp. 117-118. McGraw-Hill.