Algorithmic Enhancements to the VULCAN Navier-Stokes Solver

15 Pages·2003·0.63 MB·English

by  Litton, D. K.

AIAA 2003-3979
Algorithmic Enhancements to the
VULCAN Navier-Stokes Solver
D. K. Litton, J. R. Edwards
North Carolina State University, Raleigh, NC
J. A. White
NASA Langley Research Center, Hampton, VA
16th AIAA Computational Fluid Dynamics
Conference
23–26 June 2003
Orlando, Florida
Forpermissiontocopyorrepublish,contactthecopyrightownernamedonthefirstpage
ForAIAA-heldcopyright,writetoAIAAPermissionsDepartment,
1801AlexanderBellDrive,Suite500,Reston,VA20191–4344
AIAA 2003-3979
Algorithmic Enhancements to the VULCAN
Navier-Stokes Solver

D.K. Litton*, J. R. Edwards*, and J.A. White +

*Department of Mechanical and Aerospace Engineering
Campus Box 7910
North Carolina State University
Raleigh, NC 27695

+Hypersonic Airbreathing Propulsion Branch
Mail Stop 168
NASA Langley Research Center
Hampton, VA 23681-2199

VULCAN (Viscous Upwind aLgorithm for Complex flow ANalysis) is a cell centered,
finite volume code used to solve high speed flows related to hypersonic vehicles. Two
algorithms are presented for expanding the range of applications of the current Navier-Stokes
solver implemented in VULCAN. The first addition is a highly implicit approach that uses
subiterations to enhance block to block connectivity between adjacent subdomains. The addition
of this scheme allows more efficient solution of viscous flows on highly-stretched meshes. The
second algorithm addresses the shortcomings associated with density-based schemes by the
addition of a time-derivative preconditioning strategy. High speed, compressible flows are
typically solved with density based schemes, which show a high level of degradation in accuracy
( )
and convergence at low Mach numbers M £ 0.1 . With the addition of preconditioning and
associated modifications to the numerical discretization scheme, the eigenvalues will scale with
the local velocity, and the above problems will be eliminated. With these additions, VULCAN
now has improved convergence behavior for multi-block, highly-stretched meshes and also can
solve the Navier-Stokes equations for very low Mach numbers.

Introduction wall grid spacing. Convergence rates can degrade rapidly
for highly-stretched meshes. Furthermore, convergence
rates can degrade when large numbers of blocks are used,
The VULCAN (Viscous Upwind aLgorithm due to the lack of strong coupling between adjacent
for Complex flow ANalysis) Navier-Stokes solver is blocks. VULCAN is also designed for higher Mach
considered the NASA standard in simulating reacting number applications, and like most compressible flow
internal flows characteristic of high-speed propulsion solvers, can experience convergence degradation and
devices. [1] VULCAN can solve the Navier-Stokes solution inaccuracy for very low Mach number flows.
or Parabolized Navier-Stokes equations using a The present effort is concerned with improving the
variety of discretizations, integration algorithms, numerical efficiency of VULCAN for viscous flows on
turbulence models, and chemistry models and is multi-block, highly-stretched meshes and for general low
applicable to general structured grids on multi-block Mach number calculations. In the former context, planar
domains. relaxation based implicit methods are introduced, along
with sub-iterative procedures that allow for a large degree
MPI message-passing is used to adapt of implicit coupling among blocks. In the latter context, a
VULCAN to parallel architectures. The baseline time-derivative preconditioning strategy based on the “all-
integration strategy for Navier-Stokes applications
speed” flux formulae of [2] is implemented and tested for
within VULCAN is a diagonalized approximate laminar and turbulent flows of calorically- and thermally-
factorization scheme. This scheme is rather efficient perfect gas at low Mach numbers. The cases used for
on a per-iteration basis but is sensitive to the near-
validation are the West and Korkegi double wedge

1 of 14

AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS
AIAA 2003-3979
configuration [3], a 3-D channel flow, inviscid flow Ø rYW ø
over a bump along a 2-D channel, a flat plate Œ 1 œ
Œ M œ
simulation, and finally, a simulation of subsonic-to-
supersonic flow transition in a two-dimensional Œ rY W œ (4)
Œ NCS œ
nozzle. Œ rW œ
(cid:209) z Œ ruW +z pœ
Governing Equations G= Œ x œ
J Œ rvW+z pœ
VULCAN solves the three-dimensional Navier- y
Œ rwW +z pœ
Stokes equations, expanded to include separate Œ (E + p)Wz œ
transport equations for individual species and two- Œ t œ
equation turbulence model components. Written in a Œ rkW œ
Œ œ
generalized coordinate system, the Navier-Stokes set º rwW ß
may be written as
The Jacobian, J , of the transformation is defined as
dQ
+R(Q)=0, (1) ( )
¶ t ¶ x,h,z
J = (5)
( )
with the residual operator R(Q) given as ¶ x,y,z
( ) ( ) ( )
R(Q)= ¶ E - Ev + ¶ F - Fv + ¶ G- Gv - S (2) The components of the cell face unit normal and the
¶ x ¶ h ¶ z contravariant velocities are
The solution vector, Q, is given by: x = xx h = hx z = zx (6)
x (cid:209) x x (cid:209) h x (cid:209) z
Ø rY ø
Œ 1 œ
Œ M œ (cid:209) x = x2 +x2+x2
x y z
Œ rY œ
Œ NCS œ U =x u+x v+x w V =h u +h v+h w
Œ r œ x y z x y z
Q= 1ŒŒ ru œœ W =z u+z v+z w (7)
J Œ rv œ x y z
Œ œ
rw Definitions for viscous flux vectors, E, F , and G and the
Œ œ v v v
Œ Et œ source-term vectorS can be found in [1].
Œ rk œ
Œ œ In the above equations, Yi is the mass fraction of the ith
º rw ß
chemical species and N is the total number of chemical
CS
species.
(3)
For a calorically perfect gas, the pressure, total
and the inviscid flux terms are defined as
enthalpy, static enthalpy, and total energy are given by
ί rY1U Ͽ ί rY1V Ͽ p= rRT
Œ M œ Œ M œ
1( )
Œ rY U œ Œ rY V œ H =h+ u2 +v2 +w2
Œ NCS œ Œ NCS œ 2
Œ rU œ Œ rV œ
E= (cid:209) x ŒŒ rUu+ pxxœœ F = (cid:209) h ŒŒ ruV +hxpœœ h= CpT (8)
J Œ rUv+ px œ J Œ rvV +h pœ
Œ rUw+pxyœ Œ rwV +hy pœ E = rH - P
Œ (E + p)Uqœ Œ (E + p)Vz œ t
Œ t œ Œ t œ For a thermally perfect gas, the pressure, species
Œ rkU œ Œ rkV œ enthalpy, and static enthalpy are formed by the
Œ œ Œ œ
º rwU ß º rwV ß expressions:
p= rRT(cid:229)NCS Yi
u m
i=1
i

2 of 14

AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS
AIAA 2003-3979
h = h0 +(cid:242) TC dT (9) The actual Jacobian matrix is denoted as A. Note that the
i i T0 pi factorization of M is defined only over the interior grid
h = (cid:229)NCS hY points within a particular subdomain. In contrast, the
i i Jacobian matrix A contains elements that may multiply
i-1
corrections that are obtained from the solution of the
where Ri (Ru/µi) is the ith species gas constant. Ru is linear system in adjacent subdomains. Thus we may split
the universal gas constant and µi is the species the matrix A into M + N + E, where N contains
molecular weight.
elements of A that would multiply corrections in
adjacent subdomains and E is the factorization error.
Given this, a general iterative scheme for improving the
Planar Relaxation Implicit Flow Solvers
solution of the linear problem at a particular subdomain
for Multi-Block Domains.
may be defined as:

Algorithm
M(D Qn+1,l+1- D Qn+1,l)=- Rn - AD Qn+1,l
One of the major problems encountered
when solving three-dimensional problems on large for D Qn+1,l+1 (10)
numbers of blocks is a reduction in the overall
convergence rate as the number of blocks increases. Here, the index ndenotes a particular iteration level for
Typical domain-decomposition strategies used for
the solution of the nonlinear problem (for unsteady flows,
finite-volume discretizations only allow one or two
this could be part of another subiteration), and the index l
mesh cells of overlap between adjacent domains.
denotes a particular iteration level for the iterative
Typical implicit solvers, when formulated for multi-
improvement of the solution of the linear problem. With
block arrangements, do not consider matrix elements
this basic strategy in place, one can define an algorithm
that would multiply corrections generated in adjacent
for improving block-to-block coupling:
domains. Therefore, subdomain coupling is only

