Unconditional Stability in Convection Computations


Victor M. Ponce, Yung Hai Chen and Daryl B. Simons


Online version 2015

[Original version 1979]



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+1jn) + (1-X) (ζ j+1 n+1j+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 at a particular grid location.

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

In Eq. 2, solving explicitly for the unknown value ζ j+1 n+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 = οx2). 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 ∂2 ζ / ∂ x2 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 tt approaches zero. Table 1 summarizes the stability and convergence (consistency in the limit) properties of three common first-order schemes.

TABLE 1.- Numerical properties of three common first-order explicit schemes
of the pure convection equation.
Scheme description

(1)
X

(2)
Y

(3)
μn

(4)
Stability
μn

(5)
Convergence
μn 0

(6)
I. Forward-in-time
Backward-in-space
0 1 (1/2) u Δx (1 - C) C ≤ 1 C 1
II. Backward-in-time
Forward-in-space
1 0 (1/2) u Δx (C - 1) C ≥ 1 C 1
III. Backward-in-time
Backward-in-space
0 0 (1/2) u Δx (1 + C) All C No C


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 = the wave frequency; and βI = the 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 Eq. 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 R1 and R2 are functions of X, Y, C, and L / Δx.

Fig. 2 shows the variation of R1 and R2 as a function of X, Y, C, and L / Δx. The following conlusions can be drawn from this figure:

  1. For C = 1 the exact coincidence of numerical and analytical solutions R1 = R2 = 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 R1 = R2 = 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 > 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 R1 convergence ratio (wave amplitude) is equal to 1 for X = Y = 1/2 and for all C. However, the R2 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 ∂3 ζ / ∂x 3. 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, C2 = C3 = 0, and ζ j + 1 n + 1 = ζj n, 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 n but also on ζj n + 1 and ζj + 1 n .


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.

TABLE 2.- Numerical properties of three unconditionally stable and second-order accurate
explicit schemes of the pure convection equation.
Scheme

(1)
C

(2)
X

(3)
Y

(4)
Stability
μn ≥ 0

(5)
Convergence
μn 0

(6)
IV C ≤ 1

C > 1
(1 - C)/2

0
1

(C +1)/2C
All C All C
V C ≤ 1

C > 1
(1 + C)/2

1
0

(C - 1)/2C
All C All C
VI All C 0.5 0.5 All C All C

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 ∂3ζ / ∂x3 (or ∂3ζ / ∂t 3). 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.


8.  APPENDIX |. - REFERENCES

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

  2. Hirt, C. W., "Heuristic stability theory for finite difference equations," Journal of Computational Physics, Vol. 2, 1963, pp. 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., Hyman, M. A., and Kaplan, S., "A study of the numerical solution of partial differential equations," Journal of Mathematics and Physics, Vol. 29, No. 4, 1950, pp. 231-251.

8.  APPENDIX ||. - NOTATION

C = Courant number;

C1 = coefficient, Eq. 4;

C2 = coefficient, Eq. 5;

C3 = 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;

R1 = wave amplitude convergence ratio;

R2 = 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.


160619 09:15

Documents in Portable Document Format (PDF) require Adobe Acrobat Reader 5.0 or higher to view; download Adobe Acrobat Reader.