NASA Technical Memorandum 102538

ICOMP-90-09

..... -_ 7i_L2. Z.? 17 : . . 2 __ _T..... 2 -. 2 Z _ . . L.L . _ -2___±_.2L_7__ IL.._52 i_ _? ........ ZL.- .... : . . _i .... £_5.5__&_

A Cost-Effective Strategy for

Nonoscillatory Convection :

Without Clipping

B]P. Leonard

Institute for Computational Mechanics in Propulsion

Lewis Research Center

Cleveland, Ohio:

and

H.S. Niknafs

The Norton Company

Stow, Ohio-

.............. March 1990

.--v_

A COST-EFFECTIVE STRATEGY FOR NONOSCII,LATORY

CONVECTION WITHOUT CLIPPING

B.P. Leonard*

Institute for Computational Mechanics in Propulsion

Lewis Research Center

Cleveland, Ohio 44135

and

H.S. Niknafs**

The Norton Company

Chemical Products Division

Stow, Ohio 44224

_o

¢o

I

SUMMARY

Clipping of narrow extrema and distortion of smooth profiles is a well-known

problem associated with so-called "high-resolution" nonoscillatory convection

schemes. In this report, a strategy is presented for accurately simulating highly

convective flows containing discontinuities such as density fronts or shock waves,

without distorting smooth profiles or clipping narrow local extrema. The convection

algorithm is based on non-artificially-diffusive third-order upwinding in smooth

regions, with automatic adaptive stencil expansion to (in principle, arbitrarily)

higher order upwinding locally, in regions of rapidly changing gradients. This is

highly cost-effective because the wider stencil is used only where needed - in isolated

narrow regions. A recently developed universal limiter assures sharp monotonic

resolution of discontinuities without introducing artificial diffusion or numerical

compression. An adaptive discriminator is constructed to distinguish between

spurious overshoots and physical peaks; this automatically relaxes the limiter near

local turning points, thereby avoiding loss of resolution in narrow extrema.

Examples are given for one-dimensional pure convection of scalar profiles at

constant velocity.

*Work funded under Space Act Agreement C99066G; presently Professor of Mechanical Engineering,

The University of Akron, Akron, Ohio 44325.

**Doctoral student, College of Engineering, The University of Akron, Akron, Ohio 44325.
INTRODUCTION

Accurate simulation of highly convective flows continues to be one of the most

challenging problems in computational mechanics. The inability to adequately

simulate simple scalar profiles even in the superficially straight-forward case of onedimensional pure convection at constant velocity has been referred to as the

"ultimate embarrassment" for computational mechanics [1]. Classical (second-order

central) schemes are able to handle very smooth profiles fairly satisfactorily; i.e.,

when the convected variable has small curvature (spatial second derivatives) - or,

more specifically, when the shortest Fourier components have wavelengths of many

mesh-widths. But a sudden change in gradient (l"nvolving locally high curvature or

short wavelengths) excites spurious unphysical oscillations due to inherent spatial

(third-derivative) dispersion terms in the truncation error. By contrast, first-order

upwinding [2] can resolve a step profile monotonically, giving a spreading errorfunction in the idealized case; but the inherent spatial second-derivative (artificial

diffusion) terms in the truncation error so totally overwhelm any physical diffusion

that even smoothly varying profiles are unrecognizable after a short time.

Second-order upwinding [3] represents a significant improvement. In this case,

leading order truncation error is again a dispersive third derivative but, compared

with second-order central, fourth- (and higher even-) derivative dissipation terms

may be stronger under certain conditions, thus tending to somewhat dampen dispersive oscillations. By averaging second-order upwinding with the canonical explicit

second-order central scheme [4], Fromm [3] proposed a convection algorithm in _

which the leading phase error of the upwind scheme approximately cancels the

lagging phase error of the central scheme. However, using the same computational

stencil, it is a simple matter to eliminate the third-derivative dispersion term

entirely.

The QUICKEST scheme [5] is the canonical explicit third-order upwind convection algorithm of this type; physical diffusion terms are modelled to a consistent

order. In this case, the leading truncation-error term is a dissipative (but not diffusive) spatial fourth derivative. This is an important attribute. Numerical experience over the past decade or more has led to the following rule of thumb: leading

truncation error in a convection algorithm should be dissipative (an even spatial

derivative) rather than dispersive (odd derivative) but shouldbe of higher order than

other modelled physical terms such as diffusion (especially if these terms are supposed to be small). Thus, even-order schemes fail the first criterion (central schemes

are worse in this respect because dissipative terms are small or non-existent). Firstorder upwinding fails the second criterion. Third-order upwindingis therefore the

lowest order method satisfying both criteria for the convection-diffusion equation.

This is especially important for convection-dominated flows, as will be seen in the

