- 

NASA Contractor Report 3330 

A Rapid Implicit-Explicit Sol1 

NASA 

CR 

3330 

C.l 

to the Two-Dimensional Time-Dependent 

Incompressible Navier-Stokes Equations 

Joseph E. Davis 

CONTRACT A50807B 

OCTOBER 1980 
TECH LIBRARY KAFB, NM 

NASA Contractor Report 3330 

A Rapid Implicit-Explicit Solution 

to the Two-Dimensional Time-Dependent 

Incompressible Navier-Stokes Equations 

Joseph E. Davis 

Joseph E. Davis 

Saratoga, Califowia 

Prepared for 

Ames Research Center 

under Contract A50807B 

National Aeronautics 

and Space Administration 

Scientific and Technical 

Information Branch 

1980 
A RAPID IMPLICIT-EXPLICIT SOLUTION TO THE TWO-DIMENSIONAL 

TIME-DEPENDENT INCOMPRESSIBLE NAVIER-STOKES EQUATIONS" 

Joseph E. Davis 

20780 Sevilla Lane, Saratoga, California 

ABSTRACT 

A rapid implicit-explicit finite difference scheme has been developed 

for the numerical solution of the time-dependent, incompressible, twodimensional Navier-Stokes equations in conservation-law form using vorticity 

and stream function variables. The algorithm is second-order time-accurate 

and spatially factored. The systems of equations are solved at each time 

step by an iterative technique. Numerical results have been obtained with 

the technique for a circular cylinder at a Reynolds number of 15, and an 

NACA 0012 airfoil at zero angle of attack at Reynolds numbers of lo3 and 104. 

The results are in agreement with another numerical technique, and the computing time required to obtain the steady-state solution at the Reynolds 

number of lo4 was 49.7 set on a CDC 7600 computer using a 65 x 84 computational grid. 

INTRODUCTION 

Rapid numerical solutions to the time-dependent Navier-Stokes equations 

at moderate Reynolds numbers and low Mach numbers have been difficult to 

obtain because of the elliptic nature of the equations in this regime. 

Unlike high Mach number compressible flows, each portion of the flow field 

significantly influences every other portion of the flow field which drastically increases the amount of computation time required to obtain either 

a time-accurate or steady-state numerical solution. 

Fast numerical techniques developed for the solution of the compressible 

Navier-Stokes equations, such as the algorithms of MacCormackl and Beam and 

Warming,2 unfortunately do not perform as efficiently in the low Mach number 

regime as they do at high Mach numbers. J. L. Steger (private communication) 

has applied the latter technique to the laminar flow (Re = 104) about an 

NACA 0012 airfoil at zero angle of attack and a Mach number of 0.2 and was 

able to obtain a steady-state solution on a 71 x 33 grid in about 20 min of 

computing time on a CDC 7600 computer. 

Roache3 has presented a fast numerical technique for the incompressible 

steady Navier-Stokes equations as applied to the internal flow in flush 

*Work performed under Contract No. A50807B for Ames Research Center. 
inlets. However, his program at present does not simulate unsteady flows 

and has only been tested on a very coarse mesh, where the number of iterations required to converge to a steady-state solution may be lower than the 

number required for a fine mesh. Other excellent numerical techniques for 

solving the incompressible time-dependent Navier-Stokes equations; such as 

the methods of Mehta,lO Reddy and Thompson,4 WU,~ and Thames,6 all require 

approximately one-half hour or more of computing time on a CDC 7600 to obtain 

steady-state solutions for moderate Reynolds numbers (~10~) and sufficiently 

fine grids (-5000 mesh points). Even the very efficient method of Ghia, 

Hankey, and Hodge,7 which was developed for cavity flows, appears to require 

about 4 min of computing time at these Reynolds numbers and grid sizes to 

obtain steady-state solutions. It is therefore the purpose of this report 

to present a rapid numerical technique for solving the incompressible timedependent Navier-Stokes equations for the flow about two-dimensional airfoils 

using very fine computational grids. 

