NASA Technical Memorandum 102568

ICOMP-90-12

ULTRA-SHARP Nonoscillatory

Convection Schemes for High-Speed

Steady Multidimensional Flow

B.P. Leonard

Institute for Computational Mechanics in Propulsion

Lewis Research Center

Cleveland, Ohio

and

Simin Mokhtari

Department of Mechanical Engineering

The University of Akron

Akron, Ohio

April 1990

G._I: 4

ULTRA-SHARP NONOSCILLATORY

CONVECTION SCHEMES FOR HIGH-SPEED

MULTIDIMENSIONAL STEADY FLOW

B. P. Leonard*

Institute for Computational Mechanics in Propulsion

Lewis Research Center, Cleveland, Ohio

and

Simin Mokhtari

Department of Mechanical Engineering, The University of Akron

Akron, Ohio

O

t'q

o_

I

SUMMARY

For convection-dominated flows, classical second-order methods are notoriously

oscillatory and often unstable. For this reason, many computational-fluid-dynamicists have adopted various forms of (inherently stable) first-order upwinding over

the past few decades. Although it is now well known that first-order convection

schemes suffer from serious inaccuracies attributable to artificial viscosity or

numerical diffusion under high-convection conditions, these methods continue to

enjoy widespread popularity for numerical heat-transfer calculations, apparently

due to a perceived lack of viable high-accuracy alternatives. But alternatives are

available. For example, nonoscillatory methods used in gasdynamics, including

currently_opular "TVD" schemes, can be easily adapted to multidimensional incompressible 1low and convective transport. This, in itself, would be a major advance for

numerical convective heat transfer, for example. But, as this report shows, secondorder TVD schemes form only a small, overly restrictive, subclass of a much more

universal, and extremely simple, nonoscillatory flux-limiting strategy which can be

applied to convection schemes of arbitrarily high-order accuracy, while requiring

only a simple tridiagonal ADI line-solver, as used in the majority of general-purpose

iterative codes for incompressible flow and numerical heat transfer. The new universal limiter and associated solution procedures form the so-called ULTRA-SHARP

alternative for high-resolution nonoscillatory multidimensional steady-state highspeed convective modelling.

INTROI) UCTION

For many years the state of the art in high-speed convective modelling, especially in the field of numerical heat and mass transfer, has been dominated by firstorder upwinding, often in the guise of the "Hybrid" scheme of Spalding [1] or, more

recently, the related "Power-Law Differencing Scheme" (PLDS) of Patankar [2].

This situation has clearly evolved from an attempt to remedy the infamous problems

of unphysical oscillations and instabilities associated with "classical" central-difference methods under high-convection conditions, using practical grids. The Hybrid

scheme avoids oscillatory behaviour by switching from second-order central to firstorder upwinding for convection (and omitting modelled physical diffusion) wherever

the local component grid Peclet (or Reynolds) number exceeds a value of 2. PLDS

and the exponential differencing scheme (EDS) on which it is based [3] involve a

more subtle blending strategy, but both are also equivalent to first-order upwinding

for convection (with physical diffusion totally suppressed) for component grid Peclet

Work funded by Space Act Agreement C-99066-G.
numbers greater than about 6. For most cases of practical interest, especially in

three dimensions where very fine meshes are out of the question, the actual (laminar

or turbulent) component grid Peclet or Reynolds numbers are likely to be orders of

magnitude larger than 2 or 6 throughout the bulk of the flow domain. This means

that Hybrid and PLDS are operating as first-order upwinding almost everywhere in

the flow-field except at a very small fraction of grid points where the flow direction

happens to be parallel to grid lines or near boundaries and stagnation regions where

the convecting velocities are small. Thus, instead of solving a high-convection

problem, these methods simulate an unphysical (and anisotropic) low-convection

problem in which the effective local component grid Peclet or Reynolds number can

never be greater than 2.

As has been well established theoretically for over a decade [4-6], first-order

upwinding for convection is a highly inaccurate method because of the effects of

global artificial numerical diffusion (or artificial viscosity). Thus, in many cases of

practical interest, Hybrid or PLDS may give rather low-accuracy results; this has

been directly confirmed in several comparative benchmark studies [7-10]. Hybrid

and PLDS, as used in various evolutions of the well-known TEACH code [11], have

been most successful for flows in which the convective terms are not important; i.e.,

for essentially potential flow governed by the kinematic constraint of the continuity

equation (V • v = 0), modified by boundary regions which are usually handled by

built-in "wall-functions", and for turbulence modelling when turbulent-transport

equations are dominated by source-sink terms [12]. However, for simulating recirculating flows or strong shear flows or for predicting scalar transport, these methods

must be considered to be highly unreliable because of the inherent global artificial

viscosity or diffusivity.

Even though these defects are well documented and generally acknowledged,

the Hybrid and PLDS methods continue to enjoy widespread popularity, apparently

because available alternatives claiming better accuracy (such as higher-order

upwinding) are either not well known or have been discarded because of perceived

incompatibilities with TEACH-like solution algorithms [10]. When carefully

applied to three-dimensional laminar recirculating flows however, higher-order

upwinding methods such a QUICK [13] have been extremely successful [14], showing fine-scale details on relatively coarse grids that are simply wiped out by artificial

numerical viscosity in the corresponding Hybrid calculation [15]. In many convection-dominated flow problems, the multidimensional third-order QUICK scheme has

been shown to be stable, economical, and highly accurate; because of the wider computational stencil, the solution algorithm is most naturally cast in pentadiagonal

form [16], but can also be used with conventional tridiagonal solvers provided that

"outlying" and other compensating terms [17] are carefully transferred to the source

term.

However, in some cases (such as turbulent flows, for example) there may be

problems with higher-order upwinding because, in regions involving sudden jumps

in value, there is a tendency for the simulation to overshoot or undershoot the correct

transition value. In many cases, this presents no problem other than slight inaccuracy; however, in codes where turbulent transport variables (such as eddy viscosity,

for example) are computed as part of the solution procedure, it is possible for a

modelled transport coefficient to undershoot to a point where it takes on locally

negative values, thus resulting in violent nonlinear instability. Similar defects

occur in various forms of skew upwinding [18,8]. There have been recent attempts to

suppress the undershoot problem by blending higher-order (or skew) upwinding

methods with standard (component-wise) first-order upwinding, using sophisticated
blending strategies [19,20]. Although much more accurate than Hybrid or PLDS,

these methods still generate an unnecessary degree of numerical smearing when

attempting to simulate near-discontinuities.

The current situation in numerical gasdynamics is considerably better. The

most popular schemesare basedon second-ordercentral differencing for all spatial

derivatives, as used in the canonical second-order Lax-Wendroff method [21], for

example, together with explicitly addedglobal fourth-order artificial dissipation (a

spatial fourth-difference term) to enhance convergence in smooth regions and an

explicit locally varying nonlinear solution-dependent artificial diffusivity or viscosity (a spatial second-difference) to inhibit oscillations in regions of sharply varying

gradient [22,23]. This explicit damping of a potentially oscillatory scheme is usually

quite successful in giving second-order accuracy without excessive oscillation or

gross artificial smearing. In general terms, Hybrid and PLDS could be considered to

follow a similar philosophy: i.e., using second-order central differencing with locally

varying solution-dependent artificial diffusivity or viscosity added explicitly. But a

major difference lies in the strategy adopted for the evaluation of the magnitude of

the artificial terms. With Hybrid and PLDS it is a function of the local component

grid Peclet or Reynolds number, so that artificial diffusivity is always large

wherever physical convection dominates physical diffusion at the grid-cell level,

thus guaranteeing that modelled convection never dominates numerical diffusion.

By contrast, the magnitude of the gasdynamics codes' artificial diffusivity depends in

a nonlinear way on local behaviour of the convected field variables: very low in

smooth regions, but rapidly increasing in regions of strongly varying gradients (high

curvature), irrespective of the local Reynolds number, which could be infinite (as m

Euler-equation simulations). Thus, in these schemes, artificial viscosity is explicitly

added (automatically) only where "needed" (to suppress potential oscillations near

shock waves or contact discontinuities, for example). This technique has been developed into a sophisticated art with so-called "shock-capturing" and "total-variationdiminishing" (TVD) schemes [24,25], which may even involve local negative artificial viscosity, resulting in "artificial compression" to aid in numerically steepening

discontinuities without overshoot. The popular "Superbee" scheme is of this latter

type [26].

Apparently there has never been any attempt to explicitly adapt the highly

successful second-order gasdynamic codes or the so-called "high-resolution" (but

formally still second-order) TVD schemes to incompressible flows and convective

heat and mass transfer, even though (as will be shown) this is extremely simple to

do, using either explicit time-marching or an iterative pentadiagonal matrix

algorithm based on a variable-curvature-factor technique [27,28]. But it turns out

that the commonly used TVD flux limiters form an overly restrictive subclass of a

much more universal technique guaranteeing tight resolution of discontinuities

without overshoots or undershoots; and this "universal limiter" can be applied to

methods of arbitrarily high-order accuracy, giving very sharp resolution without

introducing artificial numerical compression (and concomitant distortion and

clipping of smooth profiles [29]). The strategy developed here is based on a simple

high-accuracy resolution program which uses third-order upwinding in smooth

regions and adaptively increases the accuracy (by locally expanding the computational stencil automatically) in high-curvature regions, while applying the

universal limiter to eliminate any possibility of overshoot or oscillation. Thus, the

Universal Limiter for Tight Resolution and Accuracy in combination with this

Simple High-Accuracy Resolution Program results in the "ULTRA-SHARP" alternative for high-speed nonoscillatory steady-state convective modelling.
In the next section, the primary short-comings of Hybrid and PLDS are briefly

reviewed. Also, the undershoot (overshoot) problems of higher-order upwinding are

demonstrated. The subsequent section describes the Normalized Variable Diagram

(NVD), the construction of the universal limiter and its implicit implementation via

the Downwind Weighting Factor (DWF). The DWF has the added advantage of converting an arbitrarily high-order scheme in N dimensions (using a correspondingly

wide stencil) into a compact (2N + 1)-point scheme involving only scalar tridiagonal

ADI solution methods, thus making the ULTRA-SHARP strategy immediately

available for incorporation into TDMA-based TEACH-type codes. Some of the more

popular TVD schemes are also described within the framework of the NVD and in

terms of the Variable Curvature Factor (VCF) which is most naturally used with a

pentadiagonal solver [27]. Finally, the performance of the ULTRA-SHARP alternative is demonstrated on a number of simple two-dimensional steady-state benchmark test problems for scalar transport in prescribed velocity fields.

SHORT-COMINGS OF HYBRID AND PLDS

Effective Grid Peclet Number

Consider the model convection-diffusion transport equation in one dimension

for a scalar (p

a_ 02¢

--= -u--+DOt c]x Ox2 (1)

for constant u and D. The forward-time, central-space (FTCS) discretization of this

equation takes the form [30]

(P_n + 1 - _ n (_i+1n - _i-In ) (_i+1 n _ 2(1) 7 + _)_-1 n )