model test problems shown later.

Unfortunately, third-order upwinding does not seem to have been widely

adopted by the computational-fluid-dynamics (CFD) community. Apparently for no

other reasons than tradition and personal inertia, many CFD researchers regard

second-order central schemes as the norm. This is appropriate for physical problems

dominated by even spatial derivatives, such as diffusion, wave-motion, and elasticity, but is singularly inappropriate for the first-derivative convection term in fluid

mechanics [6]. A similar situation occurs with higher order odd derivatives such as

the spatial third derivative in the Korteweg-deVries equation [7]. Another large

segment of the CFD community (especially in heat-transfer applications) has

adopted first-order upwinding as the "safest" form of convection modelling. This
-V

apparently stems from the influence of the Imperial College school through descendants of the TEACH code [8] using the "Hybrid" convection scheme [9] or variants of

the so-called "exponential difference scheme" [10], including the "power-law"

approximation described in a well-known textbook [11]. The difficulty with such an

approach is very simple: using these methods, one is not simulating the correct

nhvsical nroblemV A very-high-Reynolds-number (convection-dominated) problem is

arb'itrari_y replaced by an un_physical problem in which the effective component grid

Reynolds number can never be greater than 2, whereas the physical grid Reynolds

number might be hundreds or thousands (or effectively infinite).

There is also a serious logical flaw in such methods: sophisticated (and expensive) multi-equation turbulence models are usually used at each time or iteration

step simply as a diagnostic to switch off their own effects in the turbulent-transport

terms of the governing equations, replacing these terms with artificial viscosityor

diffusivity [12]. Specifically, the turbulence model is used to calculate the t_uroulent) component grid Reynolds (or Peclet) number at each control-volume face; when

this quantity exceeds 2 (in the case of Hybrid) or about 6 (for the exponential or

power-law schemes), physical turbulent transport is replaced by artificial numerical

viscosity (or diffusivity). The short-comings of first-order methods are well known,

having been shown theoretically [13-15] and through comparative benchmark

problems [16-17]. Their continued wide-spread popularity is perhaps another

manifestation of the ultimate embarrassment.

Although third-order upwinding is a more rational basis for CFD, certain

fundamental problems remain. Perhaps the most obvious is the fact that unphysical

overshoots and undershoots are excited each side of a step discontinuity in purely

convective flow [5]. And, of course, short-wavelength resolution is limited by the

third-order accuracy. Switching to higher order (upwind) schemes can increase

accuracy of resolution but cannot reduce the overshoot problem (in fact, it gets

slightly worse). Monotonic high resolution schemes can be constructed for simulating step discontinuities. This is usually achieved by using a subtle nonlinear

blending of a (traditional) second-order central base scheme with first-order

upwinding (introducing enough positive artificial diffusion locally to maintain

monotonicity) and first-order downwinding (with inherent negative artificial diffusion to locally enhance artificial "compression"). This strategy is the basis of socalled "shock-capturing" and "total-variation-diminishing (TVD)" schemes [18], and

can be related to earlier multi-step flux-correction first/second-order schemes [19]

and single-step nonlinear flux-limiting schemes such as the second-order MUSCL

scheme [20] - a monotonized version of Fromm's method - and the third-order

EULER scheme based on exponential upwinding [21,22]. Unfortunately, the best

step-resolution (supercompressive) schemes such as Superbee [23], Super-C, or

Hyper-C [24], tend to convert smooth profiles into a series of steps and plateaus. In

other words, the negative diffusivity responsible for artificial compression of discontinuities tends to concentrate moderate-curvature regions into localized unphysical

sharp changes in gradient. All such first/second-order schemes strongly "clip"

narrow extrema.

Recently, a universal limiter (UL) has been designed which can be applied to

arbitrarily high order transient interpolation modelling (TIM) of the advective

transport equations (ATE). This "ULTIMATE" CFD scheme [24] has several attractive properties Step resolution is monotonic and can be made competitive with the

A L " _ . . • •'L • • "

best supercompresslve shock-capturing schemes without introducing artlficml dlffusivity or viscosity, simply by locally increasing the order of accuracy, using adaptive stencil expansion (so that the location of the wider stencil automatically moves
along with high-curvature regions of the profile). Since numerical compression is

avoided, smooth profiles are not corrupted. However, one annoying problem

remains: very narrow extrema are slightly clipped relative to what can be achieved

with an unlimited higher order (upwind) scheme. Near sudden changes in.gradient,

the limiter needs to be "on", but near local physical extrema (but not spurious overshoots), it needs to be switched "off" automatically. The problem is one of pattern

recognition. This paper describes a simple adaptive discriminator which can identify

well-defined local narrow .physical extrema and automatically switch off the

universal limiter in such regnons, but keep it activated in regions where unphysical

