1.  INTRODUCTION

The stability and convergence properties of explicit numerical schemes of the pure convection equation have been known for some time. A unified theoretical treatment is presented herein. The von Neumann (4) and Hirt (2) analyses are used to show how unconditional stability and second-order accuracy are both possible within the framework of an explicit scheme (6).

2.  PURE CONVECTION EQUATION

The pure convection equation is of the following form:

∂ζ              ∂ζ                   _____  +  u _____  =  0       ∂t              ∂x

(1)

in which ζ = any physical quantity to be convected with the velocity u. A general explicit finite difference scheme of this equation is as follows (Fig. 1):

X (ζj n+1-ζ jn) + (1-X) (ζ j+1 n+1-ζj+1n)          Y (ζ j+1 n-ζ j n) + (1-Y)(ζj+1n+1-ζ jn+1)                    [ __________________________________ ] + u [ _________________________________ ]  =  0                            Δt                                                                 Δx

(2)

in which X and Y = weighting factors of the temporal and spatial derivatives, respectively; Δt and Δx = temporal and spatial increments, respectively; and u = the average convective velocity.

FIG. 1.- Schematic of discretization of the pure convection equation.

In Equation 2, solving explicitly for the unknown value ζ _j+1 ⁿ⁺¹ gives:

ζ j+1n+1 = C1 ζjn + C2 ζj n+1 + C3 ζ j+1n

(3)

in which:

(X + CY)         C1 =   ___________________            (1 + C) - (X + CY)

(4)

C - (X + CY)         C2 =   ____________________             (1 + C) - (X + CY)

(5)

1 - (X + CY)         C3 =  ___________________            (1 + C) - (X + CY)

(6)

and

Δt         C =   u _______               Δx

(7)

is the Courant number corresponding to the particular grid location.

3.  NUMERICAL DIFFUSION COEFFICIENT

Equation 1 is an inviscid equation, therefore, viscosity (or diffusivity, as the case may be) is absent from it. However, a numerical solution of Eq. 1 will often show a viscous (diffusive) behavior. This is to be correctly attributed to the artificial viscosity of the scheme itself. The term "numerical dispersion" is sometimes used in this connection, ostensibly to denote the fact that both amplitude and phase errors aggregate to produce a dispersive type of effect in the numerical solution.

A measure of the amount of numerical dispersion can be obtained from the approximation error of the finite difference scheme. This can be obtained by expanding the grid function ζ(jΔx, nΔt) as a Taylor series about point jΔx, nΔt, leading to (1,2):

1                    1             ∂2 ζ R = u Δx [( ____ - X) + C ( ____ - Y )] _____                      2                    2              ∂x2                                1                        1                     ∂3 ζ + uΔx2 { (1 - C ) [ ____ ( X + CY) - ____ (1 + C ) ] } _____ + ο(Δx3)                                2                        3                     ∂x3

(8)

in which R = the approximation error. From Eq. 8, for X = Y = 1/2, the scheme is of second-order accuracy. i.e., R = ο(Δx²). In addition, for C = 1, the accuracy is at least of third order. In fact, for these conditions the exact solution is obtained (1).

The coefficient of the term ∂²ζ / ∂x² in Eq. 8 is referred to as the numerical diffusion coefficient, μ_n, defined as follows:

1                         1 μn = u Δx [(____  - X )  +   C (____  - Y)]                      2                         2

(9)

Positive values of μ_n cause the wave amplitude to diffuse numerically; negative values cause the wave amplitude to amplify numerically. While in the former the consistency is impaired, in the latter stability suffers. Both stability and consistency improve as μ_n approaches zero. Table 1 summarizes the stability and convergence (consistency in the limit) properties of three common first-order schemes.

4.  VON NEUMANN LINEAR ANALYSIS

A considerable insight into the numerical properties of the various four-point explicit schemes for the computation of pure convection can be gained by a linear analysis of stability and convergence following the von Neumann approach (3,4). According to this method, a solution of Eq. 2 is sought in sinusoidal form:

ζ = ζο exp [i ( σx - βt )]

(10)

in which ζ_ο = an amplitude function; σ = wave number of the sinusoidal solution; and β = complex propagation factor, such that β = β_R + i β_I. The quantity β_R = wave frequency; and β_I = amplification (or attenuation) factor. The substitution of Eq. 10 into Eq. 2 leads to:

[ (1 - X) exp ( i σ Δx) + X ] - CY [exp (i σ Δx) - 1] exp (- i β Δt ) = _______________________________________________________                            [(1 - X) exp ( i σ Δx) + X ] + C(1 - Y) [exp ( i σ Δx ) - 1]

(11)

Substituing exp (i σ Δx) = p + i q, in which p = cos (σΔx) and q = sin (σΔx) into Eq. 11 leads to:

[p (1 - R ) + R ] + iq (1 - R ) exp (- i β Δt ) = _______________________________                               [p (1 - S ) + S] + iq (1 - S )

(12)

in which:

R = X + CY

(13)

S = R - C

(14)

Operating on Eq. 12:

M + i N exp (- i β Δt ) = 1 - __________                                     P

(15)

in which:

M = (1 - p) (1 - S ) C

(16)

N = q C

(17)

and

P = 1 - 2 ( 1- p) (1 - S ) S

(18)

Separating Eq. 15 into its real and imaginary components:

N             tan (- βR Δt ) =   __________                              P - M

(19)

[( P - M)2 + N 2 ] 1/2            exp (- βI Δt ) = _____________________                                        P

(20)

Since the analytical solution of Eq. 1 exhibits no damping, it follows that Eq. 20 is a measure of numerical damping (or amplification) of the numerical solution after one time increment Δt. A convergence ratio with regard to wave amplitude is defined as follows:

R1 = exp (βI Δt)

(21)

From Equation 19, the phase velocity of the numerical solution is:

βR           1                        N un = ______  = ______ tan -1   ( _________ )             σ          σΔt                   P - M

(22)

A convergence ratio with regard to wave phase is defined as follows:

un             1                        N R2 = ______  = _______ tan -1   ( _________)             u          σuΔt                   P - M

(23)

or likewise:

1                            N R2  = _________   tan -1   ( _______ )           C (σΔx)                   P - M

(24)

in which σΔx = the spatial resolution. Since σ = 2 π /L, in which L is the sinusoidal wavelength, it follows that:

2 π σΔx = _______                L              ____               Δx

(25)

in which L / Δx is the number of discrete space intervals per sinusoidal wavelength. Both R₁ and R₂ are functions of X, Y, C, and L / Δx.

Fig. 2 shows the variation of R₁ and R₂ as a function of X, Y, C, and L / Δx. The following conclusions can be obtained from this figure:

1. For C = 1, the exact coincidence of numerical and analytical solutions R₁ = R₂ = 1 is given for values of X and Y along the diagonal joining X = 0, Y = 1 , and X = 1,Y = 0.

2. For C ≠ 1, there is no combination of X and Y  for which R₁ = R₂ = 1.

FIG. 2.- Convergence ratios R1 and R2 as a function of X, Y and C: (a) L / Δx = 10; (b) L / Δx = 40.

The instability, R₁ > 1, of Scheme I (Table I) for Courant numbers greater than 1, and of Scheme II for Courant numbers less than 1, is confirmed. Likewise, the unconditional stability and stronger numerical dispersion of Scheme III is confirmed.

5.  STABILITY AND CONVERGENCE PROPERTIES

Linear equation with constant coefficient. In this case, the convective velocity u is constant in space and time. The numerical solution using Eqs. 3-6 with X = Y = 1/2 and C = u Δt / Δx = 1 will exactly reproduce the analytical solution.

Linear equation with variable coefficient. In this case, the convective velocity u varies in space and time. Therefore, the Courant number varies from one computational cell to the next. The von Neumann linear stability analysis shows that the R₁ convergence ratio (wave amplitude) is equal to 1 for X = Y = 1/2 and for all C. However, the R₂ convergence ratio (wave phase) is equal to 1 only for C = 1. For Courant numbers different than one a certain amount of phase error will always remain.

The stability of the scheme requires that μ ≥ 0; on the other hand. convergence of the wave amplitude improves as μ_n → 0. Using Eq. 9, several unconditionally stable yet second-order accurate explicit schemes can be formulated. The numerical properties of three of these schemes are given in Table 2. Schemes IV and V can be regarded as generalizations of Schemes I and II, respectively, for the case of varying Courant numbers. However, unlike Schemes I and II, Schemes IV and V do not diffuse the wave amplitude numerically. Scheme VI, a space-and-time centered scheme, has second-order accuracy for all C.

