Finite difference approximations of the vorticity of laminar flows at solid surfaces

evl

0

Technische Ho

NAVAL SHIP RESEARCH AND DEVELOPMENT CENTR1

C. 20007

DATUM:

eagq;

FINITE-DIFFERENCE APPROXIMATIONS OF THE VORTICITY OF LAMINAR FLOWS AT SOLID SURFACES

by

Hans J. Lugt and Yermiyahu Rimon

This document has been approved for public

release and sale; its distribution is unlimited.

DEPARTMENT OF APPLIED MATHEMATICS RESEARCH AND DEVELOPMENT REPORT

i k "v. Scheepsbouwkunde

The Naval Ship Research and Development Center is a U.S. Navy center for laboratory effort directed at achieving improved sea and air vehicles. It was formed in March 1967 by merging the David Taylor Model Basin at Carirerock, Maryland and the Marine Engineering Laboratory (now Naval Ship R & D Laboratory) at Annapolis, Maryland. The Mine Defense Laboratory (now Naval Ship R& D Laboratory) Panama City, Florida became part of the Center in November 1967.

Naval Ship Research and Development Center Washington, D.C. 20007 *REPORT ORIGINATOR SHIP CONCEPT RESEARCH orrice 01470 LDEPARTMENT OF APPLIED SCIENCE AMA . -SYSTEMS DEVELOPMENT ornce mot ANNAPOLIS caawANDIN0 OFFICER ncroncro. oirtcroo

MAJOR NSRDC ORGANIZATIONAL COMPONENTS

DEVELOPMENT PROJECT OFFICES OHM. SO, S0.10 NSRDC CARDEROCR COMMANDER TECHNICAL DIRECTOR RAWL PANNAA CITY COMANDDIO OFFICER recroircot. DIRECTOR

-1

DEPARTMENT OP OCEAN nagooLooy P7I0 DEPARTMENT OF MINE COUNTERMEASURES PI20

1

DEPARTMENT or AIROORNE MINE COUNTERMEASURES P110 DEPARTMENT OFT/MORE WARFARE AND TORFE00

DEFENSE P740 NDW-NSRDC 3960/43 (10-67) otrannaprr or tuolocAL toonou moo 11400 DEPARTMENT or HYDRINECHANICS 500 WENT OF AEROOYNAMICS GOO DEPARTMENT OF MACHINERY TECHNOLOGY A/00 DEPARTMENT OF STRUC TUNAS. MECHANKS 700

*

OCPARTMENT OF APPLIED MATHEMATICS SOO DEPARTMENT OF reoreeisos nctrooLoov AIX DEPARTMENT OF ACOUSTICS AND VIDJUIT1011

DEPARTMENTOF THE NAVY

NAVAL SHIP RESEARCH AND DEVELOPMENT CENTER

WASHINGTON, D. C. 20007

FINITE-DIFFERENCE APPROXIMATIONS OF THE VORTICITY OF LAMINAR FLOWS AT SOLID SURFACES

by

Hans J. Lugt and Yermiyahu Rimon

This document has been approved for public

release and sale; its distribution is unlimited.

TABLE OF CONTENTS,.

Page

ABSTRACT 1

ADMINISTRATIVE INFORMATION 1

INTRODUCTION 1

THE TAYLOR-SERIES ABOUT A SURFACE POINT FOR THE

VORTICITY-TRANSPORT EQUATION 3

METHODS BASED ON SERIES EXPANSIONS ALONG xi = CONST 8

DISCRETIZATION ERRORS 13 TOTAL ERRORS 15 NUMERICAL STABILITY 24 ACKNOWLEDGMENT 25 REFERENCES 26 LIST OF FIGURES

Figure 1 Oblate Spheroidal Coordinates in the Meridional Plane 3

Figure 2 Shapes of Various Oblate Spheroids in the.Meridional Plane 14

Figure 3. Error. of Surface Vorticity for the.Oblate Spheroid 77 = 0.05

Due to Various Finite-Difference Approximatiohs 16

Figure 4 Error of Surface Vorticity for Various Oblate Spheroids Due to

the Finite-Difference Approximations (20) and (21) 17

Figure 5 Error of Surface Vorticity versus 877 at the Tip 0 = 90 Degrees

Due to the Finite-Difference Approximations (20) and (21) 18

Figure 6 Error of Surface Vorticity for Various :Oblate Spheroids Due to

the Lester Formula 19

Figure 7 Error of Surface Vorticity for the Oblate Spheroid no 0.05 Due to

Various Finite-Difference Approximations 20

Figure 8 Error of the Strewn Function at = 0.1 Due to Various Finite

Difference ApprOximations for coo 22

Figure 9 Total Eiror of Surface Vorticity for the Oblate Spheroid no .= 0.05

Due to Various Finite-Difference Approximations 23

ii

NOTATION

A Cross-sectional area

P P A1,, B1, Arbitrary coefficients of the Taylor Series

a Focal distanceofoblate spheroid

CD Drag coefficient

Drag

D 2 Operator defined by Equation (6)

ds Infinitesimal distance

Relative error

h, Metric coefficients, i = 1, 2, 3

Constant defined in Equation (38)

Characteristic length

Re Reynolds number = UL/v

Equals to x1 sin + x3 cos

(r, 95, z) Cylindrical coordinate system Time

Dimensionless time U/L t

Dimensionless time D/Re

Characteristic velocity

Velocity vector with components vi

x2, x3) Orthogonal coordinate system

si Distance between two points in the xi-direction

(77, x, 0) Oblate spheroidal coordinate system

Equal to sinh

ii Kinematic viscosity

Density of the fluid

Angle between axis of symmetry and body surface

girs Wp Coefficients of stream function in the Taylor Series

Azimuthal vorticity Vorticity vector

co* Equals h2

n_{mn p} Coefficients of co* in the Taylor Series

Subscript 0 Related to body surface

Subscript i,j Related to mesh point (i,j)

ABSTRACT

Finite-difference methods for solving the Navier-Stokes equations by