achieved in an explicit manner, through the residual
evaluation at cells adjacent to block interfaces. Solve: D Qn+1,1 =- M- 1Rn
The RLX3D option in VULCAN is built For l=1, lmax:
around a planar relaxation scheme for solving the
sub-domain implicit problem, and is designed to 1: Pass appropriateD Qn+1,l elements to ghost cells of
solve the complete (not parabolized) Navier-Stokes
adjacent blocks (parallel MPI send / receive)
equations. The chosen sweep direction may be

block-specific, and the crossflow plane linear system
is approximately solved using an incomplete LU 2: Solve: D Qn+1,l+1 =D Qn+1,l +M- 1(- Rn - AD Qn+1,l)
decomposition procedure. This approach alleviates
much of the numerical stiffness associated with End loop
highly-stretched mesh cells, provided that the
crossflow plane is oriented so that to encompass the Update: Qn+1 =Qn +D Qn+1,lmax+1
coordinate direction(s) with the largest degree of

mesh stretching. Techniques such as Jacobian
This algorithm requires that an extra block diagonal
freezing are used to reduce the overall CPU load, and
implicit boundary conditions are incorporated to matrix, corresponding to the block diagonal ofA, which
further enhance stability. The RLX3D option has is normally over-written by the planar ILU factorization,
been tested for supersonic turbulent flow past a be stored in addition to M itself. The only change to the
cylinder, laminar flow between two intersecting VULCAN input deck necessary to implement this
wedges [3] and 1-D unsteady flow using a dual-time- algorithm is a flag indicating the number of subiterations
stepping approach. performed, l . If this is set to zero, then no
max

subiterations are performed and the planar relaxation
An improved implicit algorithm has been
scheme alone is used to advance the solution. Test cases
developed with better block-to-block coupling.
shown in the Results section provide indications of the
M represents the planar relaxation implicit operator
degree of improvement in convergence rates offered by
as applied over a subdomain, with its action upon a
this approach.
residual vector,R,denoted by the operation M- 1R.

3 of 14

AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS
AIAA 2003-3979
[ ]
Time-Derivative Preconditioning 1 ( ) ( )
1+M2 U – U2 1- M2 2 +4V2
Algorithm 2 ref ref ref
The time-derivative preconditioning strategy V2
M2 = ref (14)
currently implemented into VULCAN combines the ref a2
rank-one preconditioning matrix of Weiss and Smith
[4] with the “all-speed” version of the low diffusion As the incompressible limit is approached, the
flux-splitting scheme (LDFSS) of Edwards [5], eigenvalues become
developed according to a methodology presented in
U ,U ,U
Edwards and Liou [2]. The preconditioning method
[ ]
can currently be used with Runge-Kutta explicit time 1
integration and diagonalized approximate U – U2 +4V2 (15)
2 ref
factorization implicit time integration.
whereas the eigenvalues revert to their traditional values
The preconditioned Navier-Stokes equation U,U,U , and U – a as V 2 ?a 2.
set is given by: ref
dQ
P +R(Q)=0, (11) Numerical Discretization
¶ t
To ensure accuracy at all flow speeds, it is necessary that
where the preconditioning matrix, P, is defined as
the numerical discretization of the inviscid flux terms
P=I +qurvrT (12) reflect the preconditioned eigensystem. There are several
approaches for doing this, with the most rigorous being
with the use of characteristic analysis to derive preconditioned
r [ ] analogues of matrix-dissipation methods. In the
uT = Y K Y 1 u v w H k w VULCAN implementation, we instead implement a
1 NCS
r Ø ¶ P ¶ P ¶ P ¶ P ø preconditioned variant of the low diffusion flux-splitting
vT =Œ K œ (13) scheme (LDFSS) of Edwards [5], developed according to
Œº ¶ Q1 ¶ Q2 ¶ Q3 ¶ QNeqœß a methodology presented in Edwards and Liou [2]. The
interface flux E in LDFSS is split into convective and
1 1 I
q = - pressure contributions as follows:
V 2 a2 [ ] [ ]
ref E =a r C+EC - r C- EC +Ep D+p +D- p
I 1 L L R L L L R R
where N represents the total number of equations. 2
eq
(16)
The reference velocity, V , is responsible
r
ref The vector ECis the same as u in Eq. (13), while
for scaling the eigenvalues of the equation set at low
[ ]
speeds to be of the same order. V is defined as Ep = 0 K 0 0 x x x 0 0 0 (17)
ref x y z
V2 =min(cid:231)(cid:230) a2,max(cid:231)(cid:230) Vr2,KV2(cid:247)(cid:246) (cid:247)(cid:246) , (14) The “preconditioned” interface sound speed a is defined
ref Ł Ł ¥ ł ł 1
2
as
r
2
where ais the sound speed and V is the velocity Ø ( ) ø
U 2 1- M 2 2 +4V2
magnitude. V in the above equation acts as a cutoff a = μ ( ref ) ref ϧ 1/2 (18)
¥ 1 1+M 2
velocity to prevent singularities in the proximity of 2 ref 1/2
stagnation regions. In the VULCAN implementation,
where the subscript ½ represents evaluation of the
the constant K scaling the cutoff velocity is a user
quantity using flowfield information arithmetically-
input, and V is set to the inputted free-stream
¥ averaged to the cell interface. The quantities C+, C- ,
velocity.
D+, and D- are functions of left-and right-state Mach
¶ E
The eigenvalues of P- 1 are: numbers, specially defined in terms of the interface sound
¶ Q speed and the reference Mach number as follows [2]:
U , U ,U

