Convergence of four-point implicit water wave models, Dr. Victor M. Ponce,

Convergence
of four-point implicit
water wave models

Victor M. Ponce,
Horst Indlekofer,
and Daryl B. Simons
Online version 2019

[Original version 1978]

1. INTRODUCTION

The numerical solution of the equations of open channel flow (the Saint Venant equations) is an established subject in the hydraulic literature. Implicit numerical schemes are generally preferred over explicit schemes (9) because they allow larger time steps and thus reduce the computational effort. Among the implicit schemes, the four-point scheme, also referred to as the Preissmann scheme (2, 5), has received wide attention. However, the stability and convergence properties of this scheme are as yet not fully understood.

The Saint Venant equations constitute a set of quasilinear partial differential equations. At present, there exists no method for analyzing the stability and convergence of numerical solutions of equations of this type. On the other hand, there are various methods for analyzing the stabiity and convergence of linear partial differential equations (11), among them, notably, the heuristic von Neumann approach (7). Traditionally, a way out of the mathematical difficulty imposed by the quasilinear nature of the Saint Venant equations has been to linearize them, assuming that the behavior of the linear system approximates that of the more complex quasilinear system. In practice, the linear analysis serves, if not as a quantitative, at least as a qualitative description of the nonlinear phenomena.

2. STABILITY AND CONVERGENCE