- u +D

at 25x _2 (2)

using standard index notation on a uniform grid. This can be rearranged as

n+ 1 rt

#Pi = dPi - c n - 2_7 + #, 1) (*7+,

A

(3)

introducing the Courant number

c = uAt/Ax (4)

and the cell Peclet number (or Reynolds number, in the case of momentum or

vorticity transport).

Pa = 1,,1_/o (5)

Equation (1) can also be discretized using first-order upwinding for the convection term and simply ignoring the physical diffusion term, giving (for u > 0)

4
with a similar form for u < 0. Equation (6) can be rewritten in the form of Equation

(3) as

n+l n C n n ]_ ( n _ 2(_n _}_ n )

¢i = ¢i - _ (¢i +, - _,-, ) + -- ¢i + 1 ¢i-,

Pa (7)

where the effective grid Peclet number is, in this case

P_ ---2 (8)

Thus, the Hybrid scheme - defined by using Equation (3) for Pa < 2 and Equation

(6) for Pa > 2 - can be written in the form of Equation (7) with

P_ = Pa for Pa _ 2 (9)

and

P_ - 2 for Pa > 2 (I0)

Now consider an exact solution of the steady-state form of Equation (1) with a

downstream boundary condition ¢(L) = 1, and with ¢ = 0 far upstream. This can

be written

_(x) = exp[-(L-x)/M

where the length-constant * is given by

for x < L

(11)

(12) _t = D/u

The corresponding central-difference equation is

(2 - Pa)q_i+1 - 4¢i + (2+ Pzx)q_i_t = 0 (13)

and the analytical solution of this difference equation for the above boundary conditions is

2 - Pa )k

_k = _2

where k is an integer defined, for discrete values of x, by

(14)

k = (L - x)lAx for x <- L (15)

Note that the exact solution of the differential equation, Equation (11), can be

written, for discrete node values, as
_k = e- PA) k (16)

Thus, for P_ < 2, the central-difference solution is like a solution to the differential

equation, but at a larger grid Peclet number. And for P. > 2, Equation (14) shows

that the discrete solution becomes oscillatory, as is well _nown [30]. These features

are demonstrated in Figure 1, which shows the factors in Equations (14) and (16) as a

function of Pa.

Now consider Equation (13) rewritten in terms of an effective grid Peclet

number, P j-

(2 - P )¢i+l -- 4_i + (2 + Pa)(_i_l = 0 (17)

The solution of this equation is of the same form as Equation (14), with P_ replaced

by Pa*. But this can be forced to match the true solution, Equation (16), by setting

2 + PA (18)

giving, on rearrangement, an explicit formula for the effective grid Peclet number in

terms of the physical value

O

PA = 2 tanh(Pa/2) (19)

This is the basis of the exponential difference scheme [3] or so-called "optimal"

upwinding [31]. Using P^-dependent weighting factors for convection and diffusion,

Raithby and Schneider [3] develop algebraic approximations (to avoid expensive

exponential evaluations) which can easily be shown to be equivalent to an approximation to Equation (19) given by

, P_ (l + 0.005 P_) I

Pa(R&S) = , + I

2(5 + p (1 + o.05 )

-I

(20)

for all (positive) values of PA. Similarly, Patankar's power-law scheme [2] can be

interpreted as

2 P n

Pn(PLI)S) = 0 -< P_x <- 10

P6 + 2(1 - 0.1PA )5

= 2 PA> 10

(21)

In the same spirit, Spalding's Hybrid scheme [1], given by Equations (9) and (10),

could be interpreted as a (much less accurate) piece-wise linear approximation of the

hyperbolic-tangent function. Equations (9) and (10), (19), (20), and (21) are portrayed in Figure 2. As defined, Hybrid is equivalent to first-order upwinding for

convection with physical diffusion omitted (i.e., PA ° = 2) for Pa > 2. But also note

6
that the other schemesare effectively equivalent to this for P^ greater than about 6.

This is an important point which is sometimes lost in the details of the algebraic or

power-law approximations. In other words, although (as seenfrom Figure 2) Equations (20) and (21) areboth very good engineering approximations to the tanh-curve,

the exponential difference scheme itself is equivalent to first-order upwinding for

convection with physical diffusion switched off, for physical grid Peclet numbers

larger than about 6. This, of course, is entirely appropriate for the model physical

problem on which Equation (19) is based; for Pa > 6, the solution is essentially a

constant (_ = 0) at all computed points, since the effect of the downstream boundary

condition does not diffuse upstream strongly enough to significantly affect the

nearest computed point (e.g., e -6 _ 0.25%). But this (steady, source-free, one-dimensional, convection-diffusion) model problem is clearly very restrictive. For this

reason, it is dangerous to use Equation (19), or various equivalent formulations, as a

basis for a general-purpose computational-fluid-dynamics algorithm. This is clearly

demonstrated in the following section.

Source Terms, Multidimensionality and Transients

To the extent that the fluid-dynamic problem being modelled differs from

Equation (1), Hybrid and PLDS (or EDS) introduce inaccuracies due to the inherent

artificial diffusion - for PA > 2 in the case of Hybrid, and for all P. values for

PLDS. For example, the introduction of a source term into Equation (1_, representing either a real source term in one dimension or the effects of transverse transport

in a multidimensional flow, results in a serious degradation in accuracy when these

methods are used. Figure 3 shows a steady-state two-point boundary-value problem

corresponding to the model equation

/_b o2qb

u_ = D-- + S(x)

Ox _x2 (22)

using a piece-wise linear source term given by [7]

S(x) :: 10- 50 iAx for iAx <- 0.3

S(x) :: 50iAx-20 for 0.3 < lax < 0.4

S(x) :: 0 for lax _ 0.4

(23)

where i is a grid index measured from the left-hand end and x = iAx . Solutions

are shown using both Hybrid and PLDS for Pa = 1, 2, 6, and 10. Note that, for Pa =

i and 2, Hybrid (operating as second-order central differencing) gives adequate solutions; but PLDS is already clearly in error even at these small Pa values. Clearly,

for P^ > 2, the Hybrid solution does not change, since the effective P^* is frozen at

. ¢_ ,, .

• WTth PLDS, the corresponding saturated solutmn occurs at about qPa = 6, corresponding to Figure 2. Note that this is due to the shape of the hyperbolic-tangent

function and has nothing to do with any errors involved in the power-law approximation of tanh.

To demonstrate the effects of multidimensionality, Figure 4 shows the wellknown oblique-step test, corresponding to steady two-dimensional convection and

diffusion governed by the model equation
u--Ox + v--Oy = D --_x'2 + _2 (24)

with a constant convecting velocity, v = (u,v), at an angle 0 to the (uniform square)

grid, and a prescribed unit-step jump in _ on the upstream boundary. As 0 is

varied, the position of the boundary jump is adjusted so that the centerline of the

transition passes through the center-point of the grid. Figure 5 shows analytical

solutions [32] for Pa ( = Ivl Ax/D) = 100, at 0 = 30 °, 45 °, and 60 °. The analytical

solution omits streamwise diffusion and is therefore valid only at moderate-to-high

Pa values [33]. At this grid Peclet number, both Hybrid and PLDS are operating as

pure first-order upwinding for convection with physical diffusion omitted. As seen in

Figure 6, this results in gross artificial diffusion in a direction transverse to the flow

- a well-known phenomenon [4]. Note the large absolute error, given by

= ___.\",_" [qbc,,,.0 - qb xa_t[

, j (25)

summed over all computed grid points on a 25 × 25 square mesh. The effects of

varying 0 (by 1° increments) are summarized in Figure 7. Note the symmetry with

respect to 0 = 45 °, where there is a slight drop in error due to alignment across

diagonal nodes; other undulations are due to similar interactions between the

almost-discontinuous exact profile and the discrete nodal points, as the angle is

varied. Clearly, except when the direction of the convecting velocity is closely

aligned with the grid (0 _ 0 °, 90 ° .... ), the artificial diffusion of these methods at

practical (i.e., large) grid Peclet numbers results once again in highly inaccurate

predictions.

Finally, Figure 8 shows results of simulating the transient pure-convection

equation (for u = const > 0)

--+u-- =0

,3t ax (26)

for an isolated sine-squared profile 20Ax wide (representing a "smooth" function)

and a unit step, using Hybrid or PLDS (or EDS) - equivalent to pure first-order

upwinding for this infinite-P A case - and an explicit forward time step, i.e.,

_, ,1+,= _,"- c ('7- 1) (27)

At the time shown, the exact profiles have moved 45 mesh-widths to the right.

Starting from exact initial conditions, the computed results again demonstrate the

gross artificial diffusion corresponding, in the transient case, to a numerical diffusion coefficient of the form [30]

D = uAx(1 - c)/2 (28)

n U Ill

In the case shown, c = 0.45; from Equation (28), numerical diffusion is clearly worse

at smaller Courant numbers.

8
As the abovesimple test casesshow, Hybrid and PLDS are essentially equivalent to first-order upwinding of convection with physical diffusion ignored for grid

Peclet numbers larger than 2 or about 6, respectively. This must cast serious doubt

on the reliability of simulations using these methods for complex high-speed multidimensional flows using practical (i.e., relatively coarse) meshes. In addition, it

seems to make no sense to use sophisticated (and expensive) multiple-equation

turbulence models with these convection schemes- as is usually done - because

wherever the computed (turbulent) component grid Peclet (or Reynolds) number

exceeds2 (Hybrid) or 6 (PLDS), the turbulent transport terms are (actually or effectively) omitted from the momentum and scalar transport equations (and from the

turbulent transport equations themselves) - having been replaced by the artificial

numerical diffusion inherent in the truncation error of the modelled convection

terms• Apart from the questionable logic involved, one wonders how turbulence

modelling itself has been affected by such convection methods - since, except as an

expensive diagnostic to switch itself off, turbulent transport is largely ignored in

most of the flow domain (except near solid boundaries where universal wall-functions are used anyway).

PROBLEMS WITH HIGHER-ORDER METHODS

Spurious Oscillation of Symmetrical Schemes

Although the artificial diffusion associated with Hybrid and PLDS is widely

acknowledged and well documented, these methods continue to enjoy widespread

popularity [34]• This appears to be because simple alternatives, with formally

higher-order accuracy, introduce other anomalies - usually in the form of overshoots,

undershoots, or severe oscillations - which may lead to convergence problems or nonlinear instability. The primary motivation for Hybrid and PLDS seems to be associated with the suppression of deeply penetrating oscillations which occur with

"classical" second-order central differencing for P, :> 2, as discussed in relation to

• • . ._

Equation (14). Figure 9(a) shows the one-dimensional source-term boundary-value

problem, considered earlier, at PA = 10 solved using second-order central differencing. As seen, the oscillations penetrate from the downstream boundary well into the

solution domain. In fact, the penetration distance is directly proportional to PA [7].

By contrast, using second-order upwinding for convection [30], together with

second-order central differencing for diffusion, written in control-volume (CV) form

(for u :> 0) results in

where S is the CV average, giving nonoscillatory results for this one-dimensional

problem, as seen in Figure 9(b) for P_ = 10. This method also gives much better

accuracy than Hybrid or PLDS; compare with Figures 3(d) and 3(h). This may seem

to be a simple viable higher-order alternative and has, in fact, been adopted as the

basis for several codes, both incompressible [35] and compressible [36]. However,

severe undershoot and overshoot problems develop when using second-order upwinding in two-dimensional computations, as will be shown.

9
Figure 9(c) shows the same one-dimensional problem solved by the QUICK

method [37], or third-order upwinding, which, being a control-volume formulation,

takes the form (for u > 0)

Ax2 - (30)

Typically, a few oscillations are generated near the downstream boundary jump.

But the penetration distance does not continue to grow with Pa; instead, it saturates

at about 5Ax as P^ --, _ [37]. Note, the much better (third-order) accuracy in the

source-term region.-There has been some confusion in the literature concerning the

third-order accuracy of QUICK for convection; this is explained in Appendix A.

Many researchers appear to have an unsubstantiated aversion to unsymmetrical higher-order (upwind) methods, and prefer to work instead with higher-order

central differencing. For example, the fourth-order central discretization in controlvolume form is given by

A-"X u{l'2 (_i+1 + qbt) - _'(_i+2 - qb/+l - qbt + _i-1 ) l

- I'2 (_i + _z-1 ) - 16

'

(¢i+) - _, - _-1 + _,-e) )}

Ax"/[ 2 (_+i _bi) 241((Pi+2 3(Pi+l + 3_bi -(Pi-I ) q

-- l(_i -- _i-1 ) -- _-'(qbt+l - 3qbi +3qbi-1 - _-2 ) I)-_ S(x) (31)

But, as seen in Figure 9(d), although the smooth-region accuracy is good, the downstream oscillation problem is actually worse than second-order! Once again, penetration distance is proportional to P_. This is a characteristic of all higher-order

central-difference schemes; the proportionality constant actually increases with the

formal order of accuracy [38]. From this simple test, it is seen why central-differencing (of any order) is difficult to work with. At this point, second and third (and

possibly higher) order upwind methods seem to offer better alternatives in terms of

accuracy and stability.

Nonmonotonicity of Higher-Order Multidimensional Upwinding

To demonstrate multidimensional effects, Figures 10 and 11 show second-order

upwind and QUICK-2D [13] solutions, respectively, of the oblique-step test for P^ =

100, to be compared with the analytical and first-order solutions of Figures 5 and: 6.

Figure 10 clearly demonstrates the large undershoot and overshoot problems of

second-order upwinding; but note the decrease in error from that of first-order

10
u_)winding (= Hybrid = PLDS). At 45°, QUICK gives qualitatively similar results