means of the stream function-vorticity formulation require the computation of

the vorticity on solid surfaces with the aid of one-sided difference schemes. Several such schemes of second-order accuracy in the vorticity are derived for axisymmetric flows in orthogonal coordinate systems. Their accuracies are examined by comparing the appropriate solutions with the exact values of flows past oblate spheroids for vanishing Reynolds numbers (Oberbeck solution). It is found that if exact values of the solution are used at the inner points, the

various approximations of the surface vorticity differ considerably from each

other, and that the main source of error is the use of numerical derivatives. However, if the various methods of computing the surface vorticity are incor-porated in an overall numerical scheme, the distinctions among them are

re-duced considerably.

Errors in the drag coefficients computed by the various methods are an order of magnitude smaller than those of the surface vorticity. Experi-ments with numerical stability reveal that the diffusion process near the solid surface dominates the numerical stability properties of the scheme used.

ADMINISTRATIVE INFORMATION

This study was funded by the Naval Ship Systems Command (NAVSHIPS) under the

Mathematical Sciences Program, SR 003 10 01, Task 11572. The completion date was 28

August 1969.

I. INTRODUCTION

A solid body moving in a viscous and incompressible fluid under conservative forces generates vorticity only at the body surface. Subsequently, the vorticity emanates into the

fluid through diffusion and convection and disturbes the whole fluid region. The process of vorticity creation is caused by the adherence of the fluid to the solid surface. Mathematically,

this adherence is expressed by the nonslip-boundary condition which, together with proper

boundary data far away from the body; determines the entire flow field. The vorti city (7,0

produced at the surface is defined by

l

= cur v

o k )0

and can be obtained as part of the solution after the velocity field 7) has been computed. When using the primitive variables (velocity and pressure), a complete solution can be

(1)

achieved without solving for the vorticity field. However, the situation is different if

solu-tions are computed by means of the stream function-vorticity formulation in two dimensions

or with the vector potential and vorticity in three dimensions. These methods have been successfully employed in the last years (see, e.g., References 1 through 5*). Here,.the vor-ticity plays an active part. To obtain a solution, the vorvor-ticity field must be computed at each step of an iterative or a time-dependent process from a Poisson-type equation. Since the vorticity at the solid surface is the source of the flow disturbance, it is a crucial step in 'achieving as accurate a numerical solution as possible for the surface vorticity.

There are essentially two different types of difference schemes for evaluating the vor-ticity at the surface. The first utilizes series expansions with the aid of Equation (1) alone, that is, the expansion is of pure kinematical character. The second type includes the dynam-ical equations and, hence, requires less input from the flow field for the computation of the surface vorticity. The most common procedure of the first type is a second-order 'expansion.

Knowing the tangential velocity v 7, at three equally spaced grid points 03 1, 2 normal to the

surface (where (v 00 a 0 due to the nonslip condition), one finds for the vorticity at the

sur-face according to (1)

*References are listed on page 26.

" 1

= 4(v )

28 '

Here, is reduced to the scalEtr coo 'since only two-dimensional flows will be considered.

The distance between the grid points is designated by a. An example of the second type is the Lester formula6 which will be discussed in more detail later in this paper.

The objective of this study is to investigate the accuracy and stability of a number of

expansion methods. The equations are derived for axisymmetric orthogonal _{coordinate systems.}

The body contours are comprised of constant coordinate lines of one of those systems. Equa-tions for plane flow problems are not presented. They can easily be derived in an analogous

way.

Calculations have been performed for the special case of flows past oblate spheroids. This class of bodies was selected for the following reasons:

The effect of large surface curvature _{can be examined.}

An exact solution is available for vanishing,Reynolds number.7

A computer program is available for nonzero Reynolds numbers.8

The accuracies of the approximations for coo are checked by using an input both exact and approximate values of the flow field _{at the inner points. Singular body points as the}

limits of infinite surface curvature are excluded. The flow behavior at such points is not yet clarified and is still the subject of investigation.9'19

2

II. THE TAYLOR SERIES ABOUT A SURFACE POINT FOR THE VORTICITY-TRANSPORT EQUATION

Axisymmetric flows are assumed past bodies with surface contours which can be ex-pressed by a constant coordinate line in an orthogonal coordinate system (xi, x2, x3).

Sing-,

ular body points are excluded. The system (xi, x3) is fixed in such a way that the lines xi = const are normal to the surface, x2 is the azimuthal angle, and x3 = (x3)0 is the contour of the body (Figure 1). The infinitesimal distance ds is defined by

(ds).2 _(dx1)2 2

_(dx)

_(dx3)2

(3)

where the hi values are the metric coefficients.

Figure f Oblate Spheroidal Coordinates in the

Meridional Plane

z2= X, X3=77

With this convention and with no azimuthal flow, the vorticity-transport equation for

an incompressible laminar fluid flow has the form (see, e.g., Rosenhead 11)

I, 2 a(1,b,h)

(

D2 a)h, a, h at 2 h 1 h2h 3 a(xi,x3) 2 1 c)(, h20) 1h2h3 a(x1,z3) (4)

where, by definition,

D2 = h2co -- co*

Here, 0, u, t, v designate the stream function, the vorticity component in the x2-direction, the time, and the constant kinematic viscosity, respectively. The operator D2 is

/1

a20)

(i)o =

h2h ath2

.3 2 3 0

At the surface x3= (x3)0 the nonslip-boundary condition must be satisfied, that is, the

1 ao 1 ao

velocity components v1= , v3 = - .

_ _ ,,

which are used here only as auxiliary

h

h h ax

_{2 3 3 2h 1 3x1}

functions, must vanish there. This yields in terms of the stream function 0:

. ao ao

= (x3)0 : = 0, - 0

01X3

Then, the vorticity at the surface co0(x1) follows from (5) to be

Near the surface point (x)0, (x3)0 a solution of basic Equations (4) and (5) is

assumed in the form of a Taylor-Series expansion:

1 (

a a p w*[ 1)0 + 81 (x3)0 + 83 =p

I0

= +

)

