Numerical Study of the Turbulent FlowPast an Airfoil with Trailing Edge Separation
C. M. Rhie*Detroit Diesel A llison, Indianapolis, Indiana
andW. L. Chowt
University of Illinois at Urbana-Champaign, Urbana, Illinois
A finite volume numerical method is presented for the solution of the two-dimensional incompressible, steadyNavier-Stokes equations in general curvilinear coordinates. This method is applied to the turbulent flows overairfoils with and without trailing edge separation. The k-e model is utilized to describe the turbulent flowprocess. Body-fitted coordinates are generated for the computation. Instead of the staggered grid, an ordinarygrid system is employed for the computation and a specific scheme is developed to suppress the pressureoscillations. The results of calculations are compared with the available experimental data.
EGltG2
k •PP
Rec,ReD
u,v
x,yy+
Nomenclature= coefficients in the general finite
difference equations= constant in the law of the wall= convective terms normal to grid cell
boundaries= turbulence kinetic energy= nondimensionalized pressure= rate of production of the turbulence
kinetic energy= Reynolds numbers based on
characteristic lengths= source term in the finite difference
equation for the general scalar 0= x and y components of the non-
coordinate= coordinate transformation parameters= effective diffusity for the general
scalar </>= finite difference mesh spacings in £
and rj directions in the transformedplane
= cell boundary sizes in £ and 77 direc-tions in the transformed plane
= turbulence energy dissipation= constant in the law of the wall•=laminar and turbulent viscosities= nondimensional natural coordinates= nondimensional density= effective Prandtl numbers for the
turbulence kinetic energy and theturbulence energy dissipation
= tensor notation of the nondimensionalshear stress
jProfessor of Mechanical Engineering. Associate Fellow AIAA.
Introduction
NUMERICAL simulation of flow past airfoils is im-portant in the aerodynamic design of aircraft wings and
turbomachinery components. These lifting devices oftenattain optimum performance at the condition of "onset ofseparation." Thus, separation phenomena must be dealt withif the analysis is aimed at practical applications.
When airfoils operate at small angles of attack, the viscouseffects introduce only minor modifications to the inviscidflow. In this case, an ordinary boundary-layer analysis isadequate to describe the viscous flow behavior within theframework of "weak" viscid-inviscid interaction. However,when the angle of attack is increased, separation usuallyoccurs which significantly alters the inviscid flowfield and theflow exhibits a "strong" viscid-inviscid interaction. Oneeffective way to deal with this latter situation is to carry out alarge-scale numerical calculation on the basis of the turbulentNavier-Stokes equations. Although many attempts have beenmade to solve attached flows over airfoils using the Navier-Stokes equations (see, for example, Wu et al.,1 Steger,2 andGibeling et al.3), no significant progress has been made in thesolution of flows with separation. Most of the existingmethods incorporating simple turbulence models are notapplicable to turbulent separating flows. In an effort toovercome this limitation, a two-equation model (k-e model)was pursued recently by Shamroth and Gibeling,4 andnumerical difficulties were reported. The present study is apreliminary investigation into the problem of calculatingseparating turbulent flows over an airfoil using a two-equation turbulence model (k-e).
The k-e turbulence model has been used with ellipticsolution techniques to predict near-wake flows. However,these schemes are applied only to the near-wake region whereinlet boundary conditions are provided experimentally. Forexample, Pope and Whitelaw5 started the calculation behindthe trailing edge of a body in their study of near-wake flows.The inlet boundary conditions at the trailing edge had to besupplied from experimental data when available, or by in-tuition. In their work, they showed the extreme sensitivity ofthe downstream flowfield especially to the transverse velocityprofile adopted at the trailing edge. Consequently, they couldnot really assess the performance of the k-e model. Hah andLakshminarayana6 tried to predict the near wake of a singleairfoil. Sufficient experimental information was provided atthe inlet boundary located immediately behind the trailingedge. Although this technique was useful, the fact remainsthat a single method should be used for the whole flow region
Dow
nloa
ded
by U
NIV
OF
CA
LIF
OR
NIA
SA
N D
IEG
O o
n N
ovem
ber
17, 2
014
| http
://ar
c.ai
aa.o
rg |
DO
I: 1
0.25
14/3
.828
4
1526 C. M. RHIE AND W. L. CHOW AIAA JOURNAL
without relying on any experimental information. The presentstudy is aimed at examining the adequacy of the numericalcalculation with the k-e model to fulfill this need.
Computational MethodsAlthough the general scheme of the present method follows
essentially the original Navier-Stokes solution procedure byCaretto et al.7 (SIMPLE method), this procedure, togetherwith its later development (TEACH code), suffers severelyfrom geometric limitations since the equations were written inCartesian or cylindrical polar coordinates. Application of thismethod to curved surfaces such as airfoils must involve in-terpolation between grid points not coincident with theboundaries. This may adversely affect the accuracy of thesolution. In the present study, this geometric limitation isremoved by adopting a general curvilinear coordinate system.
There were other attempts to develop solution techniquesusing the curvilinear system of coordinates. Pope8 used aprocedure for flows with recirculation by presenting theconservation equations in an orthogonal coordinate system,and an orthogonal grid system is therefore required.Demirdzic et al.9 developed an arbitrary mesh finite volumemethod. They presented the equations in terms of arbitrarycontravariant curvilinear velocity components. This methodhas the flexibility in application to general geometries.However, the semistrong conservation form of the governingequation was used. The alternative of using the Cartesianvelocity components would present itself in the strong con-servation form. Hirt et al.10 developed the arbitraryLagrangian-Eulerian (ALE) finite volume procedure in thiscategory. In this scheme, velocity-pressure decoupling effectsare observed as a result of the grid arrangement where thepressure is defined in the center of the cell while u and vvelocity components are defined at the corners. This approachalso involves a very complicated solution procedure.
The present method deals with the Cartesian velocitycomponents and all the variables are located at the commongrid positions as opposed to staggered grid arrangementsemployed by other users of this approach.7'10 The velocity andpressure are coupled by a simple but effective scheme. Themethod is aimed at the future application to complex three-dimensional turbulent flow computations.
Theoretical FormulationGoverning Equations
For steady-mean flow, the time-averaged continuityequation and the Navier-Stokes equation, in conjunction withthe isotropic turbulent viscosity hypothesis, are written in aCartesian tensor form
(1)
(2)t) ] -2/3pkdij
where p is the mean density, Uj the mean velocity, and/? themean pressure. From the k-e turbulence model,11 the tur-bulent viscosity /*, is given by
(3)
where k is the turbulence kinetic energy and e is the turbulenceenergy dissipation.
The model is then composed of two equations; one for kand another for e presented as follows:
d • d / fjit dk\——pU:k=-— { —- —— ) +P—pedXj dxf \ ok d" ' (4)
(5)
where P ( = -pu-UjdUj/dXf) is the rate of production ofturbulence kinetic energy. This model contains five empiricalconstants which assume the following values:
= 0.09, 2 = l.90, (6)
The turbulence scalar transport equations [Eqs. (4) and(5)] , are only valid for fully turbulent regions. An additionalmodel must be introduced to treat the laminar sublayerregion. The wall function method12 is used in the presentstudy to eliminate the large number of grid points needed toresolve the laminar sublayer. The following functions are usedto bridge the near- wall region:
where .y =p(kpC» ) 1
dkdy
= 0
(7)
(8)
(9)
(10)
Here, y is the conventional coordinate normal to the walland TW is the wall shear stress. The subscript p refers to thegrid node next to the wall;-* and E are the constants from thelaw of the wall, with values of 0.4 and 8.8, respectively. Itshould be mentioned that this model is not capable of dealingwith laminar and transition regions. The present study paysspecific attention to the shear layer interaction and its sub-sequent relaxation in the near wake. It is expected that theleading edge transition process will not significantly affect thedownstream flow condition of the present problem.
Transformation of the Basic EquationsThe set of conservation equations typically can be written in
the Cartesian system of coordinates for a scalar variable as
— (pu<t>)+ —dx dy
d—dy
d—dy -dy (ID
where F^ is an effective diffusion coefficient and S* denotesthe source term. When new independent variables £ and r? areintroduced, Eq. (11) changes according to the general trans-formation % = %(x,y), ir) = T)(x,y). Partial derivatives of anyfunction /are transformed according to
where / is the Jacobian of the transformation given by
(12b)
Since the strong conservative form of the equation isdesirable for numerical computations, the integral form of theconservational equations over a finite volume element ispreferable. Upon introducing
(13)
04)
the following integral conservation relation is obtained fromEq. (1 1) for an arbitrary scalar dependent variable 0:
(15)
Here G7 and G2 are directly related to the contravariantvelocity components. In fact, G1/(a)1/2 and G2/(y)'/2 representthe velocity components normal to the lines of constant £ andrj, respectively, in the physical plane.
Method of ComputationGeneral Transport Equations
In terms of the notation shown in Fig. 1 for a typical gridnode P enclosed in its cell and surrounded by its neighbors N,S, E, and W, the finite difference approximation to the in-tegral conservation relation, Eq. (15), over the cell is writtenas
Ref. 13. The terms within the bracket in Eq. (18) areoriginated from cross derivatives in the diffusion terms andare the results of the nonorthogonal coordinate system. Theseterms are usually very small and can be combined with thesource term "and treated as known quantities. They areevaluated from updated values through the following finitedifference approximation
—— 09)
Pressure Correction EquationAn equation for the remaining unknown pressure is
established by combining the continuity and momentumequations. For a simple scheme, the pressure correction canbe derived from the momentum finite difference equations.To explain this procedure in some detail, it is assumed that apreliminary set of velocity components, u and v, are obtainedfrom
EWNS (20)
P=V=
where
C«=-
Bv=-
(16)
where 4> is the unknown. In the present scheme, all flowproperties are defined only at the nodes P, E, W, N, and S.Thus, in general the following approximations are made forthe above finite difference expression, Eq. (16).
pG7</>A?7 1 c j Arj) - (</>P + </>E )
(17)
Quantities such as (pG7 Ar?)e and (T/JaAr])e are obtained bylinear interpolation in the physical plane. Substituting allrelations of the type of Eq. (17) into Eq. (16), a relationbetween <£p and the neighboring values is obtained, i.e.,
A.p<l>P =
(18)
where the coefficients A involve the flow properties ofconvection, diffusion, area, etc. These coefficients are latermodified according to a "hybrid scheme,"7 which evaluatesthe convective term by the first-order upwind-differencingscheme whenever the grid cell Reynolds number (e.g.,p \Gj UdZ/otY*) is greater than 2. The details can be found in
and Sf and 'Sf are residues from SU,SV after the pressuregradient terms have been extracted from them. The super-script * for « and v denotes that they are based on theestimated pressure field p*. In general, u* and v* will notsatisfy the continuity equation. Instead a net mass source isproduced. To remove this mass source, the velocity com-ponents are assumed to be corrected by the relations
(21a)
(21b)
where p' is the pressure correction which is related to thepressure/? according to
jp=p*+p' (21c)
Subsequently, the correction equations for Gl and G2 areobtained from Eq. (13), i.e.,
= G* (22a)
(22b)
where again G J and G2 are based on u* and v*.It is noted that the last two terms of Eqs. (22a) and (22b) are
negligible if the grid system is nearly orthogonal. One may,therefore, further simplify the above expressions byneglecting these terms, and the result becomes
(23)
where B = Buy^ -Bvx^ and C= Cvx^ - Cuy^.One may expect that this highly approximate relation, Eq.
(23), leads to the establishment of the correct pressure field
Dow
nloa
ded
by U
NIV
OF
CA
LIF
OR
NIA
SA
N D
IEG
O o
n N
ovem
ber
17, 2
014
| http
://ar
c.ai
aa.o
rg |
DO
I: 1
0.25
14/3
.828
4
1528 C. M. RHIE AND W. L. CHOW AIAA JOURNAL
after consideration of the following. If the solution con-verges, the correction terms at the final converged stateshould vanish, and the converged solution, nevertheless,satisfies the basic Eq. (20). The correction equations thus donot have to be rigorously derived from Eq. (20). This form ispreferred because of its simplicity and better convergenceproperties. Upon substituting Eq. (23) into the continuityequation, one obtains a Poisson equation for the pressurecorrection, i.e.,
(24)
where mp is the integrated mass source from u* and v*established for the element under consideration.
With the second-order center-difference approximationsfor the pressure gradients on the cell boundaries, the finitedifference form of Eq. (24) is
(25)
where coefficients A involve coefficients B, C, density, etc.,and 5 represents the local imbalance of mass. Equation (25)can now be solved with the same algorithm for Eq. (18) toyield the corrective pressure. Thus, the pressure and velocitycomponents are corrected according to Eq. (21).
It has been learned from this pressure correction schemewith the ordinary grid arrangement, that the oscillatorypressure field is produced. This has been observed previouslyby others.14 The major source of this instability originatedfrom the second-order centered 2AJ£-difference ap-proximation of the pressure gradient at the grid node. Thisscheme cannot sense the lAA"-pressure oscillation. Astaggered mesh was successfully used by many7'9 to removethese oscillations since the centered lAJf-difference could beused for the pressure gradient with the staggeredarrangement. However, this technique cannot be applied tothe present method since the Cartesian velocity components uand v are not related to the grid line orientations. As analternative approach, the scheme of Hirt et al.,10 where thepressure is defined in the center of the cell, while u and vvelocity components are defined at all corners, was in-troduced. In this case, pressure oscillation remained in thediagonal direction, and this scheme is equally unsuitable forthe present purpose.
To eliminate the oscillations with the present ordinary gridarrangement, a new treatment of the locally linearizedconvective terms, Gl and G2, was introduced. It should beremembered that these convective terms are always located onthe cell boundaries and are obtained by interpolations be-tween grid nodes. These interpolations were responsible fordecoupling the velocity and pressure fields. If G* were derivedfrom Eq. (13), where u* and v* were obtained from Eq. (20),it would assume the expression
G*!=Bpl + ... (26)
where all other terms are dropped for the convenience ofdiscussion. Even though G* is related top* in the form of Eq.(26), the G* value on the cell boundary is actually obtainedfrom interpolation between the grid nodes. Thus the G; valueso determined cannot detect 16£-pressure variation. One mayremedy this situation by correcting the p\ term through a 1<5£-difference scheme on the cell boundary. The effective formwas found to be
(27)
terpolated pressure gradient term p* based on the 25£-centerdifference is replaced by the l<5£-center difference on the cellboundary. This procedure ensures strong velocity-pressurecoupling. In some test cases, this procedure was proven toshow better convergence behavior than the original SIMPLEprocedure. It is worthwhile to note that this modificationthrough Eq. (27) becomes redundant if the pressure varieslinearly within the field.
The boundary condition required to solve the pressurecorrection equation is that the normal derivative of thepressure correction p' vanishes on the boundary in thecomputational plane. The actual pressure value on theboundary is established by extrapolation from interior valueswith the conditions of vanishing normal derivatives ofpressure there.
Solution ProcedureThe overall solution procedure is similar to that of the
SIMPLE procedure.7 A field of intermediate velocitycomponents, w* and y*, is obtained by solving the momentumequations, Eqs. (20), using a preliminary pressure field.Subsequently, the continuity constraint is enforced fromsetting up an equation of the corrective pressure p' . Bysolving the resulting equation, Eq. (25), the required ad-justments to the pressure and velocity components aredetermined. Equation (18) for the remaining scalar variables,k and e, are solved in turn. This processes defined as oneiteration cycle and is repeated until a converged solution isobtained.. Convergence is reached for vanishing mass sources(within an arbitrarily small margin) throughout the com-putational domain.
Generation of the Coordinate SystemThe grid generation scheme developed by Sorenson and
Steger15 is adopted in the present study. In this method, thecurvilinear coordinates are generated by solving the ellipticequations
= 0 = 0 (28)
After the grid is constructed, a simple exponential stret-ching technique is used to cluster points near the airfoilsurface. This is done in order to resolve the high gradient flowwithin the boundary layer. In the present study, an additionalmodification is made to eliminate the coordinate slopediscontinuity in the wake region, shown in Fig. 2a. Thediscontinuity is removed by fitting fourth-order polynomialsacross the cut. The resulting grid system is illustrated in Fig.2b. With this modification, the grid points along the cut canbe treated as regular interior nodes.
Results and DiscussionThe objective of the present work is the study of turbulent
flow over an isolated airfoil and correlation of the computedresults with experimental data. Before any effort of this
where the overbar denotes the regular results obtained fromlinear interpolation between grid nodes E and P. The in-
Grid detail near body, a) After reclustering, b) After partialreconstruction.
Dow
nloa
ded
by U
NIV
OF
CA
LIF
OR
NIA
SA
N D
IEG
O o
n N
ovem
ber
17, 2
014
| http
://ar
c.ai
aa.o
rg |
DO
I: 1
0.25
14/3
.828
4
NOVEMBER 1983 TURBULENT FLOW PAST AN AIRFOIL 1529
nature was carried out, a variety of laminar flow cases,ranging from a flow past a circular cylinder to a flow past anNACA 0012 airfoil at low Reynolds number, was computedto check the numerical accuracy of the computer program sodeveloped. The results were all satisfactory when comparedwith those previously published.
The experimental data of turbulent flow past isolatedairfoils are not ^abundantly available at the present time.Recently, Hah and Lakshminarayana6 carried out ex-perimental investigations of an NACA 0012 airfoil at severalangles of incidence. Coles and Wadcock16 tested an NACA4412 airfoil at a maximum lift condition with a 13.87 degangle of attack. The present effort was aimed at the com-putation of these flow conditions for which experimental datawere available for comparison. However, it is noted that thereis no single set of experimental data covering the entireflowfield. The experimental data by Hah and Laksh-minarayana6 were restricted to the near-wake region. Theexperimental data by Coles and Wadcock16 covered only therear part of the suction surface and the near-wake region,although the pressure distribution around the whole airfoilsurface was reported. The surface pressure data by Gregoryand O'Reilly17 for an NACA 0012 airfoil were also availablefor comparison.
Results for an NACA 0012 Airfoil without Separationat/tec=2.8x!06
As a preliminary effort, the external flow past an NACA0012 airfoil was considered at Rec = 2.8 x 106. A 77 x 34 gridwas used and 49 grid points were distributed over the airfoil.The front and rear outer boundaries were located eight chordlengths away from the body, and the top and bottom outerboundaries were located at 12 chord lengths from the body.The velocity was set to freestream value over the front, top,and bottom outer boundaries. On the rear outer boundary,which was normal to the airfoil chord line, the velocitycomponents and the turbulent scalars were extrapolated fromthe inner solution by assuming that the first derivatives offlow properties along the r/ = const lines vanish. The velocityprofile normal to the exit plane was then adjusted to satisfythe principle of global conservation of mass.
0 deg A ngle of A ttackThe computation was performed in the whole domain even
for 0 deg incidence to avoid numerical errors from ex-trapolation along the line of symmetry. The surface pressuredistribution is plotted in Fig. 3 and compared with the ex-perimental data by Gregory and O'Reilly,17 and the com-putational results by Shamroth and Gibeling.4 In general, thepressure distribution shows good agreement with the ex-perimental data. Slight discrepancies are shown near the
-0.4
0.8
1.2
Present ResultRec=2.8xl06
_ __ Experiment, Gregory et. al.Rec=2.8xl06 (17)
— —— Theory, Shamroth et. al.Rec=10x 106 (4)
0.2 0.4 0.6 0.8 1.0x /C
leading and trailing edges. The minimum pressure near theleading edge was underestimated and a lower pressure waspredicted in the trailing edge region. Two possible reasons forthese discrepancies are given. First, the boundary-layergrowth was not accurately predicted because of inadequateturbulence modeling. Second, the grid spacing in the leadingand trailing edge regions was not sufficiently fine to resolvethe strong gradients in those regions.
6 deg A ngle of A ttackIn this case, all the parameters in the computation were
kept the same as the previous case except the airfoil is rotatedthrough 6 deg. The calculated surface pressure distribution isplotted in Fig. 4 along with the experimental data from Ref.17 and the computational results from Ref. 4. Again, theagreement with the data is good except in the leading andtrailing edge regions. It is anticipated that a finer grid in thoseregions would substantially improve the accuracy of thepresent results.
Flow Past an NACA 0012 Airfoil without Separationat/tec=3.8xlOs
The flow over the NACA 0012 airfoil at 6 deg incidence wascalculated for Rec = 3.8x 105 and the results were comparedwith the experimental data of Hah and Lakshminarayana.The wind tunnel blockage effect was accounted for by treatingthe computational outer boundaries as the physical locationsof the wind tunnel wall which was located at three chordlengths away from the body. The front and rear outerboundaries were located three and five chord lengths awayfrom the body. The front inlet boundary was defined by halfof the circular arc. An 87 x 30 grid was used with 59 pointsdistributed over the airfoil surface. The grid spacing in the rjdirection along the airfoil was 0.0025 chords. The treatmentof the boundary conditions was essentially the same as that inthe previous computations. The boundary-layer development
-3.0
-2.0
-i.o
Q.O
1.0
^^_ Present Result——— Rec=a.8xl06
_ _ Experiment, Gregory et. al.Rec = 2.8xl06(17)
--—— Theory, Shamroth et. al.Rec = 1.0xl06 (4)
0.2 0.8 1.00.4 0.6x /C
Fig. 4 Surface pressure distribution for NACA 0012 airfoil at 6 degincidence (Rec = 2.8 x 106).
-3.0
-2.0
-i.o
0.0
1.0
Present ResultRec=3.8xl05
Experiment, Gregory et. al.Rec=2,8xl06 (17)
0.2 0.4 0.6x/C
0.8 I.O
Fig. 3 Surface pressure distribution for NACA 0012 airfoil at zeroincidence (Rec = 2.8 X106).
Fig. 5 Surface pressure distribution for NACA 0012 airfoil at 6 degincidence (Rec= 3.8 x 105).
Dow
nloa
ded
by U
NIV
OF
CA
LIF
OR
NIA
SA
N D
IEG
O o
n N
ovem
ber
17, 2
014
| http
://ar
c.ai
aa.o
rg |
DO
I: 1
0.25
14/3
.828
4
1530 C. M. RHIE AND W. L. CHOW AIAA JOURNAL
on the wind tunnel wall was neglected and a simple slipboundary condition was applied.
Figure 5 shows the calculated surface pressure distribution.Some improvement in the accuracy of the results over thosepresented in Fig. 4 can be observed at the leading and trailingedges. This is probably the result of the finer grid used in thiscase. Some discrepancies still remain where excessively highgradients exist. It is speculated that the discrepancy near theleading edge results from inadequate resolution since inviscidflow governs the flow process there. While inadequateresolution is considered to be the dominant source of errornear the leading edge, the effect of inadequate turbulencemodeling becomes equally important near the trailing edgeregion. There, the boundary layers on both the upper andlower surfaces have considerable thickness and, therefore,significantly modify the pressure distribution. The shear-layerinteraction and the subsequent relaxation process dominatethe flow events downstream of the trailing edge which canalso alter the upstream flow conditions. Another problemwith the present approach was found in the near-wake region.Although it is well known that the mean flow behavior withinthe wake is entirely different from the wall turbulent flows,
Fig. 6 Streamlines for NACA airfoil at 6 deg incidence
0.3
y/C
0.0
x/C = 0.642 0.908 1.1747 1.9515
0. I. 0. . I. 0. I. 0. I.
a)
0.3
y/c
0.0
x/d:= o
J
642
r
0.9C)8
""""_ < ^
M
p> '
1.95
- (
15
\- 2 0 2 - 2 0 2 - 2 0 2 - 2 0 2 - xio"
b)U'V'/U ref
Fig. 7 Computed results for NACA 0012 airfoil at 6 deg incidence(Rec = 3.8 x 105). a) w-velocity profiles, b) u' v' profiles.
the effect of the upstream logarithmic profile should havepersisted within the immediate near wake and the presentapproach cannot account for this effect. The grid spacingextending from the airfoil boundary layers was not fineenough to allow an adequate description of such a near-wakeprofile, thus some error in evaluating the shear stress alongthe wake center line occurred.
Streamlines for this example problem are plotted in Fig. 6.The streamlines indicate a thicker viscous layer on the suctionsurface. The w-velocity component and Reynolds shear stressprofiles are plotted at different x locations in Fig. 7. Theseresults indicate turbulent structural behavior in the boundarylayer and the near-wake regions. Downstream of the trailingedge* the near-wake profile undergoes a fast relaxationprocess. In this region, relatively high Reynolds stresses areproduced due to high mean shear and eventually they arereduced due to the diffusive and dissipative nature of theprocess. A detailed comparison with the experimental data6
for the near wake was not attempted since the present ap-proach involves some errors in the near-wake flow aspreviously discussed. However, the evaluation of the wakecenterline velocity is presented and compared with the ex-perimental data in Fig. 8. Reasonably good correlation isshown between the data sets.
Flow Past an NACA Airfoil with Separation at Rec = 1,5 X106
A major part of this study was aimed at the prediction ofthe flow past an NACA 4412 airfoil at a maximum liftcondition with a 13.87 deg angle of attack. The blockageeffect of the confined wind tunnel test section was simulatedin the computations. The computational domain was definedas illustrated in Fig. 9. All the lengths were non-dimensiotialized by the airfoil chord. The converging anddiverging tunnel wall section were defined with a = /3=4 deg.
1.0
0.8
0.6
0.4
0.2
Present Result
Experiment, Hah andLakshminarayana (6)
O.4 0.8S/C
1.2 1.6
Fig. 8 Evolution of the wake center velocity for NACA 0012 airfoilat 6 deg incidence (Rec = 3.8 X105).
^^-^^
•
r— ̂
yjV/o
(L.ar<;/
V = 1.9308 **
13.87°
-I.45I6_J
x=-3.0 -0.5123 2.8698
Ref. (I.I466,-U209)
6.0
Fig. 9 Boundaries for computations of the flow past an NACA 4412airfoil at 13.87 deg incidence.
Dow
nloa
ded
by U
NIV
OF
CA
LIF
OR
NIA
SA
N D
IEG
O o
n N
ovem
ber
17, 2
014
| http
://ar
c.ai
aa.o
rg |
DO
I: 1
0.25
14/3
.828
4
NOVEMBER 1983 TURBULENT FLOW PAST AN AIRFOIL 1531
- 8 - -
Present Result
Experiment, Coles andWadcock (16)
0.2 0.4 0.6 0.8 1.0
x/CFig. 10 Surface pressure distribution for NACA 4412 airfoil at 13.87deg incidence (Rec = 1.5 x 106).
Fig. 11 Streamlines for NACA 4412 airfoil at 13.87 deg incidence(/tec = 1.5xl06).
The computation was performed on a 95x31 grid with 63points distribution along the airfoil surface. Uniform.ap-proaching freestream conditions were imposed on the frontouter boundary. Wake profile boundary condition simulationat the exit plane of the diverging diffuser section was aban-doned due to numerical difficulties. Instead, a uniformvelocity profile was imposed. Virtually no influence of thisrear outer boundary condition was observed on the innersolution near the airfoil. The slip boundary conditions wereimposed on the wind tunnel walls.
The calculated surface pressure distribution, plotted in Fig.10, shows very good agreement with the experimental data.The discrepancies shown near the leading edge are attributedto inadequate resolution of the strong gradient in that region.In the preliminary computations, the successive gridrefinements gave consistently improved results. Un-fortunately, further grid refinement was not possible becauseof the limitation on computer storage. Figure 11 shows thestreamlines and velocity vectors from the computation. Theseresults clearly show the trailing edge separation bubble. Thefree streamlines were modified due to the displacement effectof the bubble, which, in turn, altered the airfoil surfacepressure distribution. Detailed velocity component profileand Reynolds shear stress profiles are given in Fig. 12 alongwith the experimental data. The results are plotted vs distancenormal to the airfoil chordline. The discrepancies between thepredicted and experimentally determined mean velocity
0.5 -
y/C
o.o -
a)
x/C = 0.642
-J1 |
. i . . . . i . ,
0.908
J". I -
. i . , , . i . ,
1.1747
• 4. i . . . . i
1.9515
(
, i , , , ,o i.o o i.o o i.o o i.o
u/uref
0.5 -
y/c
o.o -
x'/C = 0
"
.642 0.908 1.17
J: 1
47
r-
1.9E>15
. . . . i .- 2 0 2 - 2 0 2 - 2 0 2 - 2 0 2
uV/u?ef
xio"
b)Fig. 12 Results for NACA 4412 airfoil at 13.87 deg incidence(Rec = 1.5 x 106): ——, presented result; n , experiment (Ref. 16).
profiles in Fig. 12a are attributed to inadequate gridrefinement and possibly inadequate turbulence modeling.
The predicted downstream near-wake profile after thetrailing edge generally shows lower minimum wake velocitiesand larger wake widths than indicated in the experimentalresults. This suggests that the mixing is less active in the nearwake, and that the rapid development of the wake is sup-pressed. This trend was also observed by Pope and Whitelaw5
in the computation of separated flows using the k-e model.It is premature to make any firm assessment on the k-e
turbulence model with the present study since the gridrefinement was not sufficient to obtain the grid independentresults. However, some remarks can be made on the tur-bulence model. The k-e model is based on simple flows such aschannel Couette flow and the decaying isotropic grid tur-bulence. The capability of the k-e model is naturally in doubtfor the complex separated flows. The qualitative picture ofthe present prediction is in accordance with what has beenreported in Ref. 5. Once again it is confirmed that a betterturbulence model is required for the accurate quantitativepredictions. The suitability of the wall function method forthe present application needs to be investigated further.
All of the computations were initialized with the freestreamflow conditions. Generally, good convergence behavior wasobserved. The calculations were made on a CYBER 175computer. The time and storage requirements of the Fortrancomputer program were 0.00083 s/grid node/iteration and23,440 words plus 37.8 words/grid node, respectively. About500 iterations were needed for the case without separation andtwice as many iterations were required for the case withseparation.
ConclusionA finite volume method for the calculation of recirculating
turbulent flows in general curvilinear coordinates wasdeveloped. Attached and separated flows over an isolated
Dow
nloa
ded
by U
NIV
OF
CA
LIF
OR
NIA
SA
N D
IEG
O o
n N
ovem
ber
17, 2
014
| http
://ar
c.ai
aa.o
rg |
DO
I: 1
0.25
14/3
.828
4
1532 C. M. RHIE AND W. L. CHOW AIAA JOURNAL
airfoil were calculated using the standard k-e turbulencemodel. Detailed comparisons on the predicted result andexperimental data were made for the wakes and separatedflow regions. It was determined that the accuracy of theresults was highly dependent on degree of resolution of thestrong leading and trailing edge gradients. Withoutseparation, the k-e turbulence model predicted values inreasonably good agreement with the experimental data. Withseparation, the k-e turbulence model predicted poor resultseven though fairly good pressure distribution on the airfoiland reasonable overall behavior were predicted. Therequirement of the better turbulence model is indicated.
AcknowledgmentsThis work was conducted at the University of Illinois at
Urbana-Champaign with partial support from U.S. ArmyResearch Office through Grant DAAG 29-79-C-184. Theauthors would also like to express their thanks to Dr. D.Sharma for sharing his knowledge on the SIMPLE method.
References1 Wu, J. C., Sampath, S., and Sankar, N. L., "Dynamic Stall of an
Oscillating Airfoil," Proceedings of AGARD Conference on Un-steady Aerodynamics, NATO, Advisory Group of Aero. Research &Development, 1977, pp. 24-1 to 24-18.
2Steger, J. L., ''Implicit Finite Difference Simulation of FlowAbout Arbitrary Two-Dimensional Geometries," AIAA Journal,Vol. 16, July 1978, pp. 679-686.
3Gibeling, H. J., Shamroth, S. J., and Eiseman, P. R., "Analysisof Strong Interaction Dynamic Stall for Laminar Flow of Airfoils,"NASA CR-2969, 1978.
4Shamroth, S. J. and Gibeling, H. J., "The Prediction of theTurbulent Flow Field about an Isolated Airfoil," AIAA Paper 79-1543,1979.
5 Pope, S. B. and Whitelaw, J. H., "The Calculation of Near-WakeFlow," Journal of Fluid Mechanics, Vol. 73, Pt. 1, Jan. 1976, pp. 9-32.
6Han, C. and Lakshminarayana, B., "Prediction of Two- andThree-Dimensional Asymmetrical Turbulent Wakes, IncludingCurvature and Rotation Effects," AIAA Journal, Vol. 18, Oct. 1980,pp. 1196-1204.
7Caretto, L. S., Gosman, A. D., Patankar, S. V., and Spalding, D.B., "Two Calculation Procedures for Steady, Three-DimensionalFlows with Recirculation," Proceedings of the Third InternationalConference on Numerical Methods in Fluid Dynamics, Springer-Verlag, 1972, pp. 60-68.
8Pope, S. B., "The Calculation of Turbulent Recirculating Flowsin General Orthogonal Coordinates," Journal of ComputationalPhysics, Vol. 26, No. 2, 1978, pp. 197-217.
9Demirdzic, L, Gosman, A. D., and Issa, R. L, "A Finite-VolumeMethod for the Prediction of Turbulent Flow in ArbitraryGeometries," Paper presented at 7th International Conference onNumerical Methods in Fluid Dynamics, Stanford University andNASA Ames, June 1980.
10Hirt, C. W., Amsden, A. A., and Cook, J. L., "An ArbitraryLagrangian-Eulerian Computing Method for All Flow Speeds,"Journal of Computational Physics, Vol. 14, No. .3, 1974, pp. 227-253.
11 Jones, W. P. and Launder, B. E., "The Prediction ofLaminarization with a Two-Equation Model of Turbulence," In-ternational Journal of Heat and Mass Transfer, Vol. 15, Feb. 1972,pp. 301-314.
12Launder, B. E. and Spalding, D. B., "The Numerical Calculationof Turbulent Flows," Computer Methods in Applied Mechanics andEngineering, Vol. 3, No. 2,1974, pp. 269-289.
13Rhie, C. M., "A Numerical Study of the Flow Past an IsolatedAirfoil with Separation," Ph.D. Thesis, Dept. of Mechanical andIndustrial Engineering, University of Illinois at Urbana-Champaign,1981.
14Rubin, S. J. et al., "Numerical Studies of Incompressible ViscousFlow in a Driven Cavity," NASA SP-378, 1975.
15Sorenson, R. L. and Steger, J. L., "Simplified Clustering ofNonorthogonal Grids Generated by Elliptic Partial DifferentialEquations," NASA TM 73252, 1977.
16Coles, D. and Wadcock, A. J., "Flying Hotwire Study of FlowPast an NACA 4412 Airfoil at Maximum Lift," AIAA Journal, Vol.17, April 1979, pp. 321-329.
17Gregory, N. and O'Reilly, C. L., "Low Speed AerodynamicCharacteristics of NACA 0012 Airfoil Section, Including the Effectsof Upper Surface Roughness Simulation Hoarfrost," NationalPhysical Laboratory, Teddington, England, Aero Rept. 1308, 1970.