to those o_second-order upwinding, but with slightly less undershoot and overshoot

and steeper resolution corresponding to the higher-order accuracy, as reflected by the

smaller absolute error. Note the inverse relationship of undershoot and overshoot

between these two methods as 0 is varied. Figure 12 shows the angular variation

(in 1 ° increments) of absolute error (at Pa = 100) for both methods; cumulative

contributions less than zero (undershoots) or greater than one (overshoots) are also

shown; note the asymmetry between overshoots and undershoots - i.e., an overshoot

at 0 becomes an undershoot at (90 ° - 0) and vice versa, in each case. Again, note the

drop in error at 45 ° and the undular behaviour at other angles. The two-dimensional

oblique-step-test results of this and the previous section have been obtained by

explicit time-marching to steady state; the conservative CV explicit flux formulation

is summarized in Appendix B.

This simple test problem clearly demonstrates the nonmonotonicity problems

associated with higher-order methods in their "unlimited" form. What is needed is a

simple limiting strategy which maintains the good resolution of higher-order upwinding while eliminating undershoots and overshoots without introducing artificial diffusion or destroying conservation. The next section shows how this can be

achieved in a straight-forward manner using an extremely simple technique which

can be extended to arbitrarily high-order accuracy.

THE ULTRA-SHARP ALTERNATIVE

Normalized Variable Diagram for the Universal Limiter

Figure 13 shows the local behaviour of the convected variable, (p(x,y,z), in the

vicinity of a control-volume (CV) face, in a direction normal to that face. Depending

on the direction of the convecting velocity (here shown to the right), label the three

indicated node values: (PD (downwind), (pu (upwind), and (Pc, (centrally located

between the other two), as shown. In part (a) of the figure, locaI behaviour is monotonic. One necessary condition of the universal limiter is sketched in the figure,

shown by the cross-hatched limits: the convected face-value, (p_, should lie between

adjacent node-values in locally monotonic regions. Note that t_s includes the limiting case, (pp = (PD (= (pf)" But if (Pc = (P-' an additional condition is necessary,

namely: (pf _ (Pc ( = (p-),'in this case, as sho_vn in part (b) of the figure. Treatment of

locally non'monotonic +behaviour is sketched in parts (c) and (d) of Figure 13. Some

flexibility is possible here, provided the interpolation is consistent, in terms of continuity, with the previous conditions.

It is more convenient to summarize the limiter constraints in terms of normalized variables. Figure 14(a) shows local behaviour near a CV face in a direction

normal to the face. Note that node-value labelling again conforms to the convecting

velocity direction - in this case, to the left. Now define the normalized variable, (l),

anywhere in the vicinity of this CV face as

_(x,y ,z) - epU

"_(x,y,z) =

(PD - CV (32)

Note, in particular, that 7, = 0 and _, = 1. Figure 14(b) gives the same information as 14(a), but in ter_ns of normali_ed variables. Now the universal limiter

constraints can be portrayed in the Normalized Variable Diagram (NVD) - i.e., the

11
c, (P) plane. This is shown in Figure 15. Symbolically, the universal limiter con

_rain_s on (pf can be written

(PC _ (Pf _ | for 0 < _C -< I (33)

with

_)f : (} at (PC = 0 (34)

and with continuous extensions beyond >__c0> 1 and _: < 0, maintaining a

monotonically rising condition: O_t/a(p C In Figure 1"5, these extensions are

• • f °

given by second-order central-differencing for (Pc > 1:

qbf = 1 + 0.5(qb C- l) for qbC> 1

and by second-order upwinding for negative (Pc"

(35)

_pf : _5 _c for ¢_c < 0 (36)

but other functional relationships are possible. The dotted lines show extensions of

Equations (35) and (36) into the monotonic range, 0 < _c < 1; note that they both

pass throug h the "second-order" point C, located at (0.5, 0.75). To avoid nonuniqueness near i_ _ 0 , the boundary OB has a steep but finite positive slope This

• C .+. .

introduces an additional constraint

dpf <- constOPC near _bC --* 0+

where const = O(100), for example.

(37)

At each stage of a pseudo-time-marching (Appendix B) or iterative solution, the

universal limiter constraints are applied as follows:

(i) For each CV face, note the direction of the (current) normal velocity component, thereby identifying (pu, (Pc, and (PD"

(ii) Explicitly compute a tentative high-order multidimensional upwindbiassed face value, (pf.

(iii) Compute the corresponding normalized face-value, (pf, and the normalized adjacent upwindnode-value, (pc-

(iv) If the point (_c, _r) falls within the triangular region of Figure 15, simply

proceed to the ne:_t CV face.

(v) If not, _)f is limited to the nearest appropriate constraint boundary at the

given (Pc- value.

(vi) The unnormalized face-value is then reconstructed by

qbf = _f(_D- qbU) + qbU

(vii) The same procedure is used for each CV face.

(38)

12
Then in a time-marching algorithm, the (limited) convective fluxes and corresponding diffusive fluxes (computed using second-order central differencing) are now

available for the explicit update step. Alternatively, iterative alternating-direction

implicit solution can be implemented by introducing the downwind weighting factor,

as now described.

The Downwind Weighting Factor

Instead of limiting the face variable directly, the Downwind Weighting Factor

(DWF) is introduced as an auxiliary variable, thereby generating a compact implicit

scheme suitable for tridiagonal solution methods. After explicitly computing the

high-order multidimensional upwind-biassed estimate, qbf, define

DWF -

(39)

Since this is the same as

@f - Oc

DWF - 1 - _'7 (40)

it is not difficult to see that, in terms of DWF, the universal limiter constraints,

corresponding to Figure 15, become

< 1 (41) 0 <_ I)WF <_ 1 for 0 <dp cor

(const - 1)_C

I)WF -< near (PC-* 0

(,- ¢_c) + (42)

in the monotonic region, with, for example,

c

1)WF - for qbC_<O

2(, - ¢c) (43)

and

I)WF ---0.5 for (PC > 1 (44)

This is shown in Figure 16. Note that the point A, given by (1,1) in Figure 15, has

been stretched out into a vertical line in Figure 16. Point C is now at (0.5,0.5).

Now rewrite the face value (in terms of the known DWF) as

_f = DWF q_u "TBC + (1 - DWF) _)cTBC (45)

13
where "TBC" stands for "to-be-computed" - in the next iteration of an implicit linesweep update. Note that the convective flux at each CV face now involves only the

adjacent upwind and downwind "TBC" node-values; but the DWF in Equation (45)

implicitly contains higher-order wide-stencil information (in addition to the universal limiter constraints). Clearly, when DWF's are computed for each CV face on a

quadrilateral grid, the implicit update stencil consists of the central CV node-value

itself, together with only the surrounding adjacent node-values - i.e., a compact

(2N + 1)-point stencil in N dimensions. This is the same compact stencil as used by a

large number of 2D and 3D commercial and research codes of TEACH-like structure,

based on tridiagonal line-sweep solution algorithms, thus immediately opening up

the possibility of incorporating (in principle, arbitrarily) high-order nonoscillatory

convection schemes into these well-established general-purpose elliptic solvers.

Using the popular compass-point notation [2], the individual convected face

values are computed as follows. Considering the west face, for example, the facevalue obtained from Equation (45) will depend on the sign of the convecting velocity

at the west face. If (Pw* is the initial higher order estimate on the west-face, the

DWF is first computed according to

DWF + - if u -> 0

_P - (Pw w (46)

or

I)WF- - if u <0

qbw - Cp _ (47)

Then (the appropriate) DWF is limited according to Figure 16. The face-value used

in the implicit update is then

or

(Pw = DWF: qbp + (1 -DWF:)qb W for u w >-0 (48)

qbw = DWF w qbw + (1 -DWF w)(pP for u w <0 (49)

In either case, _b is a linear combination of its two adjacent (upwind and downwind)

node values, with SGN(u ) as a parameter,

._f DWF[. " "SGN(u )] (50) (l)w = w q)W' (PP' w

Similarly for the other faces of this particular CV cell. This results in an update (in

three dimensions) of the form

ap_p = aw(_ W + as_ S + aB_ B + aE¢_E + aN¢_N + aT¢_T + b (51)

which is, of course, identical to the form generated by first-order upwind or secondorder central methods or combinations such as Hybrid, PLDS, or EDS [2]. Again, it

14
should bestressedthat the higher-order multidimensional information and the nonoscillatory universal limiter constraints are implicitly contained in the DWF's of

each facerather than involving "outlying" node-valuesthat are then lumped into the

explicit source term. The limiter constraints inherent in the DWF's guarantee nonoscillatory results, with stability and convergenceproperties similar to first-order

methods (DWF _ 0) - but without introducing artificial numerical diffusion.

Nonoscillatory Multidimensional Results

Figure 17 shows results of applying the universal limiter constraints to the

QUICK-2D scheme for the oblique step test at P^= 100. These ULTRA-QUICK

results should be compared to the unlimited QUICK results of Figure 11. Note that

resolution remains sharp (reflecting the third-order accuracy) even though there are

no overshoots or undershoots. Also note the concomitant reduction in error seen in

Figure 17(d). But the great power of the universal limiter resides in the fact that it

can be used with higher-order convection schemes. A dramatic increase in sharpness

in the simulation of the near-discontinuity is seen in Figure 18, showing results for

the same problem using ULTRA-fifth-upwinding for convection. As explained in

Appendix C, this scheme uses (one-dimensional) fifth-order upwinding normal to

each CV face, together with third-order (QUICK) transverse curvature terms;

second-order-central terms are adequate for diffusion. As 0 is varied from 0 ° to 90 °,

the total error is never more than about 5.2. Considering the fact that the error is a

cumulative value summed over 625 grid-points, these ULTRA-5th results are seen to

be highly accurate.

Cost-Effectiveness and Adaptive Stencil Expansion

Now that it has been established that viable nonoscillatory higher order

methods can be easily devised, a natural question to ask is: which is more costeffective in terms of overall computer usage (for a prescribed accuracy), a fine-grid

computation using a low-order method, or a coarse-grid calculation using a higherorder scheme? Gaskell and Lau [27] compared first-order upwinding with the thirdorder QUICK scheme and a nonoscillatory version (similar to ULTRA-QUICK),

using successive grid refinement on the two-dimensional oblique-step test. For a

prescribed accuracy, _st-order methods required such a fine grid, compared with

the practical-grid third-order methods, that CPU time was larger by three orders of

magnitude! This disparity would be even more impressive in three dimensions.

Another obvious question is: what order is "optimal" in the above sense of costeffectiveness? Clearly second-order upwinding (which involves the same stencil as

QUICK or ULTRA-QUICK) is not. Higher order nonoscillatory methods such as

ULTRA-5th require more operations (than ULTRA-QUICK) at each grid point; but

perhaps this can be offset by the greater accuracy available on a coarser grid.

Figure 19 gives some indication of optimality in the two-dimensional case.

Part (a) of the figure shows global absolute error versus number of grid-points for

several different methods applied to the infinite-Pa oblique step-test for 0 = 45 °.

Part.(b) shows the corresponding CPU time (without cost-adjustment for storage).

The results, corresponding to a prescribed accuracy from part (a), are cross-plotted in

part (c) of the figure. Even on a logarithmic scale, the "cost" of the first-order method

is off-scale. As expected, third-order methods are significantly better than secondorder upwinding (because of better accuracy for essentially the same cost-per-gridpoint); but fifth and higher-order methods are somewhat less cost-effective than

third. On the other hand, for a given grid-size, Figure 18 shows that ULTRA-5th, in

particular, can give dramatically sharper resolution of discontinuities than the

15
corresponding third-order scheme. This suggests a simple strategy of adaptive

stencil expansion that is proving to be extremely cost-effective in both steady-state

and time-accurate transient simulations [39]:

Use third-order ULTRA-QUICK in "smooth" regions of the flow;

then on the basis of some "non-smoothness" monitor, switch to

a higher-order (ULTRA) scheme locally, as needed.

Since "smooth" regions correspond to small values of the convected variables'

change in gradient (i.e., curvature), a suitable monitor at face (i+1⁄2), for example,

would be the absolute average "curvature" (normal to the CV face):

CURVAV = 0.5 I(qbi+ 2 - qbi+ l) -((Pi - qbi-l)] (52)

This, of course, is also a good indicator of sudden changes in gradient near discontinuities. Thus, if CURVAV is below a prescribed threshold, the algorithm uses the

optimally cost-effective ULTRA-QUICK scheme. Even when discontinuities are

present, this will account for the bulk of the flow domain - since, by definition, discontinuities (when tightly resolved) involve only a very few number of grid points in

narrow isolated regions. Wherever CURVAV exceeds the threshold, the algorithm

automatically expands to a wider stencil, such as ULTRA-5th, and uses this more

accurate scheme in the appropriate local regions. In most cases of practical interest,

the additional "cost" is imperceptible compared with a global ULTRA-QUICK

calculation; however, resolution of near-discontinuities is dramatically enhanced.

Additional thresholds can be used for further adaptive stencil expansion to invoke

(in principle, arbitrarily) higher-order resolution, if desired. In addition, it is

desirable to use a higher-order (ULTRA) scheme where the local absolute normal

"gradient"

GRAD = ]qbi+ 1 - qbil (53)

across the CV face exceeds a given threshold. This strategy has been used up to

ninth-order in inviscid transient calculations [39] in order to resolve extremely

narrow pulses in the convected variable. For steady-state calculations, it appears

that an adaptive ULTRA-QUICK/5th/7th-order convection scheme has a number of

attractive attributes, including cost-effectiveness (high coarse-grid accuracy),

reliability (excellent stability and convergence properties), and ease of coding. The

order-switching strategy is shown schematically in Figure 20. Results of the

oblique-step test (at Pa = 100, with second-order diffusion terms) are seen in Figure

21 to be essentially exact for all practical purposes. Even higher-order adaptive

stencil expansion could be developed for specialty applications.

Variable Curvature Factor

For second- and third-order methods, the normalized convected CV face-value,

_r, is a single-valued function of qb_, and the curve representing this relationship in

the NVD passes through the point U of Figure 15. Passing through C guarantees

second-order accuracy; third-order accuracy (in a control-volume formulation)

requires a slope of 3/4 at C (see Appendix A). In terms of unnormalized variables,

such schemes can always be written

qbf ----- 0.5((_D + qbC)-- VCF.(Qb D- 2_C h- qbU) (54)

16
where VCF - VCF@c) is the variable curvature factor. For third-order methods,

transverse-curvature terms [13] should be addedfor consistency. Equation (54) can

be written in terms of normalized variables, asfollows

(pf = 0.5(1 +(pC )-VCF.(l -2 ) (55)

Rearranging gives

VCF = = VCF(_c)

2(I - 2(p_C) (56)

when (p_(pc) Is known. Note that thin expression is determinate at t(])hCe= 0.5 for

curves passing through point C. It is thus a simple matter to convert universal

limiter constraints of Figure 15 into constraints on VCF; these are shown in Figure

22. For reference, the figure also shows the third-order QUICK value, VCF -= 1/8.

ULTRA-QUICK is shown by the heavy solid lines. Second-order central-differencing

corresponds to VCF - 0; whereas second-order upwinding has VCF - 1/2. Thirdorder accuracy for a given "nonlinear" scheme requires VCF = 1/8 at _c = 1/2.

In an iterative procedure, VCF is computed from Equation (56) based on

current values (of _c )" Then Equation (54) is written in terms of"to-be-computed"

values

• TBC, . TBC -. TBC . TBC,

(pf = 0"5((pTBC+OpC )-- VCF'((PD -2(PC q-q)u ) (57)

implying a pentadiagonal matrix algorithm (PDMA) for line-sweep solution. This

technique was originally developed by Gaskell and Lau [27] using a "curvature

correction", a, from the default value of 1/8; i.e., in terms of the present notation,

°(%) = vcF((p .) - ,/8 (58)

Figure 23 shows another third-order scheme based on an exponential upwinding or

linear extrapolation refinement of third-order upwinding known as EULER [28].

Note the alternative treatment in the nonmonotonic regions.

For reference, a pentadiagonal compass-point algorithm proceeds as follows:

(i) For each CV face, compute VCF + or VCF -, depending on SGN(uf).

(ii) The VCF value may be under-relaxed to enhance convergence [27].

(iii) Then for a typical, say west, face

or

% = (05 - vcF:)(p. + (05 + 2vcF:)(p w - vc . +(pw for uw- o (59)

dl)w - (0.5 - VCF2)(p w + (0.5 + 2VCF2)(p P - VCF w (pE for u < 0 (60)

In either case, (p_ is a linear combination of four adjacent node-values in a direction

normal to the face (two upwind and two downwind), with SGN(u ) as a parameter:

17
_w = "_ VCF[. w q)WW,d_W, d_p , _E ; SGN(uw)] (61)

[Inclusion of transverse curvature introduces additional terms which can be handled

explicitly.] Similarly for the other faces of this CV cell. This results in an overall

three-dimensional update algorithm of the form

apd_p = aWW_W w 4- assd2s S + aBBd_B B + aEEr_EE 4- aNNd_N N + aTT_T T

+ awd_w + as(_s + aB(_B + aE_ E + aN_ N + aT_ T + b

(62)

which is solved most effectively by alternating-direction pentadiagonal line-sweep

iteration [27]. Tridiagonal methods could, of course, be used, taking care to place

outlying - and certain other [17] - terms in the explicit source term; this, however, is

not recommended as it appears to be somewhat less felicitous than the straightforward (and highly efficient) PDMA technique.

Relationship to TVD Schemes

All of the so-called "shock-capturing" and "total-variation-diminishing" (TVD)

schemes currently used in many gasdynamic codes [24] can be immediately adapted

to incompressible convective modelling. Since these methods use the same stencil as

QUICK (without transverse curvature), they can be described within the framework

of the variable-curvature-factor technique. They can also, of course, be represented

as specific curves in the downwind-weighting-factor diagram or the normalizedvariable diagram. For example, Figures 24 and 25 show three well-known schemes

of this type in terms of the NVD and the DWF, respectively: Superbee (a supercompressive scheme), Minmod (a relatively diffusive scheme) [26], and the "harmonic"

scheme devised by van Leer [40].

All shock-capturing and TVD schemes in common use conform to an overly

restrictive limiter [25] which tends to make them more diffusive than necessary. In

terms of normalized variables, this restriction is

_f <_ 2_C for O<qb_ C-<0"5 (63)

Superbee, for example, can be described as a hybrid between second-order central

and second-order upwinding, constrained by Inequality (63) - and by ¢bf _ 1 in the

remainder of the monotonic region. Figure 26 shows the NVD for a scheme similar

in philosophy to Superbee but constrained only by the universal limiter; this might

appropriately be called ULTRA-B. In fact, the basic idea of just such a method was

originally proposed by Roe for gasdynamic methods [41] and, coincidentally, named

"Ultrabee". Figure 27 shows the 45 ° oblique-step test at P. = 100 using Minmod,

the harmonic scheme, and ULTRA-B (or Ultrabee). Note"the diffusive nature of

Minmod; the harmonic scheme gives results similar to those of ULTRA-QUICK; by

contrast, the supercompressive ULTRA-B scheme gives extremely sharp monotonic

resolution - in fact, slightly sharper than the exact solution at this finite grid Peclet

number.

18
Although the step-resolution .properties of simple supercompressive schemes

such as ULTRA-B are very impressive, it should be realized that the sharp resolution is actually obtained by the local introduction of negative numerical diffusion,

which is directly responsible for artificial steepening of discontinuities. Unfortunately, such methods have a tendency to artificially steepen all profiles - with concomitant flattening of peaks (due to mass conservation). This draw-back is seen in

Figure 28, which shows oblique convection without diffusion (at 0 = 45 °) of a semielliptical upstream boundary profile simulated using ULTRA-B; the (infinite-Pa)

exact solution is shown for reference. By contrast, Figure 29(a) shows the same test

simulated by ULTRA-QUICK/5th/7th. Note that, in this case, the sudden change in

gradient at the base of the profile is well captured without the extreme clipping

produced by the supercompressive scheme. However, there is still some loss of

resolution near the local maximum because of the action of the universal limiter.

This problem, in turn, can be largely resolved by switching off the limiter in regions

of true local extrema, as seen in Figure 29(b). The problem then becomes one of

pattern-recognition: near physical local extrema, the limiter needs to be switched

off; but it should remain active near discontinuous jumps or sudden changes in

gradient to suppress spurious unphysical overshoots or undershoots. An automatic

discriminator of this type has been constructed for one-dimensional steady-state [7]

and transient [39] cases. In Figure 29(b), the limiter is simply artificially switched

off for (1) _ 0.5. Clearly, this strategy cannot be used in general. An automatic discriminator for multidimensional steady-state flow is currently under development.

Rotating Velocity Field

To get some idea of the effect of spatially varying velocity, pure convection

(without diffusion) of boundary-specified profiles in a velocity field given by solidbody rotation represents an effective benchmark problem. In a two-dimensional

conservative, finite-volume, formulation, it is convenient to have the streamfunction values available at control-volume corners; then the CV face-average

convecting velocity is simply the difference of the corner stream-function values

across the face. The stream-function used in the proposed benchmark test is

60

o

ty(x,y) = _- [(x --xc)2 + (y -yc )2] (64)

where a) is a specified angular velocity and (x , y ) is the location of the rotation axis.

O ..... C C

The velocity field is shown schematically m Figure 30.

Figure 31 gives the exact solution for a sharp nondiffusing (Pa = oo ) complementary-error-function step in boundary value at the lower left corner. Figure 32

shows results for first-order upwinding ( - Hybrid -_ PLDS --- EDS in this infiniteP^ case). Figures 33 through 35 show unlimited forms of second-, third-, and fifthor-der upwinding. Figures 36 and 37 give Minmod and ULTRA-B results. Finally,

Figure 38 shows the ULTRA-QUICK/5th/7th results. These figures and the corresponding global errors listed in the captions are self-explanatory.

19
CONCLUSION

In terms of normalized variables, first-order upwinding is represented by 7- =

c which, m the NVD, is a straight line passing through (0,0) and (1,1). It thus c_)norms (marginally) to the universal limiter constraints (Figure 15) and, thereby,

produces nonoscillatory results. As shown bby Godunov [42], it is the only "linear"

scheme (i.e., in the present terminology, (pf is a linear function of _) which

possesses this property, as can easily be seen from Figure 15. TVD seKemes are

based on single-valued nonlinear relationships between (pf and (Pc but also require

the corresponding NVD curves to pass through (0,0) and (1',1), witE certain other (in

some respects, unnecessary) restrictions. If the NVD curve passes through (0.5,

0.75), the scheme is second-order accurate (or third-order if it passes through this

point with a slope of 3/4); if not, the scheme is only first-order accurate. Schemes

with generally larger values of _f (at a given (Pc value) tend to be more "compresi 9, • .

s ve than schemes with lower values. Thus, Superbee is more compressive than

ULTRA-QUICK, for example; whereas, Minmod is more diffusive. First-order

upwinding, of course, is seen to be the most diffusive of all nonoscillatory schemes•

Higher-order schemes cannot be represented by a single curve in the ((Pc, (Pf,)

plane; however, they can still be constrained to pass through (0,0) and (1,1), which

appears to be a basic condition for non_cillatory results. In the monotonic range

(0 ¢_ _¢. _ 1), the downwindconstraint, (1)e_ 1, guarantees uniqueness, whereas the

upwind constraint, _r _ (pc, eliminat'es stair-casing. The additional ad hoc

constraint, (p_ _%: cons_ (p-, avoids nonuniqueness near (P- _ 0. The resulting

• f .... C' 1-" v

trmngular region m the _VD, Figure 15, constitutes the steady-flow universal

limiter in the monotonic range. In the nonmonotonic regime, there is considerable

flexibility, provided the constraints are consistent with passing through (0,0) and

(1,1), and t at 0_./0(P c remains positive and finite (for uniqueness). It is not

• _ • . • •

necessary to specifjr a single-valued functmn m this regime; this is usually done

purely as a matter of convenience.

As has been shown, it is a simple matter to design higher-order nonoscillatory

convection schemes and to frame them in terms of simple explicit time-marching

techniques or (by using the downwind weighting factor) traditional TDMA iterative

algorithms. [The special class of second- and third-order schemes (including all

currently used TVD schemes) can also be constructed using straight-forward PDMA

iterative methods, using the variable curvature factor technique - essentially

equivalent to Gaskell-and-Lau's "curvature correction" method.] Existing codes

based on a TEACH-like structure involving iterative line-sweep TDMA solvers can

now be immediately upgraded to incorporate cost-effective higher-order nonoscilla-

.tory methods such as ULTRA-QUICK/5th/7th. Alternatively, because of effieieneies

revolved in vectorization and parallel processing, the simple explicit time-marching

algorithm may be more attractive. In any case, there is certainly no longer any

reason why grossly artificially diffusive methods, such as Hybrid, PLDS (or EDS),

should be used under any circumstances. The availability of ULTRA-SHARP

schemes, giving nonoscillatory high-accuracy results on coarse grids, means that

computational fluid dynamicists will at last be able to focus their attention on

physical modelling of complex fluid phenomena and on reliable practical threedimensional simulations for analysis and design, without being plagued by spurious

oscillations or nonconvergence or by the (less obvious but potentially more

dangerous) debilitating inaccuracies attributable to massive artificial diffusion.

2O
ACKNOWLEDGMENT

Portions of this work were supported by the National

under contract ECS-8904595.

Science Foundation

APPENDIX A

Formal Accuracy of the QUICK Scheme

There has been some confusion in the literature [43] regarding the formal order

of accuracy of the QUICK scheme in terms of truncation error. The problem seems to

have arisen through a confusion of the point-simulation of the convective derivative

with the corresponding control-volume formulation. It is enough to consider the

simplest case of one-dimensional constant-velocity flow using a uniform grid. At

third (and higher) order, there is a difference between

on the one hand, and the integrated control-volume (CV) conservative expression

(Dr -- (De

ax (A.2)

where the subscripts refer to the right and left CV face values of the convected

variable. A Taylor-series expansion of(A.1) results in [30]

/)(D) (Di+1 - (Di-1 (D,+l - 3(Di +3(Di-t - (Di-2 i 2Ax 6Ax + O(Ax3) (A.3)

where an upwind bias is used on the third-difference, assuming a positive convecting

velocity. This could be split into two terms, resulting in a conservative CV form

i Ax (A.4)

where

l 1

(D_ = 2 ((D,+I + (Dr)- 6 ((Di+! - 2(Di + (Di_l) (A.5)

and qbe*is obtained from Cr* by decreasing all indexes by 1.

However, Equation (A.4) does not represent the control-volume expression

given by (A.2) to third order. To see this, make Taylor expansions of the node values

about the CV face locations. For example,

21
$i+1 = cr+(P_ (t) + 2_r , "_- + 6_r 1(t)3 + _'_--r l -_- +''" (A.6)

,flAx) 1 ,,(Ax) 2 1 ,,,(Ax) 3 1 d){iv)(._) 4 (A.7)

_i-1=4r-¢) T +_(l)_,--_- -6_(T +_-_-r X -" (A.8)

So that linear interpolation across the right face would give

_rL = _(¢,+1+_ 1 )=%+% ,,( TAx 2 ) +O(_x% (A.9)

It is easily seen that the 1/6 factor in Equation (A.5) does not completely cancel the

_" term. In order to cancel this term completely, one needs the QUICK formulation

giving

1 1

(PQ = _- (_i+1 + (Pi) - _ (_,+l - 2_,+ Oi_l) (A.10)

1 (p_, Ax 3 + .

.. (A.11)

Thus, for the complete CV term,

H,

(PQ-(P? - (Pr-_e + --1 (dP:'_AxdP_) Ax3 + .-. (A. 12)

Ax Ax 16

and, of course, in the limit Ax _ 0, the term in square brackets becomes (p(iv), so that

Ax Ax

+ O(Ax 3) (A.13)

Since a conservative flux-based formulation should be consistent with the CV form,

the QUICK simulation is the appropriate one to use. Note that several other

discretization schemes can be represented in terms of a "curvature factor"

D 1

_br = 2 ((Pi+l + _bi) - CF "((P,+I - 2¢p,+ (pi_l) (A.14)

(with an analogous expression for negative convecting velocities). For example,

second-order upwinding has CF = 1/2, and Fromm's method [44] has CF = 1/4 (and,

of course, central differencing is equivalent to CF = 0). The distinction between

Equations (A.5) and (A.10) thus hinges on the factors 1/6 and 1/8, respectively. For a

steady-state control-volume formulation, the 1/8-factor gives the appropriate third22
order accuracy. Interestingly enough, for the corresponding time-accurate controlvolume formulation of unsteady convection-dominated flows, an additional factor of

1/24 (together with a Courant-number-dependent term) must be added to the 1/8-

factor, thereby reconstructing the value ofl/6. This was explained - but evidently

not universally understood - over ten years ago, when the QUICK and QUICKEST

algorithms were initially published [37].

In terms of normalized variables, Equation (A.10) becomes

_q

Or = 0.5 (1 - _C ) - 0.125 (1 - 2_C) (A.15)

or, more simply

"_Q(Pr = 0.75 + 0.75 (_'c 0.5) (A.16)

Similarly, Equation (A.14)becomes

_D

_d)r = 0.75 +(2CF + 0.5)-u (_-0'5) (A.17)

For finite CF, any scheme of this form passes through point C of Figure 15, and is

therefore at least second-order accurate, as seen from Equation (A.14). If CF = 1/8,

then the slope S is given by

S = 2CF + 0.5 = 0.75

(A.18)

so that a steady-state control-volume scheme which passes through point C with a

slope of 3/4 is third-order accurate, according to Equation (A.16). Table A.1 summarizes several second-order schemes in comparison with the third-order QUICK

algorithm.

Table A.1

Method CF S Order

QUICK 1/8 3/4 Third

Eq. (A.5) 1/6 5/6 Second

Fromm 1/4 1 Second

2nd-up 1/2 3/2 Second

Central 0 1/2 Second

23
APPENDIX B

Explicit Control-Volume Time-Marching Algorithm

Using compass-point (plus "top" and "bottom") notation, fluxes are computed

explicitly for the west, south and bottom faces. For example, for the west face

F w = Cw_w- a GRAD (B.1)

w W

where the west-face Courant number is

c : u At/Ax (B.2)

W £0

and the corresponding diffusion parameter is

a = D At/ax 9 (B.3)

W W

The face-value is first computed explicitly, using (in general) a wide upwind-biassed

multidimensional stencil; then the universal-limiter constraints are applied - again

explicitly - to give _b,,,. The gradient term is appropriately taken to be a second-order

approximation simurating the diffusive flux; i.e.,

GRADw = (_)w_ Ax = _e - _w (B.4)

[A higher order treatment of diffusion is not necessary because, for sharp changes in

gradient, diffusion is very small, whereas when diffusion is large, gradients vary

smoothly.] Then the update algorithm consists of a FORTRAN over-write assignment statement inside a multiple DO-loop:

Set:

¢_p(i,j, k) = C_p{i,j, k) + Fw(i,j, k) + Fs(i,j, k) + Fb(i,j, k)

- Fw(i+ 1 ,j, k) - Fs(i,j+ 1, k} - Fb(i,j, k + 1 ) + At S(i,j, k} (B.5)

where conservation has been used - e.g., for the east face at (i,j, k),

Fe(i,j, k) : Fw(i + l ,j, k) (B.6)

and so on. Because of nonlinearities inherent in the limiter, it is not possible to give

precise stability restrictions on the time-step. However, experience suggests that

local Courant numbers near 0.3 appear to give satisfactory results under highconvection conditions.

The steady-state limiter portrayed in Figure 15 needs to be slightly modified

for the time-marching algorithm. Rather than using the ad hoc restriction given by

OB, it is necessary to use a more precise condition corresponding to the timeaccurate limiter [29], i.e.,

24
_f _ _c./c for 0 - _c -- I (B.7)

(in addition to requiring _f =< 1 in the same range) where c is the appropriate

normal-direction Courant number at the CV face in question. Without this constraint, (very) slight undamped temporal oscillations may occur at isolated points

near strong discontinuities at finite Pa values; at such points, a steady state might

never be reached.

Although only first-order accurate in time, Equation (B.5) can give some indication of the temporal evolution of the solution from the assumed initial conditions.

From a practical point of view, this can indicate the total time necessary for reaching

a steady state. One simply observes how long it takes for initial transients near an

upstream boundary to "wash out" of a downstream boundary; then continued timemarching for 2 to 3 times this time-scale is usually enough to cause reflections and

other transients to die away. Of course, more precise stopping criteria can easily be

devised.

APPENDIX C

Higher-Order Multidimensional Fluxes

For first- and second-order methods - including combinations such as Hybrid,

EDS, PLDS, and various TVD schemes - the multidimensional fluxes at each CV

face, are identical to their normal-direction one-dimensional counterparts. For

second-order upwinding, for example, the left-face convective flux is given by upwind

linear extrapolation in the normal direction,

CFI,UXL(i,j,k) = CXL.(1.5qbi_l,j, k - 0.5(pi_2,j, k)

(C.1)

where CXL is the appropriate left-face normal (x-component) Courant number. Note

that, to a consistent order, no "transverse curvature" terms are included.

At third-order, it is appropriate to include transverse curvature; thus, for

example, the QUICK left-face convective flux is given by

CFLUXL(i,j,k) = CXL-(PeQ (C.2)

where

using the definition

1

(PQ = (PLIN -- 8-1CURVLN + 24-- CURVTL (C.3)

1

and where, ifCXL _ 0, the normal and transverse "curvature" terms are given by

CURVNL = ((Pi,j,h -- 2d_i-l,j,k + (Pi-2,j,k) (C.5)

25
and

CURVTL = (dpi_l,j+l,h - 2dpi_l,j,k + ¢i_l,j_l,_)

+ (_i-l,j,k+l

respectively; whereas, if CXL < 0,

-- 2_i - 1,j, k + (l)i - I,j, k- 1)

and

CURVNI, = (qbL+l,j, k

(C.6)

- 2¢i,i, k + ¢,__,j,_) (C.7)

CURVTL : (_i,j+l,k - 2_bi,j,k + _i-l,j-l,h ) + (_i,j,h+l -- 2_i,j,h + _bi.j,k-l) (C.8)

For higher-order convective fluxes, other transverse terms could be included.

However, numerical experimentation has shown that

(i) whereas inclusion of the third-order CURVT terms provides a significant

increase in accuracy,

(ii) higher-order transverse terms can actually reduce accuracy in some eases and

are therefore omitted, and

(iii) the coefficients of higher-order normal terms can be optimized to give better

accuracy than that obtained from values using formal interpolation formulas.

Thus, for example, the "fifth-order" convected value (at the left face) used in

this report is given by

1

m(sj 1 __3 FOURT]I + -- CURVTI,

"re = d_tjN - 6 CURVAV + 128 24 (C.9)

where CURVAV is an "average normal curvature" across the left face

1

CURVAV = 2 (qbl+l,j,k -- dPi,J,k -- _i-l,j,k + _i-2,j,k )

FOURTH is the upwind-biassed fourth-difference normal to the face

(C.10)

FOURTll = (qbi_ 1,j,k -- 4_i,j,k + 6_i-l,j,k -- 4dPl-2,j,k + qbi-a,j,k) (C.11)

(for CXL > 0), and CURVTL is as defined previously. Note the use of 1/6 (rather

than 1/8) in Equation (C.9)!

The "seventh-order" formula used in this report is given by

3 1 1

fh(7 ) 1 CURVAV + m FRTtIAV + _ SIXTIt + m CURVTI,

_'¢ = _LIN- 6 128 100 24 (C.12)

where FRTHAV is the (symmetrical) average fourth-difference across the left face,

26
FRTiiAV = (d_i+3 - 3¢i+2 + 2qbi+ I + 2qbi - 3qbi_ 1 + qbi_2) (C.13)

(suppressingj's and k's), and SIXTH is the upwind-biassed sixth-difference normal to

the face

SIXTH = (qbi+ 3 -6qbi+ 2 + 15¢i+1 + 20qbi + 15_bi_ l - 6qbi_ 2 + qbi_3) (C.14)

Formal interpolation would indicate factors of 1/8 multiplying CURVAV and 5/1024

multiplying SIXTH - together with additional transverse and cross-difference terms.

The factors used in Equation (C.12) have been found to give slightly more accurate

results.

In all cases, diffusive fluxes are modelled using the simple second-order approximation for the normal gradient. For example, for the left face,

_) (¢i,j, k --¢i-l,j,k )

e ax (C.15)

As explained previously, higher order terms are not warranted since, when diffusion

is large, profiles are smooth and Equation (C.15) is consistent with the appropriate

(QUICK) treatment of convection [13]; whereas, when steep gradients occur- necessitating higher-order convective modelling and the use of the universal limiter -

diffusion terms are very small, so that a simple diffusion model suffices.

27
REFERENCES

1. Spalding, D.B. (1972), "A Novel Finite Difference Formulation for Differential

Expressions Involving Both First and Second Derivatives," International

Journal for Numerical Methods in Engineering, 4,551.

2. Patankar, S.V. (1980), Numerical Heat Transfer and Fluid Flow, Hemisphere

Publishing Corporation, McGraw-Hill.

3. Raithby, G.D. and G.E. Schneider (1988), "Elliptic Systems: Finite-Difference

Method II," Handbook of Numerical Heat Transfer, W,J. Minkowycz, E.M.

Sparrow, G.E. Schneider and R.H. Pletcher, eds., Wiley, 241.

4. deVahl Davis, G. and G.D. Mallinson (1976), "An Evaluation of Upwind and

Central Difference Approximations by a Study of Recirculating Flow,"

Computers and Fluids, 4, 29.

5. Leonard, B.P. (1978), "A Consistency Check for Estimating Truncation Error

Due to Upstream Differencing," Applied Mathematical Modelling, 2,239-244.

6. Raithby, G.D. (1976), "A Critical Evaluation of Upstream Differencing Applied

to Problems Involving Fluid Flow," Computer Methods in Applied Mechanics

and Engineering, 9, 75.

7. Leonard, B.P. (1981), "A Survey of Finite Differences with Upwinding for

Numerical Modelling of the Incompressible Convective Diffusion Equation,"

Computational Techniques in Transient and Turbulent Flow, C. Taylor and K.

Morgan, eds., Pineridge Press, Swansea, UK, 2, 1-36.

8. Huang, P.G., B.E. Launder, and M.A. Leschziner (1985), "Discretization of

Non-Linear Convection Processes: A Broad-Range Comparison of Four

Schemes," Computer Methods in Applied Mechanics and Engineering, 48, 1.

9. Leschziner, M.A. and W. Rodi (1981), "Calculation of Annular and Twin

Parallel Jets Using Various Discretization Schemes and Turbulence-Model

Variants," ASME Journal of Fluids Engineering, 103,352-360.

10. Syed, S.A., L.M. Chiappetta, and A.D. Gosman (1985), "Assessment of

Discretization Schemes to Reduce Numerical Diffusion in the Calculation of

Complex Flows," AIAA paper 85-0441.

11. Gosman, A.D. and W.M. Pun (1974), "Calculation of Recirculating Flows,"

Heat Transfer Section Report HTS/74/2, Department of Mechanical

Engineering, Imperial College, London.

12. Leschziner, M.A. and K.P. Dimitriadis (1989), "Computation of ThreeDimensional Turbulent Flow in Non-orthogonal Junctions by a BranchCoupling Method," Computers and Fluids, 17, 371-396.

13. Leonard, B.P. (1988), "Elliptic Systems: Finite Difference Method IV,"

Handbook of Numerical Heat Transfer, W.J. Minkowicz, E.M. Sparrow, G.E.

Schneider, and R.H. Pletcher, eds., Wiley, New York, Ch. 9,347-378.

14. Perng, C.Y. and R.L. Street (1989), "Three-Dimensional Unsteady Flow

Simulations: Alternative Strategies for a Volume-Averaged Calculation,"

International Journal for Numerical Methods in Fluids, 9,341-362.

15. Kossef, J.R., R.L. Street, P.M. Gresho, C.D. Upson, J.A.C. Humphrey, and

W.M. To (1983), "A Three-Dimensional Lid-Driven Cavity Flow: Experiment

and Simulation," Proceedings of the Third International Conference on

Numerical Methods in Laminar and Turbulent Flows, Seattle, 564-581.

28
16. Gaskell, P.H., A.K.C. Lau, and N.G. Wright (1987), "Two Efficient Solution

Strategies for Use with tIigh Order Discretisation Schemes, in the Simulation

of Fluid Flow Problems," in C. Taylor et al. (eds), Numerical Methods in

Laminar and Turbulent Flow, Vol. 5, Pineridge Press, Swansea, 210.

17. Freitas, C.J. (1986), "Non-Linear Transient Phenomena in a ThreeDimensional Cavity Flow: A Numerical Investigation," Ph.D.Dissertation,

Department of Civil Engineering, Stanford University.

18. Raithby, G.D. (1976), "Skew-Upstream Differencing for Problems Involving

Fluid Flow," Computer Methods in Applied Mechanics and Engineering, 9,151.

19. Zhu, J. and M.A. Leschziner (1988), "A Local Oscillation-Damping Algorithm

for Higher-Order Convection Schemes," Computer Methods in Applied

Mechanics and Engineering, 67,355-366.

20. Turan, A. and J.P. VanDoormaal {1988), "Improved Numerical Methods for

Turbulent Viscous Recirculating Flows," NASA Contractor Report 180852,

Lewis Research Center.

21. Lax, P.D. and B. Wendroff (1964), "Difference Schemes with High Order of

Accuracy for Solving Hyperbolic Equations," Communications on Pure and

Applied Mathematics, 17,381.

22. MacCormack, R.W. (1969), "Effect of Viscosity in Hypervelocity Impact

Cratering," AIAA Paper, 69-354.

23. Jameson, A. and T.J. Baker (1986), "Euler Calculations for a Complete

Aircraft," Proceedings of the Tenth International Conference on Numerical

Methods in Fluid Dynamics, Lecture Notes in Physics, Vol. 264, SpringerVerlag.

24. Yee, H.C. (1987), "Upwind and Symmetric Shock-Capturing Schemes," NASA

TM-89464, NASA-Ames Research Center.

25. Sweby, P.K. (1984), "High Resolution Schemes Using Flux Limiters for

Hyperbolic Conservation Laws," SIAM Journal of Numerical Analysis, 2 l, 995.

26. Roe, P.L. (1986), "Characteristic-Based Schemes for the Euler Equations,"

Annual Reviews of Fluid Mechanics, ! 8, M. Van Dyke, J.V. Wehausen, and J.L.

Lumley, eds., Annual Reviews Inc., 337.

27. Gaskell, P.H. and A.K.C. Lau (1988), "Curvature Compensated Convective

Transport: SMART, a New Boundedness Preserving Transport Algorithm,"

International Journal for Numerical Methods in Fluids, 8,617.

28. Leonard, B.P. (1988), "Simple High Accuracy Resolution Program for

Convective Modelling of Discontinuities," International Journal for Numerical

Methods in Fluids, 8, 901-927.

29. Leonard, B.P. (1988), "Universal Limiter for Transient Interpolation Modeling

of the Advective Transport Equations: The ULTIMATE Conservative

Difference Scheme," NASA TM 100916 (ICOMP-88-11), NASA-Lewis

Research Center.

30. Fletcher, C.A.J. (1988), Computational Techniques for Fluid Dynamics, Vols. I

and II, Springer-Verlag.

31. Heinrich, J.C. and O.C. Zienkiewicz (1979), "The Finite Element Method and

'Upwinding' Techniques in the Numerical Solution of Convection Dominated

29
Flow Problems," in Finite Element Methods for Convection Dominated Flows,

T.J.R. Hughes, ed., AMD 34, ASME, 105-136.

32. Huget, R.G. (1985), "The Evaluation and Development of Approximation

Schemes for the Finite Volume Method," Ph.D. Thesis, University of Waterloo.

33. Leschziner, M.A. (1980), "Practical Evaluation of Three Finite Difference

Schemes for the Computation of Steady-State Recirculating Flows," Computer

Methods in Applied Mechanics and Engineering, 23,293.

34. Lewis, R.W. and K. Morgan, eds. (1989), Numerical Methods in Thermal

Problems VI, Pineridge Press.

35. Castro, I.P. and J.M. Jones (1987), "Studies in Numerical Computations of

Recirculating Flows," International Journal for Numerical Methods in Fluids,

7,793.

36. Warming, R.F. and R.M. Beam (1976), "Upwind Second-Order Difference

Schemes and Applications in Aerodynamic Flows," AIAA Journal, 14, 1241.

37. Leonard, B.P. (1979), "A Stable and Accurate Convective Modelling Procedure

Based on Quadratic Upstream Interpolation," Computer Methods in Applied

Mechanics and Engineering, ! 9, 59-98.

38. Mokhtari, S. and B.P. Leonard (1987), "Analysis of Oscillatory Convection

Schemes," Bulletin of the A merican Physical Society, 32.

39. Leonard, B.P. and H.S. Niknafs (1990), "Sharp Monotonic Resolution of

Discontinuities without Clipping of Narrow Extrema," Computers and Fluids,

(in press).

40. van Leer B. (1974), "Towards the Ultimate Conservative Difference Scheme.

H. Monotoncity and Conservation Combined in a Second-Order Scheme,"

Journal of Computational Physics, 14,361.

41. Roe, P.L. (1985), "Some Contribution to the Modelling of Discontinuous Flows,"

Large Scale Computations in Fluid Mechanics, Part 2 (Lectures in Applied

Mathematics 22), E. Engquist, S. Osher, and R.C.J. Sommerville, eds.,

American Mathematical Society, 183.

42. Godunov, S.K. (1959), "Finite Difference Method for Numerical Computation of

Discontinuous Solutions of the Equation of Fluid Dynamics," Matematik

Sbornik, 47, 271.

43. Bradley, D. and M. Missaghi (1988), "A Taylor-Series Approach to Numerical

Accuracy and Third-Order Scheme for Strong Convective Flows," Computer

Methods in Applied Mechanics and Engineering, 69, 133-151.

44. Fromm, J.E. (1968), "A Method for Reducing Dispersion in Convective

Difference Schemes," Journal of Computational Physics, 3,176.

3O
i

oo

• __ , 2_._ '

'_ 76_ ,'.oo 4.00

0 /

i

5'.O0 6'.O0

FIGURE I. - ALGEBRAIC AND EXPONENTIAL FACTORS APPEARING IN EQUATIONS (14) AND (16),

RESPECTIVELY, AS A FUNCTION OF GRID PECLET NU/'giER,P&.

o

o

4

°.

• i i

°0_, O0 I i. O0 21 • O0 3'. O0 4'. O0 5. O0 6. O0

FIGURE 2. - APPROXIMATIONS OF THE HYPERBOLIC-TANGENT FUNCTION GIVEN BY EQUATION (19).

EXACT CURVE: LIGHT CONTINUOUS LINE. HYBRID: HEAVY PIECEWISE LINEAR APPROXIMATION.

PATANKAR POWER-LAW APPROXIMATION: TRIANGLES. RAITHBY-AND-SCHNEIDER ALGEBRAIC APPROXIMATION: CIRCLES.

31
0

(a) HYBRID: P,, = 1.

X

O

/

X

(b) HYBRID: PA = 2.

_®

/_i-_-'_}'i- I _1_i 1j

X

(c) HYBRID: PA = 6.

It ' -I[

iiI-II111111 11'

X

(d) HYBRID: PA = 10.

_O

(e) PLDS: P_ = 1.

A_

(f) PLDS: PA = 2.

_0

(g) PLDS." PA = 6.

FIGURE 3. - STEADY-STATE BOUNDARY-VALUE PROBLEM WITH SOURCE TERM GIVEN BY EQUATION (23).

32
Ob = 1

®b = 0

0=1

_J

/

j" O

v J

/

I

®=0

.1.1.1.1

Sb = 0

FIGURE q. - SCHERATIC PORTRAYALOF THE OBLIOUE-STEP

TEST. AS 0 IS VARIED, THE LOCATIONOF THE BOUNDARY

DISCONTINUITY IS ADJUSTEDSO THAT THE EENTERLINE OF

THE TRANSITION PASSESTHROUGHTHE RIDPOINT OF THE

6RID.

(a) B = 30 °.

(b) O = 45°.

(c) 9 = 60°.

FIGURE 5. - THREE-DIfflENS|ONALPLOT OF¢(x,y) FOR THE ANALYTICAL

SOLUTION OF THE OBLIQUE-STEP PROBLERFOR PA (= I Vl Ax/D) = 100.

33
(a) 0 = 30° , _ = 50.6.

(b) B = q50 , _ = 6q.2.

%

60

40

20

0

Y

FIGURE 7. - VARIATION OF TOTAL ABSOLUTEERROR,_, tilTH STREAN-TOGRID ANGLE, e, IN 1o INCREflENTSFOR PLDS AT P.,. = 100 (Eg(J]VALENT

TO FIRST-ORDER UPgINDING EXCEPT NEAR 0° AND 90%. RAXIP_M _ IS

65.8 AT qqO OR 116o.

(c) e = 60° , £ = 50.6.

FIGURE 6. - CORRESPONDINGSOLUTIONS OF THE OBLIOUE-STEP PROBLEM

USING FIRST-ORDER UPWINDING--EQUIVALENT TO HYBRID, PLDS, OR EDS

AT THIS PECLET NUMBER (PA = 100). _ IS DEFINED IN EQUATION (25).

34
FIGURE 8. - FIRST-ORDER-UPIAIND SOLUTION OF THE INVISCID TRANSIENT

CONVECTIONPROBLF_FOR A UNIT STEP AND AN ISOLATED SINE-S_JARED

PROFILE, AFTER 100 TIIqE-STEPS AT c = O.qS.

(a) SECOND-ORDERCENTRALDIFFERENCING.

X X

(c) QUICK.

A_

X

(b) SECOND-ORDERUPWINDING. (d) FOURTH-ORDERCENTRAL DIFFERENCING.

FIGURE 9. - STEADYSOURCE-TERRBOUNDARY-VALUEPROBLEROF FIGURE 3 FOI_P_, = 10.

35
(b) e = 45°, # = 19.3.

(c) e = 600 • # = 21.2.

FIGURE 10. - RESULTS OF THE OBLIQUE-STEP PROBLEMFOR P& = 100 USING

SECOND-ORDER UPWINDING,

(b) 0 - 45O, _= 11.5. " ......

(c) e = 60°, _ = 13.8

FIGURE 11. - RESULTS OF THE OBLIQUE-STEP PROBLEMFOR P& = 100 USING

QUICK-2D.

36
30

20

10

0

-10

eaee

eeD_e eel%

0o q5o .................................... 90o

(a) SECOND-ORDERUPJ_INDING (NAXINIffi = 22.5 AT 250 OR 65°).

2O

10

%

O

-10 Oo q5o 900

(b) QUICK 2-D (INU(IFIUR = 13.8 AT 300 OR 60°).

FIGURE 12. - ANGULARVARIATION OF ABSOLUTEERROR, _, AT PA = 100.

UNDERSHOOTSARE SHOh'NBY f4INUS SIGNS. OVERSHOOTSARE SHOkll BY

PLUS SIGNS.

0 D

$C

v

(a) LOCALLYRONOTON[C.

0C

®U

(c) OC > SU"

0 U

v

y

I

I Uf

v

(b) ®C = ®U"

)®D

f

(d) ®C < ®U'

I® D

FIGURE 13. - CONVECTEDFACE VARIABLE, ®f, DEPENDINGON THE LOCAL BEHAVIOR OF ® NORNALTO THE CV FACE.

37
. _ _

ooY

T*!

v

(a) ORIGINAL UNNORRALIZEDVARIABLES.

cz/

I

l Uf

(b) NORNAL]ZEDVARIABLES.

FIGURE lq. - DEFINITION OF NORRALIZEDVARIABLES ACCORDINGTO EQUATION (32).

C

/

/

FIGURE 15. - UNIVERSAL LIRITER CONSTRAINTSIN THE NORNALIZEDVARIABLE

DIAGRAN (NVD). POINT C IS LOCATEDAT (0.5, 0.75).

38
c®

!

/

/

/

A

FIGURE 16. - UNIVERSAL LIRITER CONSTRAINTS(CORRESPONDINGTO FIGURE 15) IN TERRS

OF THE DONNWINDWEIGHTING FACTOR (DM:). POINT C IS AT (0.5, 0.5). NOTE THE

SINGULARITY AT ¢C = 1.

39
(a) 0 = 300 , 8 = 9.2. (a) 8 = 30°, _ = 5.0.

_10 lO

(b) e = q5°, 8 = 8.6.

1

(c) 8 = 600 , _ = 9.2.

0o q5o 900

(d)8 AS A FUNCTION OF 8 (RAXIRUR = 10.2 AT 430 OR q7°).

FIGURE 17. - ULTRA-OUICK OBLIgUE-STEP RESULTSFOR PA = 100.

(b) 0 = 45 o, # = 2.8.

(c) e = 600 , # = 5.0.

lO _-

°i

0o q5o 90o

(d) £ AS A FUNCTION OF 9 (RAXIRUR = 5.2 AT 100 OR 80°).

FIGURE 18. - ULTRA-STH-UPWINDOBLIOUE-STEP RESULTSFOR PA = 100.

4O
.10

.05

0

100

10

[]

F£PRESCRIBED

-- / N

..Z_._ 9

I I I I I_I,L I

[]

D

0

2e 8

I I I I,I,I I I I Illil

(a) ERRORAS A FUNCTION OF GRID REFINEMENT.

D

I I I I_'I'III I i l i lllll I

10 100 1000

N

(b) CPU TIME AS A FUNCTION OF GRID REFINEMENT.

UPWIND

• FIRST

rl SECOND

O THIRD

O FIFTH

SEVENTH

_ _i.._ _i

L_gl

IST 2ND 3RD 5TH 7TH

(c) cPu TIME FOR VARIOUS SCHEMES CORRESPONDING TO THE PRESCRIBED ERROR, PART (a).

FIGURE 19. - COST-EFFECTIVENESS OF HIGHER-ORDER SCHEMES.

FIGURE 20. - SCHEMATIC DIAGRAM OF ORDER-SWITCHING STRATEGY USED IN THE LOCALLY ADAPTIVE

STENCIL EXPANSION ALGORITHM.

41
Ca) e = 30°, _'= q.3.

(c) 0 = 600 , 8 = q,3,

FI6URE 21, - ULTRA-3RD/5"rH/TTH-UPklIND OBLIOUE-STEP RESULTS FOR

PA = 100.

42
-I .0

VCF

0.5

/

/

/

/

/

/

/

/

/

-0.54

Jcl1

\I

I

I

I

A

.0

_c

FIGURE 22. - UNIVERSAL LIRITER BOUNDARIES PORTRAYED IN TERMS OF THE VARIABLE CURVATURE

FACTOR. NOTE THE SINGULARITY AT_ C = 0.5. THE STRAIGHT LINES THROUGH 0 AND A CORRESPOND TO THOSE IN FIGURE 15. THE DASHED LINE SHOWS THE QUICK VALUE: VCF_I/8.

VCF

-1.0 0 1.0 _C

FIGURE 23. - VARIABLE CURVATURE FACTOR FOR THE EULER-QUICK SCHEME. THIRD-ORDER ACCURACY

REQUIRES VCF= I/8 AT_C= I/2, NOTE THE ALTERNATE TREATMENT IN THE NONMONOTONIC REGIONS

43
_f /

1

_c

FIGURE 2q. - NORMALIZEDVARIABLE DIAGRA/_FOR SUPERBEE

(HEAVY PIECEW]SE LINEAR GRAPH), THE HARMONICSCHEME

(DOTTED PARABOLICCURVE), AND MINMOD (DASHEDLINES).

ALL SCHEMESFOLLOW®f = ¢C IN THE NONI'IONOTONICREGIONS.

DgF

D

A

_c

FIGURE 25. - DOWNWINDWEIGHTING FACTORVARIATION FOR THE

SUPERBEE(HEAVY LINES), HARMONIC(DOTTED), AND MINI_)D

(DASHES) SCHEMES.

_f

I IIFIGURE 26. - NORMALIZEDVARIABLE DIAGRAMFOR THE ULTRA-B

SCHEME. COMPAREWITH SUPERBEEIN FIGURE 2q. NOTETHE

SIGNIFICANT DIFFERENCE NEAR _t - 0+.

44
(a) NINROD, _ = 26.4.

(b) HARPI)NIC SCHEME, _ = 17,1.

(c) ULTRA-B, _ = 4.0. (b) EXACT SOLUTION FOR REFERENCE.

FIGURE 27. - OBI_IQUE-STEPRESULTSOFR e = q5° AND PA = 100. FIGURE 28. - OBLIOUE CONVECTION (AT 0 = q5°, PA = 0o) OF AN

ELLIPTIC PROFILE IMPOSEDAT THE INFLOW BOUNDARY.

45
(a) ULTRA-SRD/STH/7TH, _ = 10.7.

(I)1 ULT'RA-3RD/STH/7TH WITH LIMITER-RELAXATION NEAR LOCAL

RAXIMA, _g = 10,1.

FIGURE 29. - OBLIOUE CONVECTION OF THE ELLIPTIC PROFILE.

Ob= 1

/. ,oo

/,/

0 b = 0

(x c, Yc )

FI6URE 30. - SCHEMATIC DIAfitINq OF THE ROTATIIglVELOCITY-FIELD _Y-STEP PROBLER.

FIGURE 31. - EXACT SOLUTION OF THE ROTATING-VELOCITY PROBLF_.MFOR

P,,. = _ AND A NARROM COI_LEJ'IENTARY-ERROR-FUNCTION PROFILE

AT THE UPSTREAR BO(JNORY CORNER.

J

FIGURE 32. - FIRST-OEDER UPWIND (--HYBRID- PLDS-'- EDS) RESULTS

FOR THE ROTATING PROBLER_ _ = G5.3.

46
FIGURE 33. - SECOND-ORDER-UPICIND RESULTS FOR THE ROTATING STEP FIGURE 3q. - UNLIMITED gUICKo2D RESULTS FOR THE ROTATING STEP

PROBLE.R: g = 27.3. PROBLEM- _ = 18.9.

FIGURE 35. - UNLIRITED FIFTH-ORDER UI:_IND RESULTS FOR THE ROTATING FIGURE 36. - RINROD RESULTS FOR THE ROTATING STEP PROBLF_: _ = 30.8.

STEP PROBLER: # = 16.1.

47
FIGURE 37, - ULTRA-B RESULTSFOR THE ROTATINGSTEP PROBLF_.M-

_= 7.7.

FIGURE 38. - ULTRA-3RD/5TH/7TH RESULTSFOR THE ROTATING STEP

PROBLF._ 8 = 9.1,

48
Report Documentation Page Nalional Aeronautics and

Space Administration

2. Government Accession No. 3. Recipient's Catalog No.

1. Report No. NASA TM-102568

ICOM P-90-12

4. Title and Subtitle

ULTRA-SHARP Nonoscillatory Convection Schemes for High-Speed

Steady Multidimensional Flow

7. Author(s)

B.P. Leonard and Simin Mokhtari

9. Performing Organization Name and Address

National Aeronautics and Space Administration

Lewis Research Center

Cleveland, Ohio 44135-3191

12. Sponsoring Agency Name and Address

National Aeronautics and Space Administration

Washington, D.C. 20546-0001

5. Report Date

April 1990

6. Performing Organization Code

8. Performing Organization Report No.

E-5320

10. Work Unit No,

505-62-21

11. Contract or Grant No,

13. Type of Report and Period Covered

Technical Memorandum

14. Sponsoring Agency Code

15. Supplementary Notes

B.P. Leonard, Institute for Computational Mechanics in Propulsion, Lewis Research Center (work funded by

Space Act Agreement C-99066-G). Space Act Monitor: Louis A. Povinelli, Simin Mokhtari, Dept. of Mechanical

Engineering, The University of Akron, Akron, Ohio 44325.

16. Abstract

For convection-dominated flows, classical second-order methods are notoriously oscillatory and often unstable.

For this reason, many computational-fluid-dynamicists have adopted various forms of (inherently stable) firstorder upwinding over the past few decades. Although it is now well known that first-order convection schemes

suffer from serious inaccuracies attributable to artificial viscosity or numerical diffusion under high-convection

conditions, these methods continue to enjoy widespread popularity for numerical heat-transfer calculations,

apparently due to a perceived lack of viable high-accuracy alternatives. But alternatives are available. For

example, nonoscillatory methods used in gasdynamics, including currently popular "TVD" schemes, can be

easily adapted to multidimensional incompressible flow and convective transport. This, in itself, would be a major

advance for numerical convective heat transfer, for example. But, as this-report shows, second-order TVD

schemes form only a small, overly restrictive, subclass of a much more universal, and extremely simple,

nonoscillatory flux-limiting strategy which can be applied to convection schemes of arbitrarily high-order

accuracy, while requiring only a simple tridiagonal ADI line-solver, as used in the majority of general-purpose

iterative codes for incompressible flow and numerical heat transfer. The new universal limiter and associated

solution procedures form the so-called ULTRA-SHARP alternative for high-resolution nonoscillatory

multidimensional steady-state high-speed convective modelling.

17, Key Words (Suggested by Author(s))

Upwind schemes; Higher-order upwinding; Nonoscillatory

convection; ULTRA-SHARP method

18. Distribution Statement

Unclassified- Unlimited

Subject Category 64

19. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No. of pages

Unclassified Unclassified 50

NASA FORM 1626 OCT 86 *For sale by the National Technical Information Service, Springfield, Virginia 22161

22. Price*

A03


National Aeronautics and

Space Administration

Lewis Research Center

ICOMP (M.S. 5-3)

Cleveland, Ohio 44135

Official Bu=dnmm

Penalty for Private Use $300

FOURTH CLASS MAIL

ADDRESS CORRECTION REQUESTED

!!1III

Poslage and Fees Pa_d

Nal_onal Ae,onaubcs and

Space Adm_n_stralton

NASA 45t