W

83-

* 1 1 p.

and for ik correspondingly. The following abbreviations are introduced:

[

aPw* (axl)m(ax3)n 00

Fago

_ Lazdr(dx3)s]oo 4 p = m + n (10) (7) (8) 9 1 a2 a2 D- = +

+--,

h2 a h3 a a hi (6) h2 ax 2 h32 ax/ ax h1h2

ax

h2h3 ax3 ,(x3 (9)

where the subscript 00 refers to the point (x1)0 , (x3)o . It follows from Equation (5) that

p = q 2. The boundary condition Equation (7) and the convention that the body surface is a zero-streamline imply.

The remaining 2 ± 1) coefficients Omn and tlirs can be expressed by two linearly 7)-= o

independent coefficients A and .B for each order p. In general they are functions of time. In fact, the existence of two linearly independent coefficients can be seen immediately

since for each order p, Equation (4) imposes (p 1) and Equation (5) (p + 1) conditions on ilmn and 'II,. also Reference 9). Altogether, there are (2p + 1) independent coeffi-cients Ai and Bi, which must be determined by the flow quantities at some inner points of the field.

With the subscripts x1 and x3 denoting derivatives with respect to x _{and x3, one} obtains for values of p < 2:

For p = 0: For p = 1: For p = 2: = A1₁ IY12 = B1 111 = A0 0 1 fl01 =

A +

(h2

1 3/00 1 10 = (hD00 11104 = A2 = B2₂ 1 1100 -

₂

300

Ao (from (5) )

RI-)x

_{h 2 3 hill}h2 _hh1 h3)x3 00 AO 3 2 B1

ih12)xil

00 3 5 (12) (13) (from (5) ) A (from (5) )

002 =

-[

h 2 2 h2 1 A2 2(hl )00

h2..)

(V/-

3

)xl

2

_h2)

3 1 1 h h h 31 _{x1x 1 3 xl} h _(1/2h3) h h.2 x3x3 h _hi 12