4 of 14

AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS
AIAA 2003-3979
[ ]
M~ =1[(1+M2 )M +(1- M2 )M ] (19) P- JD tSQ D Qi*,j,k =- JD tRn (25)
L 2 ref,1/2 L ref,1/2 R
This approximation allows the use of the Sherman-
~ 1[ ]
M = (1+M2 )M +(1- M2 )M (20) Morrison theorem to compute the action of P- 1on the
R 2 ref,1/2 R ref,1/2 L
residual vector in O(n) operations. The action of all other
with source Jacobian elements (those corresponding to
chemistry and turbulence source terms) on the residual
M =UL/R (21) vector is computed in a separate step, involving the use of
L,R a Householder transformations to ensure good numerical
1
2 stability.
For gas-dynamic flows, the use of the modified Mach Other additions to VULCAN required for
number definitions in conjunction with the preconditioning include the use of characteristic inflow
“preconditioned” sound speed enables the numerical boundary conditions based on the preconditioned
dissipation mechanism of LDFSS to scale properly at equations and the use of local time steps based on the
all speeds. Exceptions to this are the definitions for eigenvalues of the preconditioned system.
C+and C- , which contain a pressure-dissipation

term proportional to p - p . As shown in [2], this
L R Results
term acts to provide pressure-velocity coupling at low
speeds, and to ensure that this effect scales properly,
the term must be multiplied by 1/M2 . Planar Relaxation Results
ref,1/2

Precise definitions of the components of LDFSS
The testing of the algorithm has been focused on
may be found in [5], and a more recent extension
the West-Korkegi intersecting-wedge geometry [3] and a
valid for general fluids may be found in [6].
channel-flow analogue formed by eliminating the wedges
Time-Stepping Scheme and the clustering to the leading edge. The clustering in
the Y and Z directions remains the same, as do the length,
The solution is advanced in time by a
width, and height of the (now) rectangular geometry. The
preconditioned, diagonalized approximate
free-stream conditions for the intersecting-wedge
factorization (DAF) scheme. The preconditioned
version of the DAF scheme is written as follows: simulations are: M¥ =3, Re/m =2.11e6, T¥ =105 K. The
free-stream conditions for the channel-flow simulations
[ ]
I - JD tS D Q* =- JD tP-1Rn are: M¥ =0.5, Re/m =3.52e5, T¥ =105 K. In both cases,
Q i,j,k
[ ( )] the grid size is 65x125x125 and laminar flow is assumed.
TTxpp[II ++JJDD ttddx(llppx,c -- llx,v)]((TTxpp))--11DD QQ*i*,**j*,k ==DD QQi**,*j,k (22) Upanrlaelslesl oothne rtwheis eN moretnht ioCnaerdo,l ianlal cSauspese rwcoemrep puetirnfgo rCmeendt eirn
h [ h( h,c h,v )] h i,j,k i,j,k IBM-SP2 using a 16-block load-balanced decomposition
Tzp I +JD tdz lpz,c - lz,v (Tzp)-1D Qi,j,k =D Qi*,*j*,k of the baseline grid.