The present numerical technique evolved as the result of an investigation intended to obtain numerical solutions to the incompressible timedependent Navier-Stokes equations by utilizing the Newton-Kantorovich technique.9 This technique is a rigorous mathematical concept based on functional 

analysis that transforms a nonlinear partial differential equation into a 

sequence of linear partial differential equations. When this technique is 

applied to a second-order time-accurate differencing of the Navier-Stokes 

equations and maximum advantage is taken of previous time step information, 

it reduces to approximating the Navier-Stokes equations by a single linear 

partial differential equation. The equation becomes linear because the nonlinear term has, in effect, been replaced by a first-order Taylor series 

expansion in function space. A detailed discussion of the mathematical 

foundations of the technique as well as a summary of previous investigations 

utilizing the technique are given in Gabrielsen* and Davis, Gabrielsen, and 

Mehta.9 

As a result of that investigation, new insight was gained to improve the 

numerical technique used to solve the Navier-Stokes equations. This report 

describes these improvements and the results obtained to date using the 

improved technique. 

MATHEMATICAL FORMULATION 

The present numerical technique has been incorporated into a computer 

code that is restricted to symmetric flow fields and geometries. It is a 

special version of a more general computer code that was developed for nonsymmetric airfoils at angle of attack and is described in detail in reference 10. -The mathematical formulation contained herein describes only the 

symmetric case, although the development for the angle of attack case is 

very similar. 

The flow field and geometry under consideration are mapped into the unit 

half-circle. As a result of this transformation, the equation governing the 

2 
unsteady, incompressible flow of a Newtonian fluid may be expressed in terms 

of a transformed vorticity transport equation and disturbance stream function 

equation as 

d2p 2 

m-' 

and 

where 

&L+ a2f0 

w Tz- (1) 

+if!k 

at32 

= -H2r2, (2) 

(3) 

and 

w = vorticity 

Y = total stream function 

tJ = disturbance stream function Y-y 

(x,y) = coordinates of physical flow field 

(r,e) = coordinates of the interior of the unit half-circle where Olrll 

and TT I 8 I 2TI 

P = stretched radial coordinate where the stretching function is defined 

below 

L = dimensionless chord 

R = Reynolds number Ul/v 

U = freestream velocity 

1 = chord 

V = kinematic viscosity 

t = time 

The term J is the conservation-law form of the convective terms. 

That is, 

3 
The radial coordinate r is transformed according to a function that 

stretches the region near the body surface (r = 1) so that proper ,resolution 

of the boundary layer and separation zone may be obtained. The stretched 

radial coordinate is depicted by the variable p and is expressed as 