(

_

-ioo

A0 (from (4) and (5) ) - hih3

xi VI

_x 1 h2h3 ( 1 2 h1

[(

h h3₃ h1h2 3 1h2 x11₀₀ h 3x3 h h_{1 2}

(h

₎

h3 _h2h3 dA0

B +

2v -7)on dt h32) xi x x3₃ 2 3 x oo Ao (from (4) and (5) ) h2 3 x3x3 B 2 h 2

(h

hl chi )

1

h243. (1)

h

h)

/ h

7? \h h312 h1

2 1 x3 3 x3 2 3 _x3 + 1 dA0

2v dt

. 1 SI = _B2 11 _{(h,_. ..)3} 3 0 0

(h2

h2 hih3 3 3 1h3 00 hih3 1/2 (h,h21h3)

\

h2h3)_{x 00}

(_.)

2 h 2 + 1h 1h2 1 3 + h3 h1h3 2 2 712 2 x 3 3 X1X1 3 2 "2 ". - 2 1

+ --

1

2 (

hl 2 A

(h2 h)

2 h3h2h3 3 x3

lB

300

Ao (from. (5) ) (from (4) and (5) ) (14)

The vorticity-transport equation (Equation (4) does not contribute to the first two orders

p = 0 and p = 1. This means that the series expansion around a point of the surface under the nonslip condition does not require the knowledge of the equations of motion up to 0(82) for the vorticity. This purely kinematical behavior is also verified if one uses the Navier-Stokes equations in combination with the continuity equation instead of (4) and (5). The first two

h1h2 113 00 B

hi \

+

Fv\h

3 2 00 dA0 /133

\ h h /

cit hih2 \h2h3 _x3x3 hi )x3100 Ao

orders of such a series expansion fort are determined by the continuity equation whereas

the momentum equation enters for p = 1 only to compute the pressure.

As a consequence of this behavior, a distinction between the two types of procedures mentioned earlier can be made only if one considers orders of p > 2. In this case the number

of independent coefficients to be determined by the flow field is 2p + 1 with both Equations (4) and (5) used, whereas the computation with (5) alone makes (p + 1) (1 + p/2) coeffi-cients available. Although the first method requires less information from the flow field, a word of caution is in order. The nonlinear inertial terms do not appear in (4) for p < 2. It

can be shown (see, for instance, Reference 9) that they will occur for p > 3. This follows

from the quasilinear form of the vorticity-transport equation. Hence, the terms (12) through (14) describe the slow-motion layer adjacent to the surface. This layer must diminish with increasing Reynolds number. Thus, the radius of convergence of the series expansions must decrease too.

The local time-derivative of co occurs first in Equation (14), that is, at p = 2.

Computations have shown that if points parallel to the surface are selected, coo is not very accurate because the streamlines are almost parallel to the surface unless one oper-ates near a separation line. Therefore, best results are expected for points normal to the

surface.

METHODS BASED ON SERIES EXPANSIONS ALONG xi =CONST

If one restricts the series expansions along the coordinate line x1 = const, that is

one obtains from (12) through (14):

For p = 0:

-

lx 1 SP

p=0 P!

P\ 11

3 p=0

(P 2)!

1 2

(xi) or 2

(17)

For p = 1: For p =2: 1 -A2 2(h: )0

[2 (

(ay,

c 1; gz3 3 2 approximations are: 1

(h)

h

\

+ _(h12 3

)

h x3 2 3 3 1

h2h3 ()

h xi x o

h2)

)

1 --=A (xi) = _{A2(x )} dA0 + 2v dt (18) (hhih 2 3 _x3 (19) j) a

This expansion is accurate up to 0(833) for the vorticity and requires the determination of four coefficients. By truncating the series after p = 1, the remaining terms contain only the two coefficients A0 and A1 _{and do not include the derivative} dA0/ dt. There are various

ways to determine the two coefficients by mea4s of the flow quantities at the inner points. Because of the importance of this second-order truncation, the following five pairs of quant-ities at either the grid point 1 or the two equally spaced grid points 1 and 2are considered:

( (30

(001 (v1) ; , 02 The corresponding

0 2

-(114

--h12 )

hh2 (

12.1 ) 53 I 3 2 3 3 x 3 x 3 0 cuo 601

_(h2)t

(7,;h1

2

e

2)0 '" 3 2 0 3 1

2[4(--4)

(211

)

Ox3/ 1 3 2 0 3 1

coo = TT-

_{2 h3 8}

[I

4(v )₁

-(v1)2]

1

coo -

hh[

2az

01 -

2(---Lba

3)1

a]

3 2)0 3 1 2(h 2 h 82 3 2)0 3 10 (- aai3)21. (801

-

fr2) (20) (21) (24)

It may be mentioned that the combination col , 2 cannot be employed. Since the

approxi-mation for _(Jo is the only place in the stream function-vorticity formulation where the nonslip-boundary condition is built in, one of the flow quantities at the inner points must be either

v1, IP, or ac1i/ax3.

Further second-order approximations can be derived when utilizing velocity com-ponents which are different from vi and v3. For example, the velocity components vr and

vz in the cylindrical coordinate system (r, 0, a) are advantageous since their physical

inter-pretation is most convenient. The vorticity at the body, then is given by

hl 1 (91 laxiN

cu°

(h )

28

[(-57) 14001

(Vr)2}

t-±3

flevz)i

(V )

}]

(25)

3 0 3 0 v Z 0

This equation follows from (8) if one considers the relations and

E

-ax ax .

hi_v

h, v + v ar r äz3 3

-

_[4(v

(avr)

1 aVz 1 az

₃₀

28₃ [4(vr), r 2

(az

_{3 0}

All six equations (20) through (25) are approximations up to 0(8 ) in the vorticity.

Equations (21), (23), and (24) have the following 'special property:

If one defines the relative error E by

(c0n)_{- APPROX} E

-

_

(c)

_{- EXACT} ao)3 1 ax3 2

(az)

283;173(x1) 11 (vz (27)

Since in most practical cases, 2-0,(x1) is proportional to 7-11(xl) with sufficient accuracy, the relative error (28) becomes independent of xi. The same can be shown for Equations -(23) and (24).

Approximations of the next higher order (which include the terms of p = 2) may be divided into those which require only the kinematical relation, Equation (5), and those which in addition utilize the vorticity-transport equation, Equation (4). An example for the first group is the next higher order analogous to (21), if the derivative act,/ az3 is given at three different, equally spaced points:

1

[8

:(21-b)

- 9 --

alfr +

2(

ba az

60P2)0 83

3 1 3 2 aX3 (29)

and if one considers (a)o)_{EXACT =}

/

(hp2 )0 , p_neobtains with (17) and (21)

The second group of approximations simplifies considerably if x3 is globally normal to the surface, that is 7/3 const , and if d71-0/ dt 0. Then, only the two coefficients Z13 and Ai

remain. If one chooses, for instance, 01 and (Di as in Equation (20), one arrives at

6 (7/2)1 - CO h 2h 2'06 2 (11.2)0 3 3

9r

2(111)

_{[2 r(}

112-111 2 +

2))

h1 h,2

_{h 1h2}

2 4 21632 x3 063 x3x3 x3 3 (30)

An interesting special case arises for conical surfaces. This results in the Lester

formula.6 With h1 = h3 = 1 and h2 = x1 sin 95 + X3 cos 95 = F , where ck is the angle

be-tween the axis of symmetry and the surface, Equation (30) reduces to 6 co 0 co0 F 2 0 3 83 COS 8 2 COS2c63

42

0

Equation (31) is the formula proposed by Lester6 for axisymmetric bodies of arbitrary shapes for which the angle ck varies along the surface. Accuracy has been claimed up to 0(833). However, since the surface curvature is not considered due to the assumption

tfr axir )0 =0, r = 2, 3, ... , values for cv13 of such accuracy can be approached only if

the curvature effect becomes negligible. This will be demonstrated in Sections IV and V. The advantage of using approximations of higher orders than 83 in the vorticity (that is, p > 2) depends on the accuracy of the flow quantities at the inner points. Based on the

experience which was obtained with the IBM 360-91 and which i s recorded in Section V,

approximations higher than 0(83) in the vorticity.do not appear justified yet.

In numerical calculations, the values at the inner points obtained froth the stream

function-vorticity formulation are represented by ifr and 0, only. Thus, in order to compare the influence of exact and approximate input it is necessary to replace do/ ai,3 and v1 by the difference quotients

7.0

1

3x3) 28

(c÷1

h h ). 6

_{2 3 I 3} (0i

Oi=i)

(32)

12

-0

Equations (21), (22), and (23) then respectively assume the form

1 4(h 2h )_{3 2 0} 2

OP' + 4 02 - 03).

13 (33) (35) 1 202 co0 2(h3)0 1 (h2h3) 1 2( h2 h3)2 (34) 1 (00 _03...h2) (6 01 02) 2 83

The replacement of Equation (25) in terms of the stream function is more laborious and re-quires tfr-values also in the z1-direction. In Equation (26b) the velocity components vr and

vz must be replaced by v1 and v3 with the aid of (26a). Then,

v1 and v3 are expressed by

(

1 Vii+1. , i - 1 , i

27

3)

i,j 263 1 ' (V3)i,/ = - (Tt --.-1h,2) 01,1+1 2 81

It follows that co₀ is a function of the stream function of the form

6)0 = f(01, j 1 j- 1

02,j

(36)

;111,i+ 1 1A2,i+1

where j is the index of the variation along the lines of x = const.

IV. DISCRETIZATION ERRORS

The various approximations for wo , which were derived in the preceding section, are

now tested for flows past oblate spheroids at Re = 0. Exact values of the quantities needed at the inner points are used. Therefore, the errors computed here are due to the discreti-zation process only. The natural coordinates for this flow problem are the oblate spheroidal

coordinates with xi = 0, x2 = x, x3 = and with h12 = = a2 (cosh2 77 sin2 6), h2 = a cosh/ sin 0 where a is the focal distance (Figure 1). The surface of the oblate spheroid is 77 = no. The focal distance a is prescribed in such a way that all oblate spheroids have the radius

unity, that is, a = 1/cosh 770 (Figure 2).

For the uniform flow past an oblate spheroid with (vz). = 1 and for vanishing Reynolds number, the following exact solution due to Oberbeck7 is available:

where A2+1 1

-- sin 2 0

1-

1/

2 A2+1 A24./ 0

K=

The corresponding vorticity is

co = sinh A 2 0 A2+1 A2+1 0 0 K cosh 770 14 arc cot x)] A 2 0 A 2 0 1

arc cot A.0

2sin0

tanh

cosh 2 77 - sin 2 0

Figure 2 Shapes of Various Oblate Spheroids in the Meridional Plane

(37)

(38)

The exact solutions (37) and (39) are used to test the approximations (20) through (25) and (31) through (36). The relative error E defined by (2'7) is displayed in Figures 3 through 7. The relative error E versus 0 is plotted in Figure 3 for the app' rOximations (20) through (25) with 63 = 6 0.05 and no = 0.05. The relative error E increases toward the tip of the oblate spheroid for (20), (22), and (25) and then drops sharply to a minimum directly at the tip 0 = /2. As expected from the discussion of Equation (28), the errors due to Equations (21), (23), and (24) are independent of 0. Their accuracy is remarkably good. This is

con-firmed in Figure 4, where E is plotted versus 0 with 770 as a parameter for Equations (20) and

(21). As the oblate spheroids become spherelike, the error due to (20) approaches a uniform value. Additional calculations, which are not recorded in this paper, show that the shapes of the error curves do not change when 877 is varied. The relative error E at the tip 0 _{= al}2 is

plotted in Figure 5 versus 81 for (20) and (21).

Figure 6 shows the relative error E of the Lester approximation, Equation (31), for various no with 377= 005. The curves clearly demonstrate the good accuracy of (31) as long as the curvature is negligible. In regions of high curvature, however, values for the relative error E even surpass those obtained in methods based on the lower order accuracy (compare Figure 6 with Figure 4).

Large errors are introduced by replacing aq,/ z3 and vi by the corresponding differ-ence quotients (32); this results in the approximations (33), (34), (35). Figure 7 shows the

relative error E versus 0 for the approximations (33) through (36) with 61 = 0.05 and no = 0.05.

Again, E is independent of 0 for (33) and (35). Howeyer, the_{error is extremely high : 24.5} percent for Equation (33) and 34.4 percent for Equation (35). Therefore, it is believed that

much better accuracy would be achieved by improving the differencing methods for obtaining

ao/ ax3 and vi rather than by using higher order approximations of the surface vorticity. The reader is reminded that all the results obtained in this section are based on the knowledge of an exact solution. Such knowledge is practically nonexistent in cases of in-terest. The final appraisal of the various approximations must await a test with the com-plete numerical procedure which will be discussed in the next section.

V. TOTAL ERRORS

Numerical computations of the entire flow field do not furnish exact data for ;II and co at the inner points, as was assumed in the preceding section. Rather these data are approxi-mate, and the question arises as to how this affects the _{accuracy of the surface vorticity coo.} In order to answer this question, the flow field about the oblate _{spheroid with 7/0 = 0.05}

was computed for vanishing Reynolds number'by'rneans of_{a numerical technique.8 This is a} time-dependent stream function-vorticity formulation. Unlike_{Reference 8 where a potential}

solution was used as the initial condition, here the _{slow-motion Oberbeck solution was}

em-,

ployed. This was done to save computing time since the analytic solution does not differ

18 16 14 12 10 16 15° 300 450 600 750 900 0

Figure 3 Error of Surface Vorticity for the Oblate Spheroid

770 = 0.05 Due to Various Finite-Difference Approximations

8 = 0.05, Re =

77

16 14 12 10 4 2 0 =2 4 6 78 10 12 14 6 EQ. (20) EQ. (21) 110 = 0 135 I = 0.1 = 0.5 - 5.0 . lt Jr.841%11%11SP

Figure 4 Error of Surface Voracity for Various Oblate Spheroids Due to the Finite-Difference Approximations

77 = 0.05, Re = 0 17 0° 15° 30° 45° 6 60° 75° 90°

9

7

3

1 EQ. (20) EQ. (21)

Figure 5 Error of Surface Vorticity Versus 577 at the Tip 0 = 90 Degrees Due to the Finite-Difference Approximations (20) and (21)

Re = o

18

-5 =15 20 30 45 19 0 15° 30° 45° 60° 75° 90° 0

Figure 6 Error of Surface Vorticity for Various Oblate Spheroids Due to Lester's Formula (31)

161

24

32:

Figure 7 Error of Surface Vorticity for the Oblate Spheroid

770 = 0.05 Due to Various Finite-Difference Approximations

= o.o5, Re = 0 20 °/0 48 40 32 24 16 8

---EQ. (33) _ EQ. (34) EQ.

...

EQ. (36) 15° 30° 45° 60° 75° -90° 8

greatly from the solution of the corresponding finite-difference system. The same inner-point differencing and boundary conditions Were used in all cases. The only distinction among the cases was the approximation of we. The relative errors E computed from these results repre-sent the total error due to both the (,) - approximation and the errors inherent in the numer-ical method.

The errors in 0 and co at the inner points also depend on the choice of the

approxi-mation forco_{o' It is of interest to compare the numerical values of 0 at the inner point 1}

with the exact Oberbeck solution. The result is presented in Figure 8. Equations (20), (31), (33), and (36) were selected for approximating coe. All curves revealed large errors for 0 except in the Lester case, Equation (31). The magnitude of the relative errors is due to the

fact that the 0-values are extremely small near the solid surface which is a zero-streamline.

With the approximate values of 0 and co at the inner points, the relative error of coo

was computed for (20), (31), (33), and (36). This relative error is plotted in Figure 9 versus 0 with 770 = 0.05 and 6 ---- 0.05. _{It is of interest to note that despite the huge errors in the}

77

values of 0 at the inner points and the great difference among the various approximations, the total errors in wo were similar for all methods and were of the _{same order of magnitude.} The similarity of the curves was also in contrast to the great differences in theerrors of the

various methods due to discretization (see Figure 7). _{In fact, Equation (36), vvhich generates} the worst discretization errors among the second-order methods and also the worst errors in

cfr, had the least total errors in (Do. _{Furthermore, the neglect of the curvature in the Lester} method, Equation (31), which is of next higher order, produced larger _{errors near the tip than} did the second-order method. However, where the _{curvature effects are indeed negligible}

(0 < 0.< 72 deg), the higher order method does offer better results. The drag coefficient

CD 1

p AU2 2

computed as in Reference (9) for the various methods and the _{relative error of CD in percent}

are listed below. The quantities D, p, and A designate the drag, the density of _{the fluid, and} the cross-sectional area of the body, respectively.

21

(40)

Methods

_{C Re}

Relative Error

percent

Equation (20) 20.228 0.82

Equation (31) 20.130 1.30

Equation (33) 20. 317 0.38

20 16 12 _12 16 20 -24 32 36 -40 44 =48 EQ. (20) EQ. (31) EQ. (33) EQ. (36) 7-98r 22 8_ =

Figure 8 'Error of the Stream Function at ni - 0.1 Due to

Numerical Methods by Using Various Finite-Difference

Approximations for (Do 0.O5, R = 0 -90° 75° 0° ---15° 30° - 45° 60° 0

24, 21 18 15 12 6 3 0 3

6

9 /2 15 18 21 24 EQ. (20) EQ. (31) EQ. (33) A EQ. (36)

wJJi

wWc"Vig.'

*6A4/0( '11WeVIttr_sw 41. 27 0° 45°

Figure 9 Total Error of Surface Votticity for the Oblate Spheroid = 0.05 Due to Various Finite-Difference Approximations

8 = 0.05, Re = 0 77 23 75° 90° 306 15°

The exact value for i 0.05 was CD Re= 20.395. The results show that the errors in CD

were an order of magnitude smaller than those of coo and IA . This confirms the well-known

fact that CD,.which is a quantity obtained by integration, i.niti-Ch less sensitive than the field values of IA and co.

VI NUMERICAL STABILITY

The question of numerical stability is of great concern for constructing a solution by a time dependent or an iterative process. If such stability considerations dictate a very small time step, the numerical procedure must be repeated many times. The computer time required can become prohibitive as for example in the case of small Reynolds numbers.8'10 The inadequacy of the linearized stability analysis has already been discussed in detail in References 8 and 10. The main conclusions in these studies were: (1) the linearized sta-bility criterion fails if large vorticity gradients exist in the flow field and (2) a stasta-bility analysis which is based only on the finite-difference formulation of the inner points only is not complete because the effects of the boundary conditions must alto be considered.

Both points are also demonstrated in the present work. To obtain a numerical solution to the time-dependent flow past an oblate spheroid at Re= 0, the time t must be made

nondi-mensional differently from Reference 8. This is done by including the Reynolds numbers Re=t1L/7., in the dimensionless time as t* = Ut/LR , where L and U are the character-istic length and velocity, respectively. If this is not done, that is, if the dimensionless time t = Ut/L is used, the time dependent part of the vorticity-transport equation will vanish at

Re = 0, and the maximum stable time step A tmAx will become zero. It was observed that the solution remained stable only for a time step A V` at least ten times smaller than A/tfor Re= 10. This is explained in the following way. With the new dimensionless time t*, the vorticity-trans port equation becomes in the limit as Re-* 0 identical to the diffusion-type

equation

(a

at*

-

D2) h2 (I) = 0 (41)

/\

Where D2 , /7/2 , ^CO are the dimensionless quantities of D2, h2, , respectively.

If one examines the stability of the complete vorticity-transport equation for nonzero

Reynolds number by neglecting the nonlinear terms and by using /t one obtains equation

81

20

,\)

(42) .e

For the special case of Re= 1, Equation (41) is identical to (42). Hence, one expects that

if ADmAx of (42) is proportional to Re as shown in Reference 8, A t*mAx would be 1/10 of ADmAx for Re= 10. This is confirmed by the numerical computations.

The neglect of the nonlinear terms in the stability analysis can be justified only if the diffusion process near the so/id surface dominates the numerical stability. The fact that

A t*mAx was, indeed, found numerically- to be exactly the expected value suggests that this

conjecture is valid.

Thus, the linear dependence of A ImAx on Re reported in Reference 8 can be explained by the conclusion made in Reference 10. It is demonstrated here that the process of diffusion near a solid boundary dominates the numerical stability. This also explains Conclusion (1) on page 24 that the linearized stability criterion fails if large gradients of vorticity occur.

The above conclusion holds the Reynolds-number range for which the linear

proportion-ality was demonstrated in Reference 8. Slight deviations beyond that range do not disprove

the conclusion since the overall quality of the numerical solution becomes questionable for

such high Reynolds numbers.

At Re= 0, all solutions constructed with the approximations discussed in Section III

were computed with the same sequence of A t*, starting with A t* _{= 3- 10-6 and increasing to}

A t*mAx = 1.25 10-5. The solutions with the approximations _{given by (20), (31), (33), and} (36) remained stable and approached a steady state. However, the _{use of (24) and (35) with}

the 01 _{and 02 values gave only a very strong instability even for the minute time step of}

A t* = 3 . 10-6 . These solutions oscillated violently and diverged within ten time steps.

A similar behavior was also observed for some experiments at Re= 100 which are not

re-ported here. In this case also, some of the approximations for _{coo caused the numerical} scheme to be highly unstable-while others remained stable _{for the same time steps. This is}

in agreement with the conclusion that the boundary conditions and the diffusion process near the solid surface are of extreme importance to the numerical stability.

From the viewpoint of numerical stability, Equations (33) and (36) provided the most

favorable scheme. This was also confirmed in _{some other numerical experiments on problems}

with different geometries which will bepublished later.

ACKNOWLEDGMENT

The authors thank Mr. S. Ohring for performing _{most of the calculations and for} care-fully checking the formulations.

REFERENCES

Fromm, J. E. and Harlow, F. H., "Numerical Solution of the Problem of Vortex Street Development," Physics of Fluids, Vol. 6, p. 975 (1963).

Pearson, C. E., "A Computational Method for Time-Dependent Two-Dimensional Incompressible Viscous Flow Problems," Sperry Rand Report SRRC-RR-64-17 (1964).

Rimon, Y. and Cheng, S., "Numerical Solutions of the Time-Dependent Incompressible

Viscous Flow over a Disk or a 'Sphere," Physics of Fluids Vol 12 p. 948 (1969).

Proceedings of the IUTAM Symposium on High-Speed Computing in Fluid Dynamics, Monterey, California, Supplement II to Physics of Fluids (1969).

Aziz, K. and Hellums, J.D., "Numerical Solution of the Three-Dimensional Equations of Motion for Laminar Natural Convection," Physics of Fluids Vol. 10 p. 314 (1967).

Lester, W.G.S., "The Flow Past a Pitot Tube at Low Reynolds Numbers," University

of Oxford, Department of Engineering Sciences Report 3240 (July 1960).

Happel, J. and Brenner, H., "Low Reynolds Number Hydrodynamics," Prentice-Hall,

Inc. (1965).

Rimon, Y., "Numerical Solution of the Incompressible Time-Dependent Viscous Flow Past a Thin Oblate Spheroid," NSRDC Report 2955 (Jan 1969).

Lugt, H.J. and Schwiderski, E. W., "Flows Around Dihedral Angles," Proc. Roy. Soc.,

Series A, No 285 p. 382 (1965).

Rimon, Y. and Lugt, H. J., "Laminar Flows Past Oblate Spheroids of Various Thick-nesses," Physics of Fluids, Vol. 12 p. 2465 (1969).

Rosenhead, L., "Laminar Boundary Layers," Oxford at the Clarendon Press (1963).

INITIAL DISTRIBUTION_{_ .}

Copies _Copies

1 DNL

1 NASA, Wash D.C. , Tech Lib

CHONR 1 NASA, Marshal Space Flight Ctr

1 Dr. P. King _{Huntsville, Ala 35809, Tech Lib}

1 Dr. S. Silverman

1 NASA, Lewis Res Ctr, Cleveland, Ohio 44121 1 Mr. M. Cooper

1 Mr. R. D. Cooper Tech Lib

1 Dr. R. D. Ryan

NASA, Ames Res Ctr, Moffet Field, Calif 94305

1 Dr. B. J. Macdonald

1 Dr. W. J. McCroskey 1 Dr. L. D. Bram

1 Dr. E. D. Martin

4 NAVSHIPSYSCOM 1 Tech Lib

1 SHIPS 031

1 NASA, Langley Field, Tech Lib 1 SHIPS 0311

2 SHIPS 2052 _{1 USAEC, Tech Lib}

1 NRL, Tech Lib

1 Oak Ridge Nat Lab, Tech Lib 1 NOL, Tech Lib

4 Los Alamos Sci Lab, Los Alamos,

New Mexico 87544

1 NWC, Tech Lib

1 Dr. F. H. Harlow

7 NWL 1 Dr. C. W. Flirt

1 Mr. B. Smith 1 Dr. B. J. Daly

1 Dr. C. J. Cohen 1 Tech Lib

1 Dr. A. V. Hershey

3 Brookhaven Nat Lab, Upton,

1 Dr. E. W. Schwiderski

Long Island, N.Y.

1 Dr. P. Ugincius 1 Dr. P. Michael 1 Dr. B. Zondek 1 Dr. D. W. Lick 1 Tech Lib 1 Tech Lib 1 NUWC, Tech Lib

1 NCCCLC, Tech Lib

2 Naval PCSCHOL, Monterey

1 Dr. T.H. Gawin 1 Tech Lib 2 USNA 1 Dept of Math 1 Tech Lib 20 DDC

2 U.S. Army Math Res Ctr, Univ of Wisconsin, Madison

1 Dr. D. Greenspan

1 Tech Lib

1 U.S. Army Res Off, Durham, N.C. 27706

CRD-AA-IPL Box CM 27 3 N BS 1 Dr. H. Oser 1 Dr. W. Sadowski 1 Tech Lib

2 Nat Sci Foundation, 1520 H St.,

Wash D.C. 20550

1 Math Sci Div 1 Engin Sci Div

4 Princeton Univ

1 Prof S. I. Cheng

1 Dr. J. Smagorinsky

1 Dr. K. Bryan

1 Aerospace and Mech Eng Lib 1 Prof. M. Kruskal

3 Harvard Univ

1 Prof. S. Goldstein 1 Prof. G. Birkhoff

3 Courant Inst of Math Sci, N.Y. Univ, New York, N.Y.

1 Prof. H. B. Keller 1 Dr. A. J. Chorin 1 Dr. Burstein 4 Stanford Univ 1 Prof. A. Acrivos 1 Dr. L. G. Leal 1 Prof. M. D. Van Dyke 1 Prof. I. Flugge-Lotz

3 Brown Univ, Prov., R. I.

1 Prof. W. Prager 1 Prof. J. Kestin 1 Prof. M. Sibulkin

Poly Inst of Brooklyn

1 Prof. G. Moretti

1 Prof. R. C. Ackerberg

4 NCAR, Boulder, Colorado 80301

1 Dr. D. K. Lilly 1 Dr. A. Kasahara 1 Dr. J. W. Deardorff

1 Dr. J. Kornfeld

Bell Tele Labs, Whippany, NJ. 07981

1 Dr. G. S. Deem 1 Dr. N. J. Zabusky

3 Gen Electric Space Sci Lab, Valley Forge, P. 0. B. 8555, Phil, Pa 19101

1 Dr. S. M. Scala 1 Dr. P. Gordon

1 Dr. R. Shani

1 Prof. J. S. Allen, Penn St Univ, Dept of

Aerospace En g, Univ Pk, Pa 16802

28 Copies

1 Dr. J. E. Fromm

Dept 977, Bldg 025, Res Div, IBM Corp, Monterey and Cottle Road

San Jose, Calif 95114

1 Dr. U. Shafrir

Space Sci and Erig'Ctr, Univ of Wis, Madison, Wis 53706

1 Dr. C. W. van Atta

Univ of Calif, La Jolla, Calif 92037 1 Dr. E. R. van Driest

North American Rockwell Corp 350 South Magnolia Ave, Long Beach,

Calif 90802

1 Prof. S. Kaniel

Dept of Math, Northwestern Univ, Evanston, Illinois

1 Dr. P. J. Roache

Div 9325 Sandia Lab,

Albuquerque, N. Mexico 87115

1 Dr. S. L. Slotta

Oregon State Univ, Corvalis, Oregon 97331

Prof. W. W. Wilmarth

Univ of Mich, 1077 East Engin Bldg, Ann Arbor, Mich 48104

Dr. H. C. Kao

Northrop Norair, Orgn No 3713, Zone 31 3901 West Broadway, Hawthorne, Calif 90250

3 Boeing Sci Res Labs

P.O. Box 398, Seattle, Wash

1 Dr. A. Goldberg 1 Dr. Y. H. Pao 1 Dr. E. M. Mormon

2 Univ of Notre Dame

1 Prof. A. A. Szewczyk 1 Prof. S. A. Piacsek

2 West Virginia Univ, Morgantown, West Virginia 1 Prof. J. B. Fanucci 1 Prof. W. Squire 1 Dr. D. Thoman C23 E. Mishawaha Ave. Mishawaha, Indiana Copies

6 M.I.T., Cambridge, Mass

1 Prof. J. G. Charney

1 Prof. C. C. Lin

1 Prof. A. H. Shapiro

1 Prof. H. P. Greenspan

1 Prof. S. A. Orszag

1 Mr. M. Israeli, App I. Math Dept

Univ of Maryland, College Pk, Maryland

1 Prof. J. Weske

1 Prof. A. J. Faller 1 Prof. A. Plotkin 1 Yale Univ

Copies

1 Mr. G. N. Slinn

Battelle Memorial Inst, P.O.B. 999 Richland, Wash 99352

1 Prof. Z. Lavan

Illinois Inst of Tech, Chicago, Illinois 60616 1 Mr. J. T. Obrien

U.S. Naval Civil Engin Lab, Port Hueneme, Calif 93041

1 Prof. H. H. Chiu

Dept of Aero and Astro, Univ Heights, N.Y. New York 10453

Prof. K. E. Torrance

Cornell Univ, Ithaca, New York

1 Prof. M. Z. v. Krzywoblocki

Mich St Univ, East Lansing, Mich

1 Dr. W. E. Langlois

IBM Res Lab, San Jose, Calif

1 Prof. J. Hoppe!

Dept Chemical Engin, New York Univ, New York

tr U.S. GOVERNMEIVT

PRINTING OFFICE: 1970-397-175/71

UNCLASSIFIED

SiVirrit_

- . - DOCUMENT CONTROL DATA - R 8..D

(Security ciassilicetion_o_f_p_tle, body of 'abstract and indexing annotation must be entered when the overall report is classified)

1. ORIGINATING ACTIViTy (corpEratirmittior) ' -

-'Naval Ship Research and Development Center

Washington D.C. 20007

20. REPORT SECURITY CLASSIFICATION

_UNCLASSIFIED

25. GROUP _-.

3. REPORT TITLE

FINITE-DIFFERENCE APPROXIMATIONS OF THE VORTICITY OF LAMINAR FLOWS AT SOLID SURFACES

_

4. DESCRIPTIVE NOTES (Type of report and inclusive dates) 5. AUTHOR(S) (Fitst name, middle initial, last name)

Hans J. Lugt and Yermiyahu Riton D. REPORT DATE

April 1970

--70. TOTAL NO. OF PAGES

33

76. NO. OF REFS

11

ea. CONTRACT OR GRANT NO. b. PROJECT NO. _{SR 003 1001}

Task 11572

C.

a.

Da. ORIGINATORS REPORT NUMBER'S/ Report 3306

Db. OTHER REPORT NO(S) Any other nualbere that may be assigned this report)

10. DISTRIBUTION STATEMENT

This document has been approved for puclic release and sale; its distribution is unlimited. .11. SUPPLEMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY

NAVSHIPS

13. ABSTRACT'

Finite-difference methods for solving the Navier-Stokes equations by means of the stream function-vottidity formulation require the computation of the vorticity on solid surfaces with the aid of one-sided difference schemes. Several such schemes of second-order accuracy in the vorticity are derived for axisymmetric flows in orthog-onal coordinate systems. Their accuracies are examined by comparing the appropriate solutions with the exact values of flows past oblate spheroids for vanishing Reynolds.

numbers (OberbeCk solution): It is found that if exact values of the solution are used at the inner points, the various approximations of the surface vorticity differ consid-erably from each other, and that the main source of error is the use of numerical deriv-atives. However, if the various _methods of computing the surface vorticity are incor-porated in an overall numerical scheme, the distinctions among them are reduced

considerably. .

Errors in the drag coefficients computed by the Various methods are an order of magnitude smaller than those of the surface vorticity. Experiments with numerical stability reveal that the diffusion process near the solid surface dominates the numer-ical stability properties of the scheme used.

DD

FORM

473

(PAGE 1)

UNCLASSIFIED

UNCLASSIFIED

Se-C-urity Classification

s.gy WORDS

Finite Difference Approximations Surface Vorticity Laminar Flow LINK A ROLE WT LINK 8 nols WT ROLE LINK C .1" DD ,F,r,..1473 (BACK) UNCLASSIFIED