Figure 1 presents baseline results for the RLX3D
Qn+1 =Qn +D Q (23) planar relaxation method (l =0). The positive effect
i,j,k i,j,k i,j,k max
of using implicit boundary conditions is clearly indicated,
In this, the modal matrices T p, (Tp)- 1, etc. are as is the fact that the planar relaxation procedure allows a
x x
constructed from diagonalization transformations of much higher CFL that the baseline DAF scheme. This
the forms translates in a significant CPU savings, as the DAF
scheme at a CFL of 3.5 takes 2.9 hours to run on the
Tplp (Tp)-1 = P-1 ¶ E , (24) NCSC IBM-SP2 (16 processors) while the planar
x x,c x ¶ Q relaxation scheme at a CFL of 150 requires only 1.95
hours.
where lpx,c is a diagonal matrix containing the Figure 2 presents results from a CFL-ramping
eigenvalues of the preconditioned Euler system. exercise performed for the supersonic West-Korkegi case.
Note that the first step of the DAF procedure is an The final CFL number is reached by ramping from 0.1 to
approximation of the more exact expression 20 over the first 500 iterations, from 20 to 150 over the
next 500, and from 150 to the final value over the next

5 of 14

AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS
AIAA 2003-3979
250 iterations. In this case and in most subsequent order planar relaxation scheme with l = 0 and l =
max max
ones, the first 500 iterations are performed on a
2. As shown, more than a three-fold improvement in the
coarse mesh using a first-order accurate inviscid flux
number of iterations required for convergence is
discretization. Jacobian freezing is initiated after 250
evidenced for the first–order discretization. CPU times
iterations, with re-evaluation and factorization of the
for the first order discretizations are 107.30 minutes for
matrices performed every 5 iterations past this point.
l = 0 and 33.41 minutes for l =2. These times are
The controlling parameter l is set to one for this max max
max
also nearly a three-fold improvement. It is possible that
study. Figure 2 shows that there is little advantage to
the system load may have been different for each of these
choosing a CFL higher than 150 for this case. Figure
runs, as the fact that the CPU speedup is nearly in accord
3 shows the effect of the choice of sweep direction on
with the iteration count is somewhat surprising, given the
the performance of the iteratively improved planar
extra expense of the subiteration procedure. It is
relaxation algorithm with l = 0 and l =1.
max max noteworthy that the use of subiterations stabilizes the
Sweeping in the “i” direction (the direction of the second-order case to the point that its convergence rate is
dominant movement of the supersonic flow) is very similar to the first-order case. Otherwise, as
clearly preferable to sweeping in the “j” direction. evidenced by the results, the calculation eventually
Performing one subiteration to improve the solution diverges. In comparison with the supersonic West and
of the linear system improves the performance in Korkegi case, these results indicate that the benefits of
both cases, at least in terms of the number of subiterative improvement of the linear system solution
iterations. may be much larger for subsonic flows. The technique
Figure 4 illustrates the effect of varying appears to aid in damping and/or expelling pressure
l on the number of iterations required for disturbances that otherwise tend to reflect from physical
max
and interface boundaries.
convergence. As shown, the number of iterations

required for convergence drops as the number of
Preconditioning Results
subiterations performed increases. Also shown for

