Computation of viscous transonic flow about a lifting airfoil

J-

r

--'"

_ (NASA-CR-151999) COMPUTATION OF VISCOUS N77-23054 TRANSONIC FLOW ABOUT 1 LIFTING AIRFOIL (Numerical Continuum Mechanics, Inc.) 93 p He A05/MF AO 1 CSCL 0111, Unclas G3/02 28870 )

P

MAY 1977 RECEWEO NASA STI FACILITY ~ INPUT BRANCH ~

.. ~

.. -:~.~y

l--r

NASA CR-151999 NCMR-76-100 COMPUTATION OF VISCOUS TRANSONIC FLOW ABOUT A LIFTING AIRFOIL Prepared by: L. Walitt C.Y. Liu for AMES RESEARCH CENTER National Aeronautics and Space Administration Moffett Field, California 94035

/'

The work reported herein was perfoL~ed under Contract NAS2-9052 November 1976 Numerical Continuum Mechanics, Inc.

6269 Varie1 Avenue, Suite 200 Woodland Hills, California 91367

-r----[-~r--r-l

J-'-

!

ACKNOWLEDGEMENTS The authors are pleased to acknowledge the contributions of Mr. L. S. King of NASA AMES Research Center and Ms. J. A. Enos of Computer Sciences Corporation to the work reported herein. Mr. King developed the original computer program, which was the starting point for this study, derived the characteristic boundary condition equations employed at the downstream boundary, and provided many helpful sug- gestions and discussions in the course of this work. Ms. Enos vecto- ized the computer code employed in this study, generated some of the plots used, and provided programming support throughout this develop- ment effort.

i

J~ .. ·1

I

/ TABLE OF CONTENTS SECTION PAGE ABSTRACT.

. . . . . . . . .

. •. iii · ..

NOMENCLA TITRE. iv . . . . . . . . . . . .- . . . .

INTRODUCTION •.

1. . . . . . . . . . . 2

2. THE STOKES CODE • 6 . . . . . . . . . . . .

· .

2.1 Swi tche d Axis Computation • •

. . . 7

· , .

2.2 Vector Do Loops • . . • • . . • • • . . .

· . .

2.3 Multi-regional Timesteps •..••• · . .

TURBULENCE AND TRANSITION MODELS •.

3. 12 · .

3.1 Turbulence Modelling. • • • .

· .

3.2 Algebraic Turbulence Models • I • • • Transition to Turbulence .•.

3.3 ..

19· 4. BOUNDARY CONDITIONS FOR NON-LIFTING AND LIFTING AIRFOIlS. . • • • . • . • ....

4.1 Non-Lifting Airfoil .

2l · . . .

4.2 Lifting Airfoil • . . . • • . • • . .. .,..

RESULTS OF NON-LIFTING CALCULATION.

. . . 27

.

5.1 Me~h Used and Initial Conditions •••••••••• 27 5.2 Time-Histories of Pressure Coefficients • • • • • • 5.3 Computational Time Reduction Factor ..

· .

5.4 Velocity Vectors of Steady Flow Field • ~ .

5 . .5. Surface Pressure. DistrlbutionCom:pa.risons •• 5.~ Mach Number ContourS: '.... • • • • •. ;. ••••• 6. RESULTS OF.:LIFTING CAWULATION. .• • • • • • • • • • J4 6.1 Finite Difference Mesh ••••..••..

J4 6.2 Initial Conditions ...•.•..•......•.

6.3 Computational Time Reduction Factor . )6 6.4 General Flow Field Structure. . .. . ...

6.5 Calculated Lift and Drag .•.•••..• 40 6.6 Integrated Inviscid/Viscous Calculatio~ 42 6.6.1 General Flow Field Structure. •. • •.

6 .6 .2 Calculated Lift and Drag. • .. .•• • 44 6.6.3 Surface Pressure Distribution Comparison. 46 6.6.4 Comparison of Mach Number Contours ..• 48 '49 6.6.5 Turbulence Characteristics for Airfoil.

CONCLUSIONS AND RECOMMENDATION . . . . . . . . . . . . .

7.

REFERENCES •.•..........

8.

· ..

ABSTRACT Results are presented from a numerical investigation of the; viscous transonic flow about a stationary body in free air. The geometry chosen was a symmetric NACA 64AOIO airfoil at a freestream Mach number of 0.8, a Reynolds number of 4 million based on chord, and angles of attack -.

of O· and 2 degrees. These conditions were such that, at 2 degrees incidence unsteady periodic motion was calculated along the aft portion of the airfoil and in its wake. This unsteady phenomenon has been experimentally observed for other airfoil geometr~sat transonic speeds.

Although no unsteady measurements were made for the NACA 64AOIO air- foil at these flow conditions, interpolated steady measurements of lift, drag, and surface static pressures compared favorably with corresponding computed time-averaged lift, drag, and surface static pressures.

, \ ii1

I NOMENCLATURE a Local Sound Speed

A Surface Area

~

C)r Particle chord length Chord C d~ ~oree Drag coefficient, CD lift force Lift coefficient, CL,

'tc

Pressure coefficient, (P- PI1O)/~ C, Specific heat at constant pressure C, Specific heat arconstant volume ~ Area of airfoil section D E, Specific internal energy Stagnation enthalpy H· we/Vol;) Reduced frequency, ~ (I-M~)IM.od~ Transonic Parameter,

K

J1.

Mixing length Mass Wl Mach number

M

N Cycle number

Pressure A t ~Ul..

Dynamic Pressure t

'&

&. t -~)t.

....

Transformed. radius, (2!" + ~ yo d iv "

R4 Reynolds number ~e R~ynolds number based on momentum thickness Rlll'. Reynolds number based on x R~" Reynolds number based on chord S Entropy -t T1!Jle At T1mestep 1f Temperature U Velocity component in x direction Uoa Freestream Velocity " Veloci ty component in y direction V Volume X Streamwise coordinate Non-dimensional coordinate, x/t.C/z..)

" ~ Normal coordinate ¥. Non-Dimensional coordinate.) , - :;t/(~ ) ...

~ Scaled transonic parameter, c5'i M1 ~ VI eX Angle of attack i Ratio of specific heats

r Circulation

J ~oundary layer thickness J Airfoil thickness to chord ratio,~c € Eddy viscosity ~: Small disturbance parameter, ~~/MOtC)

e Transformed angle.) reu;' (Ii J(-i:::j ~ )

~ Momentum thickness v

'I --I

r

\' ...

~ Molecular viscosity ~ Kinematic viscosity

e Density

r: Shear Stress

..

r Maximum thickness of airfoil

7:" Characteristic Time, t:.~/c.

q> V eloci ty potential functi on

J!..'

~ Transformed velocity potential function d\~ ~v Doublet potential ( )00 Freestream condition ) Non-dimensional coordinate C

"'"'

"'" ( ) Scaled coordinate parameter Derivative with respect to X ( )" ...

-

parameter Derivative with respect to ( ))1 W .......

.....

- with respect to parameter !:t

Derivative (

)i

-

