Multigrid Relaxation of a Factorizable, Conservative Discretization of the Compressible Flow Equations

14 Pages·2000·0.73 MB·English

by  Roberts, Thomas W.

AIAA OO-2252
Multigrid Relaxation of a Factorizable,
Conservative
Discretization of the Compressible
Flow Equations
Thomas W. Roberts
NASA Langley Research Center
David Sidilkover
Institute for Computer Applications in Science and Engineering
and J. L. Thomas
NASA Langley Research Center
Hampton, VA
Fluids 2000
%
19-22 June 2000/Denver, CO
For permission to copy or republish, contact the American Institute of Aeronautics and Astronautics
370 L'Enfant Promenade, S.W., Washington, D.C. 20024
/
i
l
MULTIGRID RELAXATION
OF A FACTORIZABLE, CONSERVATIVE
DISCRETIZATION OF THE COMPRESSIBLE FLOW EQUATIONS
Thomas W. Roberts*
NASA Langleg Research Center
D. Sidilkover t
Institute for Computer Applications in Science and Engineering
and J. L. Thomas t
NASA Langley Research Center
Hampton, VA
Abstract error corresponding to advection factor. This same diffi-
The second-order factorizable discretization of culty is true for high Reynolds-number viscous flows as
well.
the compressible Euler equations developed by
Sidilkover is extended to conservation form on gen- As has been pointed out by Brandt 1, to obtain ideal
era(cid:0) curvilinear body-fitted grids. The discrete multigrid convergence rates for subsonic, inviscid flows,
equations are solved by symmetric collective Gauss- the discretization nmst distinguish those parts of the dif-
Seidel relaxation and FAS multigrid. Solutions for ferential operator which correspond to advection, and
flow in a channel with Mach numbers ranging from those which correspond to elliptic behavior. The advec-
0.0001 to a supercritical Mach number are shown, tion terms are treated efficiently by marching and the el-
demonstrating uniform convergence rates and no liptic terms are rapidly solved with a multigrid iteration.
loss of accuracy in the incompressible limit. A so- Brandt 1argams that by spfitting the system into its ellip-
lution for the flow around the leading edge of a tic and advection parts, the convergence rate of the full
semi-infinite parabolic body demonstrates that the system ought to be equal to the slowest of the two sub-
scheme maintains rapid convergence for a flow con- systems. Using this approach, Brandt and Yavneh have
tiiining a stagnation point. demonstrated textbook multigrid for the incompressible
Navier-Stokes equations 2. Their results are for a simple
Introduction
geometry and a Cartesian grid, using a staggered-grid
Steady inviscid flow is described by the Euler equa- discretization of the equations. In a closely related ap-
tions, which can be thought, of as two snbsystems. One proach, Ta'asan 3 presented a fast multigrid solver for the
subsystem corresponds to the equations governing en- compressible Euler equations. This method is based on
tropy and vorticity advection. This subsystem is hyper- a set of "canonical variables" which express the stead),
bolic in space. The other subsystem corresponds to a full Euler equations in terms of an elliptic and a hyperbolic
potential operator, which is elliptic for subsonic flow and partition 4. In Reference 3 it. is shown that ideal multi-
hyperbolic for supersonic flow. For a purely supersonic grid efficiency can be achieved for the compressible Euler
flow, space marching is the most efficient way of solving equations for two-dimensional subsonic flow using body-
the Euler equations. For subsonic flow, the ellipticity of fitted grids.
the full potential factor should be effectively handled by
multigrid. However multigrid is not effective for advec- In an earlier work 5' 6, the authors presented a multi-
tion operators, as the coarse grid only gives part of the grid scheme for the steady, incompressible Euler equa-
correction for certain smooth components of the error. tions based on a pressure Poisson discretization which
Existing multigrid methods for subsonic and transonic distinguishes the advection operator from the elliptic
flow rely on the coarse grid to smooth the entire system. part of the system of equations. Both structured grid
As such, they are fundamentally limited by the ineffec- and unstructured grid flow solvers were written using
tiveness of the coarse grid in correcting the part of the the discretization. The results in Ref. 5 demonstrated
that the scheme can achieve ideal multigrid convergence
Copyright Q2000 by the American Institute of Aeronau- rates for internal flows. The scheme has been extended
tics and Astronautics, Inc. No copyright is asserted in the to three dimensions by Sanchez 7, who has also demon-
United States under Title 17, U.S. Code. The U.S. Govern- strated ideal multigrid convergence rates for internal
ment has a royalty-free license to exercise all rights under the flows. This approach has the advantage of non-staggered
copyright claimed herein for government purposes. All other
grids--a collocated, vertex-based discretization of the
rights are reserved by the copyright owner.
*Research Scientist equations is used, simplifying the restriction and prolon-
t Senior Staff Scientist gation operations and allowing the use of simple point
t AIAA Fellow collective Gauss-Seidel relaxation. Although the pres-
surePoissoanpproacmhaybeextendetodcompressible
series about each grid vertex and considering the lead-
flowsi,t.isnotconservatiavnedisunsuitedforsuper-
ing terms of the truncation error. These terms are the
criticalflowswithshocksF.urthermorteh,eextension
artificial dissipation of the scheme.
toviscoufslowsislimitedtotileincompressiNbalevier-
The starting point fox" the scheme is the two-
Stokeesquations.
dimensional Euler equations in non-conservation form.
RecentlyS,idilkov8ehrasdeviseaddiscretizatioonf Let p be the density, ff = _u + jr, be the velocity, and p
thecompressibflloewequationtshat.overcomethsese be the pressure. The entropy s is defined as
limitationsT. hisdiscretizatiomnaybeappliedtothe
Eulerequationinsconservatfivoerm,usingthemultidi-
mensionu;dpwindschemoefSidilkov9e'r10.InRef.10
itisshownthatthisapproacshhouldleadtoascheme
thatis factorizablie.e,.,theschemdeistinguishebse- where p0 and po are a reference pressure and density,
tweenthosepartsoftheoperatotrhatrepresenatd,vec- respectively, and "_ is the ratio of the specific heats.
tion,andthatpart Then the Euler equations may be written in the vari-
of the operator that. represents po-
ables (s, u, v,p):
tential flow. In Ref. 8, such a factorizable scheme is
constructed for Cartesi,'m grids. Because the diseretiza-
a.Vs = 0, (2a)
tion is stable for Gauss-Seidel relaxation, it the conver-
gence rate does not depend upon a Courasat-Friedrich- ff-Vz7 + 1Vp = 0, (2b)
Lewy number restriction, unlike standard discretizations P
of the Euler equations which must use time marching. pc2V.a + ff.Vp = 0. (2c)
For this reason, the same convergence rate is obtained
for subsonic Mach numbers all the way to the incom- The factorizability of the scheme depends on the form
pressible limit. The scheme may be written in the form of the artificial dissipation added to this system of equa-
of a central-difference part plus an artificial viscosity. tions.
As such, it is very similar in formulation to standard The entropy is only weakly coupled to the momen-
upwind, finite-vohune discretizations for the Euler and tum and the pressure equations through the equation
Reynolds-Averaged Navier-Stokes equations, and can be of state (1). In fact, the entropy equation corresponds
written as a conventional upwind scheme with a modi- to one of the advection factors of the Euler equations.
fied artificial dissipation. Sidilkover 8 shows that for the Therefore, it may be discretized independently of the
discrete scheme to preserve factorizability the dissipation momentum and pressure equations in any appropriate
terms must be discretized in a specific way. In addition, way without affecting the factorizability of the scheme.
he shows that the artificial viscosity may be rescaled with The advection operator ff.V uses simple upwind differ-
the Mach number such that the factorizability, and thus encing in Eq. (2a). Let (,_,r/) be a general curvilinear
the accuracy, as well as the h-ellipticity of the operator coordinate system and define the contravariant compo-
is preserved ix, the incompressible limit. nents of the velocity, (U, F') by the transformation
In the present work, the generalization of Sidilkover's
(;:;,9
factorlzable scheme to curviIinear, body-fitted coordi-
nates is presented. First., the governing equations are
presented, including the form of the multidimensional
upwind artificial dissipation. It. is shown how the dissi- In this coordinate system g.Vt7 = (?O_ + f'O,. The equa-
pation terms must be discretized to maintain factorlz- tions are discretized on a uniform grid in (¢,t/) space,
ability. A discussion of the point collective Gauss-Seidel with a grid spacing A¢ = ATI = 1. The FDA of the
relaxation is presented next. Solutions for flow in a chan- first-order upwind difference approximation to tT.V is
nel_ with infetl _Iach numbers ranging from 0.0001 to
1 1 .=
a supercritical Mach number are shown, demonstrating (4)
the accuracy and convergence rates of the scheme. An
additional computation for the flow around the leading
and the entropy equation is discretized as
edge of a semi-infinite parabola demonstrates that the
current scheme does not suffer from the convergence dif-
qs = 0. (5)
ficulties near a stagnation point that were observed for
the authors' previous scheme 11
A second-order upwind discretization of the adveetion
operator has also been used in Eq. (5). However, for
Mathematical Formulation
the particular cases shown below the use of a second-
The artificial dissipation of the factorizable schemc order advection operator has an insignificant attect on
can be described hy first presenting the modified equa- the computed results. This point, is discussed in the Re-
tion, or first differential approximation (FDA), of the sults section.
discrete scheme. This is the differential equation which The dissipation for the momentum and pressure equa-
is found by expmlding the difference equation in aTaylor tions, Eqs. (2b) and (2c), is the umltidimensional upwind
dissipation
of Sidilkover 9. In vector notation, this dissi- adding the following corrections:
pation is written as
a-_"/O = qL/-4- 7l [(--?[05/at -'t- l [_='l0_0._', (9a)
a.va+ 1-vp- "'ev (p 2v.a+a.vp) =0, (6) 1 1
2 p 2pc
pc2_7.ff+ff.Vp --pc-_.(_._ff+ _Vp) =0, (7)
The cross-derivative terms in Eq. (9) are necessary for
second-order accuracy and to retain factorizability at the
where c is the speed of sound, am and ap are scaling same time. Using Eq. (8) to find /at and _', the antidis-
coefficients, and ( is a length scale proportional to the sipative terms in Eq. (9a) and (9b) are evaluated, and
grid spacing. Note that the dissipation of the momen- then rotated back into the physical components. This
tum equation is tile gradient of the pressure equation
gives a discretization of the form
residual, and the dissipation of the pressure equation is
the divergence of the momentum equation residual. It /7.Vff = qa + VD (10)
is this property of the artificial dissipation that leads to where
a factorizable scheme. Note that ttle curl of Eq. (6) is
(iCrlao,+
identical to the curl of Eq. (2b), i.e., vorticity equation of D= GI . (11)
the governing Euler equations is unaffected by the arti-
ficial dissipation. The vortieity equation corresponds to The above expression is the generalization to curvilinear
the second adveetion factor of the Euler equations. Sim- coordinates of the one given in Ref. 8.
ilarly, because the dissipation of Eq. (7) is the divergence To gain more insight about, the nature of the correc-
of the momentum equation, it is affected only by tile Jr- tion terms, the FDA of Eq. (10) may be rewritten as
rotational part of the velocity field but not the solenoidal J
+vD= +.7(c' I '1 e"lC'l
part. The irrotational part may be written as the gradi-
ent of a potential. As the pressure equation Eq. (2e) is
where _o_ 0_t, - Ouu is the vorticity, J = av_y,7- y{x,
a form of the continuity equation, it corresponds to the
is the Jacobian of the coordinate transformation, mM
full potentiM factor, and the FDA described by Eq. (7)
(4e, gn) are the contravariant basis vectors in the ((, _/)
preserves this property.
coordinate directions.* Note that. if the flow is irrota-
If the advection operator, pressure gradient, and di-
t.ional, then the first-order truncation error terms vanish
vergence terms in Eqs. (6) and (7) are discretized us-
identically and the approximation becomes second-order
ing central differences, the scheme is second-order accu-
accurate.
rate and factorizable. However, such a scheme is not
It was pointed out above thai factorizalfility depends
h-elliptic. This lack of h-ellipticity is a result of the
upon the dissipation of the momentum equation being
central difference approximation to the advection term
written as the gradient of the pressure equation, and the
in the momentum equation, Eq. (6). Sidilkover shows
dissipation of the pressure equation being written as the
that this advection operator corresponds to vorticity
divergence of the momentum equation. Likewise, the
advection 8. Replacing this operator by the first-order second-order correction to the adveetion terms must be
upwind approximation in Eq. (4) restores h-ellipticity
in the form of a gradient so that the vorticity equation
and the scheme remains factorizable, but it.now becomes
is is unaffected. With the definition in Eq. (10), taking
only first-order accurate. Note that the pressure advec-
the curl of Eq. (6) yields q.0 = 0. Thus the advection
tion term in Eq. (7) continues to be approximated by operator acting on the vorticity is unchanged from the
second-order central differences.
first-order scheme.
To obtain second-order accuracy while maintaining To get. full second-order accuracy it. is necessary to
h-ellipticity, appropriate antidissipative term_ nmst be use second-order accurate discretization of the advection
added to the momentum equation in such away that fac-
operator in the momentum equation. The FDA of an h-
torizability is preserved. The form of these terms is de-
elliptic, fully second-order upwind advection operator
pendent upon the computational coordinates (_, Ii), and
they must be written in terms of the contravariant and quo = C0{ + _'0 n
covariant components of the velocity vector. The covari-
ant. components are related to the physical components
by the transformation 1'2 IglO_+2sgn_'f'OeO,+_O,,] (12a)
when1 '1> lv'l
and
*, y,_ v = " (8)
qHO ----_0_ "k-VO,7
W¥iting the advection operator in Eq. (6) in terms of
the covariant velocity components and ignoring terms - 21 (kC_20_ + 2sg,,t 7-c-o{o,, + lvl o,g) (a2b)
containing the higher-order geometric derivatives, the
scheme may be upgraded to nearly second-order by *dg_ = _Yn- dxn, ggn = __y¢ + jx_
when l_?[ < I_::[. With t.his advection operator in the The following difference operators satisfy this condition:
momentum equation, the quantity D in Eq. (11) must ' '[i' i]
be modified in order to preserve the factorizability of tile
scheme. This modified D is given by _ = g 0o , On =g - -20 - ,
, ,[!,!] '[! !]
Dm_ -=D + lsgn_" _" (O,L? + c3¢9) o.=_ -4 , _,=_ - -4 - ,
-2 2
(13at
,,henIVl> I1"1and
0_. _ 0-
DHC) = D -4-lsgn_> f-7(O,7_ -4-O<V) Towrite the complete discrete scheme, the subscript h
is used to denote a standaxd difference to the corre-
(13h)
sponding operator, and the addition of the overbar ( )
denotes the "wide" differences given above. The second-
when IUI < I% Taking tit(" curl of the momentum difference expressions may be expressed in flux form by
equation now gives qno w = 0, i. e., the vorticity equa- taking a six-point difference centered on an edge between
tion is now second-order accurate. two vertices, and then taking a two-point difference of
To see how the multidimensional upwind dissipation of those expressions centered on a vertex. The subscripts e
and v are used to denote difference operators centered
Eqs. (6), (7) and (9) is related to the standard first-order
upwind difference scheme, consider a uniform Cartesian on an edge or a vertex, respectively. The fully discrete
grid with equal spacing in the x and g directions. If the scheme is then written as
scaling coefficients a,, and #v are taken to be one, f is
qhs = 0, (14at
the grid spacing, and the cross-derivative terms in the
dissipation terms of Eqs. (6) and (7) are ignored, the
standard first-order upwind sdmme is obtained. The
first-order upwind scheme on a nonmfiform grid further
p
replaces (_" by g_.0_ + lyon, where f_ and (u are the grid
spacings in the two coordinate directions. The multidi- = O, (14b)
mensional upwind scheme uses a single length scale, and
retains the cross-derivative terms. Consider the advec-
tion operator of Eq. (9) on a Cartesian grid. Dropping
pc2 _Zh.ff ___ff ._ hp
the cross-derivative terms yields the standard first-order
upwind discretization. . h- 1 (V__F,_)p) =0 ' (1@)
Discretization
The derivatives in that part of DH-_ given in Eq. (ll)
The factorizability of the FDA is a necessary condi- are discretized using the wide, six-point difference stencil
tion for the factorizability of the difference scheme, but. for 0_/_ and 0_V. The additional terms of DHO and
it is not sufficient. For the difference scheme to be fac- the corresponding dissipation terms for the second-order
torizat,le, the difference operators must commute in the advection operator qHo must be discretized in a very
same way as the differential opera!ors 8. Introducing the precise way, and depend upon the flow direction and
difference approximations to the partial derivative oper- the relative magnitudes of Cr and V. The details are
at,ors, described in Ref. 12 for uniform Cartesian grids. The
stencils for the general curvilinear grids are presented in
=0¢+..-, 0,7 =0,_+---, the appendix to this paper.
4', =0_+..., 0_,,=0_+..., The scaling coefficients are
1
0-,_ = max(._/, Me), _rv- max(M, M_) (15)
where _I is the local Mach number, and /ff_ is a cutoff
the following conditions must hold Mach number to prevent division by zero. The cutoff is
chosen to be O(h), and essentially becomes active near
•c_Jh_cJ.-,,hm = O,hCn•Oh_n stagnation points. The purpose of the rescaling of the
h h ,h ,h pressure equation dissipation, crp, is to prevent the ellip-
tic factor in the discrete equations from becoming the
skewed Laplacian operator in the incompressible limit.
4
Currently, we take Mc = v IAM[, where u = 5 operators in Eq. (14). The artificial dissipation in
and AM is the two-point difference in Mach number Eqs. (14b) and (14c) can be rewritten as
on an edge. For tile channel flow cases shox_m below, Me
never becomes larger than M. For the leading-edge flow,
the cutoff does become active at the stagnation point. 1 _ h
+V,_ h . ,
When Me > 3/I, it is necessary to add additional dissi-
patton to the advection operator qHho. Tile form of this
dissipation is a five-point pseudo-Laplacian,
1 c( (0_ + 0_), (16)
d_p = _max(0, Mc-M) 7 This is a conservative form of the dissipation cast in
terms of the primitive variables. The terms inside the
where J is the Jacobian of tile coordinate transforma-
square bracket in Eqs. (20) and (21) are now interpreted
tion. This operator is added to both the entropy and the
as dissipative fluxes. These terms are discretized on cell
momentum equation, as_d is cast in flux form in order to
faces according to the appropriate wide or narrow differ-
maintain conservation. No at tenq_t ha._ been made to
ences. The first-order upwind advection operator @ and
optimize either the coefficient v or the form of the oper-
the antidissipative terms /_h(/_,9) in Eq. (20) may be
ator d_p.
recast in conservation form and evaluated on cell faces.
As wilt.ten, the Eq. (14) is valid only for subsonic
The length scale ( is evaluated on the cell face by tak-
flows. This is because tile pressure differences in the
ing min((e¢, g,1), the shortest length in each grid direction.
artificial dissipation terms are not. fully upwinded in a
The scaling coefficients a,,, and o-v are evaluated using
supersonic zone. A simple modification to Eq. (14) can
the Mach number on the face.
be made by rescaling the gradients of the pressure when
In addition to the dissipation terms, examination of
the flow becomes sonic. Introducing the parameter n,
Eqs. (14b) and (14e) shows that. the pressure gradi-
defined as
ent terms are discretized using the wide stencils. The
=max(l, M2), (17) central-difference part of the conservation equations (19)
must be corrected to account for these wide differences.
the final form of the scheme is This is done by adding a term to the dissipative fluxes
for the momentum and pressure equations as follows:
qhs - d_ps = O, (18a)
p
qhl(,a -- d_p_ + V,,D_,, + _lvhp
P
Once the dissipative fluxes (Ss, 5u, Jr, 5p) have been eval-
= o, (lSb) uated, they are converted to the conservation variables
O'm(_Th" (pc2_72"_-I- _"u'_rep)2pc
by the transformation
I -P"I" p o ,,/d] ,_.
pC2_Th. ([ + ff._rhp
L_o_'i = I -p,,Is o ,,lal _,, "
) 0,, h/cV\41
-peT,,,.._u. o_+_ +v_ p =o. (18e)
(24)
where h = e+ p/p is the total ent.halpy.
Because the difference equations (14a), (14b)
and (14c), can be written as a central difference part Solution Procedure
plus a dissipation, it is straightforward to obtain a con-
The solution of the discrete equations is performed us-
servative discretization. The conservation form of the
ing a symmetric point collective Gauss-Seidel iteration,
Euler equations are diseretized using a central-difference
which has been found t.o be a very effective smoother.
finite-volume approximation,
The grid point.s are ordered such that. the forward Gauss-
Vh-(pa) = 0, (19a) Seidel sweep is in the streamwise direction. This means
that the entropy equation is marched in space during the
l_rh.(pff_) + vhp = 0.... (19b) forward sweep. The reverse Gauss-Seidel sweep is under-
v_.[(p_+v)a] = 0 (19e) relaxed with a factor of 0.5 for stability. The residuals
of the conservation equations (r¢,, ro,, ro,,, rp_) are com-
puted at each vertex using the discretization presented
where e = a.ff/2 + P/(O(_ - 1)) is the total energy. in the previous section. These residuals are then trans-
The dissipative fluxes on each cell face are computed formed to the residuals of the primitive variables by the
in terms of (s, u,v,p) using the appropriate difference inverse of the transformation in Eq. (24).
flow [_ ff p .-ighl = I
1
1
Ic,_gth = :_
Figure 1. Channel geometry. \
At each point, the update to tile solution is given by
Figure 2. /vlach number contours for a M = 0.5
M AAuv = - rr,,, (25) inlet, contour increment AM = 0.0I, for a 385 × 129
grid.
Ap ,,_ rv i,3
where M is the matrix of coefficients of the primitive
variables at. vertex (i,3). The coefficients are found by At the inflow boundary the entropy, total enthalpy
collecting the contributions to vertex (i, j) from the dissi- and flow angle are specified, and the pressure-is extrapo-
pation terms on the four surrounding faces. Because the lated from the interior. The outflow boundary condition
entropy equation decouples froln the rest of the system, is a specified pressure and s, u and v are extrapolated
this is a block diagonal matrix where the upper block is from the interior. At the upper and lower walls of the
the upper left-hand entry, and the lower block is a 3 × 3 channell the flow tangency condition tT-fi = 0 is enforced
matrix of coefficients multiplying u,.j, v,,j, and p,,j. This by setting the residual of the momentum equation nor-
matrix is easily inverted. mal t.o the wall to zero. Because of the wide stencils the
The relaxation is accelerated using a standard Full- dissipation terms on the wall are evaluated using a row
Approximation Scheme (FAS) multigrid. A sequence of of ghost vertices. The pressure is extrapolated to those
grids G_c, Gt_-1, .. • ,Go is used, where G1c is the finest vertices using the normal momentum equation. The en-
and Go the coarsest. Let Lk-l be the coarse grid oper- tropy and total enthalpy at the ghost, vertices are set to
ator, u be the vector of the conservation variables, I__, the values on the wall, and the normal component of the
be the fine-to-coarse grid restriction operator, and I_-' velocity is reflected from the interior.
be the coarse-to-fine grid prolongation operator. If fik is Solutions are obtained using a FMG cycle to initial-
the current solution on grid k, the residual on this grid ize the solution. A solution is computed on the coarsest
is rk = fk - Lkfik. This is the residual of the conserva- grid and prolonged to the next grid to ohtaln a starting
tion equations, not the prinfitive equations. This leads solution for that grid. This procedure is continued re-
to the coarse-grid equation cursively until the finest grid is reached. On each grid,
a IS(l, 1) multigrid cycle is used to solve the system.
= = + (2G)
Five multigrid cycles are run on each of the coarse grids,
and fifteen cycles are run on the finest grid. After each
After solving the coarse-grid equation for Uk-!, the fine-
symmetric Gauss-Seidel relaxation sweep, an additional
grid solution is corrected by
three streamwise sweeps are done on the wall and its
+ (27) neighboring row" of vertices.
Mach contours are seen in Fig. 2 on a 385 × 129 fine
Equation (26) is solved by applying the same relaxation grid for an inlet Mach number of 0.5. The solution
procedure that is used to solve the fine-grid equation. is symmetric except for a glitch at the outflow bound-
Multigrid is applied recursively to the coarse-grid equa- ary, which is caused by specifying uniform pressure at
tion. On the coarsest grid, many relaxation sweeps are the outflow boundary. Convergence rates are shown in
performed to insure that the equation is solved com- Fig. 3. A total of 7 grids is used, the coarsest being 7 ×3.
pletely. A conventional Iz cycle or H-cycle is used. The residual levels are renormalized aft.er each prolonga-
tion to the next finest, grid, and on the coarsest grid
Results
several relaxation sweeps are performed, giving essen-
Solutions are shown for compressible flow- in a two- tially a direct solve on that grid. The convergence rate
dimensional channel, the geometry of which is shown in on the finest grid is initially 0.25 per 1"(1, 1) cycle, and
Pig. 1. The shape of the lower wall between 0 < z < 1 the asymptotic rate is 0.38 per cycle. The residual re-
is y(z) = r sin2 rrx. For the computations shown here, duction over 5 cycles is seen to be the same on each grid,
the thickness ratio r is 0.05. The grid spacing is uniform showing that the convergence rate is O(n). The drag on
in the x-direction, and the coordinates in/.he y-direction the lower wall and the rms error in the entropy over the
are found using a simple shearing transformation. computational domain are shown in Fig. 4. On the finest
3,0E-04 0
_oeoc_oeooeo_ooeoo_o
F'Og'O
-0.0002
2.5E-04
°e°/ --_ Entropy -0.0004
"_,,\ 1 -.o - Drag
e -0.0006
2.0E-04
B
-0.0008
_. .4
1.5E-04 -0.001
vo Q
..r -0.0012
___ L,trho-u)i \ 1,0E-04 -0.0014
-0.0016
5.0E-05
-0.0018
10 20 30 40 50 O.OE+O0(_ 10 20 30 40 50
V(1,1)cycles V(1,1) cycles
Figure 3. Convergence rates for a M = 0.5 inlet, flow Figure 4. Convergence of the [-2 entropy error and
using a FMG cycle with a 385 x 129 fine grid. drag on the lower wall for a M = 0.5 inlet flow using
a FMG cycle with a 385 × 129 fine grid.
grids both are converged after essentially two multigrid
cycles. drag and entropy errors converge within two cycles on
The drag and entropy appear to be exhibiting bet- the finest grids.
ter than first-order accuracy but not quite second-order Mach contours, residual convergence and the drag and
accuracy. On the finer grids, the error decreases by ap- entropy are shown in Figs. 8, 9 and 10 for an inlet Mach
proximately one-third with the halving of the grid spac- number of 0.73, which corresponds to a supercritical inlet
ing. The entropy error decreases by about 0.28 on the Mad_ number. The peak Mach number before the shock
finest grids, while the drag is reduced by about 0.35. is approximately 1.2. The rate of residual reduction is
These values are somewhat sensitive to the choice of the initially about 0.53 per cycle on the 385 x 129 fine grid,
wall boundary condition. Curiously, when simple reflec- and the asymptotic rate is 0.68. The drag and entropy
tion boundary conditions are used, the drag converges take longer to converge than for the subcritical c_es,
better than second order, dropping by a factor of six and in fact they have not converged on any of the coarse
to eight with each grid refinement. Neither the drag grids before the solution is interpolated to the next. finer
nor the entropy is very sensitive to whether a first-order grid. On the finest grid, the drag and entropy converge
or a second-order advection operator is used in either in about eight cycles.
Eq. (18a) or Eq. (181)). This is because the flow is a All the solutions shown so far do not have a stagna-
potential flow, and the factorizability of the scheme de- tion point. In earlier work, the pressure Poisson scheme
couples the entropy and vorticity advection from the po- of Robert.s, Sidilkover and Swanson 5' 6 exhibited a de-
tential factor of the operator. Given that s and _0are terioration in the convergence rate for flows containing
identically zero analytically, first-order advection is per- a stagnation point 11. The current scheme does not stif-
fectly adequate. fer from this difficultly. A solution was obtained for the
To demonstrate that the factorizable scheme requires flow about the leading edge of a semi-infinite parabolic
no special preconditioning to handle very low Mach num- body, shown in Fig. 11. For this case, Dirichlet boundary
bers, a solution for an inlet Mach number of 10-4 is conditions were used at the far-field, where the velocity
shown in Figs. 5, 6 and 7. Because the ratio of the field was taken from the incompressible solution given by
dynamic pressure to the static pressure is 0(_I2), the a conformal mapping around the body. The freest.ream
relative roundoff error is much larger for this case than entropy was specified, and a uniform total enthalpy was
for the higher Mach number cases, so the residuals bot- computed by taking the fi'eestream Mach number equal
t.om out at a higher level. The initial convergence rate to 0.1. A fine grid of 129 x 129 grid points was computed
is 0.24 per cycle, very close to that for the Mach num- fi'om the conformal mapping, and 6 grid levels were used.
ber 0.5 case. As with the higher Mach mnnber case, In the symmetric Gauss-Seide] relaxation, it was neces-
the convergence rate is O(n). The solution is seen to sary to underrelax the streamwise sweep with a factor
be very symmetrical at the bump, with some boundary of 0.9 for stability. The wall boundary conditions were
layer beha_ior at the inflow and the same behavior at simple reflection of the velocity, entropy and pressure.
the outflow as for the the higher Mach number. The Pressure coefficient contours are shown Fig. 12, and
are seen to be symmetric about the axis of the parabola.
Convergence rates are shown in Fig. 13. Unlike tile chan-
nel flow, the convergence rates are not quite O(n), be-
coming slightly slower on finer grids. This is because
of the underrelaxation of the streamwise Gauss-Seidel
sweep, which prevents the advection factor from being
solved in that sweep. The asymptotic convergence rate
on the finest grid is approximately 0.50 per cycle. The
t
entropy error, Fig. 14, is seen to converge as well as for
the channel flow, with a reduction of around one-third
with each grid doubling.
Conclusions
A factorizable discretization of the compressible Eu-
ler equations on general curvilinear body-fitted grids has
been presented. The discretization is based on the mul-
tidimensional upwind scheme of Sidilkover, and can be
written in a form that is closely related to conventional
Figure5. Machnumbecrontoursfora
M = 10-4 upwind-differenced finite-volume schemes. Unlike con-
inlet, contour increment AM = 2 x 10-6, for a 385 x ventional schemes, this discretization lends itself to an
129 grid. extremely efficient multigrid solution procedure. The
factorizability of the scheme has the property that the
correct incompressible limit of the compressible equa-
tions is obtained without the need for any special pre-
conditioning. Because the scheme is stable for Gauss-
Seidel relaxation, fast convergence may be obtained for
essentially incompressible flow as well as high subsonic
Mach numbers. This discretization unifies compressible
and incompressible flow algorithms.
References
1. Brandt, A., "Multigrid Techniques: 1984 Guide
with Applications to Fluid Dynamics," GMD-
Studie 85, GMD-FIT, 1985.
2. BrasMt, A., and Yavneh, I., "Accelerated Multi-
grid Convergence and High-Reynolds Recirculating
Flows," SL4M J. Sci. Statist. Comput., vol. 14,
no. 3, pp. 607-626, 1993.
3. Ta'asan, S., "Canonical-Variables MuMgrid
Method for Steady-State Euler Equations," ICASE
Report 94 14, 1994.
-4
4. Ta'asan, S., "Canonical Forms of Multidimensional
Steady Inviscid Flow," ICASE Report 93--3.t, 1993.
L,(mo-u)
....... L,(n".o-v) 5. Roberts, T. W., Sidilkover, D., and Swanson, R. C.,
-6 l t.,(rho)
t.,(rho-.e)
"Textbook Multigrid Efficiency for the Steady Euler
Equations," AL4A Paper 97-1949, 1997.
-%.... ?0.... 2'o.... .... 10.... 20 6. Roberts, T. W., Swanson, R. C., and Sidilkover,
7(1,1)cycles D., "An Algorithm for Ideal Multigrid Convergence
for the Steady Euler Equations," Computers and
Fluids, vol. 28, nos. 4-5, pp. 427 442, 1999.
Figure 6. Convergence rates for a 2_'[= 10-4 inlet
flow using a FMG cycle with a 385 × 129 fine grid. 7. Sanchez, J., "An Optimally Converging Multigrid
Algorithm for the Three-Dimensional Euler Equa-
tions," M.Sc. thesis, George x,¥ashington University,
1998.
8. Sidilkover, D., "Factorizable Schemes for the Equa-
tions of Fluid Flow," ICASE Report 99-20, 1999.