NASA Technical Memorandum 106203 ICOMP-93-14 j c:_-: "1 _i:" f Large Time-Step Stability of Explicit One-Dimensional Advection Schemes INI M_ ,.4 ,._ m 0 ! .-- 0 _- U O0 Z _ 0 _) L_ B.P. Leonard Institute for Computational Mechanics in Propulsion Lewis Research Center Cleveland, Ohio and The University of Akron Akron, Ohio May 1993 w Z 0 Z _1_ Z_ZZ LARGE TIME-STEP STABILITY OF EXPLICIT ONE-DIMENSIONAL ADVECTION SCHEMES B.P. Leonard* Institute for Computational Mechanics in Propulsion Lewis Research Center Cleveland, Ohio 44135 and The University of Akron Akron, Ohio 44325 ABSTRACT There is a wide-spread belief that most explicit one-dimensional advection schemes need to satisfy the so-called "CFL condition" -- that the Courant number, c = uAt/Ax, must be less than or equal to one, for stability in the von Neumann sense. This puts severe limitations on the time-step in high-speed, fine-grid calculations and is an impetus for the development of implicit schemes, which often require less restrictive time-step conditions for stability, but are more expensive per time-step. However, it turns out that, at least in one dimension, if explicit schemes are formulated in a consistent flux-based conservative finite-volume form, von Neumann stability analysis does not place any restriction on the allowable Courant number. Any explicit scheme that is stable for c < 1, with a complex amplitude ratio, G(c), can be easily extended to arbitrarily large c. The complex amplitude ratio is then given by exp(-tN0) G(Ac), where N is the integer part of c, and Ac = c -N (< 1); this is clearly stable. The CFL condition is, in fact, not a stability condition at all, but, rather, a "range restriction" on the "pieces" in a piece-wise polynomial interpolation. When a global view is taken of the interpolation, the need for a CFL condition evaporates. A number of well- known explicit advection schemes are considered and thus extended to large At. The analysis also includes a simple interpretation of (large-At) TVD constraints. INTRODUCTION Consider the model one-dimensional pure advection equation for a scalar _b(x,t) 0____¢= _ u a_.____ (1) at #x where u is a constant advecting velocity. Take a uniform space-time grid (Ax, At) and integrate Equation (1) over Ax and At, giving afinite-volume formulation i where the bars indicate spatial averages over cell i at time-levels n and (n+ 1), and left and *Work funded under Space Act Agreement NCC 3-233. 2 right time-averaged face-values have been introduced. The Courant number is given by c = u&t/&x (3) The notation is defined in Figure 1, which also shows advective characteristics (along any one of which _b is constant) entering the left face (for u > 0). Given the set of current m spatial-average values, _, one needs to estimate the corresponding time-averaged face- values in order to explicitly update according to Equation (2). In avon Neumann stability analysis, ¢(x,O is written as a wave [Fletcher, 1990] ¢(x,0 = a(t) exp(, ) where k is the wave number and t represents the imaginary unit, substituted into Equation (1), the exact solution is [Leonard, 1980] ¢,(x,0 = A(0) exp[, k(x-uO] corresponding to a travelling wave, with _ = const along characteristics (x = ut). exact complex amplitude ratio is then Go_¢ t _ d_(x,t+At) = exp(-Lc0) ¢(x,0 where 0 is the nondimensional wave-number Note that if c > 1, Gc_ct can be written 0 = k&x ¢_utct where N is the integer part of c (5) The (6) (7) = exp(-t NO) exp(-, &cO) (8) N = INT(c) (9) (12) and _"" = A(O) exp{,k[x,-u(n+l)At]} (sinO-/21 i 0/2 / -_ = A(O) exp{,k[x,-unAt]} (sin 0/2 / (11) and AC the non-integer part (necessarily less than one) AC = c-N (10) By integrating Equation (4) from (xi - Ax/2) to (xi + Ax/2), it is not difficult to show that the exact cell-average values are 3 Hence, the exact complex amplitude ratio of cell-average values is the same as that given by Equation (8) n+l Gcx_e, - -i - Gera= t = exp(-t NO) exp(-t acO) (13) thus, giving a reference value against which to compare numerical G's. In the following sections, the von Neumann stability of some well-known explicit advection schemes is studied- first in terms of the "conventional" time-step ranges, and then for larger Courant number values. The latter extension stems naturally from identifying the sub-grid interpolation used in estimating the time-averaged face-values in Equation (2). It will be shown, in particular, that if G,s_ (c) is the (conventional, small time-step) complex amplitude ratio for a given scheme for c < 1, then the natural finite-volume extension to arbitrarily large time-step has a complex amplitude ratio of the form n÷l LT$ _i - = G._(Ac) for c >1 (14) G°um (c) exp (-L NO ) srs which should be compared with Equation (13). To the extent that GS_(Ac) is a good approximation to exp(-t AC 0), then the large-At G.L_ (C) is an even better approximation to the exact G. Viewed differently, a numerical simulation (over a given distance) will be more accurate for larger At because, as will become evident, numerical distortion depends only on the total number of time steps and arc-- not on c itself. EXPLICIT ADVECTION SCHEMES The time-averaged face-values appearing in Equation (2) can be rewritten in terms of spatial averages. For example, the time averaged left-face value is given by ¢bt = A'_l joX' _bt(r) dr = c_1 Jx:',-ca,_b"(_)d_ (15) where be(z) is the instantaneous face-value and _"(_) is the local sub-grid spatial behaviour in the region upstream of the face in question, at time-level n. The relationship should be clear from Figure 1; note that uAt = CZ_C. A similar formula holds for the right face. But it is not necessary to write this out explicitly because advective flux conservation guarantees that _br(i) = _t(i+l) (16) Different numerical schemes result from different choices of _b'(_) in estimating the local behaviour. Note in particular that for consistency, _b"(_) should obey the integral constraint ___1 [*.._,24_,(0d_ = _7 for all i (17) First-Order Upwinding One of the simplest advection schemes results from assuming that ¢n(_) is piece-wise constant, for each i, _n(_) = _: for (x,---_-) <_< (x,+-_-) (18) with discontinuities at the cell faces, as shown in Figure 2, the hatched region depicting the integral in Equation (15). This, of course, trivially satisfies Equation (17), and the left-face value given by Equation (15) for a positive Courant number less than (or equal to) one is simply ¢_,(i) = _-t for 0 <c <1 (19) and, similarly _r(i) -- _7 for 0 < c <1 (20) The update equation, Equation (2), thus becomes _,,.1 = _7 - c(_7-__1) for 0 <c <1 (21) i The corresponding complex amplitude ratio is GIu(C) = 1 - c[1 - exp(-_0)] (22) or GIu(c ) = 1 - c(1-cos0) - Lcsin0 (23) For the full range of numerical wave-numbers (0 < 0 < a-), this represents a semicircle of radius c in the lower half of the complex plane, passing through the point (1,0). The scheme is thus stable (I GI - 1)provided the CFL condition (c <_ 1) is satisfied. Figure 3 shows a polar plot of G_u for c = 0.75. Making a Taylor expansion of real and imaginary parts gives c 02 + 0(04 ) _ _ [cO + 0(03)] (24) Glv(c ) = 1 -._ whereas a Taylor expansion of the exact G gives G o_,ct(c) = cos(c0) - _ sin(c0) c202 +0(04 )-_[cO- c303 = 1- T -g- ÷ o(o')] (25) Thus Gtu indeed matches Gcx_c t only through first-order terms in 0, with discrepancies in the second-order term (unless c - 1, in which case exact cell-to-cell transfer occurs across one mesh-width: _.1 = __1). --i Second-Order Methods It is convenient to define the "order" of an advection scheme as the power of 0 up to and including which the Taylor expansion of G,m matches that of G_ct given by Equation (25). One of the best known second-order methods is that due to Lax and Wendroff [1960], equivalent in the scalar case to Leith's method [1965]. For this scheme, the sub-grid ffn(_) is assumed to be piece-wise linear satisfying Equation (17). Upstream of the left face, for example, a straight line is drawn (for u > 0) between _7-i and cbn , considered to be located at the centers of the respective finite-volume cells. A similar construction xs made across each of the other cells, resulting in discontinuities at cell faces. This is shown in Figure 4, which also indicates the integral in Equation (15) for 0 < c < 1. In this case, i.e., for a Courant number less than (or equal to) one, -- ,or - upstream of the left face, across cell (i-1). Substitution into Equation (15) results in 1 (_7 + TM c . (27) - - _,-1)- (_7- 4,_1) 6t(i)- 2 According to Equation (16), _r is obtained by increasing each index by 1. update equation as 7..1 77- c(7,.,- 7;o * c2 This gives the (28) Because of the symmetry about cell i, this formula is actually valid for positive or negative Courant number, provided Icl -< 1 (29) The negative-c case is shown in Figure 5, where the sub-grid _"(_) for computing 4t is the same formula as that in Equation (26)- but across cell i rather than (i-1). The numerical complex amplitude ratio for the Lax-Wendroff scheme can be obtained directly from Equation (28) as C 2 c [exp(L0) - exp(-L 0)1 + [exp(_ 0) - 2 + exp(-L 0)1 GLw(C ) = 1 - -_ -_- (30) or GLw(C ) = 1 - C2(1- COS0) - tC sin 0 (31) with a corresponding Taylor expansion c2 02 + 0(04) - ,[cO - c 03 + O(05)] (32) G_w(c) = 1-T 6 thus matching the expansion of Ge..ot, Equation (25), through 02 terms, as expected for a second-order scheme. Once again, exact point-to-point transfer across one mesh-width occurs when c _ 1 (and also, in this case, when c - -1). Figure 6 shows a polar plot of GLw, a semi-ellipse with vertical semi-axis equal to c and curvature near (1,0) equal to that of the unit circle. The scheme is stable provided the CFL condition is satisfied--and unstable otherwise. The Lax-Wendroff scheme is symmetrical in the sense that a given face-value (assuming lcl <--1) depends only on the two immediately adjacent cell-average values on either side of the face in question. Figure 7 shows an alternate piece-wise linear sub-grid reconstruction (for positive c) using the two immediate upstream cell-average values. Although the reconstruction of _'(_) looks superficially the same as that in Figure 5, this represents an entirely different advection scheme in that the integral of Equation (15) is computed from the left (rather than the right, as in Figure 5). This will be found to generate a second-order advection scheme. Because of the upwind bias, it is now commonly known as "second-order upwinding" [Fletcher, 1990]. Upstream of the left face, across cell (i -1), _bn(_) is given by _n(_) = _;1 + _- ;2 (__x,_l) for i-,--_- < _ < i-1 + 2) (33) Equation (15) then gives the left-face value as _bt = gi-t - _bi_2 for 0 < c < 1 (34) In fact, this equation is valid for 0 < c < 2. This is an important point (usually ignored in CFD literature) because, as described below, it suggests a natural and computationally efficient way of extending explicit advection schemes to arbitrarily large Courant numbers. Using the corresponding _,, the update equation for second-order upwinding becomes i = - - _'i-1 4'i-2 (35) and the complex amplitude ratio is G2u(C) = 1-c(_-_)+c (2-c)cos0-c (_-_-_) cos 20 , in0 (36) This has a Taylor expansion given by the following 7 C2 02+ O(04) - t [C0--(.3C262c)03+ 0(05)] G2v(c ) = 1 --_- (37) Note in particular that, for c = 1, G2u(1) = cos0 - t sin0 = G_.._t(1) (38) corresponding to exact cell-to-cell transfer across one mesh-width. that for c = 2, G2o(2) = cos20 - _ sin20 = Go,.,ct(2) with exact cell-to-cell transfer across two mesh-widths. stable throughout the entire range But, in addition, note (39) In fact, second-order upwinding is IG2oI <--1 for 0 <_ c <_ 2 (40) Two typical polar plots (for c = 0.75 and 1.75) are shown in Figures 8 and 9, respectively. In order to prove (40), it is convenient to rewrite G2v(c) in the following form, starting with Equation (36), G2v (c) = (cos0 -,sin 0) [1 - (c-1)2 (l_cos 0) -,(c-X) sin0] (41) Note, by comparing Equation (31), that this is of the form G2u(c ) = exp(-t0) Gl.w(AC) (42) where Ac = c - 1. based on Ac, Of course, l exp( - _0) I -= 1, and for the "Lax-Wendroff" component IGLw(ZXc)l _< 1 for -1 < Ac < 1 (43) thus proving (40). It should be stressed that second-order upwinding is one well-known explicit advection scheme for which the CFL condition does not apply. In order to gain some insight into the operation of second-order upwinding for a Courant number range of 1 < c < 2, Figure 10 shows the situation for c = 1.25. In this case, it is not hard to show, using Equation (15), that = ac _ (4_i-l ÷ 4_i-2) - -- (4_i-1 - 4_i-2) ÷ 4_i-t (44) with a similar formula for c4_, (with all indexes increased by 1). The term in square brackets is clearly the I._-Wendroff face-value based on AC --but for cell (i-1) rather than cell i. This means that the update equation takes the form 8 (45) with Ac = c - 1. This is algebraically identical to Equation (35). But, by comparison with Equation (28), it can be seen that the right-hand side is the Lax-Wendroff update of cell (i-1) using Ac rather than c itself. This gives a hint as to how flux-based conservative finite-volume explicit advection schemes might be extended to arbitrarily large At. The details will be examined in the next section. As is wen known, the conventional Lax-Wendroff scheme suffers from severe phase- lag dispersion [Fletcher, 1990]. On the other hand, second-order upwinding introduces phase-lead dispersion for 0 < c < 1. Fromm's method [1968] of "zero-average-phase- error" aRempts to minimize the dispersion by taking the arithmetic-mean of the respective sub-grid reconstructions, giving, from Equations (26) and (33) +(gT-gL/(_-x,_,)for(x,_,_ (x,_,T) 4_"(_) = -"_i-_ _ 2_ '/ --T ") < _ < + Ax (46) i.e., the slope across any cell is parallel to the chord joining g-values through the centers of the two adjoining cells on either side. This is shown in Figure 11. To compute the face- value, first rewrite Equation (34) for second-order upwinding as 1(_7 (_7_ _ _ _ ÷ _,__) _,L)- c 4_t(2U) = _ (47) which is seen to be the Lax-Wendroff value together with a correction term proportional to the upwind-biased second-difference. Then the left face-value for Fromm's method is seen to be, for Courant number less than (or equal to) one, _ 1 [_kt(LW) + _bt(2U)] _bt(Fromm ) - *_,"-0 c -,,-0 - ( )(_7-2_7-x_,"__) = _ and the corresponding update equation becomes (48) mg 1 _,g"÷a = ok,-"- --2c(_k,., - g,"_,) + Tc2 (g,".l -297 + _b,_l)-" c (2_) -3g7+397_1-" + (gT., -*,-2) (49) with a complex amplitude ratio GFr(c ) = GLw(C) + C(_-£)[exp(,0) - 3 + 3exp(-,0) - exp(-,20)] (50) Or