comparison is a calculation performed on the single-
The four test cases for the preconditioned system
block grid using l = 0. This calculation was
max are flow over a flat-plate, flow through a two-dimensional
performed on a Compaq ES-40, which has enough UTRC nozzle, flow between intersecting wedges, and
shared memory to store all of the matrix elements in inviscid flow over a bump in a channel. These
core. Interestingly, the single-block grid correspond to variations on standard test cases included in
performance at l = 0 is slightly worse than the 16- the VULCAN package. In all cases the maximum CFL is
max
set to 2.5, and most cases involve both turbulent and
block performance. Otherwise, the trends are what
laminar flow as well as multi-component gases. In all
one might expect. As the work increases
turbulent cases, the Wilcox (1998) k-w model is used,
significantly with the increase in the number of
while for all multi-component cases, a mixture of nitrogen
subiterations, it is instructive to examine wall clock
and oxygen is used.
time. Figure 5 shows that for this predominately
supersonic flow, there is little benefit to performing Two-dimensional flow over flat plate
the subiterations, with only about a 15%
Figure 7 shows the 65x129 grid for the flat-plate
improvement in overall execution time for the best
simulation. The first run was to compare the results of
case of lmax = 2. ramping the Mach number down from Mach 0.5 to Mach
The next test case corresponds to Mach 0.5 0.005. The free-stream conditions for the Mach 0.5
flow through a channel similar in dimension to the simulation are: Re/m = 1.11e7, T¥ = 300K. In each
West-Korkegi geometry. The single-block successive case the only parameter changed was the Mach
65x125x125 grid is decomposed into 16 blocks along number, which decreased the Reynolds number by a
the “i” coordinate. The CFL is ramped from a factor of 10 for each succeeding run. Figure 8 presents
starting value of 0.1 to a final value of 20 for the convergence histories for preconditioned and non-
planar relaxation variants, and again, 500 iterations preconditioned cases at each Mach number. Strikingly, it
are performed to first-order spatial accuracy on the can be seen that the convergence of the non-
coarser mesh before interpolating the solution to the preconditioned system is not altered by lowering the
finest mesh. Figure 6 portrays convergence histories Mach number, whereas the preconditioned system
for four runs: the first-order planar relaxation converges faster for only Mach 0.5 and 0.05. A closer
scheme with l = 0 and l = 2, and the second- investigation into the solution produced by the non-
max max
preconditioned system verified that the solution was in

6 of 14

AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS
AIAA 2003-3979
fact incorrect. This was verified when comparing shows the grid used in the calculation, while Figures 12
with the Blasius solution at distances of 0.3, 0.40, and and 13 show contours of Mach number and eddy-
0.45 meters from the front edge of the flat plate. By viscosity ratio (referenced to the laminar value),
using the Blasius solution for flow over a flat plate, respectively. In the other calculations presented in this
the boundary layer thickness can be calculated as: paper, the constant K scaling the cutoff velocity in Eq.
(14) is set to one, since the velocity everywhere is near
( ) 5.0x
d x = (26) the free-stream velocity. In the nozzle calculation,
Re however, the free-stream sound speed is chosen as the
x
reference velocity. Thus, a choice of K=1 will not activate
Table 1 shows how the non-preconditioned system preconditioning in the elliptic region. Figure 14
behaves in comparison to the preconditioned system. illustrates the effect of lowering K (equal to “qctoff” in
From this table it can be inferred that the non- the figure) on the convergence rates. The best results (for
preconditioned system is in fact converging to an a four order-of-magnitude reduction) are obtained for
incorrect solution. The preconditioned system has a K~0.1. Lower values result in a significant degradation in
more realistic value for the thickness of the boundary convergence rate, though the slope appears to be more
layer. For the Mach 0.05 run, the error for the consistent. These results indicate that the proper
preconditioned system is no larger than 6.7%. For reference-velocity choice for strongly-mixed subsonic /
the Mach 0.005 simulation, the error is less than supersonic flowfields may not be obvious, and that trial-
2.3%. Both of these numbers are in stark contrast to and-error procedures may have to be used to obtain the
the 85-95% error found in the results obtained best results. Even with this ambiguity, the use of
without preconditioning. preconditioning results in a factor of 2 savings in iteration
count. This translates into nearly a factor of 2 savings in
Figures 9 and 10 show the results for a
CPU time, as the modifications necessary for
turbulent, calorically perfect and a turbulent, two-
implementing preconditioning require very little
species air flow over a flat plate, respectively. The
additional CPU time.
convergences for each case are very similar to the
laminar flow in Figure 9. In all three examples, the Three-dimensional flow through intersecting wedges
preconditioned Mach 0.5 case showed a marked
Figures 15-18 present results from simulations of
improvement in convergence over the non-
subsonic viscous flow through the West-Korkegi
preconditioned system. This is somewhat in contrast
intersecting wedge geometry (shown in Figure 15) Free-
with the results from the eigenvalue analysis, which
stream conditions are chosen to be the same as in the two-
indicate that the non-preconditioned system should
dimensional flat-plate case. Figure 16 illustrates the
have a better overall condition number at this Mach
convergence histories for laminar flow through the
number. One reason may be the presence, in the
wedges. The three non-preconditioned solutions do not
preconditioned flux-splitting, of pressure-diffusion
converge at the same rate and show worse convergences
terms that tend to smooth out variations in the
than the preconditioned scheme at all Mach numbers,
pressure field. As will be seen in the next few
except Mach 0.005. As the Mach number decreases in
examples, this result is consistently independent of
magnitude beyond the Mach 0.5 case, the convergence
geometry. The convergence degradation indicated
rate is shown to decrease quickly. The preconditioned
for the Mach 0.005 calculations could be associated
system contains a few oscillations but maintains its
with the decrease in Reynolds number. Stiffness due
downward trend towards convergence. Figure 17 shows
to low Reynolds numbers will not be alleviated by
results for the turbulent, calorically perfect case. This
the inviscid preconditioning techniques currently
case shows trends that are almost identical to the laminar
employed in VULCAN.
flow.
Two-dimensional flow through a nozzle
The disparity in the results given by the non-
The next case considered is flow through the preconditioned system is magnified when running the two
two-dimensional UTRC nozzle. This is a standard species, turbulent simulation. As can be seen in Figure
test case for the VULCAN solver that involves 18, the non-preconditioned residual at Mach 0.005 does
decomposing the nozzle flow into two regions: an not go down, but rather oscillates wildly around 10-1.
elliptic region upstream of the nozzle throat and a Although its preconditioned counterpart does not display
parabolic region downstream of the throat. The this behavior, the convergence rate is noticeably slowed
elliptic region is solved using the DAF scheme down. The residual of the preconditioned system
(preconditioned and non-preconditioned), while the continues to go down towards convergence with minimal
parabolic region is solved using space-marching once oscillations (in comparison to non-preconditioned
the elliptic solution has been obtained. Figure 11 system).

7 of 14

AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS
AIAA 2003-3979
To verify that the results produced by the calorically-perfect gases. The second algorithm is a time-
preconditioning formulation were physically derivative preconditioning strategy that is intended to
consistent, calculations of laminar flow through the expand the range of applicability of VULCAN toward
channel analogue of the West-Korkegi geometry low-speed, nearly incompressible flows. This addition is
were also performed. At locations far enough away valid for calorically and thermally-perfect gas and is
from the corner, it was anticipated that the boundary designed for use with the baseline diagonalized
layer would develop according to the Blasius scaling approximate factorization algorithm in VULCAN. Test
shown in Eq. (26). Table 2 compares predictions cases show that the preconditioning framework greatly
from the preconditioned and non-preconditioned improves solution accuracy for low Mach numbers.
formulations versus the Blasius result. The Stiffness due to low-Reynolds number effects, is not,
boundary layer thickness obtained from the non- however, eliminated in the present formulation, leading to
preconditioned formulation turned out to have an some convergence degradation for low-speed, low
error of no less than 85%, while the preconditioned Reynolds number flows. This may require the use of
formulation provided results within a reasonable 6% “viscous” preconditioning strategies or a more implicit
of the theoretical values. Convergence trends were treatment of the viscous terms. Future work will focus on
similar to those corresponding to the intersecting combining the time-derivative preconditioning techniques
wedges and are thus not shown. with the fully-implicit formulations to arrive at a
framework capable of alleviating most sources of
Inviscid flow over a bump in a channel
numerical stiffness present in large-scale flow
The last example of the validation of the calculations.
preconditioning strategy is inviscid flow over a bump Acknowledgments
in a channel. For this particular case, the Euler
This research was supported by NASA Langley
equations are solved, thus Reynolds-number effects
under grant NAG-1-02052. IBM SP-2 computer time was
illustrated in the earlier calculations are not present.
obtained from a grant from the North Carolina
Figure 19 portrays the grid used for the calculations,
Supercomputing Center.
while Figure 20 shows the convergence histories for
Mach numbers of 0.5, 0.05, and 0.005. These
calculations reveal the expected trend of (nearly) References
Mach-number independent convergence rates when
using the preconditioning technique. In contrast, the [1] White, J.A. and Morrison, J.H., “A Psuedo-Temporal
non-preconditioned formulation displays a significant Multi-Grid Relaxation Scheme for Solving the
degradation in convergence rate as the Mach number Parabolized Navier-Stokes Equations,” AIAA 99-3360,
is lowered. July, 1999.
[2] Edwards, J.R. and Liou, M.-S. “Low-Diffusion Flux-
Finally, an important property of the
Splitting Methods for Flows at all Speeds”, AIAA Journal,
preconditioned DAF scheme is verified in Figure 21.
Vol. 36, No. 9, 1998, pp. 1610-1617.
The preconditioned scheme is designed to revert back
[3] West, J.E., and Korkegi, R.H., “Supersonic Interaction
to the compressible Navier Stokes equations for
in the Corner of Intersecting Wedges at High Reynolds
higher Mach-number flows. It is apparent from this
Numbers,” AIAA Journal, Vol. 10, No. 5, 1972.
figure that the two schemes are in fact identical for
[4]Weiss, J.M. and Smith, W.A., “Preconditioning
supersonic flow over the bump, thus demonstrating
Applied to Variable and Constant Density Flows,” AIAA
the validity of the preconditioned diagonalized
Journal, Vol. 33, 1995, p. 2050.
approximate factorization scheme for all flow speeds.
[5] Edwards, J.R. “A Low-Diffusion Flux-Splitting