P = (kl + kz)-l[tanh'l(rkg - k4) + k2] (5) 

with 

k, = I[ r 

+ tanh(k2) 1 - (r 

0 

y 1) 1 

and 

h = tanh(k2) - ro[tanh(kl) + tanh( (r y 1) 

0 

The constants r0, kl) and k2 (all positive) determine the value of p. As 

r varies from r. to 1, P varies from 0 to 1. 

The mapping into the unit circle is accomplished by the transformation 

1 KC2 Z=K+Y+l+yK (6) 

where 

z=x+iy 

K = rele 

A proper choice of the constants y and c invokes the solution for the flow 

over any one of a number of shapes, including a flat plate, a circular or 

elliptical cylinder, or thick airfoils. The trailing edge of an airfoil 

shape may be rounded off by defining 

c = [ 5 + (1 - l-w2 1 (1 - 6) , where 0 < 6 < 1 (7) 

The components of velocity u1 and u2 are defined in terms of the disturbance stream function as 

u1 +s+g, u2 =-$($$f+$) (8) 

For Reynolds numbers much larger than unity, the vorticity in the flow 

field exists only near the body and in the wake. Away from this region, the 

flow is essentially irrotational; therefore, the region of calculation is 

4 
subdivided into two parts: a small viscous region and a large irrotational 

region bounded by p, and p, with 71 2 8 < 82 (see fig. 1). 

The boundary conditions applied along the boundaries of the region of 

calculation are shown in figure 1 as follows: 

1. On the surface (p = l), the constraint of no slip is applied in the 

form 

Y=O (or J,=-Y) (9) 

and 

ay 0 

ap= 

Condition (10) is used to calculate the surface vorticity by the method of 

Israelill described in detail later. 

2. Along the line of symmetry, the vorticity w and the disturbance 

stream function $ are specified to be zero. 

3. The flow at the far boundary is constrained with first-order differential relations obtained from the Navier-Stokes equations by dropping the 

tangential derivative of the pressure and viscous terms (i.e., at the outer 

boundary); the inertia terms are dominant: 

H2r e = -2 ~(5) 

, 

and u2 is obtained from the 9 component of the Navier-Stokes equation. 

That is, 

ab12 + u22) + u w 

a0 1 1 

(11) 

(12) 

(13) 

Note that, in equation (13), w = 0 when n 5 0 I e2. At t = 0, the flow is 

irrotational (without circulation), that is, w = 0 and I/J = -yr". 

These boundary conditions are believed to be superior to specifying 

either potential flow or uniform velocity since eddies or vortices can pass 

through the downstream boundary. Also, because the velocity far away is not 

defined, the circulation there can change with time. Therefore, equation (11) 

correctly represents the vorticity transport through the downstream boundary. 

In equation (13), the absence of the tangential pressure derivative will not 

significantly affect the motion of a vortex through the boundary. 

5 
The surface pressure distribution is obtained by integrating the tangential component of the Navier-Stokes equation. That is, 

(14) 

where p(0) = 0. The pressure coefficient CP is, therefore, equal to 2p. 

On the surface, the tangential stress is given by 012 = (L/R)w. Both p and 

012 are made dimensionless with pU2. The coefficients of lift and moment 

are zero for symmetric geometries and flow fields; however, the coefficients 

of drag are computed as 

CD = 

drag 

(1/2)pU21 = 'DP + 'DF (15) 

where 

2lT 

?Y 

p a0 de 

and 

4 

I 

2lr 

C DF =-R w i!Jc de 

ae 71 

(16) 

(17) 

where the subscripts P and F, respectively, represent the contributions due 

to pressure and viscous forces. 

Equations (1) through (4) represent a nonlinear coupled set of partial 

differential equations. To reduce the difficulty in solving these equations, 

the nonlinearity may be removed by applying the Newton-Kantorovich transformation to the equations9 or by performing a Taylor series expansion in function space to the nonlinear terms. Employing such a procedure alters the 

convection terms in the-vorticity transport equation, and the linear approximation to this equation becomes 

d2p 

iizr 

2 ati : a2w 

ap a02 

- r 2 f [~rl,‘~*) + J(S) - J(s)] 

+ O[AwAY] 

6 

(18) 
where 

Aw=w-w 0 

and 

(19) 

AY=Y-Y 0 

The functions w. and Y, are initial guesses to the solutions w and Y. 

These initial guesses must be sufficiently accurate such that the truncation 

error associated with approximating the nonlinear terms is no worse than the 

truncation error of the overall numerical technique. 

NUMERICAL FORMULATION 

The Navier-Stokes equations were solved by a three-point backward time 

differencing and implicit-explicit factored central space differencing scheme. 

The spatially factored, time-differenced expression for the linearized approximation to the vorticity transport equation9 is 

2At 

- x 6p(yo) 1 tin = FAL Q(wn,Yo) + F(w ,'i' ,Yn) + G(w 0 0 0 .YoY) 1 

- $ (Tpt-- + T2w 

n-2 

> 

+ 4(At) 2 Q(Yo) Bp(Yo) n 

T2 A A w 

+ 0(AwAYAt) + O(At)3 

where the operators 60(Yo) and Ap(Yo) are given as 

(20) 

(21) 

and the cross-derivative terms Q(w 

given as 

YYo), F(wo,Yo,Yn), and G(wo,Yo,Yn) are 

Qbn ,y,> =-rdp!! 

3 dr L [& bo $)- & (Yo 6)] (23) 

7 
r dp R F(wo,Yo,Yn) = ---- 3 dr L 

aye au0 

75 ap + 

( a? ayO 

)- ae w. - ae w. 

+ aP 

ae 

a n awo 

(- 

awO 

aS ' ai - '0 --)I ap 

awO 

ap 

(24) 

n G(u~,Y~,Y > 

alu n ati ay aw alu -- apoW 'F$-$* 

(25) 

The values of the three-point backward time difference parameters are T = 3, 

Tl = -4, and T2 = 1. The transformation parameter A is 

Let the initial guesses at the solution w. and Y. be given by 

w 0 

= 2un-’ - un-2 = wn + o(At)2 

yO 

= 2yn-1 _ yn-2 = yn + O(AQ2 (27) 

except at the first two time steps where they are given as 

w =w n-1 = 

0 wn + O(At> 

y. = yn-1 = Yn + O(At) 

(26) 

(29 

(29) 

Hence, 

AWAY = O(At)4 (30) 

F(wo,Yo,Yn) = O(At)2 (31) 

G(wo,Yo,Yn) = ohW2 (32) 

except at the first two time steps. Therefore, except at the first two time 

steps, the terms F and G may be neglected because they are the same order 

as the overall truncation error. In addition, the second-order term on the 

right-hand side of equation (20) is made explicit by approximating it by 

2 6 4 (At> 6e(Yo) (Yo> n 

T2 A PA w= 

4 (y22 6e(y0) ‘~(‘0) un-l 

A A + O(At)3 (33) 

a 
The final second-order accurate time-differenced expression for the vorticity 

transport equation is then 

- $ Ge(Yo))(l - 2 6p(Yo) wn 

> 

= E Q&Y,) - $(T$-' + T,w n-2 

> 

+ 4(At) 2 $po) 6p (Yo) n-l 

T2 A A w 

+ O(At)3 (34) 

Consistent with the spatial factoring concept, equation (34) is split into 

two equations whose finite difference analogies each produce tridiagonal systems of equations. That is, 

- s Ge(Yo))~* = 2 Q(wn,Yo) - $(T$-' + T2ti n-2) 

+ 4(At)2 6e(Yo) 6 (Y',> 

T2 A A w 

n-1 + O(At) 3 (35) 

( 1-s 6p(Yo) ) wn =w * (36) 

The spatial derivatives associated with equations (35) and (36), as well as 

equation (19), were approximated by central differencing formulas everywhere 

except at the boundaries of the flow field. The truncation error for the 

vorticity transport equation is O[(AP)~ + (A8>2 + (At)2], except at the 

first two time steps where the temporal error is O(At), because a two-point 

backward difference formula is required such that T = 2, T1 = -2, and 

T2 = 0. 

The disturbance stream function equation is solved by the method of reference 12. The details of the procedure are not described here, but the 

truncation error of the disturbance stream function equation is 

~[(AP)~ + (A0)2]. 

At the outer boundary, the finite difference forms of equations (35) 

and (36) are used to describe the vorticity condition, except that the diffusion terms are deleted and forward space differencing is used in the p 

direction. This formulation results in a first-order truncation error in the 

space variable p. The disturbance stream function at this boundary is 

obtained from the finite-difference form of equations (12) and (13): 

(37) 
where 

(u2) 

n 

i,l + 

( 

n-1 

- TlU2 

n-2 

+ ‘2~2 i 1 )I , (38) 

Along the lines of symmetry the boundary conditions w=o and $ = 0 are 

specified. The no-slip boundary condition represented by equations (9) 

and (10) is satisfied by specifying the condition $=O along the boundary 

and by iteratively determining the surface vorticity that numerically 

enforces a+/ap = 0. The latter condition is accomplished by employing the 

method of Israeli.lI That is, the surface vorticity at a particular iteration "k" is related to a$/ap and the surface vorticity at the previous 

iteration "k-l" by 

k 

wJL 

a+ k-1 

=gapJL (39) 

The coefficient g is a relaxation parameter whose optimum value is determined in reference 11 for boundary-layer type flows. For the fine grids and 

relatively large time steps used in the present study, the equations defining 

the optimum value of g as given in reference 11 may be approximated to the 

first order in 

(40) 

where ArJL is the normal distance from the surface to the next grid point 

and v, R, and At are the kinematic viscosity, Reynolds number, and time 

step size, respectively. The value of g for a given geometry and computational mesh must be determined manually once for a particular Reynolds number and time step size. Then, to determine the value of g at any other 

Reynolds number or time step size, equation (40) may be used with reasonable 

accuracy. In fact, figure 2 is a log-log plot of g versus At for a NACA 

0012 airfoil at a Reynolds number of lo4 using a 65 x 84 computational mesh. 

Equation (40) is represented in the figure by the straight line and optimum 

values of g determined manually are shown by the symbols. Note that the 

agreement is reasonable. In fact, the same calculations using the manually 

determined g required only a few percent less compute time than the value 

of g determined by equation (40). Furthermore, the savings in engineering 

and compute time by not manually searching for the optimum value of g more 

than makes up for the difference. More accurate formulations of the optimum 

value of g could be determined. However, for the purposes of the current 

study, this did not appear to be warranted. 

10 
In solving the systems of equations described by equations (2), (35), 

and (36), together with the boundary conditions (37), (38), and (39), a mixed 

iterativelnoniterative procedure is used. Iteration is required because w" 

is contained in the term Q(wn,Yo> on the right-hand side of equation (35). 

Hence, equations (2), (35), and (36) must be solved at each time step by 

iterative improvement. In addition, the boundary conditions (37), (38), 

and (39) are iterative in nature, hence they must also be satisfied by iterative improvement at each time step. However, the coefficient matrices of 

equations (35) and (36) do not need to be iteratively improved at each time 

step. A considerable amount of compute time, therefore, is saved at each 

time step by calculating these matrices and performing the forward elimination upon them outside the iteration loop. Then, in iteratively improving 

the values of wn in Q(wn,Yo) and solving for w* and a new wn, only back 

substitutions need be performed. Convergence of the above-described iteration procedure was considered achieved when the change in surface vorticity 

between successive iterations was less than a prescribed value, typically 

10B3 or 10s4. 

Computations of pressure coefficients on the surface and a determination 

of the loads are made at each time step with a finite-difference integration 

formula derived by combining two four-point expressions (in order to have a 

lower effective truncation error). The details of the formula may be found 

in reference 13. 

DISCUSSION OF RESULTS 

The present numerical technique was used to perform flow-field calculations for four test cases. These test cases included a steady-state solution 

for a circular cylinder at Reynolds number 15, a steady-state solution for 

an NACA 0012 airfoil at zero angle of attack at a Reynolds number of 103, and 

a steady-state solution and a time-accurate solution for the same airfoil at 

a Reynolds number of 104. Results were also obtained for the same airfoil at 

a Reynolds number of 105. The results for that test case, however, indicated 

that-the flow field in the separation zone was unsteady at this Reynolds 

number. At this time it is not known whether this unsteadiness is due to 

numerical problems or whether it is a true flow-field phenomena; therefore, 

the results at this Reynolds number are not presented here. 

A summary of the results of all test cases is presented in table 1. 

Included in the table are the compute times required for a steady-state solution on a CDC 7600 computer as well as the number and size of the time steps 

used and the size of the computational grid. All test cases were initiated 

from an impulsive start and computed to a steady-state solution. A detailed, 

discussion of the results of each test case is given below. 

Circular Cylinder at Re = 15 

The circular cylinder test case was performed on a 33 x 40 grid in the 

half-plane. Plots of the steady-state surface vorticity and pressure 

11 
coefficient distributions are shown in figure 3. Comparisons of the results 

are made with Mehta's' numerical technique. As seen in the figure the 

results are in excellent agreement. Since only a steady-state solution was 

desired for this test case, no attempt was made to obtain time accuracy. As 

seen in table 1 the total compute time required to obtain the solution was 

2.5 sec. A-total of eight time steps were taken and the time step sizes (nondimensionalized as UAt/l) were 2.5 diam each. The cylinder traveled a distance of 20 diam during this time. The number of iterations required to 

obtain a converged solution at each time step varied from 21 at the first 

time step to 7 at the last time step. 

NACA 0012 Airfoil at Re = lo3 

A computational grid in the half-plane of 65 x 48 was used to determine 

the flow field about the NACA 0012 airfoil at a Reynolds number of 103. 

Plots of the steady-state surface vorticity and pressure coefficient distributions are shown in figure 4. As shown in table 1, a total compute time of 

24.8 set was required to obtain this solution. The size of the time steps 

was initially 0.269 chord. However, the time step size was doubled after 

10 time steps. A total of 23 time steps was used to obtain the steady-state 

solution, and the airfoil traveled a distance of 9.69 chords during this 

time. The number of iterations required for convergence at each time step 

varied from a maximum of 32 at the first time step to a minimum of 8 at the 

last time step. 

NACA 0012 Airfoil at Re = lo4 

Two test cases were performed using a 65 x 84 computational mesh for the 

NACA 0012 airfoil at'a Reynolds number of 104. The results of both test 

cases are shown in figure 5 as plots of the surface vorticity and pressure 

coefficient distributions. The results of both test cases are also compared 

with computations using Mehta's' numerical technique. 

The first test case was performed to determine the steady-state solution 

only. This test case required 49.7 set of compute time. Since time accuracy 

was not a consideration, the initial time step size was set at 0.138 chord, 

then doubled every 10 time steps. The steady-state solution was obtained 

after the airfoil had traveled about 9.67 chords, hence a total of 30 time 

steps was needed in obtaining the results. The number of iterations required 

for convergence at each time step varied from 32 at the first couple of time 

steps to 10 at the final time step. 

The second test case at this Reynolds number was performed in order to 

obtain a time-accurate solution until a steady state was achieved. Therefore, 

the calculation was initiated with a smaller time step, 0.00276 chord. However, as before, the time step size was doubled every ten time steps until 

the airfoil had traveled about 9.88 chord lengths. The total number of time 

steps required to obtain this solution was 84, while the total compute time 

required was only 88 sec. The number of iterations required for convergence 

12 
of the surface vorticity varied from a maximum of 23 to a minimum of 6. The 

maximum occurred at the first time step and the minimum occurred at the tenth 

time step. Each time the time step was doubled, the number of iterations 

increased then gradually decreased from time step to time step until the time 

step doubled again. In general, fewer iterations were required for convergence for smaller time steps than for larger ones. 

As seen in figure 5, both test cases were in excellent agreement with 

each other as well as with the solutions obtained using Mehta'sy numerical 

technique. 

CONCLUSIONS 

Based on the results obtained to date, the present method appears capable 

of obtaining rapid steady-state and time-accurate solutions to the incompressible time-dependent Navier-Stokes equations. It has been successfully applied 

to a circular cylinder at a Reynolds number of 15 and an NACA 0012 airfoil at 

Reynolds numbers of up to lo4 using fine computational grids. The numerical 

technique itself may be extended to representations of the Navier-Stokes 

equations in the primitive variables u, v, and p and therefore to threedimensional flows. It may also be extended to the time-averaged Navier-Stokes 

equations with turbulence modeling. 

ACKNOWLEDGMENT 

The author would like to express gratitude to Dr. Ralph E. Gabrielsen of 

Ames Research Center and Dr. Unmeel B. Mehta of Stanford University. 

Dr. Gabrielsen originated the idea of applying the Newton-Kantorovich concept to the time-dependent Navier-Stokes equations, and it was through his 

scientific curiosity and impetus that this investigation was brought about. 

In addition, his guidance regarding the mathematical foundations of linearizing the Navier-Stokes equations was invaluable. Dr. Mehta graciously allowed 

us to use his excellent computer code during the course of this investigation. 

The numerical technique reported in this paper is an outgrowth of his numerical technique, and much of the improvement in the technique for our applications was due to his valuable suggestions. 

13 
REFERENCES 

1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

10. 

11. 

12. 

13. 

MacCormack, R. W.: An Efficient Numerical Method for Solving the TimeDependent Compressible Navier-Stokes Equations at High Reynolds Number. NASA TM X-73,129, 1976. 

Beam, R. M.; and Warming, R. F.: An Implicit Factored Scheme for the 

Compressible Navier-Stokes Equations. AIAA Paper 77-645, Albuquerque, 

New Mexico, June 1977. 

Roache, P. J.: A Semidirect Method for Internal Flows in Flush Inlets. 

AIAA Paper 77-647, Albuquerque, New Mexico, June 1977. 

Reddy, R. N.; and Thompson, J. F.: Numerical Solution of Incompressible 

Navier-Stokes Equations in the Integro-Differential Formulation Using 

Boundary-Fitted Coordinate Systems. AIAA Paper 77-650, Albuquerque, 

New Mexico, June 1977. 

Wu, J. C.; and Sampath, S.: A Numerical Study of Viscous Flow Around an 

Airfoil. AIAA Paper 76-337, San Diego, California, July 1976. 

Thames, F. C.: Numerical Solution of the Incompressible Navier-Stokes 

Equations About Arbitrary Two-Dimensional Bodies. Ph.D. Dissertation, 

Mississippi State University, 1975. 

Ghia, K; N.; Hankey, W. L.; and Hodge, J. K.: Study of Incompressible 

Navier-Stokes Equations in Primitive Variables Using Implicit Numerical Technique. AIAA Paper 77-648, Albuquerque, New Mexico, June 1977. 

Gabrielsen, R. E.: An Effective Solution to the Nonlinear, Nonstationary 

Navier-Stokes Equations for Two Dimensions. NASA TM X-62,493, 1975. 

Davis, J. E.; Gabrielsen, R. E.; and Mehta, U. B.: A Solution to the 

Navier-Stokes Equations Based Upon the Newton-Kantorovich Method. 

NASA TM-78,437, 1977. 

Mehta, U. B.: Dynamic Stall of an Oscillating Airfoil. AGARD Paper 

No. 23, Ottawa, Canada, 1977. 

Israeli, M.: On the Evaluation of Iteration Parameters for Boundary 

Vorticity. Studies in Applied Mathematics, vol. 51, no. 1, March 1972. 

Buzbee, B. L.; Golub, G. H.; and Nielson, C. W.: On Direct Methods for 

Solving Poisson's Equation. SIAM Journal of Numerical Analysis, 

vol. 7, 1970, pp. 627-656. 

Mehta, U. B.; and Lavan, Z.: Starting Vortex, Separation Bubbles, and 

Stall - A Numerical Study of Laminar Unsteady Flow Around an Airfoil. 

J. of Fluid Mech., vol. 67, 1975, pp. 227-256. 

14 
FIGURE CAPTIONS 

Figure l.- Domain of calculation, boundary conditions, and grid notation. 

Figure 2.- Comparisons of the optimum value of g as determined manually and 

by equation 40. 

Figure 3.- Surface vorticity and pressure coefficient distributions on a circular cylinder; Re = 15. 

Figure 4.- Surface vorticity and pressure coefficient distributions on a 12% 

thick symmetrical airfoil; c1 = 0; Re = 103. 

Figure 5.- Surface vorticity and pressure coefficient distributions on a 12% 

thick symmetrical airfoil; c1 = 0; Re = 104. 

15 
TABLE 1. SUMMARY OF COMPUTATION TIMES FOR VARIOUS TEST CASES 

USING PRESENT TECHNIQUE 

TEST CASE 

CIRCULAR CYLINDER 

AT Re=15 

NACA 0012 AIRFOIL 

AT Re=103 

NACA 0012 AIRFOIL 

AT Re=104 

TOTAL COMPUTE 

TIME (CDC 76001, 

seconds 

2.5 

24.8 

49.7 

88.0 

T 

I T - 

TOTAL TIME 

STEPS TO 

STEADY STATE 

8 

23 

30 

84 

TIME STEP 

2.5 

0.269 TO 0.539 

0.138 TO 0.553 

0.00276 TO 0.70; 

GRID 

(0 X r) 

(65 X 48) 

(65 X 84) 

(65 X 84) 

16 
j=JL 

(rJ P, 

+ = -Y, awap = 0 

w=o 

*=0 

j = 1 (ro) p, L 

7T 

i=l 

i = IM 82 

Cd=0 

$ =o 

8, i 

- 

27T 

i= IL 

au2 - 

- [ 

1 aW2+u22) +u, w 

at --2H ae 1 

17 
90PT vs. At 

NACA 0012 AIRFOIL 

65 X 84 COMPUTATIONAL MESH D..- “ci’ IA I”4 

OPTIMUM VALUE OF g 

I DETERMINED MANUALLY 

1011 

.Ol .l 1 10 

TIME STEP SIZE, At 

18 
0 20 40 60 80 100 120 140 160 180 

- MEHTA 

0 PRESENT 

2 

1.6 

1.2 

.8 

.4 

‘0 

-. 4 

0 20 40 60 80 100 120 140 160 180 

Cl.RCUMFERENTlAL ANGLE, 0 

19 
-50 

-45 

-40 

-35 

3 -30 

c 

G -25 

F 

5 

> -20 

-15 

-10 

-5 

0 

) 

0 

0 

0 

0 PRESENT TECHNIQUE 

0 --0 0 

0 0 

0 0 00 

p 01 007 I,, I -I ! 

-. 2 0 .2 .4 .8 1 

DISTANCE, 

if 

20 
I 

-600 

-400 

3 

- 

z 

G -300 

F 

a 

P -200 

-1 oc 

600 

or 500 

i 

G E 400 

iii 

8 

2 

3oc 

$ 20( 

n 

lot 

I 

I 

- MEHTA 

0 PRESENT: TIME ACCURATE 

+ PRESENT: STEADY STATE ONLY 

0 .I .2 .3 .4 .5 .6 .7 .8 .9 

DISTANCE, X, 

21 

-160 

-80 

1.2 

I 

,8 

86 
licit Solution to 

ent Incompressible 

Joseph E. Davis 

20780 Sevilla Lane 

National Aernautics and Space Administration 

Washington, D. C. 

15. Supplementary Notes 

Ames Technical Monitor: Warren Ahtye 

Final Report 

16. Abstract 

A rapid implicit-explicit finite difference scheme has been developed for 

the numerical solution of the time-dependent, incompressible, two-dimensiona 

Navier-Stokes equations in conservation-law form using vorticity and stream 

function variables. The algorithm is second-order time-accurate and 

spatially factored. The systems of equations are solved at each time step b 

an iterative technique. Numerical results have been obtained with the 

technique for a circular cylinder at a Reynolds number of 15, and an 

NACA 0012 airfoil at zero angle of attack at Reynolds numbers of lo3 and lo4 

The results are in agreement with another numerical technique, and the 

computing time required to obtain the steady-state solution at the Reynolds 

number of104 was 49.7 set on a CDC 7600 computer using a 65 x 84 computational grid. 

7. Key Words (Suggested by Author(s)) 16. Distribution Statement 

Airfoils, Transient Flow, Unclassified - Unlimited 

Navier-Stokes, Laminar Flow, 

Finite Differences Star Category - 34 

9. Security Classif. (of this report) 20. Security Classif. (of this page1 21. No. of Pager 22. Rice’ 

Unclassified Unclassified 21 A02 

‘For sale by the National Technical Information Service, Springfield, Virginia 22161 

NASA-Lan9ley, 1960