Int. J. Engng Sci. Vol. 33, No. 4, pp. 469-483, 1995
Pergamon
0020-7225(94)00087-5
TURBULENT
TRANSONIC
FLOW
PRESSURE-BASED
Copyright © 1995 Elsevier Science Ltd
Printed in Great Britain. All rights reserved
0020-7225/95 $9.50+ 0.00
SIMULATION
METHOD
USING
A
Y. G. LAI, R. M. C. SOt and A. J. P R Z E K W A S
CFD Research Corporation, Huntsville, A L 35802, U.S.A.
(Communicated by C. G. SPEZIALE)
Abstrae/--In the past, turbulent compressible aerodynamic flow calculations were usually carried out
with density-based numerical methods and zero-equation turbulence models. However, pressurebased numerical methods and more advanced turbulence models have been routinely used in industry
for many internal flow simulations and for incompressible flow calculations. The application of
pressure-based numerical method and advanced turbulence models to the calculations of aerodynamic
flows is still not well demonstrated and accepted. In this study, an advanced pressure-based numerical
method and a recently proposed near-wall compressible two-equation turbulence model are used to
calculate external transonic aerodynamic flows. Several TVD-type schemes are incorporated into the
pressure-based method to better capture discontinuities such as shocks. A compressible near-wall
two-equation turbulence model is implemented into the numerical scheme to calculate transonic
turbulent flows over NACA 0012 and RA12822 airfoils with and without shocks. Good correlations
are obtained with wind tunnel data as well as with results calculated using the Baldwin-Lomax model.
Therefore, this study demonstrates that advanced turbulence models capable of handling flows with
separation and recirculation could be implemented with ease to a pressure-based method with
improved capability to resolve discontinuities to calculate turbulent compressible aerodynamic flows.
INTRODUCTION
Traditionally, density-based methods, which use density as a dependent variable in the
continuity equation, have been widely used for compressible aerodynamic flow calculations.
However, pressure-based methods, which select pressure as a primitive variable, have become
very popular for incompressible flows and have been used widely for turbulent and chemically
reacting flow calculations.
When the Euler equations are solved numerically using density-based methods, the
non-linear wave properties are taken into consideration by locally solving Riemann problems
and this leads to a robust numerical algorithm. Moreover, highly accurate schemes can be
constructed with a density-based method which are capable of capturing shocks within a few
grid points. Some of these high resolution schemes include TVD [1-3], MUSCL [4], PPM [5]
and E N O [6]. Despite the success achieved with Euler equations, however, the accuracy and
robustness of a density-based method deteriorate substantially at low Mach numbers and the
method becomes even less efficient for turbulent and recirculating flow calculations.
The weakness of a density-based method is exactly the strength of a pressure-based method.
The pressure-based method was developed to solve incompressible flows where density is a
constant. As a result, pressure is chosen as a primitive variable and the continuity equation is
used to form a Poisson-like equation for the pressure or pressure-correction. In its early form,
the method is built upon a staggered grid system [7] and is primarily used to solve
incompressible flows. For a complete description of the method, readers are referred to the
book by Patankar [8]. Because of its demonstrated robustness and efficiency, the pressurebased method was quickly adopted and widely used in industry to calculate many incompressible complex engineering flows. In recent years, various improvements have been made to the
pressure-based method in order to extend it to compressible flow calculations. One of the most
important improvements is the use of non-staggered or collocated grid in place of traditional
t T o whom correspondence should be addressed at: Mechanical and Aerospace Engineering, Arizona State University,
Tempe, A Z 85287-6106, U.S.A.
469
470
Y.G. LAI et al.
staggered grid. The use of non-staggered grid greatly simplifies the code development and
many state-of-the-art numerical methods such as multigrid, multi-zone, local grid refinement,
etc. can be incorporated in a relatively simple manner. Most of the improvements made to the
pressure-based method can be found in [9-16].
From the above discussion, the pressure-based method, combined with its ability to resolve
numerical stiffness associated with turbulent separated flows, can be very competitive when
used to calculate compressible flows. More importantly, the extension of the pressure-based
method to compressible flows makes the approach valid for all velocity regimes. Thus the
approach is very convenient to use for various computational purposes. The major objective of
this study is to examine the viability of a pressure-based method for compressible turbulent
flow calculations using some of the more advanced turbulence models. In order to accomplish
this objective, most of the recent advances made in pressure-based method are adopted.
Particularly, several TVD limiters used in density-based methods are extended to the
pressure-based method. These limiters are first examined for shock capturing capability and the
results are compared with those from a density-based method. The intention here is not to
compete with density-based methods for solving Euler equations. Rather, the comparison will
give a quantitative assessment of the ability of the pressure-based method to resolve
discontinuities such as shocks. A second objective of the present study is to assess the ease with
which near-wall compressible two-equation k - e turbulence models can be incorporated into the
pressure-based numerical method to calculate turbulent compressible flows. This computation
serves not only as a test of the numerical algorithm, because the equations of the near-wall k - e
model are known to be stiff, but also as a way to validate the near-wall k - e model for
aerodynamic flow calculations.
Among the many near-wall compressible two-equation models formulated to date, the k - e
model of [17] has been found to give good predictions of boundary layer flows up to a
free-stream Mach number of 10.31 and out-performs the k-co model of [18] even at a fairly low
free-stream Mach number of 2.244. Furthermore, the model was also evaluated in [19] and was
found to be quite promising. Therefore, the k - e model of [17] is chosen for this study.
The governing equations and the near-wall k - e turbulence model of [17] are first introduced.
This is followed by a discussion of the basics of the pressure-based numerical method.
Attention is given to the extension of TVD limiters and several numerical convergence issues.
The supersonic flow over a 4% bump is used to evaluate the shock capturing capability of the
pressure-based method and the results are compared with those obtained from a density-based
method. Finally, transonic turbulent flows over NACA 0012 and RAE 2822 airfoils with and
without shocks are calculated. The results are compared with wind tunnel data and calculations
obtained from a zero-equation model [20].
G O V E R N I N G EQUATIONS AND T U R B U L E N C E MODEL
For compressible turbulent flows, the mass, momentum, and energy equations are written in
terms of Favre-averaged quantities. With the neglect of bulk viscosity and body forces, the
mean flow equations can be expressed in Cartesian tensor form for an ideal gas as:
B0__~_~+ (•a,)., = 0,
(1)
Ot
O~a~ +
at
=
(t~a,aj),j
~
- ~ , + [#(a,,j +
a~,,)lj
-
2
(~u;'u;'),~ - ~ (#aj,,),,,
(2)
Turbulent transonic flow simulation
O/SCpT ~_( p u i C p T ) i = OP
Ot
'
Ot
--
'
471
+ ug P~ + u~p,i -e crgjagj + ~r~jug,j + fie
'
'
(3)
(/sCpu;t~-'Tt").i -~ ( K t g ),i,
/5 =/5R7",
(4)
where
2
#ij = -- -3 ~Uk,k~ij Jr- ~ ( a i , j Jr- Uj,i),
G ijl~ii,j .
In the above, R is the gas constant, the overbar (-) and (') are used to denoted the
Reynolds-averaged and fluctuating parts of an instantaneous quantity while the tilde ( - ) and
(") are used to denote the corresponding Favre-averaged and fluctuating parts, ug is the ith
component of the instantaneous velocity, (.g) is used to denote partial differential with respect
to the ith component of the coordinate and t is time. Also, note that density (p) and pressure
(p) have been decomposed according to Reynolds-averaging, while other primitive variables
use the Favre decomposition. In deriving the above equations, the Einstein summation
convention is followed and the turbulent fluctuations of the dynamic fluid viscosity (/z), the
thermal conductivity (K), and the specific heat ( C p ) are assumed negligible.
In general, The Reynolds stresses, - f u ; ' u ; ' , Reynolds heat fluxes, - f u ; ' T " and other extra
terms in the energy equation need to be specified to close the above system of equations and
this constitutes the turbulence modeling problem. The detailed modeling strategy was given in
[17], therefore, it is not repeated here. It is sufficient to point out that the Reynolds stress
tensor is closed with a gradient transport assumption, or
__
2)2
(
-/sui'u;' = ~ , ag.j + aj.g - g ak.kagj
- g /skS,j,
k2
~, = f i G L - ,
(5)
g
where /,, is the eddy viscosity, k is the turbulent kinetic energy and e is the solenoidal
dissipation rate of k. The transport of k and e is governed by the following modeled equations;
- - + (/saik), i = - ( 1 - a l M , ) 2p u ---~7-_
g u j u g,,,j
Ot
- 1511 + (cx - G 2 ) M 2 ] e
_
o~e
- -
Ot
+ ~,)k.,]
(6)
+ (~aie),~ = - C e l p U-i U ),'-r~,,;e( u g , ] - ~1 ak ,k S g,j )
where Y is the specific heat ratio. The near-wall correcting and unknown functions in (6) and
(7) are defined as
= f~2fi 2 - ~ - - 1.5
g = e-
( OV~k~ 2
29\--~x~ / ,
,
M 2
g = e -
c2,
29k
y---C,
where y is the coordinate normal to the wall, c is the local sound speed and the model
E$ 33-4-8
472
Y.G. LAI et aL
constants are given by: C~, = 0.096, O'k = 0.75, O'~ = 1.45, Gel = 1.3, Cez = 1.85, ot = 0.5, a~ = 0.4,
a2 = 0.2. Finally, the damping functions are specified as
f~,=(1
345\
+
+ ~ R t ) t a n h ( y /115),
k2
]:.,2 = e -~a'/64)~,
Rt = ~ ,
where y+ = y u d ~ is the normalized coordinate and u~ is the friction velocity defined with
respect to wall variables.
PRESSURE-BASED
NUMERICAL
METHOD FOR COMPRESSIBLE FLOWS
Governing equations in curvilinear coordinates
The governing equations given above can be expressed in a strong conservation form in
general non-orthogonal coordinates as
OJp
--
3
+
(
ov. E k) = 0 ,
(8)
O[FJg,kO¢]+jSe~,
{Jp(V • ek)ck} = -~k k
0¢ij
(9)
ot
0
~
Ot (Jpck ) +
with
E l -- E 2 X E 3
j
IE2 = E3XE1
'
j
'
E3 --- I~IXE2 .
j
'
0r
g]k=e].ek;
j _ O(Xl, X2,
Ej
0~i,
X3)
- (
,xE2) •
where 4~ can stand for the Cartesian velocity components, total enthalpy, turbulent kinetic
energy, mass concentration, etc. F is the effective diffusivity, and S, is the source term. Usually,
V - e k is defined as V k, a contravariant velocity component.
Discretization
The discretization of the above differential equations are carried out using a finite-volume
approach. First, the solution domain is divided into a finite number of discrete volumes or
"cells", where all variables are stored at their geometric centers. The equations are then
integrated over all the control volumes by using the Gaussian theorem. For example, the
integration of (8) leads to
p V - p°V°
At
~- Ge - Gw + G , - G, + Gh -- G , = 0 ,
(10)
where V is the volume of the cell, and the G's represent the mass fluxes over one of the control
volume surfaces. For example, Ge = (JrV~)e is the mass flux through the East face.
The discretization of the general convection-diffusion equation (9) can be carried out in a
similar fashion. The convective term can be discretized as
C = Ce -- Cw + Cn - Cs + Ch -- Ci = Ge4Je - Gwrkw + G, rkn - Gsck, + Gh4Jh -- Gt4)t,
(11)
Turbulent transonic flow simulation
473
with dze, for example, being the dz evaluated at the center of the East face. The diffusion term
can be split into a main diffusion (j = k) part and a cross diffusion (j ~ k) part which are given
by
O~t = j
k = 1, 2 or 3,
(12)
j ¢ k.
(13)
and
j
] r g j*
,
Without loss of generality, k = 1 and j = 2 are illustrated below. Integration of the D ~ term
over a control volume leads to
fff o.
=[
FJgl, 0_~] _ [Fjg,1 adz
O~lJe
(14)
~l]w,
By defining
D~I = FJg n,
(15)
we have
fff
jjj
D1MJ d~l d~2 d~:3 = (D~I + DeW)dze+ D~ldzE + De"~ ~bw.
(16)
Therefore, D e, the main diffusion coefficient needs to be evaluated at the faces of all control
volumes. For the cross-diffusion part, the integration of the D 21 term yields
fff
By
defining
fff
Dc21J d~l d~2 d~3 = [rjg21 O(J~]- [Fjg21 0~
t
D 21 = I~
-~Jg 21 and
assuming
Dc J dsc, d~2 d~=3= D21( dzN + •NE
(17t
O~2Je
dzne= l(dze + dze + dzN + dzNe), etc.,
- - {~S - - ~ ) S E )
--
we
have
DZw'(dzN + dzNw - dzs - dzsw).
(18)
With the above discretization process, the discretized dz equation at control volume P can be
assembled together in a linear equation form to give
Apc~p = Awdzw + AEdzE + As~bs + ANON "q-ALdzL + AHdztt
+ Aswdzsw + AseqbsE + ANWdzNW + ANEdzNE
"q-ALSdPLS q- ALNdzLN q- AHSq~HS -}- AHNdPHN
+ AWLdzWL + AWH~bWH + AEL¢EL + AEHdzEN + S ¢.
(19)
Mass flux evaluation
Since a non-staggered grid approach is used, a special procedure is required to prevent the
well-known checkerboarder problem from occurring. Following the practice of [12] and [13],
this is achieved by averaging the momentum equation at the cell faces and by relating the cell
face velocity directly to the local pressure gradient. For example, the discretized x-momentum
equation can be expressed at cell P as
ap +
p Up
anbUnb ~P \ O X / p
Y.G. LAI et al.
474
where nb refers to all neighboring cells. Dividing the above equations by a~, yields
(1 + cd;)u u = E
aub
nb ap
Unb -
d \OX/p + caput, + _,,
ap
(21)
Vp
dp. =--ft.
ap
(22)
where
c = ___p
At'
The above momentum equation is written for the P-cell, but it applies to every cell in the
calculation domain. Therefore, at cell face f, for example, we have
(l +cd~)uf = (~nb anubunb f-- d~,ox/f(OT-P-P]+ c d ~ u y - (S-~)/
(23)
If d~, (•,b (a~,b/a~,)Unb)y, (S~/a~,)I are obtained by averaging the same quantity from two
neighboring nodal points and if (21) is substituted into (23), the result is
(1 + cdy)uf = (1 + cd~)up + d~(OP~ -cd~,u~--~p(Op~ + cd~u~,
\3X/p
\3x/f
(24)
where the overbar denotes averaging from two neighboring nodal points. The above equation is
further reduced to
-
p
f
+ c d " ( u 7 - u,,),
(25)
with
d"-
dp__.
1 -
(26)
ca;
It is clear from the above equations that the cell face velocity is obtained from an average of
the two neighboring point velocity plus a second-order pressure correction and a second-order
previous time velocity correction. The pressure term serves as the mechanism to avoid the
checkerboarder problem and the previous time velocity serves as a mechanism to obtain a
time-independent steady solution.
The y-component and z-component velocities can be similarly evaluated and the contravariant velocity component at a cell face can thus be evaluated as
V~= uf . F~_~+ D . Fky + Wr . Fkz,
k=1,2,3,
(27)
with
Fk, = e k" i,
Fky = e ~" j,
Fkz = e ~" k,
k = 1, 2, 3.
(28)
E N H A N C E M E N T S OF T H E N U M E R I C A L M E T H O D
Extension of TVD schemes
In the previous section, the evaluation of fw....... i.h has not been discussed. Therefore, its
evaluation needs to be addressed. The focus here is on how to extend the TVD schemes used in
a density-based method to improve shock capturing. Without loss of generality, the evaluation
of ~e is discussed and a uniform grid is assumed. Traditionally, the pressure-based method has
used such numerical schemes as first-order upwind, central, 2nd order upwind, QUICK, etc.
Turbulent transonicflowsimulation
475
For example, the first-order upwind approach evaluated the using either the or the depending on
the flow direction at the East face. Mathematically, it can be expressed as
C~v*"= Ge the + the
2
IGelthe - th,,
2
(29)
This is one of the most diffusive schemes and, when used for supersonic flows, too much
smearing of the shocks will result.
A natural choice of a second-order scheme is the central difference scheme. Pure central
difference approach evaluates the by averaging the values at point P and E as
CCeN = Ge the
+ d,
P
2
(30)
It is widely known that pure central difference schemes could cause unphysical oscillations
despite the fact that it lowers numerical false diffusion. For most problems, some damping (or
artificial viscosity) is needed for stability. In practice, a central difference scheme with damping
is constructed as
C e = (1 - de)C Up + deCeC'V,
(31)
with 1 - de representing the fraction of upwind scheme used. That is,
Ce
Ge thP + (~E
2
=
(1 -- de)IGel tth______
he e
2
(32)
Comparing (32) with (30), it is clear why the scheme is called central difference plus damping.
de = 1 yields the pure central difference scheme, while de -- 0 results in an upwind scheme.
It is observed that the damping coefficient is constant for a central difference scheme. This
will find limited usefulness when compressible flows are calculated. Quite often, very large
damping ( [ 1 - d e ] - 1) is needed near shocks while very small damping is used in smooth
regions. Therefore, traditional schemes used in pressure-based method are not good for shock
capturing. However, it is widely known that shocks have been captured correctly with TVD
schemes using a density-based method. A close examination of these density-based TVD
schemes reveal that one of the key successes is the use of adaptive limiters. A TVD limiter has
the capability of adding more damping near shocks while using less or no damping in smooth
regions adaptively. This adaptive control of damping term is exactly what a pressure-based
method needs. In this study, three TVD-type limiters are used to evaluate the damping
coefficient de in (32). These are the MinMod limiter, the van Leer MUSCL scheme and the
Osher-Chakravarthy scheme. The MimMod limiter evaluates the damping coefficient de as
de = MinMod(1, ~e),
(33)
with
S e - v + th~,-v
~//e
--
and
U = sign(Ge).
(34)
6e - 6,,
The MinMod function is defined as: MinMod(a,/3) = sign(oe)max[0, min(lal,/3)]. This scheme
reduces to an upwind scheme if there exists local extrema such as found across shocks (Ye < 0),
or to a central difference scheme when ye -> 1, or to a second-order upwind scheme for
0<ye<l.
On the other hand, the Osher-Chakravarthy scheme yields a damping coefficient
-
1+77
de= 1 2 77MinMod(fl, "~e)q- T
MinMod(1,/3ye),
(35)
476
Y . G . LAI et al.
Upwind
Central difference
MinMod limiter
I
van Leer limiter
Osher-Chakravarthy limiter
r
Fig. 1. Mach number contours using pressure-based method and different TVD limiter schemes.
Turbulent transonic flowsimulation
477
Fig. 2. Mach number contour using density-based method with Roe's flux differencesplitting and
Osher-Chakravarthylimiter.
with n = ~ and /3 = ( 3 - 77)/(1- 77), while the van Leer MUSCL scheme gives the following
damping coefficient
]Tel + "ye
de - Ye----~
(1 +
"
(36)
In order to evaluate the performance of these TVD type schemes, a supersonic flow with a
free-stream Mach number of 1.65 over a 4% circular bump is calculated as an example. Figure
1 shows the calculated Mach number contours obtained from five different TVD schemes
incorporated into the pressure-based method. For comparison purpose, the same problem is
also calculated with a density-based method using Roe's flux difference splitting and the
Osher-Chakravarthy scheme. This result is shown in Fig. 2. It can be seen that the upwind
scheme is too diffusive and the central difference scheme causes oscillations. Overall, when
TVD schemes are used with a pressure-based method, they can also capture shocks correctly
compared to density-based methods. However, it should be pointed out that the use of TVD
limiters in the above does not mean that the schemes have TVD property. In essence, the TVD
concept is limited in the current pressure-based method relying on the assumption that all the
waves are traveling with the fluid convection speed.
On the pressure correction equation
Experience suggests that the pressure correction equation, resulted from the continuity
equation, is the most difficult equation to solve and is the dominant factor contributing to the
inefficiency of the pressure-based method. Two improvements are adopted in this study for
non-orthogonal grids. These are the SIMPLEC [21] algorithm and the implicit treatment of
non-orthogonal contributions. The first improvement is not discussed here because interested
readers can refer to [22] for details. A discussion of the second improvement is given below.
As discussed in the previous section, the velocity at the control volume faces involves local
pressure gradient. For a non-orthogonal grid this local pressure gradient has two sources of
contributions. One is derived from the main grid line direction (main term) and another is from
the cross direction (cross term). This means that the East face contravariant velocity increment,
for example, can be written as:
(Vie)'
= - D [Op"~
D (Op'~ + DfOP"~
~ O~ )e+ "0\ On ]e
C~ O~ )e'
(37)
where De, D,7, and D c are the main term and cross term coefficients, respectively. Most of the
previous studies have assumed that the D,~ and Dc terms are small and can be neglected. As a
result, a seven-point scheme (five-point for 2D) is obtained for the pressure correction
equation. However, our experience showed that this assumption is the major cause of slow
convergence for a non-orthogonal grid. In reality, the ratio, D J D e, is proportional to cot(a)
with a being the angle between ~ and 77. Therefore, the D n term can be bigger than D e for
478
Y.G. LAI et al.
non-orthogonal grid. In this paper, both D n and D e terms are retained in the pressure
correction equation which results in a 19-point scheme (nine-point for 2D).
S I M U L A T I O N S OF T R A N S O N I C T U R B U L E N T
AIRFOIL FLOWS
The pressure-based m e t h o d and the compressible two-equation k - e turbulence model [17]
are used to calculate turbulent transonic flows over N A C A 0012 and R A E 2822 airfoils with
and without shocks. At the same time, the B a l d w i n - L o m a x model [20] is also used to calculate
the same flows for comparison purposes. As a result, the relative merits and drawbacks of the
numerical scheme and the turbulence models can be assessed. In the following, the results for
the N A C A 0012 airfoil are presented and discussed first. This is followed by the R A E 2822
airfoil results.
N A C A 0012 airfoil
Two flow cases have been chosen for the N A C A 0 0 1 2 airfoil: one without shock and
separation; and another with a very strong shock and shock-induced flow separation. For the
first case, the following free-stream conditions are specified: Mach number, M~ = 0.7; Reynolds
number, Re~ = 9 x 106; geometric angle of attack, ag = 1.86. For the second case, Moo = 0.799,
Re= = 9 x 106, and ag = 2.86. The wind tunnel measurements of Harris [23] are available for
comparison. Since the calculations are carried out in free air, the angle of attack must be
corrected from the geometric angle of attack due to the presence of wind tunnel interference.
The corrections suggested by Harris [23] are ac = 1.49 and ac = 2.26 for the two cases,
respectively, and are used in our computations. The same corrected angles of attack were
widely used by other researchers in the past [24, 25].
The computational grid used for the calculation has been chosen very carefully. In
preliminary tests, both O- and C-grids have been tried and they essentially give the same
solutions using mixing-length models. However, the O-grid is much slower to converge when
the k - e model is used. Therefore, only the C-grid is chosen. With some grid refinement tests
and varying locations of the far-field boundary, a 270 x 70 grid is finally used with far-field
boundary placed at about 15 chord lengths away from the airfoil surface.
Figure 3 displays part of the C-grid used for the present computation. This grid is generated
by a combination of algebraic and elliptic procedures with some limited degree of orthogonality
control. Basically, there are 49 grids in the wake region and 172 grids on the airfoil surface.
Finer grids have been used in the trailing and leading edges of the airfoil. The mesh distribution
in the normal-to-the-wall direction is generated with a heavy power stacking to achieve the
required fine grid near a wall for a low-Reynolds-number turbulence model. As a result, the
grid spacing at the airfoil surface is around Y / C = 8 x 10 - 6 which results in Y+ < 2.0. Here, X is
taken to be the stream coordinate, Y is the wall normal coordinate, C is the airfoil chord
Fig. 3. Part of the grid layout used for the NACA 0012 airfoil calculation.
Turbulent transonic flow simulation
10-
479
V~
~
~
(a)
05-
- Cp
00-
V
~' ~
~
--:
,i
~
k-e model
. . . . : B-L model
I
I
02
0
I
04
06
I
08
I
10
x/e
15-
Co)
VV
10-
VV
"i
05-
~,VV
VV
- Cp
00-
A:
Experimentaldata [23]
--:
k-e model
. . . . : B-L model
-i
00
I
02
I
l
04
06
I
~8
I
10
~c
Fig. 4. Comparisons of the calculated NACA 0012 surface pressure distributions with measurements
for two different conditions: (a) M~ = 0.7; Re~ = 9 × 106; O~g = 1.86; (b) M~ = 0.799, Re~ = 9× 106,
ag = 2.86.
length, Y+ = Yurlvw and ur is the wall friction velocity. The above mesh system is similar to the
choices of m a n y other investigators [24, 25].
Figure 4 shows the comparison of the pressure coefficient distribution between computations
and m e a s u r e m e n t s for the two flow cases. Figure 5 shows the corresponding pressure contours
obtained f r o m the k - e model computations. It is clear that there are no flow separation and
shock waves for the first case while a strong shock is present for the second case. A close
examination also reveals that there is a flow separation right behind the shock. It is obvious
that the model predictions of the pressure coefficients are quite satisfactory for the first case,
but deviate f r o m the m e a s u r e m e n t s in the neighborhood of strong shocks for the second case.
Actually, Hoist [26] found that it is a c o m m o n deficiency of a full N a v i e r - S t o k e s m e t h o d used
to calculate the second case. In general, both the k - e and the B a l d w i n - L o m a x model yield
essentially similar results for the prediction of Cp.
480
Y.G.
L A I et al.
(a)
(b)
Fig. 5. P r e s s u r e c o n t o u r p l o t s for t h e N A C A 0 0 1 2 airfoil for t w o d i f f e r e n t c o n d i t i o n s : (a) M~ = 0.7;
R e ~ = 9 x 106; O~g = 1.86; (b) M~ = 0.799, Re~ = 9 × 106, a~8 = 2.86.
R A E 2822 airfoil
The second airfoil analyzed is the RAE2822, which has been extensively investigated
experimentally by Cook et al. [27] and computationally by many aerodynamists. This airfoil
flow, together with the wind tunnel measurements of Cook et al. [27], were recommended as a
standard test case in the 1980-1981 Stanford conference [28]. The RAE2822 airfoil is
moderately rear-loaded with a 12.1% thick supercritical section. Case 6 has been chosen for the
present computation which involves supercritical flow with shock waves. The chosen case has a
free-stream Mach number of 0.725, a geometric angle of attack of 2.92 and a Reynolds number
of 6.5 million. It is recommended that the corrected angle of attack (o~c) be chosen such that the
lift coefficient matches the measurements. Due to the limited time scope of the present study,
ac = 2.6 is picked for the computation. It is interesting to note that a corrected angle of attack
range from 2.3 to 2.79 was found by different computational groups in the 1987 Viscous
Transonic Airfoil Workshop (results were presented at the A I A A 25th Aerospace Sciences
Meeting).
A 270 × 70 C-grid is used to perform this calculation. This grid is similarly generated and
distributed as the NACA0012 airfoil grid. Near-wall Y+ is around 1 and 172 grids are
distributed on the airfoil surface with finer clustering near the leading and trailing edges. Part
Turbulent transonic flow simulation
481
(/////I//////.~/~
I//~rl////////
/ /
I I I II1111IIIIIII I // ~\ \ \ ~ ~",,\~\\\\\" I~\\X\\\\\"
x"--,"x
Fig. 6. Part of the grid layout used for the RAE 2822 airfoil calculation.
15-
10-
05-
e 0~
A:
-0
Experimentaldata [271
-1
. . . . : B-L model
-1
i
®
02
I
04
I
06
i
08
I
10
X/C
Fig. 7. Comparisons of the calculated R A E 2822 surface pressure distributions with measurements.
Fig. 8. Pressure contour plot for the RAE 2822 airfiol.
of the grid is depicted in Fig. 6. Figure 7 shows the pressure coefficient distributions on the
lower and u p p e r airfoil surfaces f r o m the k - e model, the B a l d w i n - L o m a x model, and the
wind-tunnel data. The pressure contours are shown in Fig. 8. It is seen that there is a shock
around X / C = 0.55 and a quite strong compression at the leading edge. Again the calculated
differences between the two turbulence models are small. It is obvious that the computations
underpredict Cp on the u p p e r surface of the airfoil. This discrepancy is due partly to the small
angle of attack used in the computation.
482
Y . G . LAI et al.
CONCLUSIONS
A pressure-based numerical method and a recently developed near-wall turbulence model
are used to calculate turbulent transonic flows over airfoils. Several TVD schemes are first
extended to the pressure-based method and their shock capturing capability is examined using
a supersonic flow over a 4% bump. It is found that the van Leer MUSCL scheme and the
Osher-Chakravarthy scheme are quite effective in suppressing oscillations near shocks and
provide results similar to those obtained using a density-based method.
Two airfoils (NACA0012 and RAE2822) have been chosen for the turbulence model
validation studies. The calculated results are compared with wind-tunnel measurements and
with the Baldwin-Lomax model calculations. For the cases tested, both the Baldwin-Lomax
[20] model and the near-wall, two-equation k - e model are capable of predicting turbulent
transonic flows over the airfoils, particularly in the prediction of the pressure distributions on
the surface. However, both models fail to predict the location of the shock correctly when a
strong shock and a shock-induced flow separation are present. Further investigations will be
carried out to identify the reasons for this deficiency. For separated flows, it is widely known
that the mixing-length type models are deficient while the two-equation model could provide a
more accurate answer.
Even though the k - e model does not display an apparent advantage over that of the
Baldwin-Lomax model, these calculations fully demonstrate that more advanced turbulence
models could be incorporated into a well developed pressure-based numerical method to
calculate turbulent compressible external aerodynamic flows. Furthermore, the calculations
show that the ability of the pressure-based numerical method to resolve shocks and other
aerodynamic discontinuities is quite good and its performance rivals those of density-based
methods. Consequently, an alternative to density-based numerical methods is now available.
This means that advanced turbulence models capable of handling flows with separation and
recirculation could be implemented with ease into the numerical scheme and a robust
numerical code is now available for the calculation of external as well as internal turbulent
aerodynamic flows with and without shocks.
Acknowledgements--The authors (YGL and AJP) wish to acknowledge the help they have received from their
colleagues at CFDRC; in particular, the help from Drs H. Q. Yang, M. M. Athavale, Y. Jiang and Z. J. Wang, and the
overall guidance of Dr A. K. Singhal. Typing of the manuscript by Ms J. Swann is very much appreciated. The second
author (RMCS) wishes to acknowledge the support he received from NASA Grant NAG-l-1080, which was monitored
by Dr T. B. Gatski of NASA Langley, Hampton, VA 23665.
REFERENCES
[1] P. L. ROE, A. Rev. Fluid Mech. 18, 337 (1986).
[2] H. C. YEE, J. Computl Phys. 68, 151 (1987).
[3] S. R. C H A K R A V A R T H Y and S. OSHER, A new class of high accuracy TVD schemes for hyperbolic
conservation laws, AIAA-85-0363 (1985).
[4] B. VAN LEER, J. Computl Phys. 32, 101 (1979).
[5] P. COLELLA and P. R. WOODWARD, J. Computl Phys. 54, 174 (1984).
[6] A. HARTEN and S. OSHER, SIAM J. Numer. Anal. 24, 297 (1987).
[7] F. H. HARLOW and J. E. WELCH, Phys. Fluids 8, 2182 (1965).
[8] S. V. PATANKAR, Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York (1980).
[9] H. Q. YANG and A. J. PRZEKWAS, Pressure-based high-order TVD methodology for dynamic stall simulation,
AIAA-93-0680 (1993).
[10] J. J. McGUIRK and G. J. PAGE, A I A A J. 2,8, 1752 (1990).
[11] B. R. HUTCHINSON and G. D. RAITHBY, Numer. Heat Transfer 9, 511 (1986).
[12] C. M. RHIE and W. L. CHOW, A I A A J. 21, 151 (1983).
[13] M. PERIC, R. KESSLER and G. SCHRUERER, Computers and Fluids 16, 389 (1988).
[14] R. I. ISSA, J. Computl Phys. 62, 40 (1986).
[15] W. SHYY and M. B. BRAATEN, Application of a generalized pressure correction algorithm for flows in
complicated geometries. In Advances and Applications in Computational Fluid Dynamics (Edited by O.
BAYSAL), pp. 109-119. ASME, New York (1988).
Turbulent transonic flow simulation
483
[16] Y. G. LAI, Y. JIANG and A. J. PRZEKWAS, An implicit multi-domain approach for the solution of
Navier-Stokes equations in body-fitted coordinate grids, AIAA-93-0541 (1993).
[17] H. S. ZHANG, R. M. C. SO, C. G. SPEZIALE and Y. G. LAI, A1AA J. 31, 196 (1993).
[18] D. C. WILCOX, AIAA J. 26, 1299 (1988).
[19] T. J. COAKELY and P. G. HUANG, Turbulence modeling for high speed flows, AIAA-92-0436 (1992).
[20] B. S. BALDWIN and H. LOMAX, Thin layer approximation and algebraic model for separated turbulent flows,
AIAA-78-257 (1978).
[21] J. P. VAN DOORMAL and G. D. RAITHBY, Numer. Heat Transfer 7) 147 (1987).
[22] A. J. PRZEKWAS, Y. G. LAI, A. KRISHNAN, R. K. AVVA and M. G. GIRIDHARAN, Combustion chamber
analysis code: volume 1: theory manual, NASA CR-1993 (1993).
[23] C. D. HARRIS, Two-dimensional aerodynamic characteristic of the NACA0012 airfoil in the Langley 8-foot
transonic pressure tunnel, NASA TM-81927 (1981).
[24] T. J. COAKELY, Numerical simulation of viscous transonic airfoil flows, AIAA-87-0416 (1987).
[25] C. M. MAKSYMIUK and T. H. PULLIAM, Viscous transonic airfoil workshop results using ARC2D,
AIAA-87-0415 (1987).
[26] T. L. HOLST, J. Aircraft 25, 1073 (1988).
[27] P. H. COOK, M. A. McDONALD and M. C. P. FIRMIN, Aerofoil RAE2822 pressure distributions, and
boundary layer and wake measurements, AGARD Advisory Report No. 138 (1979).
[28] S. J. KLINE, B. J. CANTWELL and G. M. LILLEY, The 1980-81 AFOSR-HTTM-Stanford Conference on
Complex Turbulent Flows: Comparison of Computation and Experiment 2, Thermoscience Division, Stanford
Univ. Stanford, Calif. (1981).
(Received 27 June 1994; accepted 18 July 1994)