For fixed values of X and Y the second-order error of Schemes IV to VI will be directly proportional to | 1 - C |. Likewise, for fixed values of X, Y and C (C ≠ 1), the second-order error will be directly proportional to ∂³ ζ / ∂x ³. The equivalence of Schemes IV to VI can be demonstrated by substituting the values of X and Y of Table 2 into Eqs. 4-6. This leads to:

C1 = 1

(26)

1 - C C2 =     ________               1 + C

(27)

and

1 - C C3 =     ________               1 + C

(28)

for all three schemes. Furthermore, for C = 1, C₂ = C₃ = 0, and ζ _{j + 1} ^{n + 1} = ζ_j ⁿ, i.e., the quantity ζ is being subjected to pure convection. For C ≠ 1 a dispersive effect is introduced because the calculation of ζ _{j + 1} ^{n + 1} is being based not only on ζ_j ⁿ but also on ζ_j ^{n + 1} and ζ_{j + 1} ⁿ .

6.  ILLUSTRATIVE EXAMPLE

A hypothetical example was set up to demonstrate the numerical properties of the schemes described in Table 2. A channel 20,000 ft long was divided into 20 reaches of Δx = 1,000 ft (304.8 m) each. A sinusoidal ζ distribution varying from 0 to 100 units in specified in the time axis as an upstream boundary condition. The distribution is specified in 20 intervals of Δt = 1,000 sec each, and the simulation is carried out for a total of 50 time intervals, or 50,000 sec. A convective velocity varying linearly in space between u = 0.5 fps (0.15 m/s) at the upstream boundary and u = 1.5 fps (0.45 m/s) at the downstream boundary is specified.

The results of the simulation are shown in Fig. 3. Stability, convergence of the wave amplitude, and mass conservation are maintained for all Courant numbers. Predictably, a small amount of numerical noise was detected.

FIG. 3. Results of hypothetical example.

This is caused by the second-order errors, and tends to be concentrated in the regions of sharp changes in ∂³ζ /∂x³ (or ∂³ζ /∂t ³). The numerical error is a function of | 1 - C | and vanishes as C → 1.

7.  SUMMARY AND CONCLUSIONS

A unified theoretical treatment of the numerical properties of a class of explicit schemes of the pure convection equation is presented. For the cases of slowly varying Courant numbers, the theory is used to show that unconditional stability and second-order accuracy can be achieved within the framework of an explicit formulation.

APPENDIX |. REFERENCES

1. Cunge, J. A. 1969. "On the subject of a flood propagation computation method (Muskingum method)," Journal of Hydraulic Research, Vol. 7, No. 2, 205-230.

2. Hirt, C. W. 1963. "Heuristic stability theory for finite difference equations," Journal of Computational Physics, Vol. 2, 339-355.

3. Leendertse, J. J., "Aspects of a computational model for long-period water-wave propagation," Memorandum RM-5294-PR, May, 1967, The Rand Corporation, Santa Monica, Calif.

4. O'Brien, G. G., M. A. Hyman, and S. Kaplan. 1950. "A study of the numerical solution of partial differential equations," Journal of Mathematics and Physics, Vol. 29, No. 4, 231-251.

APPENDIX ||. NOTATION

C = Courant number;

C₁ = coefficient, Eq. 4;

C₂ = coefficient, Eq. 5;

C₃ = coefficient, Eq. 6;

j = space discretization index;

L = wavelength;

M = parameter, Eq. 16;

N = parameter, Eq. 17;

n = time discretization index;

P = parameter, Eq. 18;

p = cos (σΔx);

q = sin (σΔx);

R = approximation error; also, parameter, Eq. 13;

R₁ = wave amplitude convergence ratio;

R₂ = wave phase convergence ratio;

S = parameter, Eq. 14;

u = mean velocity;

μ_n = phase velocity of numerical solution;

X = weighting factor of temporal derivative;

Y = weighting factor of spatial derivative;

β_I = amplification (or attenuation) factor;

β_R = wave frequency;

Δt = temporal increment;

Δx = spatial increment;

ζ = any physical quantity to be convected;

μ_n = numerical diffusion coefficient, Eq. 9; and

σ = wave number.