Scheme for Navier-Stokes Calculations,” Computers &

Fluids, Vol. 26, No. 6, 1997, pp. 635-637.
Conclusions
[6] Edwards, J.R., “Towards Unified CFD Simulations of
Real Fluid Flows,” (invited), AIAA 2001-2524CP, 2001.
Two algorithms for enhancing the

capabilities of the VULCAN Navier-Stokes solver

have been presented. The purpose of the first

algorithm is to improve convergence rates for viscous

flows on highly-stretched, multi-block meshes. The

algorithm decreases not only the number of iterations

to convergence, but also the time to convergence, an

important factor in weighing the importance of this

new addition. So far, this addition is specialized for

8 of 14

AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS
AIAA 2003-3979

0
PlanarrelaxationwithoutimplicitBC's(CFL=20) 0
- 0.5 PlanarrelaxationwithimplicitBC's(CFL=20)
PlanarrelaxationwithimplicitBC's(CFL=150)
-0.5
-1 DiagonalizedAF(CFL=3.5)

-1.5 500itsonthecoarsemesh;5000itsonthefinemesh -1
Jacobianfreezingactivatedafter250itsoneachmesh
-2 CFL=125
RelL2-2.5 ImplicitBC'sinJandKdirections -1.5 CCCFFFLLL===112565050
-3 elL2-2 CFL=280
R
-3.5
-2.5
-4
-3
-4.5
-5
1000 2000 3000 4000 5000 -3.5
CYCLE
Figure 1: Convergence of planar relaxation -4
500 1000 1500 2000
scheme CYCLE
Figure2:EffectofCFLnumberonconvergence
0
0
-0.5
-0.5
relaxationsweepin"i"direction;l =0 -1
-1 relaxationsweepin"i"direction;lmax=1
relaxationsweepin"j"direction;lmax=0
relaxationsweepin"j"direction;lmax=1
max -1.5
-1.5 l =0;singleblockgrid
lmax=0;16blockgrid
RelL2-2 RelL2-2 lllmmmmaaaaxxxx===123;;;111666bbbllloooccckkkgggrrriiiddd
l =4;16blockgrid
max
-2.5
-2.5
-3 -3
-3.5 -3.5
-4 -4
1000 2000 3000 4000 5000 0 1000 2000 3000
CYCLE CYCLE
Figure3:Effectofsweepdirectiononconvergence Figure4: Effectofsubiterationnumberonconvergence

9 of 14

AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS