- Email: [email protected]
Turbulent transonic flow simulation using a pressure-based method
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
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)
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
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)