NASA Technical Reports Server (NTRS) 20040086843: A Pseudo-Temporal Multi-Grid Relaxation Scheme for Solving the Parabolized Navier-Stokes Equations

16 Pages·1.3 MB·English

by  NASA Technical Reports Server (NTRS)

AIAA 99-3360
A Pseudo-Temporal Multi-Grid
Relaxation Scheme for Solving the
Parabolized Navier-Stokes
Equations
J. A. White
NASALangleyResearchCenter
Hampton,VA23681
J. H. Morrison
AnalyticalServicesandMaterials,Inc.
107ResearchDrive
Hampton,VA23666
14th AIAA Computational Fluid
Dynamics Conference
June 28–July 1, 1999/Norfolk, VA
Forpermissiontocopyorrepublish,contacttheAmericanInstituteofAeronauticsandAstronautics
1801AlexanderBellDrive,Suite500,Reston,VA22091
A Pseudo-Temporal Multi-Grid Relaxation
Scheme for Solving the Parabolized
Navier-Stokes Equations
(cid:3)
J. A. White
NASA Langley Research Center, Hampton, VA 23681
y
J. H. Morrison
Analytical Services and Materials, Inc., Hampton, VA 23666
A multi-grid, (cid:13)ux-di(cid:11)erence-split, (cid:12)nite-volume code, VULCAN, is presented for
solving the elliptic and parabolized form of the equations governing three-dimensional,
turbulent, calorically perfect and non-equilibrium chemically reacting (cid:13)ows. The space
marching algorithms developed to improve convergence rate and or reduce computa-
tional cost are emphasized. The algorithms presented are extensions to the class of
implicitpseudo-timeiterative,upwindspace-marchingschemes. Afullapproximatestor-
age, full multi-grid schemeis also describedwhich is used to accelerate the convergence
of a Gauss-Seidelrelaxation method. The multi-grid algorithm is shown to signi(cid:12)cantly
improveconvergence on high aspectratio grids.
Introduction schemesthatsolvetheunsteady formofthe Eulerand
2,18,19
PNS equations have also been employed using
Scramjet engines and nozzles employed by hy-
relaxation schemes that iterate at each spatialstep in
personic vehicles operate in a viscous super-
pseudo-time to reach a steady-state solution.
sonic/hypersonic (cid:13)ow regime. Such conditions allow
the useofspace-marchingsolutionalgorithmswithout
A large number of the problems of interest can in-
adversely a(cid:11)ecting the (cid:12)delity ofthe calculationwhen
clude subsonic and/or separated (cid:13)ow regions, which
the (cid:13)ow is steady and unseparated in the marching
requires that the space marching algorithmbe imple-
direction. Numerous algorithms have been developed
mented within the broader context of a code capa-
for solving the Euler, Navier-Stokes, and Parabolized
ble of solving the Navier-Stokes equations. VULCAN
Navier-Stokes (PNS) equations (e.g. Refs. 1{7) for a
(Viscous Upwind aLgorithm for Complex (cid:13)ow ANal-
thermally and calorically perfect gas. Extensions of
ysis) is a cell-centered (cid:12)nite-volume, structured grid,
space marching algorithms have also been made to
multi-blockcode whichsolves the equations governing
handle both equilibrium and non-equilibrium e(cid:11)ects
inviscid and viscous (cid:13)ow of a calorically perfect gas
(e.g. Refs. 8{15). In addition,the PNS equations can
or of an arbitrary mixture of thermally perfect gases
besolvedmuchmoree(cid:14)cientlythantheNavier-Stokes
undergoingnon-equilibriumchemicalreactions. VUL-
equations, which makes the PNS equations very at-
16 CAN allows the (cid:13)ow domain to be decomposed into
tractive foruse indesign andoptimizationmethods.
regions in which the most suitable algorithm (ellip-
Space marching methods can be broken down al-
tic or marching) can be utilized. The inviscid (cid:13)uxes
gorithmically into two groups, steady-state equation 20
are computed using the MUSCL scheme with either
solvers and pseudo-time iterative solvers. Steady- 21
the approximate Riemann solver of Roe or the low
state equation solvers, as the name implies, solve 22
dissipation (cid:13)ux splitting scheme of Edwards. The
the steady-state form of the Euler and PNS equa-
3,5,10,17 4,11,14 viscous (cid:13)uxes can be evaluated either with or with-
tions. Explicit and implicit steady-state
out cross derivative contributions. VULCAN solves
space-marching schemes have been used with consid-
elliptic (cid:13)ows by marching the unsteady form of gov-
erable success. Pseudo-time iterative space-marching
erning equations in time. Hyperbolic (cid:13)ows (e.g. the
(cid:3) Eulerequationsinsupersonic(cid:13)ow)andparabolic(cid:13)ows
Aerospace Engineer, Hypersonic Airbreathing Propulsion
Branch, Aero- and Gas-DynamicsDivision,Research Technol- (e.g. the parabolized Navier-Stokes equations) are
ogyGroup. solvedbyspace marchingwhileiteratingthe unsteady
y
SeniorResearchScientist,SeniorMemberAIAA. equations in pseudo-time to a steady state solution.
Copyright(cid:13)c 1999bytheAmericanInstituteofAeronauticsand
Astronautics, Inc. No copyright is asserted in the United States The pseudo-time iterative approach was chosen be-
under Title 17, U.S. Code. The U.S. Government has a royalty- cause it allows the time marching solution algorithms
freelicensetoexerciseallrightsunderthecopyrightclaimedherein
and convergence acceleration techniques developed to
for Governmental Purposes. All other rights are reserved by the
copyrightowner. solve the elliptic Euler equation and FNS equations
1 of15
AmericanInstituteof Aeronautics and Astronautics Paper99{3360
to be used to reduce the computational cost of the averaged variables (indicated as (cid:30)) as
space marching scheme. The explicit PNS algorithm
described in a previous paper,23 while extremely ef- (cid:26)f~i !_i
(cid:12)cient for inviscid and moderate Reynolds number 2 ... 3 2 ... 3
(cid:13)ows, was found to be too computationallyexpensive
wwlbaiahrligsetenydnescusooimlgnvnsbinteerdgraiohntfoitgssh.palaRlcoTeewhymentoahclreodcdsheeilnnliugsptmtsritubceceptarsun(cid:13)rrdeoewqsoupsfiardcVeuedUembLtyoCarsAtcthhaNe-- Q = J1 6666666 (cid:26)f(cid:26)(cid:26)(cid:26)~Nwuv~~~cs 7777777 S = 6666666 !_N000cs 7777777 (4)
ing schemes to share numerical algorithms developed 66 (cid:26)E~ 77 66 0 77
to solve and/or accelerate the convergence of the gov- 66 (cid:26)k~ 77 66 Sk~ 77
erning equations to steady state. The code has in- 666 (cid:26)~l 777 66 S~l 77
corporated within it four time marching schemes: a 66 (cid:26)g~ 77 66 Sg~ 77
multi-stage Runge-Kutta scheme, a diagonalized ap- 66 (cid:26)F~ 77 664 SF~ 775
proximate factorization scheme, a block approximate 4 5
factorizationscheme,andadiagonallydominantalter- (cid:26)f~iU
nating direction implicit scheme with a factorization ..
2 . 3
error correction sub-iteration. Each of these schemes
(cid:26)f~NcsU
can be used as the smoothing algorithm in a full ap- 6 7
6 (cid:26)u~U +(cid:24)^xp 7
proximatestorage(FAS)fullmulti-grid(FMG)scheme 6 7
andtherebyacceleratetheconvergencetosteadystate. (cid:24) 66 (cid:26)v~U +(cid:24)^yp 77
E = jrJ j 66 (cid:26)w~U +(cid:24)^zp 77 (5)
66 (cid:26)H~U 77
6 7
6 (cid:26)k~U 7
6 7
6 (cid:26)~lU 7
Governing Equations 6 7
6 (cid:26)g~U 7
6 7
The unsteady three-dimensional form of the Favre 6 (cid:26)~U 7
6 F 7
averaged Navier-Stokes and species transport equa- 4 5
(cid:26)f~iV
tionscanbewritteninnonorthogonalcurvilinearform
.
as 2 .. 3
(cid:26)f~NcsV
6 7
6 (cid:26)u~V +(cid:17)^xp 7
Qt + (E(cid:0)Ev)(cid:24) + (F(cid:0)Fv)(cid:17) + (G(cid:0)Gv)(cid:16) = S (1) (cid:17) 666 (cid:26)v~V +(cid:17)^yp 777
F = jr j 6 (cid:26)w~V +(cid:17)^zp 7 (6)
J 6 7
6 (cid:26)H~V 7
wherethesubscripts(cid:24),(cid:17),and(cid:16) denotepartialdi(cid:11)eren- 6 7
tiation with respect to the nonorthogonal curvilinear 66 (cid:26)k~V 77
coordinates and the subscript t denotes partial dif- 66 (cid:26)~lV 77
6 7
ferentiation with respect to time. The contravariant 6 (cid:26)g~V 7
6 7
components of the curvilinear nonorthogonal coordi- 6 (cid:26)~V 7
6 F 7
nate system are 4 5
(cid:26)f~iW
.
.
2 . 3
(cid:24) =(cid:24)(x;y;z); (cid:17) =(cid:17)(x;y;z); (cid:16) =(cid:16)(x;y;z) (2) (cid:26)f~NcsW
6 7
6 (cid:26)u~W +(cid:16)^xp 7
6 7
and (cid:16) 66 (cid:26)v~W +(cid:16)^yp 77
G = jrJ j 66 (cid:26)w~W +(cid:16)^zp 77 (7)
@((cid:24);(cid:17);(cid:16)) 66 (cid:26)H~W 77
J = (3) 66 (cid:26)k~W 77
(cid:12)(cid:12)@(x;y;z)(cid:12)(cid:12) 66 (cid:26)~lW 77
(cid:12) (cid:12) 6 7
(cid:12) (cid:12) 6 (cid:26)g~W 7
6 7
where (cid:24)(cid:11), (cid:17)(cid:11), and (cid:16)(cid:11) ((cid:11) = x, y, or z) are the compo- 6 (cid:26)~W 7
6 F 7
nents ofthe cell face area normalvector, while (cid:24)^(cid:11), (cid:24)^(cid:11), 4 5
where (cid:26) is the mixturedensity; f~i is the massfraction
and(cid:24)^(cid:11) arethecomponentsofthecellfaceunitnormal. th
and!_iisthespeciesproductionrateofthei chemical
The Q, E, F, G, and S vectors are de(cid:12)ned using species; Ncs is the total number of chemical species;
Favreaveragedvariables(indicatedas(cid:30)~)andReynolds p is the static pressure; u~, v~ and w~ are the Cartesian
2 of15
AmericanInstituteof Aeronautics and Astronautics Paper99{3360
00 00
velocitycomponents;andE~andH~ arethetotalenergy (cid:26) udjifi+udjifi (cid:16)^xj
and total enthalpy respectively. The variables k~ and ~l 2 h .. i 3
aretheturbulentkineticenergyandaturbulentlength f e . g00 00
scchaelmeiesqturyatiinotnervaacrtiiaobnleva(eriiathbelers!~g~oarn(cid:15)~d).~Tahreettuhrebeunleenrgcye 666 (cid:26) udjNcsfNcs+udjNcsfNcs (cid:16)^00xj00 777
00 00 F 6 (cid:16)^x(cid:28)hxx+(cid:16)^y(cid:28)xy+(cid:16)^z(cid:28)xz (cid:16)^xj(cid:26)iu uj 7
varNi=iac1snfcie00fei00.e T,haendcotnhterasvuamrianotf tvheelocsipteiecsiesarveadriea(cid:12)nnceeds Gv = 6666 (cid:16)^y(cid:28)xfy+(cid:16)^y(cid:28)fyy+(cid:16)^z(cid:28)yzg(cid:0)(cid:0)(cid:16)^xj(cid:26)vg0000u0j000 7777 (11)
as g 66 (cid:16)^x(cid:28)xz+(cid:16)^y(cid:28)yz+(cid:16)^z(cid:28)zz (cid:16)^xj(cid:26)w uj 77
P g 66 (cid:16)^xev+(cid:16)^yfv+(cid:0)(cid:16)^zgv g 77
66 (cid:16)^xDk~x+(cid:16)^yDk~y+(cid:16)^zDk~zg 77
6 7
U =(cid:24)^xu~+(cid:24)^yv~+(cid:24)^zw~ 66 (cid:16)^xD~lx+(cid:16)^yD~ly+(cid:16)^zD~lz 77
66 (cid:16)^xDg~x+(cid:16)^yDg~y+(cid:16)^zDg~z 77
V =(cid:17)^xu~+(cid:17)^yv~+(cid:17)^zw~ (8) 6 7
6 (cid:16)^xDF~x+(cid:16)^yDF~y+(cid:16)^zDF~z 7
W =(cid:16)^xu~+(cid:16)^yv~+(cid:16)^zw~ 6 7
4 5
The components of the species moleculardi(cid:11)usion ve-
th th
locity of the i specie in the j direction, udjifi, are
expressed using a Fickian model and their turbulent
00 00
counterparts, udjifi , are closed using a grafdienet dif-
Theviscous(cid:13)uxvectors Ev,Fv,andGv arewritten fusion model as
as g
00 00 (cid:22) (cid:22)t @f~i
(cid:26) udjifi+udjifi = + (12)
(cid:0) Sc Sct @xj
(cid:18) (cid:19)
h i
(cid:26) udjifi+u0d0jifi00 (cid:24)^xj where f(cid:22) aned (cid:22)t agre the molecular and turbulent vis-
cosity and Sc and Sct are the laminar and turbulent
2 h .. i 3
f e . g00 00 Schmidt numbers. The Cartesian components of the
66 (cid:26) udjNcsfNcs+udjNcsfNcs (cid:24)^xj 77 viscous energy (cid:13)ux are
6 00 00 7
6666 (cid:24)(cid:24)^^xy(cid:28)(cid:28)hxxfxy++(cid:24)(cid:24)^^yy(cid:28)(cid:28)fxyyy++(cid:24)(cid:24)^^zz(cid:28)(cid:28)xyzzg(cid:0)(cid:24)(cid:24)^^xxjj(cid:26)(cid:26)ivu00uu0jj0 7777 fevv == uu~~(cid:28)(cid:28)xxxy++vv~~(cid:28)(cid:28)xyyy++ww~~(cid:28)(cid:28)yxzz(cid:0)qqyx++DDk~k~yx (13)
Ev = 6 (cid:0) g00 00 7 (9) (cid:0)
66 (cid:24)^x(cid:28)xz +(cid:24)^y(cid:28)yz +(cid:24)^z(cid:28)zz (cid:24)^xj(cid:26)w uj 77 gv = u~(cid:28)xz +v~(cid:28)yz+w~(cid:28)zz qz+Dk~z
66 (cid:24)^xev+(cid:24)^yfv+(cid:0)(cid:24)^zgv g 77 (cid:0)
66 (cid:24)^xDk~x+(cid:24)^yDk~y +(cid:24)^zDk~zg 77 where(cid:27)kistheturbulentPrandtlnumberforthedi(cid:11)u-
6 7 sionofturbulentkineticenergy. TheCartesianviscous
66 (cid:24)^xD~lx+(cid:24)^yD~ly+(cid:24)^zD~lz 77 stresses are de(cid:12)ned as
66 (cid:24)^xDg~x+(cid:24)^yDg~y+(cid:24)^zDg~z 77
6 7 @u~ 2 @u~ @v~ @w~
6 (cid:24)^xDF~x+(cid:24)^yDF~y+(cid:24)^zDF~z 7 (cid:28)xx = (cid:22) 2 + +
6 7 @x (cid:0) 3 @x @y @z
4 5 (cid:20) (cid:18) (cid:19)(cid:21)
@v~ 2 @u~ @v~ @w~
(cid:28)yy = (cid:22) 2 + +
@y (cid:0) 3 @x @y @z
(cid:20) (cid:18) (cid:19)(cid:21)
@w~ 2 @u~ @v~ @w~
(cid:28)zz = (cid:22) 2 + + (14)
@z (cid:0) 3 @x @y @z
(cid:20) (cid:18) (cid:19)(cid:21)
@u~ @v~
00 00
(cid:28)xy = (cid:22) +
(cid:26) udjifi+udjifi (cid:17)^xj @y @x
(cid:18) (cid:19)
2 h .. i 3 @u~ @w~
6 f e . g00 00 7 (cid:28)xz = (cid:22) @z + @x
6 (cid:26) udjNcsfNcs+udjNcsfNcs (cid:17)^xj 7 (cid:18) (cid:19)
6 00 00 7 @v~ @w~
666 (cid:17)^x(cid:28)hxfx+(cid:17)^y(cid:28)fxy+(cid:17)^z(cid:28)xzg(cid:0)(cid:17)^xj(cid:26)iu00u0j0 777 (cid:28)yz = (cid:22)(cid:18)@z + @y(cid:19)
Fv = 666 (cid:17)^y(cid:28)xy+(cid:17)^y(cid:28)yy+(cid:17)^z(cid:28)yz (cid:0)(cid:17)^xj(cid:26)vg00uj00 777 (10) and the total heat (cid:13)uxes are given as
6 (cid:17)^x(cid:28)xz+(cid:17)^y(cid:28)yz +(cid:17)^z(cid:28)zz (cid:17)^xj(cid:26)w uj 7
6 (cid:0) g 7 Ncs
6 (cid:17)^xev+(cid:17)^yfv+(cid:17)^zgv 7 (cid:22)t @T~ (cid:22) (cid:22)t
66666 (cid:17)^(cid:17)^xxDDk~~lxx++(cid:17)^(cid:17)^yyDDk~~lyy++(cid:17)^(cid:17)^zzDD~lk~zzg 77777 qx =(cid:0)(cid:18)(cid:21)+ Prt(cid:19) @x (cid:0) (cid:26)(cid:18)Sc + Sct(cid:19)XiN=c1s~hif~ix
66666 (cid:17)^(cid:17)^xxDDF~g~xx++(cid:17)^(cid:17)^yyDDFg~~yy++(cid:17)^(cid:17)^zzDDg~F~zz 77777 qy =(cid:0)(cid:18)(cid:21)+ P(cid:22)rtt(cid:19)@@Ty~ (cid:0) (cid:26)(cid:18)S(cid:22)c + S(cid:22)ctt(cid:19)Xi=1~hif~iy (15)
4 5 3 of15
AmericanInstituteof Aeronautics and Astronautics Paper99{3360
Ncs 24
(cid:22)t @T~ (cid:22) (cid:22)t Menter areavailableintwoforms;thebaselinemodel
qz =(cid:0) (cid:21)+ Prt @z (cid:0) (cid:26) Sc + Sct h~if~iz and the Shear Stress Transport (SST) model. A cor-
(cid:18) (cid:19) (cid:18) (cid:19)Xi=1 rection to each model for compressibility e(cid:11)ects has
where (cid:21) is the gas mixturemolecularthermalconduc- alsobeen implementedusingthe compressibledissipa-
tivity and Prt is the turbulent Prandtl number. The tionmodelofSarkaretal.29forthek~ (cid:15)~basedmodels
(cid:0) 28
turbulence modeling equation di(cid:11)usion terms D(cid:30)~x, and the compressibilitymodelsuggested by Wilcox,
D(cid:30)~y, and D(cid:30)~z, where (cid:30) is either k~, (cid:15)~, !~, g~, or F~, wlahgigcihnginfcuonrcptoiorant,e30s Sfoarrkaalrl’osfmthodeek~las!w~ebllaasesdZmemoadne’lss
are modeled as (cid:0)
(including Menter’s).
(cid:22) (cid:22)t
D(cid:30)~x =(cid:18)Prl(cid:30) + Prt(cid:30)(cid:19)(cid:30)~x NeTahreWnaelalrTwuarlblubleehnacveioMroodfetlhse k~ ~(cid:15) two equation
D(cid:30)~y = P(cid:22)rl(cid:30) + P(cid:22)rtt(cid:30) (cid:30)~y (16) mlowodRelseyinsocldonstnroulmlebdetrhmrooudgih(cid:12)ctahteioinns(cid:0)troofdAucbtiido.n31ofTthhee
(cid:18) (cid:19)
(cid:22) (cid:22)t near wall behavior of the k~ !~ two equation models
D(cid:30)~z = Prl(cid:30) + Prt(cid:30) (cid:30)~z is treated in either a high R(cid:0)eynolds number manner
(cid:18) (cid:19) by integrating the high Reynolds number form of the
Turbulence Models equations to the wall without any special wall treat-
00 00
TheReynoldsstress term((cid:26)ui uj )modelisdepen- ment or by integrating the low Reynolds form of the
dent on the type of turbulence model selected. Cur- Wilcoxk~ !~ two equation model28 to the wall.
(cid:0)
rently,onlytwoequationmodelsgareavailableinVUL- If integration to the wall is not feasible for a
CAN. However, the two equation model implemen- givenproblem,thewallmatchingfunctionapproachof
28
tation allows (cid:13)exibility when modeling the Reynolds Wilcox can be used. This approach models the wall
stresses. If one of the eddy viscosity based two equa- shearstress andheattransferbyenforcingacompress-
tion turbulence models is selected, the Boussinesq ap- ible law of the wall which includes additional terms
32
proximation is made where the Reynolds stresses are that model streamwise pressure gradients. Wilcox
modeled as showed that inclusion of the pressure gradient terms
in the wall matching function dramatically improved
00 00 @u~i @u~j 2@u~k 2
(cid:26)ui uj =(cid:22)t + (cid:14)ij (cid:26)k~(cid:14)ij predictions of shear stress. The wall matching func-
(cid:0) @xj @xi (cid:0) 3@xk (cid:0) 3 +
(cid:20)(cid:18) (cid:19) (cid:21) tion has been successfully used on grids where the y
(17)
g (=p(cid:28)w=(cid:26); (cid:28)w is the wallshear stress) varied between
If an explicit algebraic Reynolds stress two equa-
0.5and150.0andproduced wallskinfrictionandheat
tion turbulence model is selected, then the Reynolds
transfer predictions that agreed well with results ob-
stresses are modeled as
tained on grids that were integrated to the wall for
+
(cid:26)ui00uj00 = 2(cid:26)k~(cid:14)ij +2(cid:22)(cid:3)t Sij 1Skk(cid:14)ij y <1:0. Thisisaveryusefulcharacteristicwhichsig-
(cid:0) (cid:0)3 (cid:0) 3 ni(cid:12)cantly reduces (cid:13)ow solution sensitivity to the grid
(cid:18) (cid:19)
(cid:3)(cid:3) (cid:11)4 and greatly reduces the required amount of a priori
g +2(cid:22)t (SikWkj+SjkWki) +
!~ knowledgeabout the range ofy when generating the
(cid:3)(cid:3) (cid:11)5h 1 i computationalgrid.
2(cid:22)t SikSkj SklSkl(cid:14)ij (18)
(cid:0) !~ (cid:0) 3 Pressure and Thermodynamic Models
(cid:20) (cid:18) (cid:19)(cid:21)
where VULCAN was written to modelthe working gas ei-
1 @u~i @u~j ther as a single component calorically perfect gas or
Sij = + (19)
2 @xj @xi as an arbitrary mixture of thermally perfect gases.
(cid:18) (cid:19)
For a single component calorically perfect gas, the
1 @u~i @u~j
Wij = (20) solutionvectorofconserved variablesinEq.(4)iscon-
2 @xj (cid:0) @xi
(cid:18) (cid:19) tracted to contain only a single continuity equation
and !~ = ~(cid:15)=k~. The de(cid:12)nitions for (cid:22)(cid:3)t, (cid:22)(cid:3)t(cid:3), (cid:11)4, and (cid:11)5 (cid:26) (with !_ = 0), three momentum equations (cid:26)u~, (cid:26)v~,
can be found in Abid et al.25 and (cid:26)w~, the energy equation (cid:26)E~ and the turbulence
The two equation turbulence models in VULCAN equationsfor(cid:26)k~ and(cid:26)~l. The pressure, staticenthalpy,
can be categorized as k~ ~(cid:15) based models, Wilcox’s28 totalenthalpy,andtotalenergyforasinglecomponent
k~ !~ based modelsandM(cid:0)enter’s24 k~ !~ models. The caloricallyperfect gas are
k~(cid:0)(cid:0)(cid:15)~basedmodelscanbesolvedwithe(cid:0)ithertheBoussi- p=(cid:26)RT~ (21)
nesq model or an explicit algebraic Reynolds stress
model.25{27 The Wilcoxk~ !~ based modelsare avail- (cid:13)R
(cid:0) 28 h~ = T~ (22)
able inthree forms;ahigh Reynoldsnumber model, (cid:13) 1
28 (cid:18) (cid:0) (cid:19)
a lowReynolds numbermodel, and an explicit alge-
braic Reynolds stress model.27 The k~ !~ models of H~ =h~+ 1 u~2+v~2+w~2 +k~ (23)
(cid:0) 2
4 of15 (cid:0) (cid:1)
AmericanInstituteof Aeronautics and Astronautics Paper99{3360
p 35
E~ =H~ (24) VULCAN. Wilke’sformula is then applied tomodel
(cid:0) (cid:26) the mixture molecularviscosity as
whereRisthegasconstantand(cid:13) istheratioofspeci(cid:12)c
Ncs
heats. (cid:31)i(cid:22)i
(cid:22)= Ncs (32)
Thepressure, speciesenthalpy,staticenthalpy,total i=1 j=1(cid:31)j(cid:30)ij
enthalpy and total energy for an arbitrary mixture of X
thermallyperfect gases are modeled as where P
1 1 2
(cid:22)i 2 Mj 4
Ncs f~i 1+ (cid:22)j Mi
p=(cid:26)RuT~Xi=1 Mi (25) (cid:30)ij = (cid:20) p(cid:16)8 1(cid:17)+(cid:16)MMji (cid:17)12 (cid:21) (33)
T~ (cid:16) (cid:17)
h~i =h~oi + CpidT (26) and (cid:31)i is the ith chemicalspecies molefraction.
To The molecular thermal conductivity, (cid:21), of a single-
Z
component gas is modeled using the Prandtl number
Ncs
such that
h~ = h~if~i (27) Cp(cid:22)
i=1 (cid:21)= (34)
X Pr
where u is the universal gas constant, i and Ri The molecular thermal conductivity of a multi-
R M
(= u= i)arethemolecularweightandgasconstant component gas is modeled using either using
R Mth o
of the i species, respectively, and hi is the refer- Eqs. (32,33,34) or by using Sutherland’s formula for
th
ence enthalpy at a reference temperature, To, of i the chemical species molecular thermal conductivity
chemical species. The speci(cid:12)c heat at constant pres- (cid:21)i such that
sure,enthalpy,andGibbsfreeenergyforeachchemical
3
species are modeled with a polynomialof static tem- 2
T~ o;i+ i
perature. Forexample,the Cpi polynomialismodeled (cid:21)i =(cid:21)o;i T S (35)
as either33 To;i! T~+Si
Cpi =Ai+BiT~+CiT~2+DiT~3+EiT~4 (28) where (cid:21)o;i, To;i, and Si are constants that are speci(cid:12)c
Ri to the chemicalspecies of interest and are availableas
part of the thermodynamics data base supplied with
or as34 VULCAN. Wassiljewa’s formula36 is then applied to
Cpi = Ai+Bi+Ci+DiT~+EiT~2+FiT~3+GiT~4 (29) modelthe mixture molecularviscosity
Ri T~2 T~ Ncs
(cid:31)i(cid:21)i
The coe(cid:14)cients Ai;Bi;Ci;Di;Ei;Fi; and Gi are read (cid:21)= Ncs (36)
i=1 j=1(cid:31)j(cid:8)ij
from thermodynamic curve (cid:12)t data base (cid:12)les com- X
piled for each polynomial from the data in McBride where P
33,34 1 1 2
et al. (cid:21)i 2 Mj 4
1+ (cid:21)j Mi
Molecular Transport Models
(cid:8)ij = (cid:20) (cid:16) (cid:17) (cid:16) (cid:17)21 (cid:21) (37)
The molecular viscosity, (cid:22), of a single-component p8 1+ MMji
gas is modeled using Sutherland’s formulawhere
(cid:16) (cid:17)
3 Chemical Source Term Model
2
T~ To+S The chemical source term, !_i, represents the mean
(cid:22)=(cid:22)o (30) th
To! T~+S of the net rate ofproduction of the i species in each
chemicalreaction. Thissourcetermcanbemodeledei-
where (cid:22)o, To, and S are constants that are speci(cid:12)ed ther withorwithoutturbulence/chemistry interaction
for the gas of interest. The molecular viscosity, (cid:22)i, of e(cid:11)ects. The chemical source term, without turbu-
each component of a multicomponent gas is modeled lence/chemistry interaction e(cid:11)ects, is modeled in the
39
using Sutherland’s formula manner of, such that
(cid:22)i =(cid:22)o;i TTo~;i!32 TTo~;i++SSii (31) !_i =MiXNj=cr1(cid:16)(cid:23)j00i(cid:0)(cid:23)j0i(cid:17)"Kfj mNY=cs1Cm(cid:23)j0m (cid:0)KbjnNY=cs1Cn(cid:23)j00n#
(38)
where (cid:22)o;i, To;i, and Si are constants that are speci(cid:12)c Eq.(38)istheproductionrateoftheithspecies asde-
to the chemicalspecies of interest and are availableas termined by the lawof massaction where Cm a0nd Cn
part of the thermodynamics data base supplied with are the participating species concentrations, (cid:23)jm and
5 of15
AmericanInstituteof Aeronautics and Astronautics Paper99{3360
00
(cid:23)jn their stoichiometric coe(cid:14)cients, and fj and bj where I represents the cell interface, L and R repre-
the reaction rate coe(cid:14)cients of the forwarKd and baKck- sent the left and right states, and A^~ is due to Roe’s
ward reactions. In addition, the constraint imposed 21
Ncs approximateRiemannsolver. A (cid:12)low(cid:12)dissipation(cid:13)ux
by mass conservation exists such that i=1!_i = 0. (cid:12) (cid:12) 22
splitting scheme based on the me(cid:12)tho(cid:12)d of Edwards
The forward and backward reaction rate coe(cid:14)cients
P that addresses some of the entropy violation proper-
are calculated based on the Arrhenius law and the
ties of Roe’s scheme is also implemented. The low
equilibriumrate coe(cid:14)cient, which are de(cid:12)ned as
dissipation (cid:13)ux split (cid:13)ux construction is de(cid:12)ned as
Kfj =AjT~Bjexp(cid:18)(cid:0)REuaT~j(cid:19) (39) EI =a21 (cid:26)LC+EcL(cid:0)(cid:26)RC(cid:0)EcR +Ep DL+pL+DR(cid:0)(p4R3)
bj = Kfj (40) where Ec(cid:2) is the convection co(cid:3)ntribu(cid:2)tion to the (cid:13)ux(cid:3)
p
K eqj and E is the pressure contribution to the (cid:13)ux. The
K (cid:6) (cid:6)
The constants in the Arrhenius law, j; j; and aj, parameters C and D (see Ref. 22) are functions
A B E
fora hydrogen-airreaction system and acomplete de- of the left and right contravariant Mach number and
scriptionofthemethodused toobtaintheequilibrium a21 =(aL+aR)=2withathefrozenspeedofsound. In
rate coe(cid:14)cient can be found in Drummondand Hus- both (cid:13)ux construction methods, the MUSCL scheme
saini37 or Carpenter.40 is used to reconstruct values ofthe primitivevariables
The chemicalsource term modelwhich includes the at the left and right states. The primitive variables
e(cid:11)ects of turbulence/chemistry interaction uses as- used are the species mass fraction f~i, the density (cid:26),
sumedprobabilitydensityfunctions(PDF)tocompute the Cartesian velocity components u~, v~, and w~, the
the mean species production rate. In this approach, static pressure p, and the turbulence variables k~, ~l, g~,
the law of mass action is written as and ~.
F
The viscous (cid:13)uxes are discretized using a (cid:12)nite-
Ncr Ncs 0 Ncs 00
00 0 (cid:23)jm (cid:23)jn volume representation of a central di(cid:11)erence opera-
!_i = i (cid:23)ji (cid:23)ji fj m bj n 46
M (cid:0) 2K C (cid:0)K C 3 tor which can be con(cid:12)gured to evaluate the cell
j=1 m=1 n=1
X(cid:16) (cid:17) Y Y interface gradients using either an approximategradi-
4 (415)
ent construction where the cross-derivative terms are
The mean reaction rate coe(cid:14)cients, fj and bj, are
K K neglected or a full gradient construction in which all
computedusingatruncated Taylorseries expansionof 46
components of the derivative are included.
theratecoe(cid:14)cient equationswithrespect tothemean
and (cid:13)uctuating static tem00perature, where the static EllipticSolution Procedures
temperature (cid:13)uctuation T is obtained fromthe solu- TheNavier-Stokesequations,Eq.(1),arehyperbolic
00 00
tion of the energy variance, e e ,00transport equation. intimeandellipticinspaceand,therefore,canbeinte-
The higher order moments of T introduced by the gratedintimetoprovideasteadystatesolutionthatis
Taylor series expansion are cgomputed from either an elliptic in space. VULCAN has several time integra-
assumed Gaussian or (cid:12) PDF as described in Ga(cid:11)ney tion procedures from which to choose, depending on
41,43
et al. The mean forward and backward species the type of problem being solved. A temporally sec-
0 00
Ncs (cid:23)jm Ncs (cid:23)jn ond order accurate low storage, explicit, multi-stage
production terms m=1Cm and n=1Cn are com- Runge-Kuttascheme47 can be used tointegrate either
puted from a multivariate (cid:12) PDF where the sum of
Q Ncs 00 00 Q to a steady state solution or in a time accurate man-
the species variance i=1fi fi transport equation is ner. Alternately, due to the sti(cid:11) nature ofthe govern-
solvedandthevarianceisusedinthemannerdescribed ing equations and the expense involved in computing
42 P 43
in Girimaji and Ga(cid:11)ney egt al. the chemical source terms, implicit time integration
Spatial Discretization schemes can be selected to integrate the equations
to steady state. An Euler implicit time integration
The governing equations, Eq. (1), are discretized
scheme is written in delta formas:
withacellcentered, (cid:12)nitevolumerepresentation. The
convection terms are treated in a second order ac- I n
[ +(cid:14)(cid:24)EQ+(cid:14)(cid:17)FQ+(cid:14)(cid:16)GQ SQ](cid:1)Q= R (44)
curate manner with the primitive variable extrapo- J(cid:1)t (cid:0) (cid:0)
20
lation/interpolation scheme (MUSCL) of van Leer. @E @Ev @F @Fv
Discontinuitiessuchasshocksarecapturedbyemploy- where EQ = @Q @Q , FQ = @Q @Q, GQ =
@G @Gv (cid:0)@S n+1 (cid:0)n n
ing(cid:13)uxlimitersor byvariantson the MUSCL scheme @Q @Q , SQ = @Q, (cid:1)Q = Q Q , and R is
44,45 (cid:0) (cid:0)
which employdi(cid:11)erentiable limiters. The numeri- thesteadystateresidualattimen. Thedirectsolution
cal (cid:13)ux at the cell interface maybe constructed using of Eq.(44) results in a large banded N N system of
(cid:2)
a (cid:13)ux-di(cid:11)erence-split scheme where the (cid:13)ux construc- equations(where N isthenumberofequations)which
tion is de(cid:12)ned as must be inverted at each timestep.
Solving Eq. (44) for even a medium sized three-
EI = 1 EL+ER 1 A^~ QR QL (42) dimensional problem is impractical due to storage
2 (cid:0) 2 (cid:0)
(cid:0) (cid:1) (cid:12)(cid:12) (cid:12)(cid:12)(cid:0) (cid:1) 6 of15
(cid:12) (cid:12)
AmericanInstituteof Aeronautics and Astronautics Paper99{3360
requirements. However, Eq. (44) may be factored
to reduce the storage. Three approaches to the (cid:3) n
[M+(A(cid:24)+B(cid:24)+C(cid:24))](cid:1)Qi;j;k= R
approximate factorization of Eq. (44) are currently (cid:0)
(cid:3)(cid:3) (cid:3)
implemented in VULCAN. The approach with the [M+(A(cid:17) +B(cid:17) +C(cid:17))](cid:1)Qi;j;k=M(cid:1)Qi;j;k
least factorization error and highest storage require- (cid:3)(cid:3)
[M+A(cid:16) +B(cid:16) +C(cid:16)](cid:1)Qi;j;k=M(cid:1)Qi;j;k (49)
ments is a planar Gauss-Seidel relaxation which uses
n+1 n
aDiagonallyDominantAlternatingDirectionImplicit Qi;j;k =Qi;j;k+(cid:1)Qi;j;k
49
(DDADI) scheme to factor the cross stream plane. where
The forward sweep in i of a nonlinear planar Gauss- I
M= SQ (50)
Seidel relaxation of Eq. (44) maybe written J(cid:1)t (cid:0)
Both factorized schemes above result in a tridiagonal
I
[ +B(cid:24)+B(cid:17)+B(cid:16) SQ](cid:1)Qi;j;k+ systemofblockN N matriceswhichmustbeinverted
J(cid:1)t (cid:0) ineither the two((cid:2)(cid:17);(cid:16)) orthe three ((cid:24);(cid:17);(cid:16))directions.
A(cid:17)(cid:1)Qi;j(cid:0)1;k+C(cid:17)(cid:1)Qi;j+1;k+ Note that Eq. (50) contains only the source term Ja-
cobian.
A(cid:16)(cid:1)Qi;j;k(cid:0)1+C(cid:16)(cid:1)Qi;j;k+1= The schemes described above are storage and op-
n n+1 n+1 n n n eration intensive. A third scheme commonly referred
R (Qi(cid:0)2;j;k;Qi(cid:0)1;j;k;Qi;j;k;Qi+1;j;k;Qi+2;j;k)
(cid:0) to as a diagonalizedapproximatefactorization (DAF)
(45)
50
scheme signi(cid:12)cantly reduces the storage and num-
and the backward sweep as
ber of operations by diagonalizing the Jacobians and
I reducingtheN N blocktridiagonalstoasystemofN
[ +B(cid:24)+B(cid:17)+B(cid:16) SQ](cid:1)Qi;j;k+ (cid:2)
J(cid:1)t (cid:0) scalar tridiagonals. The DAF scheme may be written
as
A(cid:17)(cid:1)Qi;j(cid:0)1;k+C(cid:17)(cid:1)Qi;j+1;k+ (cid:3) n
[I J(cid:1)tSQ](cid:1)Qi;j;k = J(cid:1)tR
(cid:0) (cid:0)
A(cid:16)(cid:1)Qi;j;k(cid:0)1+C(cid:16)(cid:1)Qi;j;k+1= (cid:0)1 (cid:3)(cid:3) (cid:3)
T(cid:24)[I+J(cid:1)t(cid:14)(cid:24)((cid:21)(cid:24);c (cid:21)(cid:24);v)]T(cid:24) (cid:1)Qi;j;k =(cid:1)Qi;j;k
n n n n n+1 n+1 (cid:0)
(cid:0)R (Qi(cid:0)2;j;k;Qi(cid:0)1;j;k;Qi;j;k;Qi+1;j;k;Qi+2;j;k) T(cid:17)[I+J(cid:1)t(cid:14)(cid:17)((cid:21)(cid:17);c (cid:21)(cid:17);v)]T(cid:0)(cid:17)1(cid:1)Q(cid:3)i;(cid:3)j(cid:3);k =(cid:1)Q(cid:3)i;(cid:3)j;k
(46) (cid:0)
(cid:0)1 (cid:3)(cid:3)(cid:3)
where A(cid:11), B(cid:11), and C(cid:11) are block N N matrices T(cid:16)[I+J(cid:1)t(cid:14)(cid:16)((cid:21)(cid:16);c (cid:21)(cid:16);v)]T(cid:16) (cid:1)Qi;j;k =(cid:1)Qi;j;k
(cid:2) (cid:0)
associated with replacing the (cid:14)(cid:11) of Eq. (44) with a n+1 n
Qi;j;k =Qi;j;k+(cid:1)Qi;j;k (51)
(cid:12)rst order di(cid:11)erence in the (cid:11)=(cid:24) ,(cid:17), and (cid:16) directions.
(cid:0)1
Each sweep is factored in the (cid:17);(cid:16) plane as where T(cid:11) and T(cid:11) are matrices of the left and right
eigenvectors,(cid:21)(cid:11);caretheinviscideigenvalues,and(cid:21)(cid:11);v
(cid:3) n
[M+(cid:18)(A(cid:17)+C(cid:17))](cid:1)Qi;j;k = R are a diagonalizedformof the viscous eigenvalues.
(cid:0)
(cid:3)
[M+(cid:18)(A(cid:16) +C(cid:16))](cid:1)Qi;j;k=M(cid:1)Qi;j;k (47) Space Marched Equations
n+1 n The Navier Stokes equation, Eq. (1), can be made
Qi;j;k =Qi;j;k+(cid:1)Qi;j;k
parabolic in space while remaininghyperbolic in time
where by dropping all di(cid:11)usive terms in the streamwise ((cid:24))
1
(cid:13)uxvectorandbyemployingVigneron’s treatmentof
I
M= +B(cid:24) +(cid:18)(B(cid:17)+B(cid:16)) SQ (48) the streamwise pressure gradient terms. This results
J(cid:1)t (cid:0)
in an equation of the form
and (cid:18) is a parameter that determines the order of ac- Qt+E(cid:24)+(!v 1)E(cid:20)!(cid:24)v +(F Fv)(cid:17)+(G Gv)(cid:16) =S
curacy of the time derivative and also can be used (cid:0) (cid:0) (cid:0)
(52)
to control the stability of the scheme. A value of 2
where
is suggested in Ref. 51 and Ref. 52 and is used here 0
as well. Note that the diagonal term M contains the .
.
diagonal contributions of (cid:14)(cid:24)EQ, (cid:14)(cid:17)FQ, (cid:14)(cid:16)FQ as well 2 . 3
0
as SQ and is shared by both factors. This relaxation
schemecouldalsobewrittenasaconventionaltwofac- 66 (cid:24)^xp 77
6 7
torscheme18 where the B(cid:17) andB(cid:16) termsare removed 6 (cid:24)^yp 7
from M and combined with their A and C contribu- E(cid:20)!v = jr(cid:24)j 66 (cid:24)^zp 77 (53)
6 7
tions. However,thiswouldincrease factorizationerror J 6 0 7
6 7
without reducing cost or storage appreciably. 6 0 7
6 7
A three factor approximate factorization method 6 0 7
6 7
that has lower storage requirements maybe written 6 0 7
6 7
6 0 7
6 7
6 0 7
6 7
4 5
7 of15
AmericanInstituteof Aeronautics and Astronautics Paper99{3360
The (!v 1)E(cid:20)(cid:24)!v term that appears in Eq. (52) is and
a correcti(cid:0)on to the streamwise (cid:13)ux di(cid:11)erence, E^(cid:24), M= I +B!(cid:24)c SQ (58)
J(cid:1)t (cid:0)
that parabolizes the streamwise pressure gradient.
This correction term contains Vigneron’s streamwise However, anapproachthat hasless factorizationerror
pressure-gradient coe(cid:14)cient, !v, which is constrained is to factor Eq. (56) in the cross stream plane using
such that the hyperbolic-parabolic nature of the gov- DDADI
erning equations is preserved with respect to the (cid:3) n
streamwise direction. The constraint equation that [M+(cid:18)(A(cid:17)+C(cid:17))](cid:1)Qi;j;k = R
(cid:0)
Vigneron obtained from an eigenvalue analysis is (cid:3)
[M+(cid:18)(A(cid:16) +C(cid:16))](cid:1)Qi;j;k=M(cid:1)Qi;j;k (59)
2
(cid:27)(cid:13)M(cid:24) n+1 n
!v =min 1; 2 (54) Qi;j;k =Qi;j;k+(cid:1)Qi;j;k
" 1+((cid:13) 1)M(cid:24) #
(cid:0)
and
2
M(cid:24)2 = (cid:24)x2+(cid:24)y2U+(cid:24)z2 (cid:13)p=(cid:26) (55) M= I +B!(cid:24)c +(cid:18)(B(cid:17) +B(cid:16)) SQ (60)
J(cid:1)t (cid:0)
where (cid:13) is the frozen(cid:0)ratio of spec(cid:1)i(cid:12)c heats and (cid:27) is a
The DDADI factorization is an expensive algorithm,
safety factor (usually = 0.95) included to account for
especially when the number of equations is large.
nonlinearitiesthatarenotconsideredintheeigenvalue
1 Therefore, a modi(cid:12)ed planar Gauss-Seidel relaxation
analysis.
whichuses the diagonalizedapproximatefactorization
Eq. (52)is written in acorrection formina manner
19 scheme is proposed that has lower storage and oper-
similarto Morrison and Korte. Morrison and Korte
ation requirements for use on moderate aspect ratio
found that the implementation of Vigneron’s stream-
grids. This scheme will converge poorly when the
wisepressuregradientapproximationinacellcentered
streamwise grid aspect ratio becomes large; however,
(cid:12)nite volume context required special treatment to
by combiningthe DAF factorizationwith wallmatch-
avoid the introduction of arti(cid:12)cial pressure gradient
ing functions, the wall normal grid spacing (and thus
termswhichinturncauseerrors thatarestronglygrid
19 aspect ratio) can be relaxed and acceptable conver-
dependent. As in Morrison and Korte, !v is evalu-
gence behavior recovered. A planar relaxation based
atedatthecellcenterbyaveragingthecellfacemetrics
on DAF maybe written as
and cell face primitive variables obtained from either
a (cid:12)rst or second order fullyupwind reconstruction. !c (cid:3) n
I+J(cid:1)t B(cid:24) SQ (cid:1)Qi;j;k = J(cid:1)tR
(cid:0) (cid:0)
Space Marching Solution Procedures
h (cid:16) (cid:17)i
Following Newsome et al.,18 Eq. (52) is solved by T(cid:17)[I+J(cid:1)t(cid:14)(cid:17)((cid:21)(cid:17);c (cid:21)(cid:17);v)]T(cid:0)(cid:17)1(cid:1)Q(cid:3)i;(cid:3)j;k =(cid:1)Q(cid:3)i;j;k
constructingthe(cid:13)uxesandsourcetermsforaconstant (cid:0)
(cid:0)1 (cid:3)(cid:3)
(cid:24) plane (three-dimensional) or line (two-dimensional T(cid:16)[I+J(cid:1)t(cid:14)(cid:16)((cid:21)(cid:16);c (cid:21)(cid:16);v)]T(cid:16) (cid:1)Qi;j;k =(cid:1)Qi;j;k
(cid:0)
and axisymmetric)of cells and then applying the for- (61)
n+1 n
ward sweep step ofthe planar Gauss-Seidel relaxation Qi;j;k =Qi;j;k+(cid:1)Qi;j;k
scheme presented in Eq. (45). Eq. (45) is modi(cid:12)ed
Eitheroftherelaxationprocedures givenaboveisthen
slightlyfor the space marching procedure to become
used to solve each (cid:24) plane of cells beginning at the
I !c in(cid:13)ow plane and sweeping in positive (cid:24) through the
[J(cid:1)t +B(cid:24) +B(cid:17)+B(cid:16) (cid:0)SQ](cid:1)Qi;j;k+ computational domain. Each (cid:24) plane of cells is iter-
ated in pseudo-time until convergence is reached (a
A(cid:17)(cid:1)Qi;j(cid:0)1;k+C(cid:17)(cid:1)Qi;j+1;k+
four order reduction in the sumof the L2 normof the
A(cid:16)(cid:1)Qi;j;k(cid:0)1+C(cid:16)(cid:1)Qi;j;k+1= continuity, momentum,and energy residuals is gener-
n n+1 n+1 n allysu(cid:14)cient). The solutionisthen projected intothe
R (Qi(cid:0)2;j;k;Qi(cid:0)1;j;k;Qi;j;k) (56)
(cid:0) next plane of cells as an initialcondition. This proce-
!c
where B(cid:24) is the Jacobian of the parabolized stream- dureisrepeatedforeachplaneofcellsuntiltheout(cid:13)ow
wise (cid:13)ux, E(cid:24) +(!v 1)E(cid:20)!(cid:24)v. plane is reached.
(cid:0)
Extending the approach of Newsome et al. to in-
clude source terms, the cross-stream plane is factored Convergence Acceleration
using a two factor scheme VULCAN has three techniques for accelerating the
(cid:3) n convergence of Eqs. (1) and (52). The (cid:12)rst technique
[M+(A(cid:17)+B(cid:17)+C(cid:17))](cid:1)Qi;j;k= R
(cid:0) is local time stepping, the second technique is coarse-
(cid:3) to-(cid:12)ne grid sequencing, and the third technique is the
[M+(A(cid:16) +B(cid:16) +C(cid:16))](cid:1)Qi;j;k=M(cid:1)Qi;j;k (57)
FullApproximateStorage (FAS) multi-gridscheme of
n+1 n 53
Qi;j;k =Qi;j;k+(cid:1)Qi;j;k Brandt.
8 of15
AmericanInstituteof Aeronautics and Astronautics Paper99{3360
The FAS grid transfer operators for the solution, constructed using second order, fully upwind extrap-
residual, and coarse grid corrections are those intro- olation and the cross-stream (cid:13)uxes were constructed
54
duced by Jameson. The multi-grid strategy can be using MUSCL with (cid:20) = 1=3 and the smooth lim-
48 45
either a V(1,0) or W(1,0) type cycle. The correc- iterofVenkatakrishnan. Agridsub-steppingfeature
tions after each prolongation are smoothed using a available in VULCAN was employed to control oscil-
scalarcoe(cid:14)cientimplicitresidualsmoothingtoremove lations near the shock. Space marching results in a
high frequency errors introduced by the prolongation lack of downstream information during the stream-
operator. To provide a su(cid:14)ciently smoothinitialcon- wise (cid:13)ux construction makingit impossibleto employ
dition for the (cid:12)ne grid, the Full Multi-grid (FMG) (cid:13)ux limiters in the \classical" sense. However, VUL-
method can be invoked. The use of FMG is similar CAN can be made to \sub-step" on the grid. Sub-
to coarse-to-(cid:12)ne grid sequencing with the exception stepping is a procedure that linearly sub-divides the
that multi-grid cycles are performed on each level of input streamwise cells, thereby reducing the gradients
the coarse-to-(cid:12)ne grid sequence. in the streamwise direction and increasing accuracy
The multi-gridschemealsohasseveralfeatures that while also reducing oscillations. Figure 3 shows the
arenotcommonlyemployed. Thegoverningequations typicalconvergence behaviorof several computational
are allsolved inafullycoupled fashiononallgridlev- planes from the solution. Each plane was initialized
els. Additionally,it was found useful for higher Mach withthesolutionfromthepreviousplaneandtheCFL
number(cid:13)ows toconstruct the viscous (cid:13)uxes aswellas wasrampedfromaninitialvalueof100toa(cid:12)nalvalue
the source terms on all coarse grids. The turbulence of 1000 over 5 iterations. All planes converged 4 or-
modelequations are solved onallgridlevels; however, ders intheL2 normin15iterationsorless. Multi-grid
some special treatment for wall boundary conditions was not found to appreciably improvethe already ex-
are employed. The wall boundary condition on the cellent convergence rate ofthisproblem. The problem
dissipation equation is either restricted from the (cid:12)ne wasalsorunusingtheDAFbasedGauss-Seidelscheme
grid or computedon the coarse grid by evaluatingthe using a CFL ramped from 1 to 10 over 5 iterations.
wallmatchingfunction on the coarse grid. The DAF scheme required 20 iterations on average
Implementation of the FAS and FMG methods is to reduce the L2 four orders of magnitude. However,
identicalforthehyperbolic-ellipticandthehyperbolic- theDAFschemehadapproximatelythesamecost per
parabolic forms of the equations with one exception. plane due to its lower cost per iteration.
Whensolvingthe hyperbolic-ellipticformofthe equa- Figure 4 presents the Mach contours from a su-
tion, the grid is coarsened in all directions uniformly. personic, calorically perfect ((cid:13) = 1:4) calculation of
However, when solving the hyperbolic-parabolic form turbulent (cid:13)ow over a 0.5 meter long (cid:13)at plate. The
of the equations, the grid is coarsened in the cross unit free stream Reynolds number is 2:64 107 and
(cid:2)
(cid:13)ow plane only. For two-dimensionaland axisymmet- the wall temperature ratio Tw=T1 = 0:65. The base-
ric problems this results in grid coarsening in the (cid:17) line k~ !~ turbulence modelwas solved to the wallon
(cid:0)
direction only; while for three-dimensional problems, acomputationalgridof200cellsinthe(cid:24) directionand
thisresultsincoarseninginthe(cid:17)and(cid:16) directionsonly. 153cellsinthe(cid:17)direction. Thestreamwise(cid:13)uxeswere
constructed using second order, fully upwind extrap-
Results olation and the cross-stream (cid:13)uxes were constructed
Four example problems are presented to demon- usingMUSCL with(cid:20)=1=3andthe smoothlimiterof
45 +
strate the capabilities of the multi-grid,space march- Venkatakrishnan. The y of the (cid:12)rst cell center o(cid:11)
ing, relaxation algorithm. The (cid:12)rst test case is in- thewallwaslessthan0.5fortheentireplate. Figure5
viscid, calorically perfect, Mach 2 (cid:13)ow over a two- provides the computed out(cid:13)ow velocity and eddy vis-
dimensional bump in a channel; the second is calor- cositypro(cid:12)ledemonstratingthattheboundarylayeris
ically perfect, turbulent, viscous, Mach 6 (cid:13)ow over well resolved. Figure 6 presents several typical \post-
a (cid:13)at plate; the third is calorically perfect, laminar, transition" convergence histories at a location where
Mach3(cid:13)owoverathree dimensionalcompressioncor- the grid aspect ratio was 300. The line Gauss-Seidel
(cid:25)
ner; and the fourth is a two-dimensional supersonic relaxation space marching scheme was used with a
turbulent di(cid:11)usion (cid:13)ame. W(1,0)cycle FMGwith 3grid levels. The L2 normof
Figure 1 shows the two-dimensional computational theresidualisreduced 12orders ofmagnitudeoneach
grid for a Mach 2 simulation of inviscid (cid:13)ow over a streamwise step in 32 multi-gridcycles. Although not
bumpinachannel. The gridconsists of64cellsinthe shown, it is importantto note that the residual ofthe
streamwise ((cid:24)) direction and 32 cells in the cross (cid:13)ow k~ and !~ equations both converged at the same rate
((cid:17)) direction. The grid aspect ratio can be seen to be and to the same level as shown here. The problem
nearlyoneforallcellsandthegridisorthogonaltothe was run without multi-grid and was found to require
lower boundary. Figure 2 presents density contours more than 130 cycles to converge to the same level as
from a space marched solution obtained using line shown in (cid:12)gure 6.
Gauss-Seidel relaxation. The streamwise (cid:13)uxes were The third case is the supersonic laminar three-
9 of15
AmericanInstituteof Aeronautics and Astronautics Paper99{3360