The essential feature of a numerical model is the replacement of a derivative such as ∂ƒ/∂s by a ratio of finite differences such as Δƒ/Δs. In doing so, the numerical model must satisfy certain requirements of stability and convergence. Stability refers to the ability of the numerical scheme to march in time without generating unbounded error growth. Convergence refers to the ability of the scheme to reproduce the terms of the differential equation with sufficient accuracy. Stability [in the sense of O'Brien, et al. (7)] is governed by round-off errors; convergence is governed by discretization errors. In general, stability is impaired if the discretization is made finer, while the converse is true from the convergence viewpoint.

In practice, stability is a necessary condition for model operation, since an unstable model is of little or no use. In implicit models, stability is usually achieved at the expense of convergence. Therefore, convergence is the focal issue: Once the model is shown to be stable, then it becomes necessary to assess its convergence properties.

The numerical properties of four-point implicit schemes have been studied by Liggett and Cunge (5) and Fread (3). Liggett and Cunge used a simple system of linear partial differential equations to derive conclusions on the numerical properties of the four-point implicit scheme. While they fully realized the limitations of their analysis, they nevertheless pointed out that the Saint Venant equations would be expected to behave in a similar manner. Actual computations have shown, however, that the 'theoretically better" value of the weighting factor θ = 0.5 is impractical from the stability viewpoint, and that in practice a θ ≥ 0.6 may be required (2).

Fread (3) has analyzed the numerical properties of the four-point implicit scheme. He used a linearized version of the Saint Venant equations, but simplified the equation of motion by neglecting the convective inertia and bed slope terms, and the equation of continuity by neglecting the wedge storage term. The complete linear analysis has been made possible by the study of two of the writers (8), who formulated the analytical solution for the propagation of sinusoidal perturbations in open channel flow on the basis of a linearized version of the Saint Venant equations (6).

The present study aims to determine the convergence properties of the four-point implicit numerical models of unsteady open channel flow. The celerity and attenuation of the numerical solution will be compared with those of the analytical solution (8). The von Neumann technique (7,11) will be used in this regard. The conclusions relate to the assessment of convergence as a function of the relevant physical and numerical parameters.

3. GOVERNING EQUATIONS

One-dimensional unsteady open channel flow, per unit of channel width, can be expressed by a set of two partial differential equations (4):

Water Continuity Equation:

∂d ∂u ∂d
^____ + d ^____ + u ^____ = 0
∂t ∂x ∂x (1)

Equation of Motion:

1 ∂u u ∂u ∂d
^____ ^_____ + ^____ ^_____ + ^_____ + (S_f - S_o ) = 0
g ∂t g ∂x ∂x
(2)

in which u = mean velocity; d = flow depth; g = acceleration of gravity; S_ƒ = friction slope; S_o = bottom slope; x = space variable; and t = time variable.

Equations 1 and 2 must satisfy the unperturbed flow (equilibrium or steady uniform flow) for which u = u_o and d = d_o as well as the perturbed flow for which u = u_o + u' and d = d_o + d'. The terms with a superscript represent small perturbations (unsteady flow components) to the steady uniform flow. Introducing the perturbed variables in Eqs. 1 and 2 yields, after linearization (neglecting second-order terms) (6):

  ∂d'             ∂u'              ∂d'
^____ + d_o ^_____ + u_o ^_____ = 0
  ∂t               ∂x               ∂x (3)

  1       ∂u'          u_o      ∂u'           ∂d'                    u'         d'
^____ ^_____ + ^____ ^_____ + ^_____ + S_o ( 2 ^____ - ^____ ) = 0
  g       ∂t            g       ∂x            ∂x                     u_o        d_o
(4)

4. ANALYTICAL SOLUTION

Two of the writers have shown (8) that the analytical solution of Eqs. 3 and 4 is characterized by the propagation celerity and attenuation of sinusoidal perturbations. For the celerity, they chose the dimensionless celerity c_{_*} defined as follows:

                         L
                     (^____ )
           c            T
c_{_*} = ^____  =  ^_______
          u_o           u_o
(5)

in which L = wavelength of the sinusoidal perturbation; and T = wave period. For the attenuation, they chose the logarithmic decrement δ defined as follows:

a_t_{_o}₊_T = a_t_{_o}e^{^δ} (6)

in which a_t_{_o} and a_t_{_o}₊_T are the wave amplitudes at the beginning and at the end of a period T.

Both c_{_*} and δ are functions of the Froude number of the equilibrium flow F_o and the dimensionless wave number σ_{_*}, in which F_o = u_o / ( gd_o )^1/2 and σ_{_*} = (2π / L) d_o /S_o (Appendix Ι).

5. NUMERICAL SOLUTION

The four-point implicit scheme is based on the following formulas (5), Fig. 1:

Fig. 1 Definition sketch for four-point implicit scheme.

θ (1 - θ)
ƒ (x,t ) = ^___ ( ƒ_j+1 ^{n + 1} + ƒ_j ^{n + 1} ) + ^_______ ( ƒ_j+1ⁿ + ƒ_j ⁿ )
2 2

(7)

∂ 1
^_____ ƒ (x,t ) = ^_____ ( ƒ_j+1 ^{n + 1} - ƒ_{j + 1}ⁿ + ƒ_j ⁿ⁺¹ - ƒ_j ⁿ )
∂t 2∆t

(8)

∂ θ 1 - θ
^_____ ƒ (x,t ) = ^_____ ( ƒ_j+1 ^{n + 1} - ƒ_j ⁿ⁺¹ ) + ^______ ( ƒ_j+1 ⁿ - ƒ_j ⁿ )
∂x ∆x Δx

(9)

in which ƒ(x,t) is any dependent variable; Δx and Δt are the space and time increments; and θ is the weighting factor of the implicit scheme. The weighting factor θ is introduced in the function and its space derivative to allow greater flexibility in the actual operation of the model. While the stable range is 0.5 ≤ θ ≤ 1.0, a strongly stable range is 0.6 ≤ θ ≤ 1.0.

The substitution of Eqs. 7-9 into 3 and 4 leads to discretized equations in u' and d'. The solution is postulated in the following exponential form:

u'
^____ = u_{_*} exp {i (σ_{_*} x_{_*} - β_{_*} t_{_*} )}
u_o (10)

d'
^____ = d_{_*} exp {i (σ_{_*} x_{_*} - β_{_*} t_{_*} )}
d_o (11)

in which x_{_*} = xS_o / d_o; t_{_*} = tu_oS_o / d_o and β_{_*}_R = 2πd_o / Tu_oS_o, and u_{_*} and d_{_*} are dimensionless amplitude functions.

The substitution of Eqs. 10 and 11 into the discretized equations in u' and d' yields:

Δt_{_*} Δt_{_*} ξ
u_{_*} [ i ^_____ (θξ + 1) tan α ] + d_{_*} [ i ^_____ (θξ + 1) tan α + ^____ ] = 0
Δx_{_*} Δx_{_*} 2 (12)

                ξ                   Δt_{_*}
u_{_*} [ F_o^² ^_____ + i F_o^² ^______ (θξ + 1) tan α + Δt_{_*} (θξ + 1) ]
                2                  Δx_{_*}
              Δt_{_*}                               Δt_{_*}
+ d_{_*} [ i ^______ (θξ + 1) tan α - ^_____ (θξ + 1) ] = 0
              Δx_{_*}                               2
(13)

in which i = √-1; ξ = [exp (-iβ_{_*} Δt_{_*}) - 1 ]; α = σ_{_*}Δx_{_*} / 2; Δt_{_*} = Δt u_o S_o / d_o and Δx_{_*} = Δx S_o / d_o.

Equations 12 and 13 constitute a set of two homogeneous equations in u_{_*} and d_{_*}. For the solution to be nontrivial, the determinant of the coefficient matrix must vanish. Accordingly:

                                    C                                          C           α          ξ^² F_o^²
[ (θξ + 1)^² (1 - F_o^²) ( ^____ ) ^² tan^²α + ξ (θξ + 1) ( ^____ ) ( ^____ ) + ^_______ ]
                                    c_{_*}                                          c_{_*}          σ_{_*}             4
                                           C                               C             α
+ i tanα [ ξ (θξ + 1) F_o^² ( ^____ ) + 3 (θξ + 1)^² ( ^____ ) ^² ( ^____ ) ] = 0
                                           c_{_*}                               c_{_*}             σ_{_*}
(14)

in which C / c_{_*} = Δt_{_*} / Δx_{_*}. The value C is referred to as the Courant number, and it is defined as the ratio of the physical celerity c to the "grid celerity" Δx / Δt. Equation 14 is a complex quadratic in ξ. It can be reduced to:

ξ² (M₂ + iN₂) + ξ (M₁ + iN₁) + (M_o + iN_o) = 0 (15)

in which:

M_o = (1 - F_o²)P² (16)

N_o = 3PQ (17)

M₁ = 2θM_o + Q (18)

N₁ = F_o^²P + 2θN_o (19)

F_o^²
M₂ = ^_____ + θQ + θ^²M_o
4 (20)

N₂ = θF_o^²P + θ^²N_o (21)

C
P = ( ^____ ) tanα
c_{_*} (22)

and

C α
Q = ( ^____ ) ( ^____ )
c_{_*} σ_{_*} (23)

Calling:

A = Q^² - P^² F_o^² (24)

B = -PQ F_o^² (25)

C = (A^² + B^²)^{^1/2} (26)

A + C
D = ( ^________ ) ^{^1/2}
2 (27)

A - C
E = ( ^________ ) ^{^1/2}
2 (28)

and solving the quadratic Eq. 15, the two roots are:

ξ_{_1,2} = M + iN (29)

in which:

-M₂ ( M₁ ± E ) + N₂ ( -N₁ ± D )
M = ( ^{_______________________________} )
2( M_₂^² + N_₂^² ) (30)

and

M₂ ( -N₁ ± D ) + N₂ ( M₁ ± E )
N = ( ^{_______________________________} )
2( M_₂^² + N_₂^² ) (31)

Since ξ = [exp ( -i β_{_*} Δt_{_*} ) - 1], it follows that:

N
( β_{_*}_R Δt_{_*} ) = - ^________
1 + M (32)

and

exp ( β_{_*}_I Δt_{_*} ) = [( 1+ M )^² + N ^² ]^{^1/2} (33)

The dimensioness celerity of the numerical solution c_{_*}_n is defined as follows (8):

β_{_*}_R 1 -N
c_{_*}_n = ^_______ = ^________ tan ^-1 ( ^_________ )
σ_{_*} σ_{_*} Δt_{_*} 1 + M (34)

or likewise:

                     1                            -N
c_{_*}_n = ^{_____________} tan ^-1 ( ^_________ )
                      C                         1 + M
            2α ( ^_____ )
                      c_{_*} (35)

The logarithmic decrement of the numerical solution δ_n is defined as follows (8):

                 β_{_*}_I            π ln [( 1 + M )² + N ²]
δ_n = 2π ^______ = ^{_______________________}
                 β_{_*}_R                              -N
                                   tan ^-1 ( ^_________ )
                                                  1 + M (36)

Since c_{_*} is a function of F_o and σ_{_*}, and M and N are functions of F_o, σ_{_*}, α, C, and θ, it follows that c_{_*}_n and δ_n are also functions of F_o, σ_{_*}, α, C and θ. F_o and σ_{_*} are the physical parameters, F_o being the Froude number and σ_{_*} a dimensionless wave number (σ_{_*} = 2πd_o / LS_o). The values α, C, and θ are the numerical parameters, in which α is the spatial resolution (α = π Δx / L); C is the Courant number (C = c Δt / Δx); and θ is the weighting factor of the implicit scheme.

5. CONVERGENCE ANALYSIS

The convergence analysis is based on the following ratios:

R₁ = exp ( δ_n - δ ) (37)

c_{_*}_n
R₂ = ^______
c_{_*} (38)

The value R₁ is the attenuation ratio and R₂ is the translation ratio. For R₁ > 1, the numerical damping is smaller than the physical damping; for R₁ < 1, the numerical damping is greater than the physical damping. For R₂ > 1, the numerical translation is greater than the physical translation; and for R₂ < 1, the numerical translation is smaller than the physical translation. For R₁ = R₂ = 1, there is an exact coincidence between analytical and numerical solutions.

The sensitivity of the convergence ratios to the physical and numerical parameters has been studied by varying the parameters within a specified range. For the Froude number F_o, the following range has been studied: 0.01 ≤ F_o ≤ 0.99. This excludes the critical and supercritical regimes. For the wave number σ_{_*}, the following range has been studied: 0.01 < σ_{_*} < 1,000. This range includes both the kinematic and diffusion waves (σ_{_*} ≅ 0.01), and the inertia-pressure waves (σ_{_*} ≅ 1,000). The spatial resolution has been varied between π / 100 ≤ α ≤ π / 10. This range includes both a very poor resolution (α ≅ π / 10) and a very high resolution (α ≅ π / 100). The Courant number has been varied in the range 0.25 ≤ C ≤ 4.0. The weighting factor θ has been varied in the range 0 ≤ θ ≤ 1.

The analysis has been performed for the primary wave only (8), by taking the minus sign associated with D and E in Eqs. 30 and 31. Space limitations preclude a complete analysis of both primary and secondary components of the wave propagation.

Fig. 2 Attenuation ratio R₁ as function of θ, σ_{_*}, F_o and C for: (a) α = π / 10; (b) α = π / 40; (c) α = π / 100.

Figures 2(a)-2(c) portray the attenuation ratio R₁ as a function of F_o, σ_{_*}, α, C, and θ. The examination of theses figures leads to the following conclusions:

For kinematic / diffusion waves (σ_{_*} ≅ 0.1) and inertia-pressure waves (σ_{_*} ≅ 1,000), the Froude number has a negligible effect on R₁. For dynamic waves (σ_{_*} ≅ 10), however, the effect of F_o on R₁ is very marked.
For a given α, the attenuation ratios for kinematic/diffusion waves are not unlike those for inertia-pressure waves. In fact, Fig. 2 shows that for a given α, the R₁, values for σ_{_*} = 0.1 and σ_{_*} = 1,000 are indistinguishable.
The higher the spatial resolution (the lower the α value) and the lower the Courant number, the closer is R₁ to 1 for all θ.
For kinematic / diffusion and inertia-pressure waves, the convergence of the wave amplitude (R₁ → 1) is obtained for values of θ equal to or slightly greater than 0.5. For dynamic waves, a value of θ much greater than 0.5 may be required to avoid an insufficient amount of numerical damping (R₁ > 1) which is tied to numerical amplification (instability).

The following general conclusions are obtained from Fig. 2:

The simulation is reasonably good for kinematic / diffusion waves and inertia-pressure waves, provided sufficient care is taken to minimize the amount of numerical damping. Large values of θ (θ approaching 1 from the left) can be used in conjunction with high spatial resolution. For poor spatial resolution, smaller θ values (θ approaching 0.5 from the right) may be required in order to avoid an excessive amount of numerical damping.
For dynamic waves (waves with both resistance and inertia), the accuracy of the simulation is highly dependent on the value of θ. The optimum value of θ is a function of σ_{_*} and α, and may be difficult to determine in practice.
Values of θ < 0.5 always cause numerical amplification; values of 0.5 ≤ θ < 1 may cause numerical amplification or attenuation; a value of θ = 1 always causes numerical attenuation (R₁ < 1). The existence and amount of numerical attenuation in the range 0.5 ≤ θ < 1 is a function of F_o, σ_{_*}, α, and C.

The analysis of the translation ratio R₂ shows that this ratio is nearly unity for a wide range of parameters F_o, σ_{_*}, α, and C. Predictably, for midrange dimensionless wavenumbers (σ_{_*} = 10) and poor spatial resolution, the ratios R₂ depart from unity. As an illustration of this trend, Fig. 3 shows the values of R₂ for varying F_o, σ_{_*}, C, and θ, and α = π / 10.

Fig. 3 Translation ratio R₂ as function of θ, σ_{_*}, F_o and C for: α = π / 10.

6. SUMMARY AND CONCLUSIONS

A comprehensive theoretical treatment of the convergence of the four-point implicit numerical model of shallow water waves is presented. The propagation celerity and attenuation factor of the numerical analog are derived, and convergence is tested by establishing the ratio of attenuation and translation given by the numerical and analytical solutions. Convergence is shown to be a function of the complex interaction of five parameters: (1) The Froude number of the equilibrium flow F_o; (2) the dimensionless wave number σ_{_*} (σ_{_*} = (2π / L)d_o /S_o); (3) the spatial resolution α (α = π Δx / L); (4) the Courant number C (C = c Δt / Δx); and (5) the weighting factor θ of the implicit scheme.

The convergence analysis is carried out for Froude numbers in the subcritical regime. Furthermore, due to space limitations, only the primary wave convergence is analyzed.

The following conclusions are warranted by this study:

The simulation is reasonably good for small σ_{_*} (kinematic / diffusion) waves and for large σ_{_*} (inertia-pressure) waves, provided sufficient care is exercised to minimize the amount of numerical damping.
For intermediate waves σ_{_*} (dynamic waves, in which both friction and inertia predominate), the accuracy of the simulation is highly dependent on the correct value of the weighting factor θ. In practice, an optimum value of θ that will assure both stability and convergence may be difficult to determine.
Values of θ < 0.5 always cause numerical amplification; values of 0.5 ≤ θ < 1 may cause numerical amplification or attenuation; a value of θ = 1 always causes numerical attenuation.

7. LIMITATIONS AND APPLICATIONS

The findings of this study have been obtained by using a linearized version of the Saint Venant equations. Therefore, emphasis should be placed on the qualitative nature of the results. At present, the performance of a nonlinear model can only be assessed through carefully designed numerical experiments. The linear analysis presented herein provides an appropriate framework for the design of these experiments.

8. ACKNOWLEDGMENTS

The research subject of this paper was supported in part by the Colorado State University Experiment Station. Horst Indlekofer was awarded a NATO fellowship to make possible his visit to Colorado State University.

APPENDIX Ι - EQUATIONS FOR PROPAGATION CELERITY c_{_*} AND LOGARITHMIC DECREMENT δ OF SINUSOIDAL PERTURBATIONS IN OPEN CHANNEL FLOW (8)

The equations are: c_{_*} = 1 ± D; δ = -2π (B ∓ E)/(|1 ± D |), in which A = (1/F_o^²) - B^²; B = 1/(σ_{_*} F_o^²); C = (A^² + B^²)^{^1/2}; D = [(C + A)/2]^{^1/2}; E = [(C - A)/2]^{^1/2}; F_o = u_o / (gd_o)^{^1/2}; and σ_{_*} = (2π / L)(d_o / S_o).

APPENDIX ΙΙ. REFERENCES

Chaudhry, Y. M., and D. N. Contractor. 1973. "Application of the Implicit Method to Surges in Open Channels," Water Resources Research, Vol. 9, No. 6, Dec., 1605-1612.

Cunge, J. A., and M. Wegner. 1964. "Numerical Integration of the Saint Venant Equations by an Implicit Finite Difference Scheme," La Houille Blanche, Grenoble, France, Vol. 19, No. 1, Jan., 33-39.

Fread, D. L. 1974. "Numerical Properties of Implicit Four-Point Finite Dfference Equations of Unsteady Flow," NOAA Technical Memorandum NWS HYDRO-18, National Oceanic and Atmospheric Administration, Mar.

Liggett, J. A. 1975. "Basic Equations of Unsteady Flow," Unsteady Flow in Open Channels, K. Mahmood and V. Yevjevich, eds., VoL. I, Water Resources Publications, Fort Collins, Colo.

Liggett, J. A., and J. A. Cunge. 1975. "Numerical Methods of Solution of the Unsteady Flow Equations," Unsteady Flow in Open Channels, K. Mahmood and V. Yevjevich, eds., Vol. I, Water Resources Publications, Fort Collins, Colo.

Lighthill, M. J., and G. B. Whitham. 1955. "On Kinematic Waves I, Flood Movements in Long Rivers," Proceedings, Royal Society of London, Vol. A 229, May, 281-316.

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, 223-251.

Ponce, V. M., and D. B. Simons. 1977. "Shallow Wave Propagation in Open Channel Flow," Journal of the Hydraulics Division, ASCE, Vol. 103, No. HY12, Proc. Paper 13392, Dec., 1461-1476.

Price, R. K. 1974. "Comparison of Four Numerical Methods of Flood Routing," Journal of the Hydraulics Division, ASCE. Vol. 100, No. HY7, Proc. Paper 10659, July, 879-899.

Quinn, F. H., and E. B. Wylie. 1972. "Transient Analysis of the Detroit River by the Implicit Method," Water Resources Research, Vol. 8, No. 6, Dec., 1461-1469.

Roache, P. J. 1972. Computational Fluid Dynamics, Hermosa Publishers, Albuquerque, NM.

APPENDIX ΙΙΙ. NOTATION

The following symbols are used in this paper:

A, B, C, D, E = parameters, Eqs. 24-28;

a = wave amplitude;

C = Courant number;

c = wave celerity;

d = flow depth;

F_o = Froude number;

ƒ = function;

g = acceleration of gravity;

i = √-1;

j = space discretization index;

L = wavelength;

M, M_o, M₁, M₂ = parameters;

N, N_o, N₁, N₂ = parameters;

n = time discretization index;

P, Q = parameters;

R₁ = attenuation convergence ratio;

R₂ = translation convergence ratio;

S_ƒ = friction slope;

S_o = bed slope;

T = wave period;

t = time variable;

t_o = initial time;

u = mean velocity;

x = space variable;

α = spatial resolution parameter;

β_{_*} = complex propagation factor;

β_{_*}_{_i} = imaginary part of β_{_*};

β_{_*}_{_R} = real part of β_{_*};

Δx = space step;

Δt = time step;

δ = logarithmic decrement;

θ = weighting factor; and

σ_{_*} = dimensioinless wavenumber.

Subscripts

o = equilibrium flow; and

n = numerical.

Superscripts

' = perturbed component; and

* = dimensionless parameter / variable.

220701 18:30

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