Time-averaged property ( , Wall quanitiy J..v vi

I

~

I SECTION' 1 INTRODUCTION The principal objective of this research effort is to numerically calculate the viscous transonic flow field about a stationary lifting airfoil in free air. To our knowledge, this calculation represents

_

the first-of-a-kind in transonic aerodynamics.

Previous transonic lifting airfoil computations concerned steady inviscid flow ,wi th varying degrees of approximation to the Euler's equations. ,For example, Murman and Cole (Ref. 1), Krupp and Murman (Ref. 2), and more recently Murman, Bailey, and Johnson (Ref. 3) solved the transonic small disturbance equation for flow past thin airfoils including imbedded shock waves. Transonic small disturbance theory solves the small perturbation equation (Ref. 4) with an approximate tangency condition at the airfoil surface, and under the assumption of isentropic flow. Jameson solved the transonic full potential equation (Ref. 5) using an exact tangency condition at the airfoil sur- face, and under the isentropic assumption (Ref. 6). Furthermore, Jameson's method permits computation of.imbedded shocks.

Magnus and Yoshihara (Ref. 7) solved the unsteady form of the Euler equations through an explicit finite difference scheme. A numer- ical solution of Euler's equations, including the possibility of embed- ded shocks (through shock capturing), embodies no other assumptions but the inviscid flow assumption. Coincident with the inviscid assumption, all three methods cited above employ a Kutta condition to determine the magnitude of the circulation for lifting bodies.

For a lifting transonic airfoil with little or no separation in the neighborhood of the shock, combination of turbulent boundary layer theory (Refs. 8 and 9) with any of the above inviscid methods yields useful predictions of the shock location, lift coeffiCient, and drag coefficient, provided tha:c a displacement thickness c.:orrection is made to the airfoil shape. However, when the shock wave boundary layer interaction induces considerable flow separation, P'RECEDING PAGETILANK NOT FItME»

/ '-./ integrated inviscid-boundary layer theories are of dubious value in prediction of the shock location, lift, and drag of the system.

Liepman and Roshko (Ref. 10) point out that transonic shock-wave boundary layer interactions involve an interplay between strength and posi tion of shock waves on the body and boundary layer: character"~ i. e. I separation. The complete e~uations of motion are thus needed to properly

-

account for these nonlinear interactions.

-

Numerical solution of the time-dependent Navier-Stokes e~uations, including turbulence and transition via models, yields predictions of the shock locatioD"lift coefficient,and drag coefficient of an air- foil under conditions of considerable flow separation.

Thus, for thick transonic airfoils at moderate angles of attack, or for thin transonic airfoils at high angles of attack, numerical solution to the time-dependent Navier-Stokes e~uations is the only means for flow field prediction.

A time-dependent Navier-Stokes calculation does not invoke a Kutta condition at the airfoil trailing edge, and includes viscous effects; hence, the possibility of vortex shedding, with attendent periodic flow is within the scope of this computation. For example, !~cDevitt, Levy, and Deiwert (Ref. "11) and Finke (Ref. 12) measured unsteady periodic flow about transonic airfoils for a limited range of Mach numbers and angles-of-attack. Finke's experiments on a symmetric NACA 63 -012 airfoil at Mach .70~Beynolds number 1.25 x 10 based on chord, and two degrees incidence, show periodic flow at the airfoil trailing edge with a reduced fre~uency* K of 2.5 for t0~ oscillating shock. The amplitide of shock oscillation is about one percent"of the chord at two degrees angle of attack and increases to 10% of the chord at eight degrees incidence.

* The reduced fre~uency K is defined as the product of the angular

fre~uency and the chord divided by the free stream velocity.

(

r-

J

/ Finke (Ref. 12) also took interferogram pictures over one full pexiod of shock oscillation for a quasi-elliptical NRL-profile, at a Mach number of .71, a Reynolds number based on chord of 1 x 10 , and at 5 degrees angle of attack. These interferogram pictures, which represent lines of constant density, are shown in Figure 1. The pictures are depicted in intervals of 700~s. In picture 3, the shock - is positioned at 40 percent chord, in picture 5 the shock degenerates to a Mach wave and leading-edge sepaxation determines the flow. Finally, picture 8 is similar to picture 2 indicating periodic flow. Shedding of vorticity is also indicated in the pictures of Figure 1.

The geometry for thic investiGation was a symmetTic NACA 64AOIO airfoil at a freestream Mach number of .80, a Reynolds number of 4 x 10 based on chord, and at angles of attack of 0 and '2 degrees. At 2 degrees incidence unsteady periodic motion was calculated along the aft portion of the airfoil and in its wake. Although no unsteady measurements were made for the NACA 64AOIO airfoil at these flow condi tions, interpolated steady measurements of lift', drag, and surface static pressures compared favorably with corresponding computed time-averaged lift, drag, and surface static pressures.

In order to solve for the unsteady viscous transonic flow field about a lifting body in free air, research was conducted in three principal areas. The are as follows: 1. A computer code, called "STOKES", was developed which solved the Navier-Stokes equations with Multi-Regional Timesteps (lffiT) and Vector-Do-Loops (VDL). The MRT/VDL logic signifi- ca~tly reduced computational time.

2. Turbulence and transition models were incorporated into the STOKES compute~ code.

3. Boundary conditions were generated along the perimeter of the region of calculation which simulated the lifting airfoil immersed in a free airstream.

A discussion of the STOKES computer code is presented in Section 2, the turbulence and transition models embodied in the STOKES code are described in Section 3, boundary conditions are presented in Section 4, while results of the non-lifting (ol. = 0) and lifting (oi = 1) airfoil calculations are discussed in Sections 5 and 6, respectively. Section 7 presents the conclusions reached in this re:search effort.

(

SECTION 2 THE STOKES CODE A computer code has been developed, called "STOKES", for calcu- lating viscous, compressible, time-dependent flow fields about two- dimensional aerodynamic bodies. This code solves the Navier-Stokes e~uations including turbulence, transition, and free air boundary conditions along the perimeter of the domain of calculation. The ..

I finite difference e~uations embodied in the STOKES code were originally developed by Trulio (Ref. 13). These finite difference analogs of the e~uations of motion are such that their self consistency property is maintained. That is, the finite difference e~uations for continuity, momenta, and internal energy imply an exact finite difference e~uation for total energy. The Trulio finite difference relations are second order aCGurate in space and first order accurate in time according to Taylor's series analysis.

Turbulence and transition are discussed in Section 3, while the free air boundary conditions are described in Section 4. This Section is concerned with the numerical method employed to solve the finite difference analogs of the e~uations of motion embodied in the STOKES code.

The principal innovation embodied in this computer code concerns computational time reduction. A method has been developed which produced computational time reduction factors between five and ten. For a non- lifting airfoil the flow field can be computed in less than 30 minutes on a CDC 7600 computer, while for a ,lifting airfoil' the calculation re~uires less than two hours.

The STOKES computer code employs three prinCipal time-reduction methods.

1. Writing the program to proceed along ~uasi-streamlines - rather than the conventional processing along- potential-like lines.

( )

I --r- .. '[

r

r

/ 2. Incorporation of special assembly language subroutines to take maximal advantage of the "pillelining" capa- bility of the more advanced comlluters.

). Setting up a procedure whereby the computational domain is automatically divided into several regions in which different timesteps govern numerical stability.

2.1 Switched Computational Axis Let us consider computation of the flow field about a two~dimen sional body whose axis is parallel to the x coordinate axis, and immersed in an airstream moving along the x-axis •• The y-coordinate is considered the normal direction. The computation takes place on a finite difference mesh comprised of the intersection of quasi-streamlines, which are nearly parallel to the x-axis and extend upstream to downstream, and potential-like lines, which are nearly paralled to tne y-axis, and initiate below the body and end above the body.

Due to computer storage limitations, computations were previously conducted along the potential-like lines of the mesh, starting from the line farthest upstream of the body and ending at the line farthest downstream of the body. This mode of computation is referred to as unswitched. Usually the potential-like lines are shorter than the quasi-streamlines; hence, they are comprised of fewer points. With fewer points to a calculational line the do-loops in the program are executed fewer times, and there are correspondingly more interruptions for data transfer between central and peripheral memory banks. Further- more, most potential-like lines of the finite difference mesh pass through the body; hence, special branch point logic is required to account for the body aerodynamics.

r·

I

/ The switched axis computer code performs, cOIll1'uta"cions along the streamline-like lines of the finite difference mesh rather than along the potential-like lines. For each c~lculational line the do":loops in the 'program are' execut'ed ~ gr~eater number of times ar:d'the number of data transfer interruptionS are reduced. Furthermore, the body geometry becomes two special streamlines of the mesh in the !

switched mode ,of.. computa.tion,' Through' most of, the'- 'streamline computational sweeps there is no branching·of .logic. For the above reasons, the switched axis mode of computation is more efficient than the unswitched mode (Ref. 14).

2.2 .The Vector Do-Loop Recently at NASA-Ames Research Center, Dr. H. Lomax started devel- opment of CDC 7600 assembly language subroutines to replace operations done in regular FORTRAN do-loops. These subroutine~ were designated "Vector Do-Lo'ops," 1. e., VDL*. The Vector-Do-Loop takes advantage of the instruction stack of the CDC 7600 computer and employs the concept of pipe ling in the computational sequence. In Reference 14, a simple vector do-loop is explained and an actual example is given showing the speed-up factor between a vector do-loop subroutine and the equivalent FORTRAN program.

The Vector-Do-Loop substitutes computation during the waiting periods that normally. occur in FORTRAN do-loops. The case of the vector routine is a ,tight loop with computations going on at different stages of development, while the :FORTRAN Do-Loop is essentially a series calculation. The Vector-Do-Loop fits into the instruction stack and thus executes quickly. By making optimal use of the CPt! registers, the Vector-Do-Loop can do the work of a FORTRAN-Do-Loop faster and more efficiently.

* The vector do-loops have also been referred to as the. SSPEEDY Codes.

/

2.3 Multi~Regional Timesteps The Multi-Regional-Timestep (11RT) logic divides the domain of computation into a set of concentric regions about the airfoil. In each concentric region a different timestep is employed, which is govel."Tled by the sta bili ty of that particular region. The time step

_

generally increases by a factor of two as one moves from the inner concentric region towards the outer one. Thus, the inner regions are calculated frequently, while the outer regions a~e calculated in- frequently.

The MRT logic embodied in' the STOKES computer code is based on concepts developed by Magnus and Yoshihara (Refs. 7 and 15) and Trulio (Refs. 16 and 17). In solving for the inviscid field about an oscillating airfoil (Ref. 15), Magnus and Yosq~hara employed four con- centric mesh regions about the airfoil. Three are employed in the vicinity of the airfoil nose and the fourth is a coarse cartesian mesh which contains the other three and covers the remainder of the domain of calculation. The three in the no.se region are a fine s~ewed mesh adjacent to the airfoil, a fine cartesian mesh containing the skewed mesh, and a medium cartesian mesh containing the other two. In their explicit differencing scheme the allowable timestep is limited by stability criteria to an amount which is directly proportional to the spatial mesh increments. For each explicit timestep advancing the solution in the coarsest cartesian mesh of the system, 8 are made in the medium cartesian mesh, 64 are made in the ~ine cartesian mesh, and 256 are made in the s?ewed elliptic mesh around the nose. Thus, Magnus and Yoshiha"ra employed £o.ur. concentric regions in their compu- tation of flow about an oscillating airfoil.

The multi-regional timestip formulation employed herein is patterned after that of Magnus and Yoshir:ara.. A set of concentric regions is developed each of which extends from the upstream to the downstream boundary and contains the airfoil. In this formulation, which takes

(

r \

I _ 1

I place in the switched axis computer code , the variables of motion are calculated for a given layer of zones as often as required by considerations of stability of that layer of zones and its immediate neighbors, rather than the former method of updating dependent variables for all zones-as often as required for the least stable zone of the finite difference mesh. Since the finite difference mesh is composed of streamline-like lines finely spaced in the neighborhood of a surface and coarsely spaced as one moves away from the surface, there are many different timestep levels in a typical aerodynamic problem.

Therefore,' a high percentage of the zones, i.e., those in the large timestep levels, will be calculated very infrequently and great savings is computational time should occur.

The finite difference mesh is broken up into a series of regions comprised of a group of streamline-like lines. The regions are determined by stability for the least stable zone along a given streamline. The timestep levels monotonically decrease as one approaches theaero~ dynamic body from below, and the timestep levels monotonically increase as one moves away from the top surface of the aerodynamkbody. In the multi-regional timestep al'proach, one-dimensional variable timestep techniques are employed to sove a two-dimensional problem. Some of the techniques published by TrJllio (Refs.16,/7) are incorporated in the multi-regional timestep approach considered herein. Particularly the concept of having adjacent concentric regions differ in timestep by a factor of two. Trulio has stated that under these conditions the conservation properties of the finite difference equations discussed in the introduction to Section 2 can be preserved after the solution is advanced in the c'oarsest region of the system.

Let us consider an airfoil finite difference mesh with four concentric regions, each differing by a factor of two in timestep, and let At be the minimum timestep, a "microstep". It is required to update the variables of motion on this mesh through the total time

T

/ interval of" 8 At, designated as a macro-ste:p. The finest region will be u:pdated through 8 increments in At, or through 8 "microste:ps", the next finest region will be u:pdated through 4 timeste:ps of 2 bt, the medium regi.on will be updated through 2 timesteps of 4 At and the coarsest region through the macroste:p.

Complexity is minimized relative to the generalized variable

-

timeste:p :procedure of Trulio (Refs. 14 and 16), since all :points along a given streamline belong to the same timeste:p level. Hence I since the streamline spacing controls the timeste:ps in most two-dimensional :problems, the sinl:ple multi-regional timeste:p a:p:proach should yield most of the benefits of a more com:plex general two-dimensional variable timeste:p method. Development of the MRT logic is :presented in much more detail in Reference 18 • • <:

/ SECTION 3 TURBULENCE AND TRANSITION MODELS An algebraic turbulence model, originally formulated by Baldwin and Rose (Ref. 19 ), is selected for the computation of turbulent flow over an airfoil at transonic speed. The model is, in effect, the

_

mixing length theory to which relaxation along a streamline-like trajectory is incorporated. The simple algebraic model has been recently examined by Baldwin and Rose (Ref. 19 ), Deiwert (Ref. 20 ), and others to be useful in shock boundary-layer interaction problems.

All versions of available algebraic models are discussed. A criterion for boundary layer transition is also presented.

3.1 Turbulence Modelling Numerical modelling of turbulence has become quite practical in the past decade with the advancement of high-speed computers. Though a universal model with wide range of applicability is far from reality, there is ample evidence that existing models have served well even in complex situations such as shock-wave boundary-layer interaction. All models of turbulence are supposed to be general in scope, and until recently, cross-comparisons between models (mainly studies done at NASA- Ames', Refs. 19 , 20', and. 21) are few. For transonic flow, there is no definitive .conclusion as to the best turbulence model to employ.

The usage of numerical models naturally bypasses the more funda- mental approach to turbulence studies via statistical theory, which might be at times academically pl~ing but unrealistic in engineering applications. In general, turbulence modelling is divided into two categories: the algebraic models such as mixing length theory, and the transport models which are described by one~or more differential equa- tions governing some quantity like turbulence energy, turbulence vorticity or shearing stress. The original work of Prandtl and its

-I

I

r

/ subsequent extension by Cebeci and Smith (CS) (Refs. 22 and 23), Rose (Ref. 19), etc. are examples of the first class; the classical Kolmogorov model (Ref. 24) (1942) and the Saffman model (1970) (Ref. 25 ) fall into the latter category.! In adopting a transport model, one must solve, in addition to the basic conservation laws, other differential equations from which turbulence stresses are determined. Transport models try to depict the physics of turbUlence transpor,t, generation, dissipation .

and diffusion. In addition, some models (such as Saffman's) show the correct analytical behavior near the wall (as demanded by the law of wall). The predictive capabilities for incompressible boundary layer flows by those models are convincingly established. Turbulent flows in more than two spatial dimensions, including separation, compressi-' bility, rotational effects, and containing boundary layers interacting with shock waves have not been subject to examination by those models*.

In short, the transport models, as promising as they are, have yet to

be " thoroughly tes ted 'by pro blems"more complex than plane-boundary

layer flows.

In view of the existing complexities in the unsteady transonic flow problem, the desired economy in computation, and the added degree of complication in the nonlinear equations, we must see an alternative to the formulation by turbulence model equations. Thp alternative should be. Cible to render a reasonably good description of the turbulent boundary layer development without a 'd~,sproportional amount of c;,)lI1I>ut:;I'\:.:..i..o"::'ll ti1..t3.

' ... ".

Wilcox (Ref. 26), applying Saffman's model, has shown good results

*

in the study of turbulent boundary separat~on and reattachment at moderate (2.96) Mach number.

/ 3.2 Algebraic Turbulence Models The mixing length theory, originated by Prandtl (Ref. 27), pro- vides the foundation to all algebraic models. Modifications introduced by van Driest (Ref. 28), Cebeci and Sciith (Refs. 22 and 23), and recently Baldwin and Rose (Ref.19), Shang and Hankey (Ref. 29), and Deiwert (Ref. 20) all direct to improve theapplicabili ty of the model. Alge- braic models bypass the necessity of solving additional differential equations. From a computational standpoint, the eddy viscosity based on an algebraic model is post processed from mea~-flow information. Our past application of the CS mixing length theory to internal flow problems in an impeller has shown good qualitative results (Ref. ·30).

Quantitative· comparison is not possible due to the complete lack of experim~ntal data. Hence, some version of an algebraic turbulence model is prefered to the m6re complex transport model. Despite the mixing-length common ingredient, there are variations in each indi- vidual formulation. The variations range from the unmodified theory to a relaxation model incorporating special treatment for the separated regions. The relaxation model was found significantly better than the unmodified algebraic model. According to Shang and Hankey (Ref. 29), it was significantly better than the Saffman's transport model for flow over a flat plate. Since separation on the transonic airfoil is a real possibili ty, incorporation of thtJ relaxation effect becomes qui::'e desirable. For a detailed comparison of various formulations, we list them in the following table.

i I I

t ! I I I

, ...

r :0' ....

I/,.{, { (f-·I).

(f ., , }f,lrt ),,:II-I.r ".M@

s.·, (I-~·rel&.)]

(£.,(IHCl-I)} Ell)' .J f-=--~--- V"..,.,t' X,tIr{(,,·drl \)' .

S_-"J J4"

I

I

I ...

J "I'" .e .

@ ~ &. Ii";' J '.U" ......

I( lL), .... No 1.,./ ) Jr.-~

_

/If.", !'ri,l@ I

a;- ...

,."

,

.. """,l."

~ ...tUt ....

" "J

•.

~.,-· ~ 1,,/ • ••• ... ~ J"itl""J ... •• Pi.".rt , ~U., l: s;.

"t So E. 5ot" t="'I"ltII.~·[!J·(I-e·.{1 ""':'=j[( ......

' t)'" ,.

•• L" 1(~ ....

Ie ~"i ~ ~YJ

f.. t

U,S/

. ••

,I b••

, rOo rill

"r'

t(1.:.'1,t , I , _jf~~.!..), - 7 t:M

..el"

.... I.' , ' E:l.~H;·fF;;)~m.

,.~~.·tI_'"

("

,,4

i"'..... I<:.,wot.!

I ~t

7 U

('II!

.'''.) i , C) (r .... .. - ./.

;f ,..

8o.Jt oj ,,1 It

I tlo

..

"",'7 .... ." n'''';J.• :PiI!"",t- Mo~.' J"" ..

s.~t ~"win f'.~il;ftNI !I,=v'.I"" ~l ~.=O.""7,(J"-!J.).

~

IIJ..J I., r'"

•

/

/.)yr ••• ' • • .. i..

''''' i" 1L

Y U",),,)'

8.r.{,.·. •• &..1.,/1,..

r. \(.-.£.) ,.

tit.

~

U.... \'-

•• .-..

i.

..

.41 +(f'-~I)e""'-

1~ ~i. ~ '" _u, ..

...

,.

' ..

£ 1'.-=.

S~, .... t -I"r.

E.",t/!j

i/.,,-8.f.i,

,-_._- .\-

:"

, E"'~

.J..e f.:: __ WI,..

I

• ...

].

1Ute -

.... ... r

lJ ("", '."";'..

< ..

· ...

+/f/

~{~!'.t.. r t.)

_____ 8.1 1 .. ...

l' . . (t-.,,(-olJ.'1~"'~J il((<)., .~.::::-:---"': ="',f.

(€ "" /i•.,- r--"-·I::.:'~f.lJC' ....

.... I/ ...

M.·" r

i"

.1._ €.=i:'./i; _~_-_'_._.I~!, e.

,

"': ~~:7:':~.T'.)

S,.",

[".,,. 'ff.>=(fE-)'''((ft)'f."(ff

I __

., -'- <r'" ~YI ili f.)

/t//;

o,.,r

(li, .,

I~I{· =

, 0:'" /Jo .

"Mt.

....

.f""" 1- A."J ( s;.

_____ e.= .J.re "..,) I ...1...

]

'J

_ J

)'J ..

....

"-1" .

r' /,..

,l..) 7177,.. •

......

.•

f·

·

!L -r:!$itu" 1.

&._~_/u._e_+--E.-~-c-IC-J,-'.-r-_11_&_.

r 0.0118

",(-l)

'1 % __ ....

1/,'1> w (/)'

•. ____

..

I • ND

=

/t

!--i!'

,,"''7' (1.( .

i I:. 1.. P""/-

':26):,//f

., n.

y= th,;[I(/-.f1 A t.", wlwtt

S-'",i}" [Z~;_~[~~:r

•.

1c-:-:r~K 1~~~~".JiI-t-~~_

-JI_N...:"_~:..·_a._I"'_i.

).t, f.)

r

___ .

I 'l!J. .....

";or u (f.-.

.,. ~·/r

au iF ...

"M"" I-e~d-J) 0.01"

f 5od-/-l

.... (1_!i NO = = Mi., U

I~

____ '/ • :: \ I.

'/=I,·,'./J' I. ])=: CoS P..

EO ___

=y+('=)(:;-+/!J It''I1'!:!'

• r=

.. = ,."=1

.. -A'1t f:.=~.

-,....,.'

E wAfI~ .. A= ~, ..

1;\;.

1-' q:.

h.:'i (:i;: ...

t.'J

oc I:rj

"t:I ot=:.:; O~I ~F; re!j." /[yi ~ ~~ ~r,,,~ /0. t"10 Inn.,. O"~r Jo<i 1?el.~.t,'."

......

\J\ ~; ,,' r ..

r

I The formulation we shall adopt in the transonic flow problem is \ _ basically a hybrid relation primarily based on Rose's relaxation model and the modification suggested by Deiwert. Ingredients of the present algebraic model are blocked in heavy-lined rectangles in the preceeding table. The model for the turbulence stress "'C' ..

~J can be summarized as follows.

( I ) _ ..L ".. •.

with .. - ~ " ......

( .21

The eddy viscosity E is estimated by the mixing length theory which-subdivides the shear layer into an inner and an outer region.

Inner region t 2

(~)"+(~)

e:

=

I I 1) %J lJ%,;

= klyD

tI 0.4

=

kl normal distance from the nearest wall

=

Y

D = 1- exp (-y/A)

A = 26v / J l--rwJ /p i

w

v = kinematic viscosity coefficient at the

w nearest wall

T = shearing stress at the nearest wall

w i;

r

I Outer Region

€o = lmax

= boundary layer thickness Selection

e:: = Min (€I'€O)

eq.

(3)

Relaxation along a streamline-like traj~ctory ~

(4)

fCn = c(5-"~) -t- [~.cP -tlJ-APJ[ i - exp (- "'))

if

cei, < E (~-.6~)

where othenlise And ~ is a parameter defined along a streamline.

-

Two major components, due to ~eiwert., are introduced into Rose's - 2. - 1-

dU,i) + (" ()~) in place formulation. One is the adoption of

( "% . ~%./" J -

of I () 7lZ + d Uj' I

to avoid the complete vanishins of, € in a ";;~I ;;tt,; recirculating zone. Another one is the modification of the relaxation process in which E(~-~~) is used in place of a fixed € evaluated at" some reference station. Moreover, Deiwert found that relaxation over a streamline-like contour was more appropriate particularly for flows over a curved boundary, such as airfoil or turbine blade. Both modifications, indeed minor in nature, are convenient to implement with our computer code in which the scanning is done along streamline-like trajectories. The necessity of incorporating the relaxation effect has been substantiated by )3aldwin and Rose (Ref. 19), Deiwert (Ref. 20), Shang and Hankey (Re'f. 29 ).

Its usefulness for flows in an unsteady transonic flow will be borne out in our forthcoming computation.

PRECEDING PAGE ELANK NOT Flt"MEPJ

1-

. ..... l

Transition to Turbulence Laminar flow at large Reynolds numbers becomes unstable, then the growth of disturbance in the boundary layer builds up until transition to turbulence occurs. The pOint of transition is strongly affected by the streamwise pressure gradient and the turbulence level of the free stream. To account for these factors, several empirical methods are available (for example, van Driest and Blumer, Crabtree, Granville, Smith and Gamberoni, van Ingen, Michel). It is not possible to give a thorQ~gh comparison for those methods. In o\;n,,· blade-to-blad.:.: computation for flows in an i:ilpeller passage (Ref' 30), ooth Granville (Ref' 31)., and M:ichel's (Ref 32) formulations were examined. We found that Michel's simple algorithm provided a clear-cut prediction of transition point and it was extremely if easy to implement. Since- boundary layer transition is such a .( dubious subject in numerical computation, our guideline in the selection of a criterion is again "the simpler the better." Unless future experiments contradict our selection, we shall adhere to Michel's criterion for the present application. The criterion gives a transition Reynolds number, (Re~)t ,based on the local 'Q rans Reynolds number Rex'

.z.z4oa ) 1).#

(5)

(J(eo) == /.17 I -I- ~ - R~,.(

(

t~n~ ~~' The local Re can be estimated from the incompressible momentum

e

thickness e ..

~

(6J

I

If the local Re is larger than (Ree)t transition to e rans, (~ . turbulent flow has taken place. Michel's criterion, resulted from correlation of experimental data, is supposingly valid for the range of Rex between 0.1 x 105 and .60 x 10 • Criterion of this nature signifies that transition to turbulence occurs at a point, rather than in a region, and relaminarization is not

_

possible. A.M.O. Smith had compared Michel's algorithm against Granville's, Smith found the simple criterion of Michel quite satisfactory in the description of transition to turbulence.

ff ' ..

-r

SECTION 4 BOUNDARY CONDITIONS FOR SIMULATION OF FREE AIR FLOW ABOUT AN AIRFOIL This research effort concerns flow about the symmetric NACA 64A010 airfoil, which has a thickness to chord ratio of 10%, at two different angles of attack. At zero incidence the airfoil does not develop lift; hence, the circulation of the system is zero. For the 2 degrees incidence ..

case the transonic field about this airfoil does contain circulation, since lift is being developed about the body. In this section boundary conditions are prescribed which simulate far field free air flow conditions for the non-lifting and lifting cases above.

4.1 Boundary Conditions fo~ Non-Lifting Airfoil Figure 2a shows a schematic of the domain of calculation of a symmetric airfoil having a chord C, a thickness r, and at zero degrees incidence. Let us consider inviscid, isentropic flow about this airfoil. For the case of a thin airfoil, the flow is governed by the transonic small disturbance potential equation (Ref. 4 ).

where: (8)

x=

.....

(9)

(10)

U=

(11)

v:::: (

l

vt Let us define a scaled parameter y such that

. 11 t-

if = J"l M(){) ~

(12) where: " Combination of Equations (7) and (12) yields (13) where:

(1- M';)/ MoO J~

(14) 1(=

- (lS)

4>! / UI10 €:

-

c:P '" ....

f Ii '; -:::: (16)

dJ~ /UoIE

Q>~

(17)

d~MD6

€:::: ;Based on equations (10), (11), (lS), (16), and (17), the streamwise and normal velocity components are, respectively: , (18)

lA =. U (J+ € d> ~ )

ao

V =' Uo(J/~) 4J~

For the non-lifting symmetric airfoil of Figure 2a the effective far field potential solution to Equation (13) is a q.oub1et.

<t>;, == Zlr ~ f ~ .. ! K!j." j! (20)

I where D is the area of the airfoil section. Equation (20) indicates that the doublet potential strength diminishes inversely with the square of the transformed radius ~ = (x + Ky2)"~ frOlll the origin of ,..

the coordinate system, which is located at the 50 percent chord station (Figure 2a). Thus, the effects of the non-lifting airfoil on the perimeter of the system are small provided this perimeter is selected far from the airfoil.

As a result of the above, the following boundary conditions were imposed on the perimeter of the domain of calculation shown in Figure 2a.

(1) Along the upstream boundary the Murman-Cole inviscid solution (Ref. 1) is prescribed for all time.

(2) Frictionless flow is imposed along the lower horizontal lateral boundary of Figure 2a.

(3) Frictionless flow is imposed on the dividing streamline upstream and downstream of the airfoil surface.: (4) No slip flow is imposed along the airfoil surface.

(5) A two-dimensional unsteady method of characteristics, including dissipation, is employed to determine the velocity field exiting from the downstream boundary (Ref. 14).

At the upstream boundary the Murman-Cole solution is essentially the doublet of Equation (20). The lateral boundary is far enough away from the airfoil, so that any reasonable boundary condition will work including frictionless flow. Finally, at the downstream boundary the viscous equations of motion are solved, via the method of charact- eristics, to provide a very realistic model of the exiting flow.

4.2 Boundary Conditions for Lifting Airfoil A lifting airfoil develops circulation; hence, the boundary condi-t,ions enforced on the perimeter of the ,system greatly influence the flow in the neighborhoo of the body. t.et us consider the airfoil shown in Figure (2b) ,at an angle of attack ot, and having the circu-

lation r. The coordinate system (x, ) is located at the 50 percent

/ I~ cho:rd station, the, airfoil thickness is 7:'., -and its cho:rd is C.

Krupp and Murman (Ref. 2) and Small (Ref. 33) have determined the far field potential solution to Equation (12) for the Iifting case.

The resultant equation is as follows: d>': ~' -J! e t rt-(H) ~ ~)(uc 9 - r«-HJ Co.(3~ + ...

( 21)

D 2l1' I6T~ f (~Tfl'K. f

-

wher .i' ~ is the (20) and, doublet potential given by Equation

e = ta~-I(~ k: lilt)

(22) -r ~ {t+ 't(. ,,).1;.

(23) In Figure (2b) the point P is shown 'on the upper lateral boundary ....

haYing the scaled radius r and the angle Q. The second term of Equation (21) ia directly dependent on the angle Q and the circulation

r ; hence, the effects of airfoil circulation are felt throughout the

domain of calculation for the lifting case. Therefore, circulation must be properly accounted for far from the body to simUlate free air boundary conditions.

In order to preserve circulation throughout the domain of calcu- lation, two sets of boundary conditions were investigated in this study.

In the remainder of this section these boundary condition sets are described and their value ascertained.

In the first set of boundary conditions, the velocity field on

the upstream and lateral boundaries of the system shown in Figure (2b) was calculated at each timestep from the far field small disturbance Equations (18), (19), (20), and (21). The method of' characteristics, which included disstpation, was utilized to compute velocities at the downstream boundary (Ref. 14).

To obtain a more accurate velocityf'ield on the upstream and lateral boundaries of the system, the far field small disturbance equations

/ (18), (19), (20), (21) were solved in conjunction with the TSFOIL (Ref. 3) inviscid fi~ld, which served as the initial conditions for the lifting calculation. The TSFOIL. velocity fi~ld was used to numerically evaluate the time-independent doublet contribution to the far field potential given by Equation (21). The numerically determined doublet contribution to the far field velocities, which resulted from the initial' conditions, was then saved and employed in Equation (21) at later times.

At each timestep of the caloulation the lift force was computed about the airfoil and the local lift coefficient Ci was determined.

Based on the local lift coefficient the circulation about the airfoil, which was preserved along the perimeter of the system, was calculated from the following equation.

(24) Equation (24) is consistent with the small disturbance theoryapprox- imations.

The boundary conditions defined in this first set did not work.

Vortices shed from the airfoil, whose rotation produced local velocities in the neighborhood of the 'upper lateral boundary (Figure (2b) ) opposite in sense to these prescribed at the lateral boundary from Equation (21). Material could not flow out of the system; therefore, the pressure built up along the upper la.teral boundary eventually wiping out the wake of the airfoil.

As a result ,of the above, a second set of boundary conditions was prescribed in the far field. The velocity field on the upstream boundary was computed at. each timestep from the far field small dist- urbance Equations (18), (19), (20), and (21) in a manner identical to that described above. However, along the lateral boundaries of the system (Figure 2b), far field small disturbance theory, via Equations (18), (19), (20), and (21), and the TSFOIL solution, were used

1-"-"

1-

_.

I l r

/ to compute only the external static pressure field along these boundaries. Based on the external far field theory pressures and the interior stresses along the lateral boundaries, the velocity field was computed from the equations of motion. Finally, the method of charact- eristiCS, described previously, was employed to compute velocities at the downstream boundary.

The boundary conditions defined in the second set worked. During the shedding process the pressure boundary condition permitted material to flow out of the system at the upper lateral boundary. This per- mitted the shed vortices to pass through the system and resulted in ?-' reas()nable descp:pti9Il;of .. the flQw: field.for the wake .of the airfoil.

T

/ SECTION CALCULA'L':rON :OF NON-LIFTING RESULTS: 64AOIO airfoil the NACA for first conducted Computations were was number the Reynolds Although of attack.

angle zero degrees at the turbulence made with were originally computations the 4 x 106 , model with the in progress is now computation off; a turned model ..

was designated problem The laminar mesh.

with a finer on and turned 0" 101.

"Problem are Problem 101.0 for results numerical airfoil The non-lifting format: the following in presented Conditions used and Initial (1) Mesh Coefficients of Pressure Time-Histories (2) Factor Time Reduction (3) Computational Field Flow of Steady Vectors (4) Velocity Comparisons Distribution Pressure (5) Surface Contours Number (6) Mach Conditions and Initial Mesh Used 5.1 3.

Figure in is presented Problem 101.0 to solve The mesh employed of 3~quasi-streamlines intersection of the is comprised mesh This spacing is The streamline-like-line potential-like-lines.

and 130 layer downstream boundary in the five points about to provide designed shock.

airfoil of the leeward and were. investigated Figure 3 mesh of of the levels The timestep con- run and a timestep multi-regional factor between a speed-up the were conditions freestream Uniform determined.

run was timestep stant is the ratio factor The speed-up investigation.

numerical in the assumed time the total (cover a macrostep to do would take the time it of time- the multi-regional without region) the coarsest of increment with the same computation to do the would take time it to the steps speed-up and their timesteps Multi-regional timesteps.

multi-regional minimum that the was found 14. It Reference in are described factors

r

/ factor speed-up and the 3 was 1.33~ Figure the mesh of for timestep spacing the ~otential-like-line that it was found 2.15. Furthermore, was the time- controlled the airfoil edge of leading of the the vicinity in levels.

step the having investigated, mesh was difference finite A previous spacing potential-like-line with finer but mesh points, of same number ..

finer potential-like- With edge.

airfoil leading of the vicinity in the timestep a minimum mesh had difference finite this line spacing is this mesh 1.69. Furthermore, factor of speed-up and a of .684~s in the spacing the potential-like-line limited by severely still streamline-like- than the rather edge, airfoil leading of the vicinity spacing.

line even factor the speed-up to increase warranted study is Further edge.

leading the airfoil near potential-like-lines by spreading further was generated 101.0 of Problem the solution field for initial The (Ref. 3).

from TSFOIL solution small disturbance inviscid from the interpolated method was from Murman's generated (i initial field The \ written 3 by a computer code of Figure mesh difference the finite onto purpose.

for this especially a to interpolate developed "INTER", was called code, A computer the mesh, termed finite difference on a generalized flow field mesh, difference finite generalized onto another points" "unprimed the interpol- events in of The sequence points".

the "primed termed as follows: are ation process contains that unprimed zone in the is located primed point 1. A it.

is employed expansion series double Taylor's order 2. A first unprimed of the corners at the four the data to interpolate primed point.

onto the quadrilateral internal specific density, and components, The velocity 3.

interpolated.

energy are (

J-

j

r r

/ The INTER computer code was employed to interpolate the inviscid flow field onto the 64AOIO airfoil mesh of Figure 3. The inviscid field developed by TSFOIL for the 64AOIO airfoil was generated for zero angle of attack at llJach .80. The inviscid field is developed on a rectangular mesh comprised of 50 horizontal lines and 89 vertical lines and considers the airfoil as a flat plate through use of the ..

approximate tangency condition. The airfoil mesh (see Figure 3) is comprised of 34 streamline-like-lines and 130 potential-like-lines.

The interpolated pressure field on the 64AOIO airfoil surface is presented in Figure 4. Measured pressure coefficients, (Gross and 6 ,

Steinle) Ref 34) at Reynolds number 4 x 10 are also included for

comparison purposes. It is seen from Figure 4 that the inviscid pres- sure distribution is a good approximation to these data, except in the vicinity of the shock-wave-boundary-layer interaction, and at the trailing edge of the airfoil.

5.2 Time-Histories of Pressure Coefficients Starting from an initial field generateo. from the inviscid small 'disturoance ) solution . :from'TSFOIL, the zero incidence case was run 1100 computational cycles to a characteristic time* ~ OI 3.62. The 1100 cycles required 28 minutes on the CDC 7600 computer.

Pressure-time histories indicated that the field was near steady- state at this characteristic time.

In order to demonstrate that the airfoil flow field is approaching a steady-state, time histories of the pressure coefficient were monitored at five points along the surface of the airfoil. The pressure coeffi- cient is plotted as a function of the characteristic time 1: in Figure 5. It is seen Irom Figure 5 that at a characteristic time z: * Unit characteristic time corresponds to the time it takes a free stream particle to travel one chord length.

r

/ [ ... , of 3.62 the time-histories are nearly horizontal. Although time- histories were not generated in other regions of the flow field, it

is believed that at ,= 3.62 the remainder of. the field is also near

steady-state.

5.3 Computational Time Reduction Factor

_

The STOKES Computer code employed to solve Problem 101.0 embodies three computational time reduction methods, namely: 1. Switched axis 2. Vector-Do-Loops 3. Multi-Regiona·l Timesteps It is the purpose of this section to determine the speed-up factor afforded by each method, and the total speed-up factor.

Let us consider the multi-regional timestep levels (MRT) first.

Table I presents a tabulation of the timestep levels associated with the streamline-like-lines of Figure 3 at a characteristic time ~ of .6600.

TABLE I Tabulation of Timestep Levels at a Characteristic Time r of .6606 Timestep Timestep Number of X-Line Range Level us MicJ:'osteps per Macrostep 1 8 .775655 30-33 2 4 1.55131 25-29 3.10262 2 16-24 4 6.20524 1 1-15 There are four timestep levels in Table Ii hence, the computation requires eight microsteps per macrostep of computation. Each streaw~ine-

r

l

/ the to corresponds K, where K=l the index with is numbered like-line four I that from Table is seen 3. It of Figure boundary lower lateral are K-lines 2, nine in Level K-lines are five in Levell, K-lines are speed- basis a 4. On this in Level are fifteen K-lines 3, and in Level number ratio of the is the which be calculated, can factor, SMRT' up the compu- gove:rns timestep the smallest occur if that of point-steps to MRT.

occur due actually that of point-steps number to the tation (33)(130)(8) :: 3.Jo.,58 (2S) '-')(""""2"'-) -""+--'l--'S""-) -="5 ),...".( ':"""'4 )~+ ""'(~9 4:'""'1i)""( 8:::""1)~+ -r( = ~1~3~O .,-( -r( SMRT be could maximum that not the is computed above factor The speed-up higher limit the 1be pot~ntial-like-lines 101.0.

for Problem achieved in the dropper potential-line mesh and a A revised levels.

timestep of two.

another factor give code should STOKES and (Ref. 14)

of 1.S

a factor gives axes of STOKES switched The

1.S

factor of an additional provides the Vector-Do-Loops employed factor is speed-up the total Therefore, 14).

(Ref.

(26) 7.00 (3.1058) ~ (1.S) S = (loS) run through code was STOKES S.2,the in Section discussed As was The 1100 at r=3.62.

near steady-state to reach a 1100 macrocycles Based 7600 computer.

Ames CDC on the 28 minutes required macrocycles of STOKES version timestep constant an unswitched (26), on Equation the minutes on fifteen hours and in three 101.0 solve Problem would computer.

CDC 7600 Field Flow of Steady Vectors

s.4 Velocity

flow of the steady plots vector velocity and 7 present Figures 6 to proportional are The vectors airfoil.

NACA 64AOIO about the field J.

of Figure positions the mesh and emanate from speeds the particle in the the field field and edge flow leading the Figure 6 shows

{ )

r

r vicinity of the shock wave, while Figure 7 shows the trailing edge flow field.

A deceleration along the dividing streamline of the airfoil is shown in the flow field of Figure 6. Furthermore, subsequent expansion about the forward portion of the airfoil is also indicated in Fig~e 6.

_

The boundary layer must be thin on the forward portion of the airfoil since it is not descernible in the vector plot.

Aft of the shock wave, which occurs at x=.74 ft, a thickening of the boundary layer is indicated in Figure 6. A boundary layer flow can be clearly seen at the axial station x=l.OI ft.

A recirculation region is clearly indicated in the neighborhood of the airfoil trailing edge in Figure 7. This recirculation region is only calculable through the Navier-Stokes equations; such a region is not within the scope of inviscid methods which employ a Kutta condition.

, The calculated surface pressures are compared to measurements in the next section.- 5.5 Surface Pressure Distribution Comparisons Surface pressure distributions obtained from the STOKES numerical solution, the .TSFOIL inviscid solution (Ref. J), and the measurements of Gross and Steinle (Ref. 34) are shown in Figure 8. It Is seen that the STOKES pressure distribution is nearly identical to these data, except in the vicinity of the shock and at the trailing edge. The calculated shock position, i.e., x/c = .475, appear correct; however, the shock transition is smeared out relative to these data; it is believed that the finer mesh computation now in progress will produce a sharper shock transition. The TSFOIL prediction produces a shock location aft of the experimental value, i.e., x/c = .525, but transition is sharper. Finally, at the airfoil trailing edge the TSFOIL prediction .

--T

r--

r

I l shows an increase in pressur 3, while both the STOKES and experimental results indicate no such increase. The flattening of the trailing edge pressure distribution is dUEl to viscous effects. Differences between the calculated trailing edgf3 pressures and corresponding data may be attributed to an early numerical separation caused by the absence of turbulence. It is anticipated that the turbulent computation now in ..

progress will greatly narrow differences between the calculated and measured pressure fields.

5.6 Mach Number Contours

Calculated local Mac:h number contours are shown in Figure 9.

The freestream flow at a Mach number of .80 decelerates to a stagna- tion point at the leading edge of the system, then accelerates past sonic flow to Mach 1.08, shocks down to subsonic conditions, and decelerates subsonically throughout the remainder of the airfoil chord.

Furthermore, a low speed recirculation region exists in the nei'ghbor- hood of the trailing edge of the airfoil.

The gradual transition through the shock wave is indicated in the Mach number contours of Figure 9. The airfoil accelerates to Mach 1. 08 at x = .74 ft, then gradually goes through shock transition to x = .94 ft where subsonic: flow occurs. This. smearing of the shock will be reduced in the finer mesh calculation now in progress.

JJ ;;.L _

I SECTION 6 RESULTS OF LIFTING CALCULATION Computations for the NACA 64AOIO. airfoil at two-degrees angle of attack were conducted with the turbulence and transition models operational. The turbulent lifting airfoil problem was designated "Problem 102.0". During the process of running Problem 102.0, it was found that an integrated inviscid/viscous calculation provided more realistic results t,han a complete viscous calculation. The principal reason for this was that the finite difference mesh employed was too coarse to define the thin boundary at the airfoil leading edge. The integrated inviscid/viscous calculation was designated as "Problem 102.i".

The lifting airfoil numerical results of Problem 102.0 are presented in the following format: 1. Finite Difference Mesh 2. InitialConditions 3. Computational Time Reduction Factor 4. General Flow Field Structure

5. Calculated Lift and Drag

After the results of Problem 102.0 are described a discussion of Problem 102.1 follows.

6.1 Finite Difference Mesh The mesh employed to solve Problem 102.0 is presented in Figure 10.

This mesh is comprised of the intersection of 68 streamline-like-lines and 130 potential-like-lines. The streamline-like-line spacing is designed to provide about six points in the leeward boundary layer aft of the airfoil shock.

The mesh of Figure 10 was developed from the half:mesh employed for the NACA 64AOIO airfoil at zero angle of attack (Figure 3). This half-mesh was reflected to produce a symmetric mesh about the full i: ~

!

asso- airfoil and The of attack.

zero angle at 64A010 airfoil a finite to produce degrees two rotated through then mesh were ciated ...

imposed were then lines Horizontal of attack.

at angle mesh difference the potential-li~e-lines system and of the boundaries lateral.

at the boundaries.

at these to terminate forced were boun-- turbulent edge the leading the airfoil of the neighborhood In of the first thickness than the smaller is much thickness dary layer mesh the Therefore, surface.

the airfoil to zones adjacent layer of to It is designed mesh.

as a medium 10 is considered of Figure to define accuracy with sufficient time running a practical provide phenomena.

mechanical airfoil fluid the important Conditions 6.2 Initial was 102.0 of Problem the solution field for initial The from TSFOIL solution small disturbance inviscid from the generated of attack, angle and 2 degrees at M =.8 64AOIO 3) for the (Ref.

Figure 10.

mesh of onto the and interpolated sides and windward the leeward on distributions pressure Initial figure from this is seen 11. It in Figure are shown the airfoil of at the airfoil rapid expansion a very flow undergoes leeward that the subsonic to then is shocked flow, and edge -to supersonic leading

85% of

of approximately the airfoil along at a station conditions chord.

the airfoil to inter- are compared drag coefficients lift and The inviscid Table II.

)4) in (Ref.

and Steinle of Gross measurements polated

J-

rl

r

-

/ TABLE II Comparison of Lift and Drag Coefficients for the NACA 64AOIO Airfoil Method of Generation Drag Coefficient Lift Coefficient C CD L Murman Inviscid .0413 .8709 Solution (Ref. 3 ) Measured Value (Ref. 34) .047 .415 It is seen from Table II that the inviscid theory overestimates the lift and underestimates the drag by a large factor. If the inviscid lift prediction were closer to measurements and the flow were steady, the Murman far field solution could be imposed on the upstream and lateral boundaries of the mesh and the problem run in this way. Due to the large discrepancy in lift and the possibility of unsteady flow, the boundary conditions described in Section 4.2 were employed on these boundaries.

6.3 Computational Time Reduction Factor Starting from the initial field of Section 6.2, the two degree incidence case was run 2615 macrocycles to a characteristic time r of 9.05. The 2615 cycles required 1.88 hours on the CDC 7600 computer at -Ames Research Center. Since Multi-Regional Timestep (MRT) logic is employed, each macrocycle is comprised of many microcycles. There- fore, the 1.88 hours of computational time is divided according to the timestep levels employed in the calculation. The timestep levels employed for Problem 102.0 are presented in Table III.

.

. :

TABLE III Computer Time Division Macrocycles Timestep Levels Computer T.ime (hrs) .85?

1.023 1.88 The results of TableIII indicate that the first 965 macrocycles occurred at four timestep levels, i.e., 8 microsteps per macrostep, and the remaining 1650 macrocycles occurred at three timestep levels, i.e., 4 microsteps permacrostep. Furthermore, the first 965 cycles required/85? pours of CDC ?600 time, while the remaining 1650 cycles took 1.023 hours.

Based on the results of TABLE III the speed-up factor afforded by the MET logic can be evaluated in a manner similar to that of Section 5.3. Through the first 965 macrocycles four timestep levels were employed. The MRT speed-up factor, SMRT' that results can be determined from Table' IV below.

TABLE IV Tabulation Qf Timestep Levels for First 965 Macrocycles Timestep Level layers of Zones in that Level

)7

= 2.392 The switched axis gives a factor of 1.5 (Ref~,,13), and the Vector-Do- Loops provide another factor of 1.5 (Ref. 13); therefore, the total speed-up factor' over the first 965 macrocycles is

~

Through the remaining 1650 cycles three timestep levels were used.

The MRT speed-up factor) 8 can be determined from Table V below.

MRT TABIE V Tabulation of Timestep Levels for Remaining 1650 Macrosteps Timestep Level Layers of Zones in the Level 1 14

3 34

- (67)(4)(ltO~ _

8 - 130 (14)(4) + (19) 2 + 34) - 2.0938

MRT Accounting for the switched axis and Vector-Do-Loops yields a total speed-up factor of The overall speed-up, 8, is the cycle weighted average of 8 and 8 , i.e., 8 = 4.96 Therefore, in the absence of a switched axis code, Vector-Do-Loops, I;;: and MRT, Problem 102.0 would require 9.3 hours on the Ames CDC 7600 computer.

{ ri

J

I

I As was discussed previously, the timestep levels are limited by the potential-like-line spacing of the finite difference mesh in the vicinity of the airfoil leading edge, not the streamline spacing.

Therefore, zone dropper logic, which drops potential-like-lines of the finite difference mesh away from the immediate neighborhood of the airfoil leading edge, would increase the number of timestep levels

,_

from four to five and decrease the number of computations required.

It is believed an additional speed-up factor of two would result; hence, lifting airfoil problems could be solved within an hour.

6.4 General Flow Field Structure In this section velocity vector plots are shown which illustrate the sequence of events as the viscous flow field develops from the inviscid initial conditions. Figures 12 through 17 present the velo- city field of Problem 102.0 at various characteristic times ranging from zero to 8.62.

Figures 12 to 14 indicate upstream movement of the leeward airfoil shock from its initial position to its farthest upstream location.

Figure 12 shows the initial inviscid field with a leeward shock clearly indicated at the 85 percent chord station. At a characteristic time

t of 2.10 (Figure 13), the leeward shock has moved to about the 40

percent chord station with a separated flow trailing the shock. In fact a. clockwise vortex is seen about to shed at the trailing edge.

A large clockwise vortex is shedding from the leeward side in Figure 14 (~=3.75). Furthermore, the leeward shock has moved to the 30 per- cent chord station, which is approximately its farthest upstream location.

Downstream motion of the shock to a near equilibrium position is indicated in the velocity vector plots of Figures 15 through 17.

Figure 15 shows the velocity field at a characteristic time ~ of 5.44.

At this characteristic time the leeward shock is at the 35 percent

/ chord station. The shock has moved to approximately the 50 percent chord station atr=S.27 (Figure 16). Figure 17 at't:=S.62 also indicates a leeward shock position at about the 50 percent chord station.

The main points of Figures 15 through 17 are that (a) the leeward shock seems to arrive at a near-equilibrium position, (b) vortex shedding appears to continue throughout the calculation, and (c) the wake of the airfoil assumes a sinusoidal pattern.

6.5 Calculated Lift and Drag The main result of the velocity vector plots is that the flow at the airfoil trailing edge and in its wake is unsteady. To obtain a more quantitative description of this unsteadiness time-histories of the lift and drag coefficients are examined.

The lift coefficient time history is presented in Figure IS.

Starting from the inviscid lift coefficient of .S709 the lift coeffi- cient generally decreases in a transient way until a characteristic time r. of 7.6 • For characteristic' times greater than 7.6, the lift-time- history.is periodic, with a complete period of oscillation shown in Figure IS. According to inviscid small disturbance theory (Equation 24), the circulation is proportional to lift; hence, the airfoil sheds vorticity in a transient way until ~ of 7.6 and then periodic shedding occurs. From Figure IS, the period of oscillation is about 1.4 in characteristic time units~ A corresponding time history of the drag coefficient is shown in Figure 19. Starting from the inviscid drag coefficient of .0413, the drag coefficient increases to about .OS at a characteristic time rof 3.0, then decreases to about .03 at~=6.0, and finally becomes periodic after a characteristic time of 7.6. The period of oscillation of the drag coefficient appealS to be similar to that of the lift coefficient.

The lift and drag coefficients of Figure 18 and 19, respectively, are time-averaged over the period of oscillation. The time-averaged

T

I lift

coefficient C is .137 and the time-averaged drag coefficient

L

is .038. A comparison of these values with the corresponding

~

l interpolated measurements of Table II indicates that the calculated values are low. The calculated drag coefficient is in the ballpark, but smaller than measured, while the lift coefficient is far less than measured.

The numerical solution was examined to determine the reasons for the low lift and drag predictions. The principal reason found was that the fin:.te difference mesh in the vicinity of the leading edge of the airfoil was too coarse to define the thin boundary layer there.

The thin leading edge boundary layer is contained within the first layer of zones adjacent to the airfoil on both its leeward and wind- ward sides. As 9. result the hiGh suction pressures near the i1050 we::e not obt:.ined numerically and the supersonic rer:ion was much slT.alle::.:' t.han that observed experimentally.

A thin leading edge boundary layer, whose thickness is smaller than the width of the first layer of zones of the finite difference mesh, has the following two numerical effects: 1. Calculated airfoil leading edge shear stresses are smaller than actual shear stresses; therefore, the predicted drag is less than the actual drag.

2. Calculated average Mach numbers in the first layer of zones about the airfoil leading edge are less than they should be. These lower Mach numbers permit the periodic flow aft of the leeward shock to erode away the expansion region on the leeside.

An alternate way of explaining the smaller calculated supersonic region is to introduce the concept of displacement thickness*. A coarse finite difference mesh at the airfoil leading edge is similar to adding a large displacement thickness to the airfoil which inhibits the leeward expansion.

* The authors are indebted to Dr. Gary T. Chapman of NASA Ames for

suggesting this concept.

I ,I - The calculated flow field contains a leeward supersonic region and a shock wave; however, the magnitude of the pressure drop du.e to the leeward expansion, and the shock strength are reduced below the actual expansion pressure drop and shock strength.

There are two methods of correcting this numerical problem.

1. Re-calculate the flow field with a finer mesh on the ..

forward portion of the airfoil solving the full Navier-Stokes equations.

2. Employ another technique, which does not require a very fine mesh, to solve for the leeward supersonic region, and integrate this solution with the Navier- Stokes solution everywhere else.

It is believed that a fine mesh Navier-Stokes solution is imprac- tical at this particular time. The finite difference mesh required is about ten times finer than what we are now using. This means that the problem will take about 20 hours to solve on the CDC 7600 computer.

This computer expenditure makes Navier-Stokes computations about lifting transonic a.i:cfoils too costly for most applications.

In this research program a numerical experiment has been conducted to investigate the results of integrating another solution with the medium mesh Navier-Stokes solution of Problem 102.0. The forthcoming section describes ,this numerical experiment.

6.6 Integrated Inviscid/Viscous Calculation The inviscid steady field from TSFOIL, employed as the initial conditions, has a much more accurate supersonic region than the medium mesh periodic Navier-Stokes solution of Problem 102.0.

Furthermore, the leeward steady supersonic regi.on compub3d in Problem 102.0 during the early stages of motion, i.e., at a characteristic time of 1. 6, is a more accurate prediction than the TSFOIL field.

- f ---r

This is because the airfoil sr~pe is properly accounted for and some viscous effects are included. The TSFOIL solution, which employs the approximate tangency condition, is interpolated to account for the body shape (see Section 5.2).

The calculated leeward supersonic field of Problem 102.0 at a characteristic time of 1.6 was integrated with the Navier-Stokes

_

solution everywhere else in two steps.

1. At characteristic times between 7.6 and 9.05 in its period of oscillation, the periodic flow field computed in Problem 102.0 was patched together with the steady leeward supersonic region computed atr =1.6.

2. The interface region between these solutions, which contained the leeward shock, was recomputed by running the STOKES code through a characteristic time interval of .025.

After the STOKES run, an integrated solution resulted ,,,hich had a proper supersonic region, a leeward shock wave- Joundary layer inter- action, and a separated trailing edge region. l'he int,egrated inviscid/ viscous calculation was designated as Problem 102.1.

In the remainder of this section the results of Problem 102.1 are presented. The format for presentation is as follows: 1. General Flow Field Structure 2. Calculated Lift and Drag 3. Surface Pressure Distri butlon Comparisons 4. Comparison of Mach Number Contours 5. Turbulence Characteristic for Airfoil 6.6.1 General Flow Field Structure Shedding of vorticity from the lee side of the airfoil is indi- cated in the velocity vector plots shown in Figures 20. 21, 22, and 23.

Figure 20 presents the velocity field at ""t:' =7.831. A leeward vortex is forming at approximately the 50 percent chord station, while a

I a vortex rotating counter-clockwise is shedding at the airfoil trailing edge. The counter-clockwise shed vortex induces clockwise circulation about the airfoil; thus, lift is enhanced. At a characteristic time of 8.204 (Figure 21), the eye of the leeward vortex of Figure 20 has,-moved to about the 65 percent chord station and the counter-clockwise shed vortex has moved further downstream. Figure 22, at a characteristic

~

time of 8.649, shows the leeward vortex of Figure 21 at the airfoil trailing edge and in the process of shedding. Furthermore, the sense of' this vortex is now clockwise. This shed clockwise vortex induces a counter-clockwise circulation about the airfoil; hence, the lift is reduced. At a characteristic time of 9.052, the velocity field shown in, Figure 23 no longer has the clockwise leeward vortex in it; however, a counter-clockwise vortex appears to be shedding from tbe trailing edge. In fact, the velocity fields of Figures 20 and 23 are somewhat similar, suggesting a periodic flow.

The velocity vector plots of Figures 20 to 23 are similar to the interferogram pictures taken by Finke (Ref. 12) and shown in Figure 1.

Shedding of vortiCity is clearly indicated in both the interferogram pictures and velocity vector plots as well as a sinusoidal wake pattern.

Finally, the alternate shedding of counter-clockwise and clockwise rotating vortices in Problem 102.1 indicates that the instantaneous lift coefficient should fluctuate with time. Lift and drag are presented in the next section.

6.6.2 Calculated Lift and Drag The periodic nature o~ this flow is quantitatively confirmed in calculated time-histories of the lift and drag coefficients shown in Figure 24. It is seen from Figure 24 that the instantaneous lift curve has a period of about 1.4 in characteristic time units.

This period converts to a frequency of )66 hertz, or to a reduced frequency* ~ of 4.49. The drag coefficient time history curve

* The reduced frequency K is defined as the product of the angular

frequency and the chord divided by the freestream velocity.

·r

/ .... aD".:" to> exhibit the same ..• frequency. As is shown in the velocity vector plots above, during the period of oscillation of the lift coefficient alternate shedding of counter-clockwiSe. and clockwise vortices occur.

Independent. experiments on a .. symmetric NAC'A 631-012 airfoil at Hach number .70, aeyno1ds number 1.25 x 10 , and two degrees incidenc~ ...

described in Section 1, show a periodic flow at the trailing edge with a reduced frequency K of 2.5 for the OSCillating shock. The dis- crepancy in the reduced frequency between the two cases can be attri- buted to differences in airfoil geometry, Reynolds number and Mach number._ The calculated instantaneous lift and drag coefficients of Figure 24 were time averaged and compared to (a) steady lift and drag measurements of Gross and Steinle (Ref. 34) and (b) lift and drag predictions of small disturbance inviscid theory (Ref. 3). The calculated time-averaged lift and drag coefficients of Problem 102.1 are compared to correspon-~ ding data. in Figure 25. The predicted time-averaged lift coefficient is close to the faired experimental curve while the predicted drag coefficient is lower.

A quantitative comparison of predicted lift and drag with inter- polated data of Gross and Steinle (Ref. 34), and inviscid predictions of TSFOIL are shown in Table VI.

TABLE VI Comparison of Lift and Drag. Coefficients for the NACA 64AOIO Airfoil Drag Coefficient Lift Coefficient Method of Generation C CD L Murman Inviscid Solution .0413 .8709 (Small Disturbance) Measured Value (Gross and' Steinre, . .047 .415 Interpolated, TMX-62468) STOKES Numerical Calculation .034- .384 (Problem 102.1)

f It is seen from Table VI that the calculated lift coefficient is within seven percent of the measured value, while the inviscid prediction is greater by a factor of 2.1. Furthermore, the calculated drag coefficient is somewhat lower than measured, while the inviscid prediction is closer to the, data point.

The lower calculated drag coefficient is a result of the coarse- ..

ness of the finite difference mesh at the leading edge of the airfoil.

The leeward and windward boundary layers on the NACA 6J~A010 airfQil are much thinner than the thickness of the first layer of zones adjacent t9 the airfoil. Thus, the numerically determined shear stress at :.the airfoil leading edge is greatly underestimated, and results in a lower drag coefficient prediction. It is believed that for thin boundary layers at the airfoil leading edge, boundary layer theory must be inte- grated with Navier-Stokes computation everywhere else to produce correct drag coefficients.

6~6.3 Surface Pressure Distribution Comparisons Instantaneous calculated pressure distributions on the leeward and windward sides of the NACA 64AOIO airfoil are shown in Figures 26 and 27, respectively. The curves of these'figures correspond to the characteristic times 7.831, 8.204, 8.64~, and 9.052, which are identi- cal to the characteristic times of the velocity vector plots of Figures 20 to 23. The leeward instantaneous pressure distributions (Figure 26) depict the periodic shedding of vortices when compared to the velocity vector plots of Figures 20 to 23.

Let us examine the leeward pressure distributions of Figure 26.

Upstream of the 40 percent chord station the pressure field is steady.

Curve 1, at a characteristic time of 7.831, has a high pressure point at the 50 percent chord station followed by a rarefaction whose minimum

is at the 58 percent chord station. This correlates with the velocity

vector plot of Figure 20 at 1:" =7 .831. Figure 20 shows a stagnation point at the 50 percent chord station followed by an expansion about the separated

I 26), of Figure (Curve 2 of 8.204 time a characteristic At region.

followed chord station percent at the 60 point is pressure the hi€;h plot of vector the velocity with correlates this by & r.arefaction; 60 percent is at the point the stagnation that which shows Figu~: 21, In region.

separated about the expansion by an followed chord. station chord .50 percent at the has formed region a separated -.

other words, 1:"=8.204.

station at chord 60 percent to the and has moved ,""7.831 stai~ion at moved point has the stagnation shows that ('t"~.649) Figure 26 Curve 3 of pressure rarefaction minimum with the station chom the 80 percent to of Figure plot vector The velocity station.

chord the 90 percent at edge of bailing to the has moved the vortex shows that (r.=8.649) time characteristic 26, at a of Figure Curve 4 Finally, airfoil.

the 'I'he' rarefaction.

followed. by a point no stagnation 9.0.52, has of at r=9.0,52.

has shed.

the vortex that Curve 4 indicate of results of Figure distributions pressure windwam The instantaneous chord For percentage airfoil.

of the edge the trailing near oscillate For side.

on the windward is steady the flow percent than 50 less are indicated.

oscillations .50 percent greater than chord percentage the periodic result of are a oscillations edge trailing These windward of the airfoil.

side on the leeward occurs process that shedding vortex the windward on distributions pressure calculated Instantaneous 27) 26 and (Figures airfoil NACA 64A010 of the leeward sides and of A comparison of oscillation.

the period over time-averaged were measure- steady With·co~sponding pressures time-averaged calculated the theory are small disturbance inviscid from predictions and with ments calculated time-averaged the On balance Figure 28.

in presented data.

these cia ta, although with corresponding well compare pressures predicts theor'J The inviscid attack.

angle of wi th interpolated are the STOKES while chord station 85 percent at the location a shock is at the shock that data indicate and the location shock calculated station.

chord 45 percent the !

.'

\

T-'

'J-'

I , 6.6.4 Comparison of Mach Number Contours Mach number ,contours about the NACA 64A010 airfoil are compared

-

in Figures 29 and 30. Figure 29 presents the TSFOIL inviscid solution, while Fig~~~ ~O represents the STOKES prediction at a charact- eristic time of 8;65. The most striking differences between the two figures a~ the extent of the leeward supersonic region and the locus 01' the leeward shock. The inviscid flow field expands about the leeward side of the airf6il to a supersonic Mach number of 1. 3 and then goes through shock t.:!'allsition a.t the 8; percent chord station. In contrast,

the viscous field at r=8.6; b~s a supersonic expansion to Mach 1.2

anc. goes through shoclt transition at about the 4; percent chord station.

The results of Fi.gures 29 and )0 are similar w;.th respect to the thickness of the shock transition region. Both fields have a leeward shock layer approximately 5 percent of a chord thick. The finite shoc:k layer of Figure 30 is C\ resul t of the STOKES numerical method and the finite difference me'sh employed in the shock region. If shock f1tt.L~g is not employed, shock transition usually takes place across four zones of the finite dlfference mesh. Hence, the spacing of these four zones defines the tldckness of the shock layer. In this particular case the layer is approximately; percent of a chord thick. A finer mesh in the shock region can reduce the thickness of this layer.

The viscous wake of the airfoil is also indicated in Figure JO. In fact, at a. station x of a.bout 2.80 the eye of a yortex can be seen in the ~Ach number contour plot. At approximately this station a closed

l

a indicating .90 contour, closed Mach a contour contains M&ch .80 the vortex.

eye of from the moves away as one in velocity decrease of plot number contour the Mach out that pointed it 1s Finally, time a characteristic field at of the a snapshot represents jO Figure shedding and vortex oscillates ~hock actually 8.65. The leeward of leeward pressure the instantaneous of an examination ...

From occurs.

of the amplitude that the is seen 26, it of Figure distributions is chord. This of a one percent i. e., about sma1~, ..

is oscillation the 1) for Section 12 and of Finke (Ref.

the data with consistent • 1.25 x 10 Mach .70 and Reynolds number airfoil at 63,1-012 NACA a chord percent of about one was that the shock amplitude found Finke Q o as the angle Thus, at ol.. =8.

of a chord to 18% and increased 0<. =2 £01: equi- the leeward shock approaches an reduced towards zero attack is of is achieved.

and steady. flow position librium for Airfoil Characteristics Turbulence 6.6.5 of the leeward shock flow aft of the ~Ature the transient Due to turbulence the Rose of the validity to assess difficult it is wave.

numerically However, investigation.

this numerical in model employed of the lee- upstream and downstream profiles eddy viscosity determined on profiles corresponding with compared qualitatively shock ward layer interaction.

plate undergoing a shock-wave boundary a flat the molecular E to eddy viscosity of the the ratio Profiles of along stations at two chord in Figure are presented p viscosity data correspond These ai:r:f'oil.

NACA 64A010 side of the leeward the hand on the left profile 8.65. The ~ of time a characteristic to is station.

percent chord the 31.4 is at figure, which of the side of the hand side on the right profile while the the shock, of upstream is downstream o£ station, chord 61.1 percent is at the which figure, nearly linear is

the E(}A. profile

the shock of Upstream the shock.

boundary of the turbulent portion and a sublayer the laminar through

~- .r

r

/ layer, constant through the remainder of the boundary layer and then rapidly decays as the local freestream is approached. Downstream of the leeward shock the velocity profile is separated in the vicinity of the wall. The separated velocity profra,a causes the eddy viscosity to approach zero near the wall and results in a nearly linear E ~ profile in the shear layer above the wall. The lineare~ region is~£o11bwed by a constant value throughout the remainder of the boundary layer and then a rapid decay to local freestream conditions. A comparison of the calcu- lated ~ ~ profiles upstream and dmvnstream of the shock indicates that the turbulence intensity increases thr.ough the shock transition.

The qualitative behavior of the calculated STOKES E~ profiles is similar to corresponding E profiles calculated by Baldwin and Rose (Ref. 19) on a flat plate. Baldwin and Rose used a relaxation eddy coef- ficient model to investigate a shock wave-boundary layer interaction on a flat plate at Mach 2.39 and Reynolds number per meter of ;.7 x 10 .

The Baldwin and Rose profiles upstream of and downst~~am of the shock are presented in the upper right corner of Figure 31. It is seen from Figure 31 that the Baldwin and Rose .E profiles upstream of and downstream of the flat plate shock are similar in shape to corresponding NACA 64AOIO ~Ip. profiles upstream of and downstream of the leeward shock.

r

, SIOOTION 7 CONCLUSIONS AND RECOMMENDATION The present investigation has 4emonstrated that the STOKES computer code is an accurate, practical tool for solving vis~ous. compressible lifting and non-lifting airfoil problems at transonic Mach numbers .- and for Reynolds numbers ranging from the laminar to turbulent regimes.

Applicability to the turbulent Reynolds number regime was achieved by employing the Baldwin-Rose mixing length theory to which relaxation , along a streamline-like trajectory is incorporated.

The numerical method embodied in the STOKl!.S code offers a complete.

accurate description of the viscous flow field about a transonic air- foil, including the important fluid mechanical effects of 1. shock wave-boundary layer interactions.

2. transonic buffet and 3. circulation for the'lifting case.

Furthermore, a Kutta condition, whose formulation for unsteady flow I i is quite uncertainJis not required in this method.

Accuracy of the numerical method was demonstrated in numerical- experimental com:p&risons for two transonic airfoil cases. The three fluid mechanical effects defined above were computed and verified experimentally. At zero degrees angle of attack a shock wave-boundary layer interaction was successfully computed and verified for the NACA 64AOIO airfoil. At two degrees incidence periodic motion (trans- onic Duffet) was calculated along the aft portion of the airfoil and in its wake; this periodic motion has been observed for other airfoil geometries at transonic speeds. Although no measurements for the un- steadiness of the flow field were made for the NACA 64A010 airfoil at these conditions, steady measurements of lift. drag and surface I?tatic pressures compared favorably with corresponding computed. time- avemged lift, drag, and surface static pressures. Finally, at two degrees incidence the flow field was computed with far field boundary .51

"!

r

!

conditions that simulated the free air case and preserved the airfoil circulation throughout the flow field.

Practicality of the STOKES computer code was demonstrated by the Ames CDC 7600 computer time requirements for the two cases solved.

For the non-lifting case the steady flow field was computed in 28 minutes, while for the lifting case the periodic flow field was ..

computed in 1.88 hours.

The major cri ticisJI of the present method of solution is perhaps the accuracy with which the shear stresses are computed at the leading edge of the airfoil. Due to the coarseness:of the leadin~~edge finite difference mesh relative to the boundary layer thickness there. the calculated leading edge shear stress distribution and hence, the drag, was lower than measured. Therefore, it is .re':'!ommended that in the neighborhood of the airfoil leading edge, boundary layer theory be integrated with Navier-Stokes computation everywhere else. This will produce correct drag coefficients with a minor increase in computational time.

--

l

J-

·r

SECTION REFERENCES Plane Steady of D., "Calculation Cole, J.

E. M. and Murman, 1.

114-121, No. I, pp Vol. 9, Journal, AIM Flows," Transonic 1971.

January, Flow Transonic of "Computation E. M., and Murman, J. A.

2. Krupp, 10, Journal, Vo1.

AIAA Bodies," and Slender Airfoils Lifting Past 1972.

July pp. 880-886, No.7, -~A "TSP'OIL ~.-L. John:sor,; and F. R. ra.iley, }1. Munnan, 3. :me :rncludinp'- Calculat.ion::., Tram:onic for Two-f.imensional CoiniJUter -Code SP-J47, NASA Evaluation", \~ave-Drag

and ~

Wall Effects vlind-Tunnel Computers, Advanced Re9.uir1ng Analyses Aerodynamic Part II, 1975.

769-788, pp.

p. 205, Gas Dynamics, of A., Elements and Roshko, H. W.

4. Liepmann, 1957.

New York, Wiley, p. 204.

5. Ibid, A.. Lecture Jameson, D., and P., Korn, Gara. bedian, F., 6 . Bower, Wing Si~tems. Sunercritica1 and Mathematical Economics Notes in 1973.

Springer-Verlag, Vol. 108, II, Section Airfoils," Flow Over Transonic Itlnviscid H.

R. and Yoshihara.

Magnus, 7.

1970.

December 2157-2162, 10, 12, PPt Vol.

Journal, AIM - AMES INIET Boundary Guide C. B., "Users Davies, J. D. and Murphy, 8.

1973.

January I," NASA TMX-62-211, MK Layer Program "On the Calcu- t W. C., a.'"ld Rose L. L., D., Presley, J.

9: Murphy.

AGARD Flow," Reattaching and Separating Supersonic lation of Germany.

1975. Gottinge.n, 27-)0, 23. May Paper Proceeding 3. p. 270.

Reference 10;' "Transonic G. S., and Deiwert, L. L., B., Levy, J.

11. McDevitt, 14, Vol.

AlAA Journal, Airfoil," Thick Circular-Arc About a Flow May 1976.

606-613.

No.5, pp.

on Layer Interaction 12. Shock Wave-Boundary IC •• "Unsteady Finke, 168, No.

Proceedings Conference Flow," AGARD Transonic in Profiles Sympoisum, Panel Fluid Dynamics FLOW SEPARATION, AGARD-CF-168, 1975.

27-30 May Germany Gottingen, Air AFTON Codes," or the Structure and J. G.. "Theory Trulio, 13.

1966.

AFWL-TR-66-l9, No.

Report Technical Labon.to~J Force Weapons Compu- of Investigation J •• "An L., Enos, S., Walitt, King, L.

14.

Navier-5tokes the for Solving Methods Reduction tational Time 1975.

102 t October Tl~-73,- t tt NASA Equations

r I

r- ··-r

[

r

/ 26. Wilcox, D. C., "Calculation of Turbulent Boundary-layer Shock-Wave Interaction," AIAA Journal Vol. II, No. 1+, p. 1592, 1973.

27. Prandtl, L., "Bericht uber Untersuchungen zur ausgebildeten Turbulenz," ZAMM, Vol. 5, p. 136, 1925.

28. Van Driest, E. R., "On Turbulent Flow Near a Wall," J. Aerospace Sci., Vol. 23, No. 11, p. 1007, 1956.

I ....

29. Shang, J. S. and Hankey, W. L., Jr., "Numerical Solution of the Navier-Stokes Equations for Supersonic Turbulent Flow over a Compression Corner," AIM Paper 75-4, 1975.

30. Walitt, L., Harp, J. L., Jr., and Liu, C. Y., "Numerical Calculation of the Internal Flow Field in a Centrifugal Compressor Impeller", NASA CR-134984, Dec. 1975.

31. Granville, P. S., "The Calculat.ion of the Viscous Drag" of Bodies of Revolution," Rep. 849, David Taylor Basin, 1953.

32. Michel, R., "Etude de la Transition sur les Profiles d'ailej Establissement d'un Critere de Determination due Point de Transi- tion et Calculde la Trainee de Profile Incompressible," O.N.E.R.A.

Report 1/1578A, July 1951.

33. Small, R. D., TAE Report 273, Technian, !srael, February 1976.

Gross, A. R. and Steinle, F. W. "Pressure Data from 64AOIO 34.

Airfoil at Transonic Speeds in Heavy Gas Media of Ratios of Specific Heats from 1.67 to 1.12", NASA TMX-62-468, August 1975.

eld Figure 1. I nte over o s econds .

G GE

(a ) No n -Lifting Symaetr1c Airfoi downs tre atl >, ~-- upstream boundary (b) Lifti ng S ymIle tric Ai rfoi l lateral ~ ~ry ...

downstream boun da ry '--- ups t reaa boundary l atera l bounda ry 2 . S ch e ti c 0 t flo lif ing nd non-lifting airfoils.

5-.00 4.00 F~Fe Trajljn~ .00 I')tre~mline 3h Aj~foil of l.OO comprincri menh .00 ; .

ft computation X , 0.00 .

il airfo -1.00 F~gn 6l~AOlO for Lc~din~ potential-liT·f'-liltf'~ mc:;h il - l.OO o rf Ai and difference 3.00 te - Fin:i like-lines 3.

4.00 - FiOlre . 00 0.

-1

l,() 116R .9 2 .

ntrlbutlon n ::>O.lutlo 1\ , a 4 in c

c

of fII' R n , 64AOIO .

NACA

=

ompariso Ot; C on

o

4.

gure n

.3

, 2 o

/e .9 116

/e = . )205

- .50 .

...

/e = . 1 1

~ c: :l

. --

(.. t::.h ,. ~ = • \J., _~

xl

- . ...

c: ...

• • • oJ

x./c - . 0C "25

: I.:re 5· ?re St'!'O c ef:':'c: n't.

tine h storie or. :'A ... A 6[' AC_O airfoi_ su r ... ace .

I ...,

'\ 1

..

r- .~ t.- , H c:; c.

-

, c.

~ :;3 C- c 't: \ ,-i c.

Co- ,

-

,....

c.- o t ..

re:, Q

\

~ J . ~

1 >~ r""\

. -\ ~~.

' 0 .,..: 'I ~"6 ~ -t

b)(

, r-

\

:r

-

": :;

\ ~

1 ;:::I

t.J

\

. ~ r=..

"!

:;

\

I

\

co

\

~

I

, , i i - .. ... ..

.. . - ,., ...

..

co co N ..

- ... ...

N co • ...

-

.: ...

. .; ...

co , .. , co co , I I

"

(

/

/

/

o-

f

'-' ~ f-; ..:. ~

"'"

to:; r-: t~ ~

·

t) 1-\

-- · (

$ ....

·

..:..

· ~

ES ?: ; ::"" <: ') C' (

-

( r- ~ ~ .

.: U ~ ..-{ ~- r --" ..: c.:.

( / / ~ r-

.-

~ c .

\

, " :.

..

.- c..

,..

( l C . ~

-

..

;:: . ., r0- ~

t

' ;-1 .

,

--'

.-

r c:

-

l .

., ~ " ...

. -{ • .J (; C ~ .

. -{ I- r;, i ~." .. :J ( .-1 C ....

.

CO r t r , r-l ....

~

p II ."

-

-

'IJ

---

lit

-----

" .

..- I

..

t ..

.

< ...

!

" .

.

~" t .: " .- " .... incidence, ..,..

•.

'::'50 zero .

l--'

at . ·: .

.'

r '

·

~

·

· Ft

,;: ' .,. airfoil - - - - t

~~~~

,.

r '

?"~·

·

64AOIO ' '

" 3·-i~

' S

~·~

NACA

~~~:-----==~'==E'

"S the

;;:~

.

::; .

..

, ~----::;: ..

;;;; X.f about .

,

~

; tI - '.

....;:!

..

x

~~

; ..

~

contou~ 4 t • 7 ..

c---.

~!!

[ ~t:;; number Re~ ~~,, ."

.

, .

Mach ---:;0 .&.

= I.', ft1.tJII Local 9.

_ ~ __ Figure ~_~ s __ _ • m • • _~ ~m u t .

.....

I rJ ~ __ ~. U '

\~

HitCH __ ...

In _____ _ loem a - '" ..

• .

1"'

..... --~ - ----- -, -

<:) <:)

·

(J) o CD

· :1'

o (J) en

·

·

N <:)

-

:t:o

..,0 >- .

o <:) N

·

-

I

· N

I o (J) <:) CD

·

-

<:) <:) w1- ----~;_----_r------r_----_.------._------r-----~------~----~-- ~ 1- .

t\Jo eg ees esh for AC 64 0 0 Figure 1 . ed'urn f'n e d'ffe ence e - 02 . 0- a ng e 0 attack Pr oble j 'nes a nd 130 pote I

I

- Leeward Surface - 1. 2 -. 8 - . 4 o U Win dward Sur f ace III ~.4 Ul CJ) III ~ IJ.4 .8 1.2 1.

6 ()

o

1 °

°

Pe rc e nt a Cho rd Figure 11 . Initia pressure i tr "b ut 'o ns on he leew rd nd "nd ar urface =2 0 Re = of the C 64 010 a'r oil " == . 8 x 10 " , =c

-

.U

airfoill

64A010

MACA

~----~----~----~----~----~----~----_t

the

____ ~ ____ abo~t ~

l d

____ ~ ---------~~~

fie

____ --...

~

.0

flow

____ ~ ____ ~

initial

____

rob1ea

~

P

lL....lo

the

,. 1 ,.

____ ,.

~ ~ > , __ ,.

" X ,.

,.

plet

), ,.

b - , ,.

~ot

,.

, Yee~r ,.

" ,.

> b .80, ,.,.

.,.

Veloe1ty

MO()

12.

Figue

tl----~----~----~----~----~----~----~-- 'OJ 1.U i. 1 .

w

·

0 ~ ....

.... -

~

~

~ Cl -ioJ (.)

~

e

~ (.)

·

.,

1 · ~

-ioJ as .

~ QJC\I ~ o fM r-f a EI ~ r-I ~

t.o-

1 11 ~~ 0)0( r-f o.:t

i

< ,,1 ~ ~.- C\I

!

-ioJ'l' '-1 0- -ioJ CX) .

r-f ~ C\I

~~

~I _'f'4 >-ioJ .

C""'I r-f ~ ~ ..

.. 0

..

~ .; ~ I! ... .. .. ..

T A

'f> 10'0 time characteristic a at .0 flow f1e1d em Pr obl airfoil , x 64AOlg IS MACA ~c~ J ~ot~ - ~ plot d r .8 - Yeetor ~ 3.7.5.

Velocity of .

....

Figtn'e ..

:a o

C

~ ~~ ~~ -< ~ \J~ ~~ .&.

~ .....

o

to

,..;- Pt H • () ~

~~i

o

o

(J r;) .....

c

co

co

~ CIl QI

co

"0 I )

..=

.!: CIl M E E .!: :::l 0 8

;... 0

CIl E .!: .D E ~ r-i :::l CO ' -4 ~

o ·CO e.c c: tOO II <: ~ p::: . ,-1 E-< CO :z -I /I CI) ..c CI)

e

""

..c '0

e

I ~ CI)

e

.s::; .0 S 0 E-< Ii.

~ __________ ~ __ ~ ______ ~ ______ --~r----------,-----------t- 0

o

CO o

o o

r-I o .

OJ '4ual~111aoJ 2~~O

ti

--

tic ~~~~~~~~~~~~~~~-----r ____ field ~ 102.1 flolf ______ blem ~"_ ol1 Pro a1.rf 6, ______ "- 4 x 10 64AQ ...

________

eo<)

"- ...HACA

fR

the " 1 -2 ________ of "- ,0/.

p~ot . 80 go ________ tor "- ~~ vec ____ t1 -" 1 .

f Ve10el ________ -" __ 20.

_"- ul Figure l.

0- -< • :z=

time -~~ ~- ristic characte a t field fiow b1ell l1 o ro irt '!'

- , MACA ~ C>

of ·

t

eX

r • ••• D4 v M - •••• , y h h .

6'5 ' • -clt o 8 .

•••• •• ••

-

tI ~ .

tlae ~---~- _ -~

-

0#

-

--

characteristic a ~~ ..

_

a + ';:IIt.,t".

x 64AOIO 4

-

CA NA

/!R'04

the 2"; of "" 0( plot .

tor - vec MD6 > , a j : ,."

...

..

a), 9.0 ),,, clocit f 'b"" - j," a.

" ..

1:.

1)," %", 19u .1 1."

....

, 06 , 0 , 0

o

8 9 o

C ract r s t c . 6

o

6 7 10 Ch cte istic Time , ~ F' ure 2 • T ' e his ories of the lift n dra co e fficien for the AC 6 AOlO airfoil at two de rees angle a attack' Problem 102.' - integrated inviscid/viscous calculation; M = . 80. Re 4 x 10 , ~

o

2.

with s ient

o

2. 0 c fi ef co drag 1.6 "9d a - 62468) egrees d lift (THX ;: 102 . 1 e o 1.

ged 0( , ra k , R", teinl ac ave S , e Problem Att .

aod tim

I

- of ;: . 8 tion ss ~ } ro ated G :

o

Ansde ul ta alc da c l . 4 Dat Llv rica al ca nt r:i me nu erime o

o

Ex

o

ur dg

-1

- . 8

o

.8 Ots

' .

- . 4 o

~ c:: (l) M () (l) o eo 1\1 1-1

0 (l) M () (l) 0 .jJ .jJ u 4-1 4-1 u o u 'M 4-1 4-1 u 4-1 'M ....:l

/ Cu

o

o u . 8 1.

1. 6 80 100

20 o 60

o

Percent of Cho d ~ . ur e he leeward gur ace 6 . Pressure Coeffic'en t "str"bu ons on R = x 10 of h N C 6 0 0 "r oil " M ~ = "80 c.

~ 0 . Problem 102 . 1 .

Cu ve Ilo U ~ ~ !=: ' .-!

-.-!

, 8 1.

1. 6

o 60 o

o

Perce t of Chord ' nd ar su a ce P essure coefficie n istr'b tions n the r e x 10 o the AC 6 A010 'rfoil" .. 00 - , 80 Re = 0< = 0 - P ob em 10 , 1 _ «-c;.,.

nd S e nl (Interpolated s urman T n on c S s r b

---

-1. 2

-

- \

-

I

I

-.8

I

I

I

\

I

uP-; -.

I

\

A U c:

\

. ~ .~

--

"""

u ).< ~ CJ) ~ ).< P-; . 4 . 8 1.2 ... . 6 4~ ----r----"'- ---,-----.---, o

o

60 80 100 Pe rcent 0 Cho r d F' ure 28. Comparison 0 calculated press es i h g xpe im n 1 d t a and an so c s 1 'stu ba nce theory ' 011 . 80 Re = x 10 0( - 2 • oO

t .

.:"

ck

..

,:

"

/

CI

':

£-

-0

\ \ 1; . '

"-

, ,t ".

an

~tI J

es

..

~ , ..

I. .

o-de

----,.

t .u

a.t

.80 .

•• .• «

-

·011

~

t;U

"

t.

1.

1.tI ~----~----~--~----~----~----~----~----~--~~--~----~----~----~----~----t

1."

--~---- 1.

....

~ • • .

....U

\

..

.

·1 ...

··

1.url----~----~----~--~~--~----~-- ....... ...

--

:u .

. .

~ 10.

.

attac - 1.

of - 102 .

ngle b Pro rees : , d a - c:.

~ .80 .

airfoil - I O ~ 4AO J t-: ----~----~--~----~----~--~----~----~--~~--~----~----~--~----~----~--~~--~----~----~--~----~----~--_t l Lt·T

r

\ ,...; ..:T M M M I- e- ...:.

,...;

~I

.-i N CJ < .c ~ \1"\ <r s::

~

0 0 .

.

~

Q) II :z: ~.

~

0 ;K '("'\ 0 N r-!

' - r-!

\tJ 0- 0 0 'X 0 . .

.

I

°

X

'e

.-i 0 0 .

.

.

~ TlJlI.l:oU p,Io aou"ilsw lU oOl oll o ., o o o o

•

o CJ o .

r-!

~ ~~ o

______ __

________ ~

======::~~ ~-- ------~ r-!

CO"I N o

o

o o

.

.

.

OU