overshoots or very short-wavelength spurious oscillations would otherwise occur.

(i)

(ii)

(iii)

The overall strategy is as follows:

In "smooth" regions (identified by small values of local absolute'curvature of

the convected variable), the unlimited third-order upwind QUICKEST scheme

is used; this accounts for the overwhelming bulk of the flow domain - especially

in multidimensional flows.

In relatively large-gradient or strong-curvature regions, automatic adaptive

stencil expansion occurs locally, using (in principle, arbitrarily) higher order

upwinding with the universal limiter activated.

Near well-defined local extrema, identified by the automatic discriminator, the

limiter is switched off, thus allowing an appropriate degree of resolution

depending on the narrowness of the extremum.

The next section summarizes the basic ideas underlying the universal limiter.

Then four challenging scalar test profiles are considered under conditions of pure

one-dimensional convection at constant velocity. Results for a number of solution

methods are shown in order to demonstrate some of the difficulties mentioned above

with respect to various manifestations of the ultimate embarrassment. Construction

of the automatic discriminator is then briefly described. Finally, results are shown

for a cost-effective variable-order scheme consisting of QUICKEST in smooth

regions, ULTIMATE seventh-order upwinding in moderate-curvature re_ions and

ULTIMATE ninth-order upwinding in high-curvature or high-gradient regmns.

UNIVERSAL LIMITER

Details of the ULTIMATE scheme can be found elsewhere [24]; however, a brief

description of the universal limiter is given here for convenience. Figure 1 shows a

one-dimensional control volume, suggesting the local behaviour of the convected

variable (p in terms of locally normalized variables [22]

= (¢ - ¢v)/(¢o" ¢v) (1)

where (l), is the (unnormalized) upstream node value and qbl) the downstream

value. NSte from Equation (1) and the figure that, in terms of normalized variables,

_v =0 and 7o =1 (2)

Let _r be the normalized face value on the downstream control-volume face;

simila/-ly, _, isthe normalized upstream face value. Figure 1 shows interpolative

constraints on _f and qb,_; i.e.,
and

(Pc = = 1 (3)

o <_-- _hi <_. _ C (4)

where _c" is the normalized central node value. The inequalities given by (3) are

two necessary conditions on _f prescribed by the universal limiter when local

behaviour is monotonic across the three nodes shown in Figure 1, i.e., the two nodes

straddling the face in question and the adjacent upstream node. However, this is not

enough to guarantee computational monotonicity. For this, consider the explicit

update step for _c under conditions of constant velocity (positive to the right in

Figure 1). In terms of normalized variables, ....

(5)

where c - uAt/Ax is the Courant number.

In order to maintain monotonicity (locally), 7o must satisfy

_'uTM -<-",'c"+' _ --<7. "+' (6)

Taking conservative worst-case conditions implies

0 < _.+l < 1 7)

---- .it. C

The right-hand inequality is assured by the interpolative constraints on _f and 7. ;

the left-hand inequality implies, using Equation (5),

_ =<_hi+ _"/c (8)

Once again, a worst-case estimate for _hi (= 0) results in an additional (time-step)

constraint on _'f:

'_f <- "_c" Ic (9)

Inequalities (3) and 9_ constitute the universal limiter constraints on _f with respect to _c n when v" is within the monotonic range:

0<'_= n< (Pc = 1 (10)

Outside of this range (i.e., _n" < 0 or ~ " (Pc > 1), various strategies can be devised.

Figure 2 shows one simple po§sibility:

_f = _c" (nonmonotonic range) (11)

together with the monotonicity constraints. This is a graphical portrayal of the

universal limiter. Using this limiter, the procedure is as follows:

(i) Construct an explicit estimate for the unnormalized control-volume face value,

(Pf, by any (in principle, arbitrarily high order) method.

(ii) Compute the corresponding normalized value, _f, according to Equation (1).

5
(iii)

(iv)

(v)

(vi)

If the point (_'c", _'r) lieswithin the shaded region of Figure 2, proceed with the

unadjusted face value, ¢bf.

If the point lies outside this region, replace _f with the nearest allowable _f

value at the same _ n C"

Reconstruct Cr = _u + _'f (¢D" CV)"

Repeat for each control-volume face.

Note that this strategy is significantly different from that used in currently

popular shock-capturing and TVD schemes [18]. In such schemes, using the present

normalized-variable terminology, ..

where _'7 is a single-valued monotonic function of _" passing through the points

(0,0), (0._, 0.75), and (1,1). This restricts such metho_ls to second (or, at best) third

order accuracy in space, and second order in time. Figure 3 shows the functional

form of _f_(_,') for three well-known methods: Roe's Minmod and Superbee [23],

and van Leer_s MUSCL [20]. As commonly used in shock-capturing schemes, the

TVD limiter [25] can be described by

_c "<: _f_ =-< min(2_c" ,1) (13)

• II

in the monotonic range, with _/ = _c n elsewhere. These constraints, together with

Equation (12), represent muck more restrictive conditions than those of the universal limiter [24]. An important point to stress is that the universal limiter can be

used with arbitrarily high order (in both space and time) explicit methods whereas

the usual TVD strategy is restricted to essentially second order.

CONVECTION OF SCALAR PROFILES

Consider one-dimensional constant-velocity pure convection (i.e., no diffusion)

of a scalar profile initially given by

_(x,0) = _o(X) (14)

After time t, the exact solution is the initial profile translated by a distance

where u is the convecting velocity; i.e.,

ut,

¢(x,t) = ¢o(x-ut) (15)

Figure 4 shows the exact solution of four test profiles under uniform convection to

the right. Discontinuity resolution is tested by a simple step rather than the square

wave used by some authors [25] (at best, a square wave simply gives twice as much

information as a step; for the more dispersive schemes, oscillations produced by the

square-wave's step-up can interfere with those produced by the step-down, thus

producing a confusing interaction pattern). Smooth behaviour is represented by an

isolated sine-squaredwave 20Ax wide; there is no discontinuity in value or

gradient, but there is a discontinuity in curvature at each side of the profile. A semiellipse 20Ax wide is very challenging because of a combination of sudden and

gradual changes in gradient. A narrow Gaussian with o = 1.94Ax has been chosen

to test narrow peak resolution.

6
In the following simulations, the Courant number for all cases is

c : 0.45 (16)

Calculations are run for 100 time-steps, so that exact solutions should have shifted

45 mesh-widths to the right. The exact solutions are shown for reference in each

case. All of the schemes have been programmed in conservative control-volume

form, so that "mass" is conserved (to machine accuracy) in all cases. Upstream

(_ = 1) and downstream (¢ = 0) boundaries are sufficiently far away (beyond the

region shown in the figures) that they do not interfere with the solutions. To get

some idea of how far the profiles have moved, note that the three maxima are each

30Ax apart. Thus, for example, the initial position of the Gaussian peak was halfway between the final positions of the sine-squared and semi-ellipse maxima.

Figure 5 shows the results of using first-order upwinding. Because of strong

artificial diffusion (spatial second derivative) in the truncation error, the computed

profiles have strongly interacted. Individually, the step simulation would look like a

spreading error-function, whereas the other profiles soon degenerate into spreading

(and amplitude-decaying) Gaussians, since all short-wavelength components in the

original profiles are quickly damped out.

If one defines the grid Peclet number as

Pa = uAx/D (17)

where D is the physical diffusion coefficient, it is clear that the physical problem

under consideration corresponds to infinite PA" AS is well known [26], first-order

upwinding for infinite P. is indistinguishabl_ from any higher order method with

an effective grid Peclet nu_mber given by

Pa* = 2/(1-c) (18)

In other words, for transient problems, first-order upwinding introduces artificial

numerical diffusion of the form

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

11l.llrl

Such artificially diffusive results should be recognized as being totally unacceptable.

It is indeed rather surprising that such inaccurate methods are still the industry

standard in many branches of computational engineering, especially numerical

convective heat transfer [27,28].

Figure 6 shows results for the explicit second-order central Lax-Wendroff [4]

scheme. In this case, there is a large phase lag due to the spatial third-derivative

term in the truncation error. Only the smooth (sine-squared) profile bears any

resemblance to the corresponding exact solution. Clearly, for highly convective

flows, the explicit Lax-Wendroff method is very disappointing for all but the

smoothest of profiles. The situation is not improved by using explicit second-order

upwinding, shown in Figure 7. In this case, the third-derivative truncation-error

terms are of opposite sign to those of the Lax-Wendroff method, resulting in phaselead dispersion. For Courant numbers near 0.5, the individual second-order upwind

results are almost a reflection of the Lax-Wendroffprofiles.
Implicit second-order convection schemes are not necessarily better than

explicit methods. The so-called "Crank-Nicolson-type" convection scheme [26] -

defined by analogy with the well-known method for the diffusion equation [29] - is

actually worse than explicit methods in the purely convective case, because there are

strong phase-lag dispersion terms and no (even-derivative) dissipation terms in the

truncation error. The highly oscillatory nature of this method is seen in Figure 8.

The second-order implicit linear-finite-element method [26], shown in Figure 9,

again has no dissipation but gives somewhat better results than other second-order

methods because the leading third-derivative dispersion term is smaller.

Clearly, it would seem somewhat more felicitous to design a high-convection

scheme in which the oscillatory third-derivative dispersion term in the truncation

error were entirely eliminated. This was achieved over ten years ago in a simple

explicit formulation known as the QUICKEST scheme [5]. In this case, leading truncation error is a small (fourth-derivative) dissipation term which strongly damps any

dispersive tendencies of the non-zero fifth-derivative term without corrupting the

modelled physics, viz.: the absence of physical (second-derivative) diffusion terms.

Figure 10 shows a dramatic improvement over the dispersive second-order methods.

Note particularly the good phase behaviour (which is relatively insensitive to

Courant number). For reference, Figure 11 shows results for an explicit fourth-order

central method [24]. Phase-lagdispersion in this case stems from the leading fifthderivative truncation error (which is only lightly damped by higher order dissipation). Results are only slightly better than those of the best second-order method.

At the risk of belabouring the point, the above results seem to indicate fairly

clearly that third-order upwinding represents a much more rational basis for CFD

[6] than first- or second-order methods or higher order central schemes. As

mentioned previously, there are still two serious short-comings, both evident in

Figure 10: (i) unphysical undershoots or overshoots near regions of rapid change in

gradient (high curvature); (ii) a limit to short-wavelength resolution. The first of

these deficiencies can be overcome by applying the universal limiter; Figure 12

shows the resulting ULTIMATE QUICKEST computation. Step resolution is

monotonic and comparable to that of MUSCL, Figure 13, a limited form of Fromm's

method. (One should note that at c = 0.5, Fromm's method is equivalent to

QUICKEST [24] ). But it is clear that MUSCL tends to clip and flatten local extrema

more strongly. This is directly linked to the more restrictive conditions imposed by

the commonly used "TVD" limiter, given by (13).

The second deficiency of third-order upwinding can, of course, be corrected by

using higher order (upwind) methods. The strategy proposed here is to use an

appropriate higher order method locally, with third-order upwinding in relatively

"smooth" (low-curvature) regions. For reference, Figure 14 shows results for

(unlimited) ninth-order upwinding [24] used globally; ninth-order appears to be

necessary to fully capture the peak of the narrow Gaussian. But, of course, this is an

unlimited scheme, and undershoots and overshoots are excited near discontinuities.

ULTIMATE ninth-order, Figure 15, eliminates the oscillations and gives very tight

step-resolution, but again introduces slight clipping at extrema. If one could

somehow arrange to automatically switch off the limiter near local extrema (while

keeping it active near discontinuities), a very desirable convection scheme would

result. This is the job of an adaptive discriminator which can distinguish between

well-defined local physical extrema and short-wavelength spurious oscillati0ns.

Note that long-wavelength numerical oscillations, such as those occurring in secondorder simulations - which might otherwise "confuse" the discriminator - present no

problem (because they are excluded from the basic scheme).

8
CONSTRUCTION OF THE DISCRIMINATOR

The following algorithm attempts to distinguish between artificial numerical

peaks, such as those potentially occurring near the step or the base of the semiellipse, and true physical extrema, such as those of the exact profiles in the test cases

studied. Artificial peaks would be associated with short-wavelength numerical

oscillations, with rapidly changing value, gradient, and curvature; as mentioned

above, long-wavelength numemcal oscillations are excluded from the third/higher

order algorithm. If a local extremum is associated with short-wavelength oscillations, the discriminator chooses the limited algorithm; however, if the curvature is of

the same sign in adjacent regions, the discriminator relaxes the limiter constraints

at the extremum and the immediate upstream and downstream adjacent nodes.

In setting up the convective fluxes (for the right-face of each CV cell), consider

an increasing-i "DO-loop" sweep. For the discriminator, choose a stencil of seven

points: .(1)_.3- (1),.2, (Pi,-1, _b_,._i+l, (1)i+2, ¢i+3" Then compute the differences between

eacn pair o_consecuuve points:

D1 = (dpi.2,-¢i.3), D2 = ((_i.i - d_,.2),..., D6 = (dpi+3- _i+2) (20)

For convenience, assume there is a local maximum; a minimum requires reversal of

some of the subsequent inequalities. Limiter constraints are active unless otherwise

stated. The algorithm proceeds as follows:.

(i) Check if D1 and D2 are both positive and D3 and D4 both negative; if not,

go to step (iii); if they are, then:

; if true, switch to an unlimited version at (ii) Check if 1921</911 and 1931<

the current /-value and skip the remaining steps; otherwise:

(iii) Check if D2 and D3 are both positive and D4 and D5 bothnegative; if not,

go to step (v); if they are, then:

(iv) Check if ID3I </D21 and ID41< ID51; if true, switch to an unlimited version at

the current /-value and skip the remaining steps, otherwise:

(v) Check if D3 and D4 are both positive and D5 and D6 both negative; if not,

keep limiter constraints active; if they are, then:

(vi) Check if ID4[ < ID31 and 1951 < ID6[ ; if true, switch to an unlimited version at

the current /-value; if not, proceed with the limiter constraints active.

Figure 16(a) shows a sketch of a case in which the discriminator would keep the

universal limiter activated; in Figure 16(b), the limiter constraints would be

removed at each of thepoints shownby a hollow circle. The discriminator routine is

by-passed unless one of(p ÷., (1)i, or (l)i. is a local extremum. As noted above, limiter

relaxation occurs (if at all) _n groups _" three points, with the discrete extremum in

the middle.

AI)APTIVE STENCIL EXPANSION

The ultimate convection scheme proposed here uses the unlimited QUICKEST

(third-order upwind) algorithm in smooth non-steep regions; i.e., wherever the

average (right-face) absolute "curvature'

CURVAV = 0.5 1(1),+2" ¢,+, - qb, + (21)
and "gradient"

GRAD = I_b,+_ - _l (22)

are both less than pre-assigned thresholds, and provided [24]

0.2 _-< _c < 0.8 (23)

If CURVAV exceeds THC1 or the larger THC2 , the algorithm automatically

switches locally (at the CV face in question) to ULTIMATE seventh- or ninth-order

Winding, respectively. If GRAD exceeds THG it locally switches directly to

TIMATE ninth-order upwinding. This order-switching strategy is sketched in

Figure 17. Also, if the limits expressed in (23) are locally exceeded, the limiter

constraints are activated for the basic third-order scheme. Near local extrema, the

discriminator automatically relaxes the limiter constraints, as described in the

previous section. Clearly, the (dimensional) threshold constants need to be chosen

carefully so as to capture the desired degree of accuracy without "overkill" in terms

of cost; this requires some experimentation for optimization; a more universal

(nondimensional) procedure is currently being explored.

Adaptive stencil expansion is a very cost-effective strategy because the widerstencil computations are performed only where needed - at a relatively small

number of grid points in narrow regions, by definition. This is even more effective in

two and three dimensions. Figure 18(a) shows results for the complete algorithm;

the small arrows show where the automatic discriminator has relaxed the limiter

constraints (at this particular time-step). Figure 18(b) shows the corresponding

distribution of GRAD, and 18(c) that of CURVAV. Finally, Figure 18(d) shows the

local order of the algorithm to be used (at the next time-step) for the threshold

constants chosen. It should be clear that the locations of the limiter-relaxation

points and the order-switches move along with the profiles automatically, as time

progresses.

CONCLUSION

A cost-effective strategy for high-convection modelling has been introduced.

The algorithm is based on the (third-order upwind) QUICKEST scheme in smooth

regions, as this is the lowest order method in which the leading truncation error is

dissipative but not diffusive. QUICKEST s phase error is relatively low, but

overshoots or undershoots can be excited by the unlimited scheme near discontinuities, and short-wavelength resolution would, of course, be limited to third order.

A more sophisticated scheme is therefore devised in which monotonicity is

Uparanteed by a universal limiter which can be applied to any order of accuracy.

propriate degrees of higher order resolution are automatically introduced locally,

using adaptive stencil expansion controlled by monitoring, local (absolute) gradient

and curvature of the convected variable. Potential clipping of narrow extrema is

avoided by automatic relaxation of the limiter constraints in local regions based on

pattern-recognition decisions of an adaptive discriminator. By judicious choice of

threshold constants, the overall method can produce extremely accurate results on a

coarse mesh at relatively little additional cost above that of the base third-order

scheme. The same principles can be applied to multidimensional flows and nonlinear systems. Work is proceeding on these applications.

10
lo

.

.

4.

5.

.

.

.

.

10.

13.

14.

15.

REFERENCES

Mitchell, A.R.: Recent Developments in the Finite Element Method. Computational Techniques and Applications: CTAC-83, J. Noye and C.A.J. Fletcher,

eds., Elsevier North-Holland, 1984, pp. 2-14.

Courant, R.; Isaacson, E.; and Rees, M.: On the Solution of Nonlinear

/_yaPerbolic Differential Equations by Finite Differences. Commun. Pure Appl.

th., Vol. 5, 1952, pp. 243-255.

Fromm, J.E.: A Method for Reducing Dispersion in Convective Difference

Schemes. J. Comput. Phys., Vol. 3, 1968, pp. 176-189.

Lax, P.D.; and Wendroff, B.: Systems of Conservation Laws. Commun. Pure

Appl. Math., Vol. 13, 1960, pp. 217-237.

Leonard, B.P.: A Stable and Accurate Convective Modelling Procedure Based

on Quadratic Upstream Interpolation. Comput. Methods Appl. Math. Eng.,

Vol. 19, No. 1,June 1979, pp. 59-98.

Leonard, B.P.: Third-Order Upwinding as a Rational Basis for Computational

Fluid Dynamics. Computational Techniques and Application: CTAC-83, J.

Noye and C.A.J. Fletcher, eds., Elsevier North-Holland, 1984, pp. 106-120.

Anderson, D.A.; Tannehill, J.C.; and Pletcher, R.H.: Computational Fluid

Mechanics and Heat Transfer, Hemisphere Publishing Corp., Washington,

D.C., 1984.

Gosman, A.D.; and Pun, W.M." Lecture Notes for a CourseEntitled:

Calculation of Recirculating Flows. Heat Transfer Section Rep'_rt I-iTS/74/2,

Imperial College Mechanical Engineering Dept., 1974.

Spalding, D.B.: A Novel Finite Difference Formulation for Differential Expressions Involving Both First and Second Derivatives. Int. J. Numer. Methods

Eng., Vol. 4, No. 4, July-Aug. 1972, pp. 551-560.

Raithby, G.D.; and Schneider, G.E.: Elliptic Systems: Finite-Difference

Method H. Handbook of Numerical Heat Transfer, W.J. Minkowycz, et al., eds.,

Wiley, 1988, pp. 241-291.

Patankar, S.V.: Numerical Heat Transfer and Fluid Flow, Hemisphere

Publishing Corp., Washington, D.C., 1980.

Syed, S.A.; Gosman, A.D.; and Peric, M.: Assessment of Discretization

Schemes to Reduce Numerical Diffusion in the Calculation of Complex Flows.

AIAA Paper 85-0441, Jan. 1985.

de Vahl Davis, G.; and Mallinson, G.D.: An Evaluatin of Upwind and Central

Difference Approximations by a Study of Recirculating Flow. Comput. Fluids,

Vol. 4, No. 1, Mar. 1976, pp. 29-43.

Raithby, G.D.: A Critical Evaluation of Upstream Differencing Applied to

Problems Involving Fluid Flow. Comput. Methods Appl. Mech. Eng., Vol. 9,

No. 1, Sept. 1976, pp. 75-103.

Leonard, B.P.: A Consistency Check for Estimating Truncation Error Due to

Upstream Differencing. Appl. Math. Modelling, Vol. 2, No. 4, Dec., 1978, pp.

239-244.

11
26.

27.

28.

29.

16. Huang, P.G.; Launder, B.E.; and Leschziner, M.A.: Discretization of NonLinear Convection Processes: A Broad-Range Comparison of Four Schemes.

Comput. Methods Appl. Mech. Eng., Vol. 48, No. 1, Feb. 1985, pp. 1-24.

17. Leschziner, M.A.: Practical Evaluation of Three Finite Difference Schemes for

the Computation of Steady-State Recirculating Flows. Comput. Methods Appl.

Mech. Eng., Vol. 23, No. 3, Sept. 1980, pp. 293-312.

18. Yee, H.C.: Upwind and Symmetric Shock-Capturing Schemes. NASA TM89464, 1987.

19. Boris, J.P.; and Book, D.L.: Flux-Corrected Transport HI. Minimal-Error FCT

Algorithms. J. Comput. Phys., Vol. 20, No. 4, Apr. 1976, pp. 397-431.

20. van Leer, B.: Towards the Ultimate Conservative Difference Scheme. V. A

Second-Order Sequel to Godunov's Method. J. Comput. Phys., Vol. 32, 1979, pp.

101-136.

21. Leonard, B.P.: Adjusted Quadratic Upstream Algorithms for Transient Incompressible Convection. AIAA 4th Computational Fluid Dynamics Conference,

AIAA, 1979, pp. 226-233.

22. Leonard, B.P.: 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, 1981, pp. 1-36.

23. Roe, P.L.: Some Contributions to the Modelling of Discontinuous Flows. Large

Scale Computations in Fluid Mechanics, Part 2 (Lectures in'Applied Mathematics, Vol. 22), B.E. Engquist, S. Osher, and R.C.J. Somerville, eds., American

Mathematical Society, Providence, RI, 1985, pp. 163-193.

24. Leonard, B.P.: Universal Limiter for Transient interpolation Modeling of the

Advective Transport Equations: The ULTIMATE Conservative Difference

Scheme. NASA TM 100916(ICOMP-88-11), 1988.

25. Sweby, P.K.: High Resolution Schemes Using Flux Limiters for Hyperbolic

Conservation Laws. SIAMJ. Numer. Anal., Vol. 21, No. 5, Oct. 1984, pp. 495-

1011.

Fletcher, C.A.J: Computational Techniques for Fluid Dynamics, SpringerVerlag, 1988.

Minkowycz, W.J.; et al.: eds., Handbook of Numerical Heat Transfer, Wiley,

1988.

Lewis, R.W.; and Morgan,K.: eds., Numerical Methods in Thermal Problems

VI, Pineridge Press, 1989.

Crank, J.; and Nicolson, P.: A Practical Method for Numerical Evaluation of

Solutions of Partial Differential Equations of the Heat-Conduction Type. Proc.

Cambridge Philos. Soc., Vol. 43, Pt. 1, Jan. 1947, pp. 50-67.

12
zl_//./

l

f //_//

I

I

I

liv v

'I

CONTROL

•-,-,I_ut

VOLUI'_ I

.I

FIGURE 1 Locally monotonic behaviour across a control-volume cell in terms

of normalized variables.

/

FIGURE 2 Universal limiter constraints in the normalized-variable diagram.

13
/

$UPERBEE

MUSCL

MINMOD

!

FIGURE 3 Normalized variable diagrams for Superbee (heavy lines), Minmod

(dashed), and MUSCL.

FIGURE 4 Exact solution for the four scalar test profiles considered.

14
FIGURE5 First-order upwind results compared with exact solution.

FIGURE6 Lax-Wendroff results compared with exact solution.

15
FIGURE7 Second-order upwind results compared with exact solution.

FIGURE8 Crank-Nicolson results compared with exact solution.

16
FIGURE9 Linear-finite-element results compared with exact solution.

FIGURE10 Third-order upwind (QUICKEST)results compared with exact

solution.

17
FIGURE11 Fourth-order central results compared with exact solution.

FIGURE12 ULTIMATEQUICKESTresults compared with exact solution.

18
FIGURE 13 Results using MUSCL compared with exact solution.

...... _______ -l_ _

FIGURE 14 Unlimited ninth-order upwinding.

19
i||

FIGURE 15 ULTIMATE ninth-order upwinding.

FIGURE 16

(a) (b)

Action of the discriminator. (a) Limiter switched on.

(b) Limiter switched off at three central points.

20
FIGURE 17

CURVAV

THC2

THC1

Schematic of order-switching strategy in terms of CURVAV

and GRAD.

21
(a)

(b)

(c)

Cd)

I _ 9th

t f th

.......... 1 1...... _rO.........

FIGURE 18 Third/higher-order adaptive scheme with discriminator.

(a) Computed results. (b) Distribution of GRAD relative toTHG.

(c) Distribution of CURVAV relative to THC1 and THC2.

(d) Order of the algorithm to be used locally in next time-step.

22
NatbonalAeronaullcs and

Space Administration

1. Report NO. NASA TM- 102538

ICOMP-90-09

Report Documentation Page

2. Government Accession No.

4. Title and Subtitle

A Cost-Effective Strategy for Nonoscillatory

Convection Without Clipping

7. Author(s)

B.P. Leonard and H.S. Niknafs

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

3. Recipient's Catalog No.

5. Report Date

March 1990

6. Performing Organization Code

8. Performing Organization Report No.

E-5336

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 C99066G); presently Professor of Mechanical Engineering, ,The University of Akron, Akron, Ohio 44325.

H.S. Niknafs, The Norton Company, Chemical Products Division, Stow, Ohio 44224; doctoral student, College of

Engineering, The University of Akron, Akron, Ohio 44325. Space Act Monitor: Louis A. Povinelli

16. Abstract

Clipping of narrow extrema and distortion of smooth profiles is a well-'_nown problem associated with so-called

'_high-_'esolution 't nonoscillatory convection schemes. In this report, _ strategy is presented for accurately

simulating highly convective flows containing discontinuities such as density fronts or shock waves, without

distorting smooth profiles or clipping narrow local extrema. The convection algorithm is based on non-artificiallydiffusive third-order upwinding in smooth regions, with automatic adaptive stencil expansion to (in principle,

arbitrarily) higher order upwinding locally, in regions of rapidly changing gradients. This is highly cost-e'ffective

because the wider stencil is used only where needed-in isolated narrow regions. A recently developed universal

limiter assures sharp monotonic resolution of discontinuities without introducing artificial diffusion or numerical

compression. An adaptive discriminator is constructed to distinguish between spurious overshoots and physical

peaks; this automatically relaxes the limiter near local turning points, thereby avoiding loss of resolution in

narrow extrema. Examples are given for one-dimensional pure convection of scalar profiles at constant velocity.

17. Key Words (Suggested by Author(s))

Upwind schemes

High-order convection

Nonoscillatory convection

Clipping eliminated

18. Distribution Statement

Unclassified- Unlimited

Subject Category 34

19, Security Classif. (of this report) 20. Security Classif, (of this page) 21. No. of pages 22. Price"

Unclassified Unclassified 24 A03

NASAFORM1626OCT86 *For sale by the National Technical Information Service, Springfield, Virginia 22161