- Email: [email protected]
Stability and accuracy of finite element schemes for the one-dimensional kinematic wave solution
Advances in Water Resources 25 (2002) 427–438 www.elsevier.com/locate/advwatres
Stability and accuracy of finite element schemes for the one-dimensional kinematic wave solution Fouad H. Jaber, Rabi H. Mohtar
*
Agricultural and Biological Engineering Department, Purdue University, West Lafayette, IN 47907, USA Received 22 May 2001; received in revised form 2 November 2001; accepted 8 January 2002
Abstract Solving the kinematic wave equations for overland flow using the conventional consistent Galerkin finite element scheme is known to result in numerical oscillations due to the non-symmetric first spatial derivative terms in the kinematic wave equations. In this paper the lumped and the upwind finite element schemes are evaluated as alternatives to the consistent Galerkin finite element scheme. Stability analysis of the upwind scheme shows that the damping effect, that could reduce the oscillations, is small for the high Courant numbers encountered in overland flow problems. The upwind scheme, using upwind factors of 0.1 and 1.0, did not provide any improvement to the stability of the lumped and the consistent schemes. The lumped scheme considerably reduces oscillations without significant reduction in the overall solution accuracy. No analytical guidelines for time-step criteria that will insure stability and accuracy were provided by the stability analysis performed for the three schemes. Problem specific accuracybased dynamic time-step criteria was developed and evaluated for the lumped scheme. These time-steps were found to be on average, double the size of the dynamic time-steps for the consistent scheme. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Finite element analysis; Numerical stability; Dynamic time-step; Upwind scheme; Lumped scheme; Consistent scheme; Fourier analysis; Kinematic wave equations; Overland flow
1. Introduction 1.1. Numerical stability Numerous studies reported on solving the kinematic wave problems using finite element analysis [11,24,30]. The numerical solutions of these equations have encountered stability and convergence problems. Gresho and Lee [10] defined the conventional Galerkin finite element method as having (a) a centered difference like treatment of advection terms and (b) a consistent scheme for coupling time derivative terms. Christie et al. [5] stated that central difference approximations of first derivatives should be avoided since they give rise to numerical oscillations for reasonable grid sizes. They added that the Galerkin finite element method using symmetric pyramid basis linear shape functions reproduce central difference formulations with their inherent oscillatory properties. Hughes [15] reported that diffi*
Corresponding author. Tel.: +1-765-494-1791; fax: +1-765-4961115. E-mail address: [email protected] (R.H. Mohtar).
culties are encountered when the advection operators are non-symmetric and result in spurious oscillations. According to Fiedler and Ramirez [7], flow depths are generally small enough that small oscillations may destroy the solution of overland flow problems. They also added that the oscillation control is a formidable challenge given the oscillatory nature of the equations. Zhang and Cundy [32] used a MacCormak finite difference scheme, which is an explicit two-step finite difference scheme, to solve overland flow problems. They reported that numerical oscillations were evident in most hydrographs and in some cases the solution was aborted. Abbott and Basco [1] reported that the shallow water equations are known to be prone to high frequency Fourier components, 2Dx oscillations, when using central difference schemes because the value of the slope estimated by the difference is independent of the value of the dependent variable at the mid-point between two nodes. Chan and Williamson [4] reported that partial differential equations, containing first derivative terms in space, exhibit numerical oscillation in the Galerkin finite element solution if the size of the mesh exceeds a critical value based on the Peclet number being
0309-1708/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 9 - 1 7 0 8 ( 0 2 ) 0 0 0 0 5 - 2
428
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
smaller than two. The authors did not discuss cases where the Peclet number is infinitely large such as the case of overland flow problems. Nevertheless, an argument might be made for a fine finite element grid and small time-step to solve the problem, however, this will require unacceptable computational cost [29]. Some literature related to the numerical oscillation of the Galerkin solution of the overland flow problem illustrate the upwind methods as a possible solution for oscillations caused by the centered difference scheme. Chan and Williamson [4] described upwind schemes as replacing the central difference scheme by backward difference. They added that this is usually at the expense of accuracy. Heinrich et al. [12] reported that upwind finite element schemes consist of using weighting functions of non-symmetric forms different from those used as shape functions. They concluded that using the appropriate shape functions would result in an accuracy loss equivalent to the finite difference upwind schemes. Gresho and Lee [10] and Leonard [18] argued that the upwinding scheme by itself is a smoothing approach that while providing a stable solution can be very misleading with regards to the accuracy of the problem. Tisdale et al. [29] developed a streamline upwind scheme for a two-dimensional overland flow solution. The scheme consists of isolating the advective terms in the equation according to a main slope and solve the problem using upwind techniques. They reported that this scheme solved the problem without oscillations with negligible effect on accuracy. However, their study was based on a modification of the original equation and is limited to simple surfaces where the downslope node of each element can be determined using a subroutine within the program. This method is yet to be tested for complex runoff surfaces. Fiedler and Ramirez [7] used an upwind MacCormack finite difference scheme for solving overland flow problems. They found that in addition to upwinding, a smoothing function might be needed to control oscillations. The smoothing function that they propose is a function that smoothes regions of large gradients such as high frequency oscillations while leaving smooth areas undisturbed using a weighted average of flow depth along the nodes. Anderson et al. [2] explained that oscillations are due to the phase errors in the Fourier components of the numerical solution. They defined a stable numerical scheme as the one for which errors from any source are not permitted to grow in the sequence of numerical procedures as the calculation proceeds from one marching step to another. They also suggested the Fourier analysis to establish stability limits for linear partial differential equations with the possibility to extend these guidelines to non-linear equations in an approximate sense. Katopodes [17] used a Fourier analysis to study the stability and the oscillations of the linear kinematic wave equation for channel flow with constant
kinematic wave celerity and no lateral inflow. He stated that in the presence of discontinuities most discrete methods generate spurious oscillations behind numerical jumps, which can only be eliminated by selective algorithmic damping. The damping need to be very sharp and localized in the area of jumps because it is usually accompanied with a strong dispersion of celerity. Singh [26] showed that linear stability analysis could be applied to the linearized kinematic wave equations and still allows us to identify obviously unsuitable schemes or appropriate step length for conditionally stable schemes. The use of the lumped matrix to replace the consistent matrix has been reported to reduce oscillations in channel flow computations [31]. 1.2. Time-step criteria The task of choosing the proper time-step has often been considered a matter of experience. Bajracharya and Barry [3] stated that reducing the time-step does not necessarily result in more accurate solutions. They added that optimal solutions that would provide stability and accuracy with minimal computational time are yet to be found. Hromadka and DeVries [14] argued, after a series of tests varying Dx and Dt, that the use of the kinematic wave method for channel routing needs evaluation for use in hydrologic models unless guidelines are developed to control the arbitrary use of the kinematic wave in design studies. They added that kinematic wave programs need internal checks to select Dx and Dt such that an accurate solution is achieved. Jaber and Mohtar [16] developed empirical equations for dynamic time-steps [19–21] that will insure accuracy and stability of the consistent scheme for the kinematic wave solution. These dynamic time-steps vary with the number of elements and the time of concentration of the watershed. Their study was limited to the conventional Galerkin method using the consistent scheme. They reported significant oscillations for time-steps larger than the dynamic time-steps without suggesting ways to deal with these oscillations. A more stable scheme needs to be identified and tested and new time-step criteria need to be developed for this scheme. 1.3. Research objectives and limitations When using the kinematic wave equation, kinematic shock may arise where faster waves overtake slower ones, which leads to breaking of the wave profile. Singh [26] presented approximate methods in shock fitting to remedy this problem. It must be emphasized that the use of numerical techniques ignoring the existence of these shocks and not using a shock fitting technique is not a valid application of the kinematic wave theory. The problems used in this study to test different schemes and to develop and test dynamic time-steps did not result in
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
kinematic shocks. The kinematic wave model in the form presented in this paper is the basis of several commercial and research models available in the market that approximate watersheds as runoff planes. It is one of the methods available to calculate excess runoff in HEC-1, the Hydrological Simulation Program-Fortran (HSPF), the Institute of Hydrology Distributed Model (IHDM), the EPA Storm Water Management Model (SWMM) and the USDA-Water Erosion Prediction Project (WEPP) among others [25]. There has been some research to solve the kinematic wave equation for converging [28] and diverging watersheds [27]. In both cases additional terms are added to the conservation of mass equation and change it from its original form. Results that apply to the original kinematic wave equation do not necessarily apply to the modified equations. The overall objective of this paper is to improve the numerical solution of the one-dimensional kinematic wave equation for overland flow by addressing the two identified causes of oscillations: (a) the consistent mass matrix for coupling the time derivatives and (b) the central difference treatment of the advective term. The specific objectives include: 1. Evaluate the effectiveness of introducing non-symmetric upwind scheme to replace the symmetric central difference treatment of the advective term. 2. Evaluate the effectiveness of replacing the consistent scheme by a lumped matrix. 3. Develop and evaluate a dynamic time-step criteria for the most stable scheme identified from the first two objectives. 2. Numerical schemes for solving the kinematic wave equation A brief description of the kinematic wave equations for overland flow follows. Details and derivation of these equations are given by Chow et al. [6], and Henderson [13] among many others. The one-dimensional kinematic wave equation is governed by the continuity equation:
429
re ðx; tÞ is the excess rainfall rate (m/s); S0 is the element bottom slope in m/m; Sf is the friction slope in m/m; and x (m) and t (s) represent the space and time domains, respectively. Franco and Chaudhry [9] stated that interpolating h and q provides the most accurate and stable solution of this kinematic wave problem. The Galerkin finite element solution of this problem consists of expressing independent variables as: n X /ðeÞ ¼ Ni ðxÞUi ; ð4Þ i¼1
where Ui are known values of the continuous function /ðeÞ of each of the nodal points of an element (e), Ni ðxÞ are shape functions that approximate the function of /ðeÞ based on the nth nodal values [30]. The shape functions at nodes j 1 and j for one-dimensional elements are shown in Fig. 1 and are represented mathematically as: Nj1 ¼
Xj x Xj Xj1
and
Nj ¼
x Xj1 : Xj Xj1
ð5Þ
A description of the three numerical schemes considered in this paper is presented in the following section. 2.1. The consistent scheme In the consistent scheme oh=ot is assumed varying linearly between nodal values and is expressed as ( ) oUj1 oh ot ¼ Nj1 Nj : ð6Þ oUj ot ot The residual integral of the Galerkin finite element formulation is
Z Xj ðeÞ oh oq ð7Þ R ¼ ½N ðxÞT þ rðx; tÞ dx ¼ 0: ot ox Xi Integrating Eq. (7) and arranging the results yield the system of ordinary differential equations (ODE) of the form:
oh ½C þ ½Bfqg ff g ¼ 0; ð8Þ ot
oh oq þ ¼ re ðx; tÞ ð1Þ ot ox and the conservation of momentum equation for the kinematic wave assumption is reduced to: S0 ¼ Sf :
ð2Þ 2
The unit width flow (q) in m =s and the vertical flow depth (h) in m in Eq. (1) are related by the following energy equation: q ¼ ahb ;
ð3Þ 1=2 ðS0 =nÞ;
b ¼ 5=3 where, using Manning’s equation, a ¼ in which n is the Manning roughness coefficient (m1=3 =s);
Fig. 1. The weighting functions for the upwind c ¼ 1, c ¼ 0:1 and c ¼ 0.
430
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
where [C] and [B] and ff g are global matrices constructed using the direct stiffness method from the following element matrices, respectively: ðeÞ l 2 1 ðeÞ 1 1 1 C ; B ¼ ¼ 6 1 2 2 1 1
ðeÞ l rðx; tÞ and f ¼ ; ð9Þ 2 rðx; tÞ where l is the length of the element in m. 2.2. Lumped scheme The lumped scheme [23] is developed by defining the variation in oh=ot with respect to x to be constant between the midpoints of adjacent elements. The scheme uses alternative shape functions N for the time derivative defined as
Nj1 ¼ 1 H s 2l and Nj ¼ H s 2l ; ð10Þ where H ðxÞ (Heavyside function) is defined as H ðxÞ ¼ 0
if x < 0;
H ðxÞ ¼ 1
if x >¼ 0;
where l is the length of the element in m. The function that was used in this work and satisfies these conditions is (Fig. 1): F ðxÞ ¼ 3
ðx Xj Þðx Xj1 Þ ðXj Xj1 Þ2
:
ð14Þ
The residuals of the advective term in the kinematic wave equation oq=ox for one element are: Z Xj oN ðeÞ q dx Wj1 Rj1 ¼ ox Xj1 Z Xj oN ðeÞ and Rj ¼ q dx; ð15Þ Wj ox Xj1 solving these integrals for c ¼ 1:0 results in the following: ðeÞ
Rj1 ¼ 0; ðeÞ
Rj ¼ qj1 þ qj or in matrix form
Rj1 qj1 0 0 ¼ : Rj qj 1 1
ð16Þ
ð17Þ
s is the distance from node j 1. The element matrix [C] for the lumped scheme is: ðeÞ l 1 0 : ð11Þ C ¼ 2 0 1
For an upwind factor c ¼ 0:1 the results becomes
qj1 Rj1 0:45 0:55 ¼ : ð18Þ Rj qj 0:55 0:55
The [B] matrix remains the same as in the consistent method.
For c ¼ 0 the matrix reduces to the non-upwinded [B] matrix of Eq. (9).
2.3. Upwind finite element scheme 3. Methodology for dynamic time-step development An ‘‘upwind’’ finite element scheme to solve onedimensional kinematic wave equations is presented in this paper as a function of Xj1 and Xj , which are, respectively, the coordinates of the nodes of an element with nodes j 1 and j based on Heinrich et al. [12] who developed an upwind finite element scheme for convective transport equation. Element matrices for this scheme are developed in this paper. The one-dimensional upwind scheme of Heinrich et al. [12] requires that the advective term weighting function of the equation oq=ox for node j 1 be Wj1 ¼ Wj1 ðx; cÞ ¼ Nj1 þ cF ðxÞ
ð12Þ
and for node j be Wj ¼ Wj ðx; cÞ ¼ Nj cF ðxÞ;
ð13Þ
where c is the upwinding coefficient (0 for no upwinding and 1 for fully upwinded). F ðxÞ is some positive function such that it is equal to 0 at the nodes and Z Xj 1 F ðxÞ dx ¼ l; 2 Xi
3.1. Overall methodology The one-dimensional kinematic wave overland flow problem governed by Eq. (1) with specified boundary and initial conditions was solved using constant and varying rainfall and Manning’s a for a certain mesh using different time-steps. For each time-step, the average percentage error between the analytical and numerical solutions was as: " !# Pm n 1 X j¼1 qaij qnij Pm error ¼ 100; ð19Þ p i¼1 j¼1 qaij where m and p are the number of sampling points in space and in time, respectively, qa and qn are the analytical and numerical solutions for flow, respectively. Singh [26] presented a general form of the analytical solution for rainfall and Manning’s a varying in space. Eqs. (20) and (21) represent the rising and constant parts of the hydrograph while Eqs. (22) and (23) are for the receding limb of the hydrograph,
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
1=b Z x 1 hðxÞ ¼ re ðxÞ dx ; aðxÞ 0 ð1bÞ=b 1=b Z x Z 1 x 1 re ðxÞ dx dx tðxÞ ¼ b 0 aðxÞ 0 for 0 6 t 6 tr
ð20Þ
ð21Þ
and 1=b Z x
1 re ðxÞ dx ; ð22Þ aðxÞ 0 Z x
ð1bÞ=b Z x 1=b 1 1
re ðxÞ dx dx tðx; x Þ ¼ tr þ b 0 aðxÞ x
hðx; x Þ ¼
for tr 6 t 6 tf ;
ð23Þ
where re is the input rainfall rate (or rainfall excess for infiltrating slopes), b ¼ 5=3, tr is the storm duration in hours, tf is the total simulation time, and x is the intersection of the characteristic curve passing through the origin with the t ¼ tr line. Using Eqs. (20) and (21), for a fixed value of x, a hydrograph may be developed to calculate depth and discharge of water as a function of time. Similarly the water profile of the slope for a given time can also be calculated from Eqs. (20) and (21). From Eq. (23), an expression for x can be generated in function of t and x and used in Eq. (22) to solve for the receding part of the hydrograph. This part requires an iterative solution because of the non-linearity of the resulting equation. A relationship between time-step and error is generated for a series of solutions using several time-steps. As the value of the time-step increased, the error is expected to grow, provided that the spatial discretization error is small. There exists a time-step value that integrates the problem within a specified error. Using a smaller timestep results in an unnecessarily longer computational time, while a larger time-step results in large errors and leads to spurious oscillations that could destroy the solution. The specified problem is solved using several mesh sizes generating a relationship between dimensionless mesh size and optimal time-step. The results are summarized in a regression equation that is used to define a time-step estimate. The time of concentration, which depends on the storm, the slope and the roughness factor of the watershed, is introduced in the equation as a problem specific parameter so the criteria can be used for other scenarios. The time-step estimate equations are tested using a different case to ensure the validity of the new dynamic time-steps. Jaber and Mohtar [16] developed and tested dynamic time-step equations for the finite element consistent scheme of the kinematic wave equation for overland flow. The development and the testing were done for different combinations of rainfall and surface characteristics constant and varying in time and space. The
431
dynamic time-step estimates Dt, for the consistent scheme for case 1: re and a constant, case 2: re varying and a constant, case 3: re constant and a varying and case 4: re and a varying, are shown in Table 1. In [16] L0 was represented as 1=N , where N is the number of elements. The same procedure is used in this paper. Dynamic time-steps will be developed for the most stable and accurate scheme and the results are compared to the dynamic time-steps of the conventional consistent Galerkin scheme of Table 1. 3.2. Characteristics of the scenarios used in the dynamic time-step development The scenarios used for the dynamic time-step development in this research are the same as the ones used by Jaber and Mohtar [16]. In order to illustrate the solution accuracy requirement in terms of time-steps and element sizes, nodal values of water depths were computed for 5, 10, 20, 25 50, and 75 elements for each case. Various time-steps were used to integrate the system of ordinary differential equations of Eq. (8). Each scenario was solved using the following conditions: contributing area of unit width L equal to 152.4 m (500 ft), average rainfall excess intensity equal to 2.74 cm/h (2.5E ) 5 ft/s), average bottom slope equal to 0.05 m/m, and average Manning’s roughness coefficient n equal to 0:035 m1=3 =s. For cases where re and a varied in space the following relations were used: re ðxÞ ¼ 0:036x þ 5:49;
ð24Þ
aðxÞ ¼ 6:015 e0:00984x ;
where re is in cm/h and a is in s=m1=3 . The time of concentration of each simulation was calculated in seconds according to the following equation [26]: ð1bÞ=b 1=b Z x Z 1 L 1 tc ¼ re ðxÞ dx dx; ð25Þ b 0 aðxÞ 0 where b ¼ 5=3. The storm duration is set to exceed the time to concentration so that pre-equilibrium, equilibrium and post-equilibrium conditions could be evaluated. The storm length during which rainfall is continuously falling is set to be 0.4 h. The total simulation time is 1 h and the time of concentration was determined to be 0.21 h. Table 1 Consistent scheme dynamic time-step equations for the four cases Constant excess rainfall Varying excess rainfall
Constant a
Varying a
Dt ¼ tc L1:16 0 Dt ¼ 0:81tc L1:08 0
Dt ¼ 0:96tc L1:27 0 Dt ¼ 1:66tc L1:43 0
L0 is the dimensionless mesh size l=L, where l is the length of element and L is the hillslope length, tc is the time of concentration; a is the Manning’s coefficient defined in Eq. (3).
432
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
3.3. Dynamic time-step evaluation scenarios
4.2. Fourier analysis
The evaluation problems used in this study were defined as: hillslope length L ¼ 144 m (472.5 ft), S0 ¼ 0:02 m/m, roughness n ¼ 0:009 m1=3 =s. The rainfall event simulated is re ¼ 19 cm/h (1.731E ) 4 ft/s) applied for 0.133 h. The above values of re and a were used in the cases where re and a were constant. For the spatially varying re and a the following equations were used:
A Fourier analysis is performed on the consistent, lumped and upwind schemes to analyze the stability of these schemes. The Fourier or the Von Neuman stability analysis is applied to Eq. (27). Singh [26] showed that in terms of errors, Eq. (27) reduces to
re ¼ 0:2736x þ 39:51;
ð26Þ
a ¼ 8:0903 e0:0079x :
This a accounts for a variation in Manning’s roughness coefficient n from a value of 0.02 to 0:04 m1=3 =s. It also allows for a variation in slope from 0.04 to 0.13 m/m. The time of concentration ðtc Þ was found to be 0.054 h for case 1 and 0.07 h for cases 2–4. The timestep was calculated using the generated equations for each scenario as appropriate.
4. Stability analysis Stability analysis is performed for the consistent, lumped and upwind schemes using eigenvalue and Fourier analysis. It is meant to provide insights into the accuracy and stability of the schemes. Also it is an attempt to provide analytical time-step criteria for the three schemes. Eq. (1) can be linearized as follows: oh oh þ l ¼ re ðx; tÞ; ot ox
ð27Þ
where l ¼ abhb1 a is defined in Eq. (3) and b ¼ 5=3. The ODE system of equations (Eq. (8)) can be rewritten as ½Afhnew g ¼ ½P fhold g þ ff g;
ð28Þ
where ½A ¼ ½C þ 0:5Dtl½B and ½P ¼ ½C 0:5Dtl½B. Multiplying by ½A1 , we obtain 1
1
fhnew g ¼ ½A ½P fhold g þ ½A ff g:
ð29Þ
4.1. Eigenvalue analysis Some numerical oscillations that may arise in the values of h from the finite element formulation are related to the eigenvalues of the matrix product ½A1 ½P in Eq. (29) [23]. The eigenvalues of the resulting product matrix were found to be positive for all x and Dt. This is a sufficient condition to conclude that no oscillations result from the spatial integration of Eq. (1) which leads us to believe that the oscillations that arise are from the finite difference integration in time.
oe oe þ l ¼ 0; ot ox
ð30Þ
where e is the error term. Using the consistent Galerkin finite element method Eq. (30) can be written as follows: Dx 2Dx Dx l l e_j1 þ e_j þ e_jþ1 ej1 þ ejþ1 ¼ 0; 6 3 6 2 2
ð31Þ
where Dx is the element length, e_ is the time derivative of flow depth error, j 1, j, and j þ 1 represent the previous, the present and the next nodes, respectively. Integrating in time using the Crank–Nicholson finite difference scheme Eq. (31) reduces to nþ1 n nþ1 nþ1 n 2 e e ej1 ej1 ejþ1 enjþ1 j j þ þ 3 6 6 h i Dt l 1 n nþ1 n enþ1 þ e e þ e : ð32Þ ¼ j1 j1 jþ1 jþ1 Dx 2 2 Let the Courant number ðCN Þ ¼ ðlDt=DxÞ and let enj ¼ eat eikm x , where a and km are wave numbers in time and space, respectively; n 1, n and n þ 1 represent the previous, present and next time-steps, respectively. Dividing Eq. (32) by enj and using the relations ia ia ia ia e þpeffiffiffiffiffiffi ffi ¼ 2 cos a and e e ¼ 2i sin a, where i ¼ 1, results in the following amplification factor for the consistent scheme: ¼ eaDt c
1 þ ð1=2Þ cos w CN i sin w ; 1 þ ð1=2Þ cos w þ CN i sin w
ð33Þ
where W ¼ km Dx. Anderson [2] defines a stable numerical scheme as a scheme where the error is not permitted to grow in the sequence of numerical procedures as the calculation proceeds from one marching step to the next. Furthermore, since ejnþ1 ¼ eaDt enj , it is clear that jeaDt j need to be less than or equal to 1 in order for the scheme to be stable. Similar calculations for the lumped scheme result in an amplification factor equal to: eaDt L ¼
1 CN i sin w : 1 þ CN i sin w
ð34Þ
The fully upwind lumped scheme amplification factor was found to be:
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438 aDt eU ¼
CN CN CN cos w þ i sin w 1þ 2 2 2
CN CN CN 1 þ cos w i sin w : 2 2 2
ð35Þ
The amplification factors can be written in the following form: eaDt ¼ Dðcos x i sin xÞ;
ð36Þ
where x is the phase angle and D is the algorithmic damping, which is equal to the absolute value of the complex amplification factor. The algorithmic damping is a real number and needs to be smaller than 1 for the scheme to be unconditionally stable. Algorithmic damping for the consistent and the lumped method are equal to unity, while for the upwind method it is 8( 2 ) < CN 2 CN 2 2 DU ¼ þ cos w ðCN sin wÞ þ 1 : 2 2 99 ,8" 2 2 #2 ==1=2 < CN CN CN cos w þ sin w : 1þ : ;; 2 2 2
433
when there is no damping in the amplification factor [17], while a strong sharp selective damping at low wavelengths can eliminate these oscillations. Fig. 2 shows that the upwind scheme results in strong sharp damping capable of reducing the oscillations for CNs equal to one but it is accompanied by deterioration in the accuracy of the solution. Figs. 5–7 show the hydrographs for the 0.1 upwind, the fully upwind and the no upwind schemes, respectively, using the time-steps that caused the first reported oscillations. It is clear from these figures that the oscillations in overland flow problems start at the time of concentration when the flow and the flow depth are at their peak. The CNs were calculated for the peak depth when the time-step results in oscillations for 75, 50, 25 and 20 elements systems for case 1 described in Section 3.3. The CNs were found to be equal to 3 for all cases. Fig. 2 shows that for CN ¼ 3 the damping for the upwind scheme is much weaker than the damping at CN ¼ 1:0. This reduction in damping effect and sharpness makes the upwind method inadequate to reduce oscillations.
Katapodes [17] showed that the behavior of the algorithmic damping of a certain scheme in terms of the dimensionless wavelength 2P=w can describe the dissipation level per wavelength. The dimensionless wavelength can be interpreted as the number of nodes per wavelength. Fig. 2 shows the behavior of algorithmic damping in terms of the dimensionless wavelength. The results are plotted for different CNs for the consistent lumped and fully upwind lumped schemes. It is important to mention that in the case of overland flow CN is a dynamic entity that is function of the depth h. Fig. 2 shows that the three methods are unconditionally stable since D is always smaller than one. The consistent and the lumped schemes have a constant algorithmic damping of one. These are known as no damping systems. In the presence of discontinuities most discrete methods generate spurious oscillations behind numerical jumps,
5. Assessment of the numerical schemes
Fig. 2. Algorithmic damping in function of the dimensionless wavelength for a Courant number equal to 1 and Courant number equal to 3.
Fig. 3. Error variation with time-step for the consistent – C, the lumped – L, the fully upwind consistent – CU, and the fully upwind lumped – LU schemes.
In order to assess the performance of the upwind and the lumped schemes performance, several preliminary runs were performed. Fractions and multiples of the dynamic time-steps developed by Jaber and Mohtar [16] for the consistent no upwind scheme, which are the upper limits for a non-oscillatory solution with an accuracy of 95%, were used for this evaluation. 5.1. Upwind schemes Upwind schemes using factors c ¼ 0:1 and c ¼ 1:0 were used to evaluate upwinding effect on the lumped and consistent finite element solutions. The same problem used to evaluate the dynamic time-steps in Section 3.3 is used to assess the upwind finite element schemes for the overland flow problem.
434
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
As shown in Fig. 3, the full upwind scheme (c ¼ 1:0) resulted in higher error than the consistent and lumped schemes for the same time-step. The solution is nonoscillatory for time-steps equal to the dynamic time-step but it became and inaccurate and oscillations were large enough to destroy the solution when the time-step used is 1.5 times the size of the dynamic time-steps for large mesh sizes. The upwind scheme did not reduce the oscillations of this problem. The c ¼ 0:1 upwind scheme had accurate results for time-steps equal or smaller than the dynamic time-step but did not result in any reduction in the oscillations of the solution for time-steps equal to 1.5 times the dynamic time-steps as shown in Fig. 4. Figs. 5–7 show that using the time-step that
Fig. 4. Error variation with time-step for the consistent – C, the lumped – L, the upwind 0.1 consistent – CU, and the upwind 0.1 lumped – LU schemes.
Fig. 5. Upwind c ¼ 0:1 hydrograph for case 1.
Fig. 7. Consistent no upwind c ¼ 0 hydrograph for case 1.
causes oscillations in the consistent scheme results in oscillations in the upwind formulations with c ¼ 0:1 and c ¼ 1:0. It also shows that there is an improvement in the oscillatory movement but it is not enough to provide an accurate solution. This confirms the results from the Fourier analysis. The damping for a CN ¼ 3 improved the general shape of the hydrograph of the upwind scheme (Fig. 6) compared to the non-upwind scheme (Fig. 7), but it was not strong enough to provide an accurate solution. These results indicate that the upwind scheme does not deal with oscillations that result in the one-dimensional kinematic wave problems for overland flow and is not recommended for such problems. 5.2. Lumped scheme Preliminary runs of the evaluation problems using the lumped scheme instead of the consistent scheme showed significant improvement in the results. The same level of accuracy compared to the consistent scheme was reached using larger time-steps. Figs. 3 and 4 show that the lumped schemes remain accurate within 5% of the analytical solution for time-steps much larger than the consistent scheme. Since the analytical methods described above do not provide any indication of how to select a time-step that will insure an accurate nonoscillating solution, the need to develop new dynamic time-steps for the lumped scheme empirically becomes evident.
6. Results of the dynamic time-step for the lumped scheme 6.1. Dynamic time-step criteria
Fig. 6. Fully upwind c ¼ 1:0 hydrograph for case 1.
Fig. 8 shows the selected time-steps that insure accuracy within 95% and stability for different mesh sizes for the four cases, respectively. The figure also shows the regression of these time-steps against the dimensionless mesh size. The time-step used in the graphs is a dimensionless time-step. It is the actual time-step divided by the time of concentration of the system. The resulting regression equations are shown in Table 2.
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
Fig. 8. Regression for developing dynamic time-steps for all cases of the lumped scheme.
435
Fig. 10. Evaluation of dynamic time-step for case 2 lumped scheme.
Table 2 Dynamic time-step equations for the four cases of the lumped scheme Constant a Constant excess rainfall Varying excess rainfal
Dt ¼ Dt ¼
0:89tc L0:96 0 0:79tc L0:86 0
Varying a Dt ¼ 1:27tc L1:16 0 Dt ¼ 0:79tc L0
L0 is the dimensionless mesh size l=L, where l is the length of element and L is the hillslope length, tc is the time of concentration; a is the Manning’s coefficient defined in Eq. (3).
6.2. Dynamic time-step evaluation Fig. 11. Evaluation of dynamic time-step for case 3 lumped scheme.
6.2.1. Steady rainfall problems The above-developed dynamic time-steps were evaluated using different problems than the development problems and proved to be adequate for solving kinematic wave problems for the lumped scheme. The use of half the size of the dynamic time-step did not result in a significant reduction in error, while the use of 1.5 times the size of the dynamic time-step resulted in highly oscillating solutions for the four cases. The results are shown in Figs. 9–12. Fig. 13 shows the change of the dynamic time-steps with the number of elements for both the consistent and lumped schemes for the four cases, respectively. It also shows that on the average the lumped scheme dynamic time-steps are double the size of the consistent scheme. Additionally, the lumped scheme solution remains stable for time-steps double the size of the ones used in the
Fig. 9. Evaluation of dynamic time-step for case 1 lumped scheme.
Fig. 12. Evaluation of dynamic time-step for case 4 lumped scheme.
Fig. 13. Comparison between dynamic time-steps of lumped and consistent schemes for all cases.
436
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
consistent scheme solution with the error remaining under 2%. The lumped scheme reduced the oscillations of the solution without affecting the accuracy. This result agrees with Yue [31] in his analysis of the threedimensional Navier–Stokes equation for open channel flow. Yue [31] explained the better performance of the lumped scheme by the fact that the lumped matrix is a sparse matrix with diagonals playing a major role thus eliminating numerical noise contributed by non-diagonal elements and thus providing a more stable solution. 6.2.2. Unsteady rainfall problems Jaber and Mohtar [16] have shown that the dynamic time-step for the consistent scheme can be applied to unsteady rainfall using the time of concentration equation presented by Parlange et al. [22] for unsteady rainfall kinematic wave problems. The dynamic timesteps were shown to be effective for an accurate nonoscillating solution. Another type of problems that occurs often is the rainfall event that is intermittent with bursts of rain followed by periods of no rain. Using the evaluation scenario of Section 3, a rainfall event consisting of 19 cm/h was simulated from 0 to 0.1 h, followed by no rain for 0.05 h, then another burst of rain equal to 9.5 cm/h was simulated from time equal to 0.15–0.3 h. The total simulation time was 0.5 h. Since the dynamic time-step is inversely proportional to the excess rainfall rate, the larger the rainfall the smaller the time-step required to provide stable and accurate solution. This problem was simulated using the time-step that will guarantee convergence and accuracy for the most critical storm, which is the first burst of rainfall of 19 cm/h (Fig. 14). Half of the dynamic time-step, the dynamic time-step and 1.5 times the time-step were used to simulate this problem. Fig. 15 shows that the error was less than 5% for half the dynamic time-step and the dynamic time-step while 1.5 times the dynamic time-step resulted in very high error and oscillations. Fig. 16 shows the hydrographs for the numerical and analytical solution of this problem for half the dynamic time-step and the dynamic time-step, respectively. The analytical solution matches the numerical solution in both cases except for a small error in
Fig. 14. Rainfall distribution for unsteady rain problem.
Fig. 15. Error percentage in the lumped formulation solution when using fractions and multiples of the dynamic time-step for unsteady rainfall.
Fig. 16. The analytical and numerical hydrographs solutions for the lumped formulation using half the dynamic time-step and the dynamic time-step for unsteady rainfall.
the rising part of the second storm. This shows that the use of the dynamic time-step will give approximately the same result as a time-step half its size. No hydrographs are shown for a time-step equal to one and a half times the dynamic time-step because the oscillations were large enough to destroy the solution.
7. Application and implementation The dynamic time-step developed in this study can be used as the upper limit to time-step size in solving the kinematic wave problem. An adaptive time-step scheme such as the one developed by Franca and Haghighi [8] can be integrated where the dynamic time-step is used for the equilibrium phase and larger time-steps can be used for the rising and the recessing limbs of the hydrograph. The dynamic time-step equations developed could also be integrated within any hydrologic computer programs utilizing the lumped scheme of the one-dimensional overland flow kinematic wave solution. These time-steps are functions of the grid size and the time to concentration, the latter being a function of rainfall, slope, hillslope length and Manning’s roughness. All these parameters are readily available and needed to run the hydrologic processes. The time-step can then be calculated automatically depending on the user-defined mesh size.
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
8. Conclusions A non-symmetric upwind scheme to replace the symmetric central difference treatment of the advective term of the kinematic wave equation was developed and evaluated. The upwind scheme did not eliminate the oscillations of both the lumped and the consistent scheme. It required smaller or equal time-steps compared to a non-upwinded scheme for the same accuracy level. Replacing the consistent scheme with a lumped matrix was also evaluated. The lumped scheme was found to significantly improve the solution without any reduction in solution accuracy. Stability analysis shows that the oscillations in overland flow problem starts at high Courant numbers when the upwind scheme damping is not effective. The improvement in the lumped scheme results is attributed to the sparse diagonal matrix that eliminates numerical noise from off-diagonal terms. The study provided a new dynamic time-step criteria for the solution of the one-dimensional overland flow problem for the lumped scheme. A different set of problems including an intermittent rain case was used to test and evaluate the dynamic time-steps. The dynamic time-steps were found to be good criteria to determine the size of the time-step that generate accurate and stable solution. The dynamic time-step equations could be integrated within hydrologic computer programs utilizing the lumped scheme of the one-dimensional overland flow kinematic wave solution. The time-step can be calculated automatically depending on userdefined mesh size. The dynamic time-step can be the upper limit for time-step size of an adaptive time-step finite element scheme for overland flow problems. The lumped scheme of the Galerkin finite element method for transient problems is the most suitable when solving the one-dimensional overland flow equation.
Acknowledgements The authors would like to thank Purdue University Agricultural Research Program office (ARP), the Environmental Sciences and Engineering Institute (ESEI) and the Purdue Research Foundation (PRF) for their financial support of this research.
References [1] Abbott MB, Basco DR. Computational fluid dynamics: an introduction for engineers. Essex, UK: Longman Scientific and Technical; 1989. [2] Anderson DA, Tannenhill JC, Pletcher RH. Computational fluid mechanics and heat transfer. New York: McGraw-Hill; 1984.
437
[3] Bajracharya K, Barry DA. Accuracy criteria for linearized diffusion wave routing. J Hydrol 1997;197:200–17. [4] Chan EKC, Williamson S. Factors influencing the need for upwinding in two-dimensional field calculation. IEEE Trans Magn 1992;28:1611–4. [5] Christie I, Griffiths DF, Mitchell AR, Zienkiewicz OC. Finite element methods for second order differential equations with significant first derivatives. Int J Numer Meth Engrg 1976;10: 1389–96. [6] Chow VT, Maidment DR, Mays LW. Applied hydrology. New York: McGraw-Hill; 1988. [7] Fiedler FR, Ramirez JA. A numerical method for simulating discontinuous shallow flow over an infiltrating surface. Int J Numer Meth Fluids 2000;32:219–40. [8] Franca AS, Haghighi K. Adaptive finite element analysis of transient thermal problems. Numer Heat Transfer, Part B 1994; 26:273–92. [9] Franco JL, Chaudhry FH. A comparison of various finite element procedures for watershed routing. In: Computational methods in water resources XII. Southampton, UK, Boston, USA: Computational Mechanics Publications; 1998. p. 537–44. [10] Gresho PM, Lee RL. Don’t suppress the wiggles – they’re telling you something!. Comput Fluids 1981;9(2):223–55. [11] Heatwole CD, Shanholtz VO, Ross BB. Finite element model to describe overland flow on an infiltrating watershed. Trans ASAE 1982;24(1):630–7. [12] Heinrich JC, Huyakom PS, Zienkiewicz OC, Mitchell AR. An upwind finite element scheme for two-dimensional convective transport equation. Int J Numer Meth Eng 1977;11:131–43. [13] Henderson FM. Open-channel flow. New York: Macmillan; 1966. [14] Hromadka TV, DeVries JJ. Kinematic wave routing and computational error. J Hydraul Eng, ASCE 1988;114(2):207–17. [15] Hughes TJR. A simple scheme for developing upwind finite elements. Int J Numer Meth Eng 1978;13:1359–65. [16] Jaber FH, Mohtar RH. Dynamic time step for the one-dimensional overland flow kinematic wave solution. J Hydrol Eng, ASCE 2002;7(1):3–11. [17] Katapodes ND. Fourier analysis of dissipative FEM channel flow model. J Hydraul Eng, ASCE 1984;110(7):927–44. [18] Leonard BP. A survey of finite differences of opinion on numerical muddling of the incomprehensible defective confusion equation. In: Hughes TJR, editor. Finite element methods for convection dominated flows, vol. AMD 39. New York: American Society of Mechanical Engineers; 1979. [19] Mohtar RH, Segerlind LJ. Dynamic time step estimates for two dimensional transient field problems using square elements. Int J Numer Meth Eng 1998;40:1–14. [20] Mohtar RH, Segerlind LJ. Dynamic time step estimates for onedimensional linear transient field problems. Trans ASAE 1999; 42(5):1477–84. [21] Mohtar RH, Segerlind LJ. Dynamic time step and stability criteria comparison. Int J Thermal Sci 1999;38:475–80. [22] Parlange JY, Rose CW, Sander G. Kinematic flow approximation of runoff on a plane: an exact analytical solution. J Hydrol 1981;52:171–6. [23] Segerlind LJ. Applied finite element analysis. 2nd ed. New York: Wiley; 1984. [24] Sharda VN, Singh SR. A finite element model for simulating runoff and soil erosion from mechanically treated agricultural lands. 1. Governing equations and solutions. Water Resour Res 1994;30(7):2287–98. [25] Singh VP, editor. Computer models in watershed hydrology. Colorado: Water Resources Publications; 1995. [26] Singh VP. Kinematic wave modeling in water resources: surfacewater hydrology. New York: Wiley; 1996. [27] Singh VP, Agiralioglu N. Diverging overland flow. Adv Water Resour 1981;4(3):117–24.
438
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
[28] Singh VP, Woolhiser DA. A nonlinear kinematic wave model for watershed surface runoff. J Hydrol 1976;31:221–43. [29] Tisdale TS, Scarlatos PD, Harick JM. Streamline upwind finite element method for overland flow. J Hydraul Eng, ASCE 1998;124:350–7. [30] Vieux BE, Segerlind LJ, Mohtar RH. Surface/subsurface flow equations interactions: identifying sources of groundwater con-
tamination. Report no. 89-G1569-02, Institute of Water Research, Michigan State University, East Lansing, MI; 1990. 45 p. [31] Yue J. Selective lumping effects on depth-integrated finite element model of channel flow. Adv Water Resour 1989;12(2): 74–8. [32] Zhang W, Cundy TW. Modeling of two-dimensional overland flow. Water Resour Res 1989;25(9):2019–35.
Stability and accuracy of finite element schemes for the one-dimensional kinematic wave solution Fouad H. Jaber, Rabi H. Mohtar
*
Agricultural and Biological Engineering Department, Purdue University, West Lafayette, IN 47907, USA Received 22 May 2001; received in revised form 2 November 2001; accepted 8 January 2002
Abstract Solving the kinematic wave equations for overland flow using the conventional consistent Galerkin finite element scheme is known to result in numerical oscillations due to the non-symmetric first spatial derivative terms in the kinematic wave equations. In this paper the lumped and the upwind finite element schemes are evaluated as alternatives to the consistent Galerkin finite element scheme. Stability analysis of the upwind scheme shows that the damping effect, that could reduce the oscillations, is small for the high Courant numbers encountered in overland flow problems. The upwind scheme, using upwind factors of 0.1 and 1.0, did not provide any improvement to the stability of the lumped and the consistent schemes. The lumped scheme considerably reduces oscillations without significant reduction in the overall solution accuracy. No analytical guidelines for time-step criteria that will insure stability and accuracy were provided by the stability analysis performed for the three schemes. Problem specific accuracybased dynamic time-step criteria was developed and evaluated for the lumped scheme. These time-steps were found to be on average, double the size of the dynamic time-steps for the consistent scheme. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Finite element analysis; Numerical stability; Dynamic time-step; Upwind scheme; Lumped scheme; Consistent scheme; Fourier analysis; Kinematic wave equations; Overland flow
1. Introduction 1.1. Numerical stability Numerous studies reported on solving the kinematic wave problems using finite element analysis [11,24,30]. The numerical solutions of these equations have encountered stability and convergence problems. Gresho and Lee [10] defined the conventional Galerkin finite element method as having (a) a centered difference like treatment of advection terms and (b) a consistent scheme for coupling time derivative terms. Christie et al. [5] stated that central difference approximations of first derivatives should be avoided since they give rise to numerical oscillations for reasonable grid sizes. They added that the Galerkin finite element method using symmetric pyramid basis linear shape functions reproduce central difference formulations with their inherent oscillatory properties. Hughes [15] reported that diffi*
Corresponding author. Tel.: +1-765-494-1791; fax: +1-765-4961115. E-mail address: [email protected] (R.H. Mohtar).
culties are encountered when the advection operators are non-symmetric and result in spurious oscillations. According to Fiedler and Ramirez [7], flow depths are generally small enough that small oscillations may destroy the solution of overland flow problems. They also added that the oscillation control is a formidable challenge given the oscillatory nature of the equations. Zhang and Cundy [32] used a MacCormak finite difference scheme, which is an explicit two-step finite difference scheme, to solve overland flow problems. They reported that numerical oscillations were evident in most hydrographs and in some cases the solution was aborted. Abbott and Basco [1] reported that the shallow water equations are known to be prone to high frequency Fourier components, 2Dx oscillations, when using central difference schemes because the value of the slope estimated by the difference is independent of the value of the dependent variable at the mid-point between two nodes. Chan and Williamson [4] reported that partial differential equations, containing first derivative terms in space, exhibit numerical oscillation in the Galerkin finite element solution if the size of the mesh exceeds a critical value based on the Peclet number being
0309-1708/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 9 - 1 7 0 8 ( 0 2 ) 0 0 0 0 5 - 2
428
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
smaller than two. The authors did not discuss cases where the Peclet number is infinitely large such as the case of overland flow problems. Nevertheless, an argument might be made for a fine finite element grid and small time-step to solve the problem, however, this will require unacceptable computational cost [29]. Some literature related to the numerical oscillation of the Galerkin solution of the overland flow problem illustrate the upwind methods as a possible solution for oscillations caused by the centered difference scheme. Chan and Williamson [4] described upwind schemes as replacing the central difference scheme by backward difference. They added that this is usually at the expense of accuracy. Heinrich et al. [12] reported that upwind finite element schemes consist of using weighting functions of non-symmetric forms different from those used as shape functions. They concluded that using the appropriate shape functions would result in an accuracy loss equivalent to the finite difference upwind schemes. Gresho and Lee [10] and Leonard [18] argued that the upwinding scheme by itself is a smoothing approach that while providing a stable solution can be very misleading with regards to the accuracy of the problem. Tisdale et al. [29] developed a streamline upwind scheme for a two-dimensional overland flow solution. The scheme consists of isolating the advective terms in the equation according to a main slope and solve the problem using upwind techniques. They reported that this scheme solved the problem without oscillations with negligible effect on accuracy. However, their study was based on a modification of the original equation and is limited to simple surfaces where the downslope node of each element can be determined using a subroutine within the program. This method is yet to be tested for complex runoff surfaces. Fiedler and Ramirez [7] used an upwind MacCormack finite difference scheme for solving overland flow problems. They found that in addition to upwinding, a smoothing function might be needed to control oscillations. The smoothing function that they propose is a function that smoothes regions of large gradients such as high frequency oscillations while leaving smooth areas undisturbed using a weighted average of flow depth along the nodes. Anderson et al. [2] explained that oscillations are due to the phase errors in the Fourier components of the numerical solution. They defined a stable numerical scheme as the one for which errors from any source are not permitted to grow in the sequence of numerical procedures as the calculation proceeds from one marching step to another. They also suggested the Fourier analysis to establish stability limits for linear partial differential equations with the possibility to extend these guidelines to non-linear equations in an approximate sense. Katopodes [17] used a Fourier analysis to study the stability and the oscillations of the linear kinematic wave equation for channel flow with constant
kinematic wave celerity and no lateral inflow. He stated that in the presence of discontinuities most discrete methods generate spurious oscillations behind numerical jumps, which can only be eliminated by selective algorithmic damping. The damping need to be very sharp and localized in the area of jumps because it is usually accompanied with a strong dispersion of celerity. Singh [26] showed that linear stability analysis could be applied to the linearized kinematic wave equations and still allows us to identify obviously unsuitable schemes or appropriate step length for conditionally stable schemes. The use of the lumped matrix to replace the consistent matrix has been reported to reduce oscillations in channel flow computations [31]. 1.2. Time-step criteria The task of choosing the proper time-step has often been considered a matter of experience. Bajracharya and Barry [3] stated that reducing the time-step does not necessarily result in more accurate solutions. They added that optimal solutions that would provide stability and accuracy with minimal computational time are yet to be found. Hromadka and DeVries [14] argued, after a series of tests varying Dx and Dt, that the use of the kinematic wave method for channel routing needs evaluation for use in hydrologic models unless guidelines are developed to control the arbitrary use of the kinematic wave in design studies. They added that kinematic wave programs need internal checks to select Dx and Dt such that an accurate solution is achieved. Jaber and Mohtar [16] developed empirical equations for dynamic time-steps [19–21] that will insure accuracy and stability of the consistent scheme for the kinematic wave solution. These dynamic time-steps vary with the number of elements and the time of concentration of the watershed. Their study was limited to the conventional Galerkin method using the consistent scheme. They reported significant oscillations for time-steps larger than the dynamic time-steps without suggesting ways to deal with these oscillations. A more stable scheme needs to be identified and tested and new time-step criteria need to be developed for this scheme. 1.3. Research objectives and limitations When using the kinematic wave equation, kinematic shock may arise where faster waves overtake slower ones, which leads to breaking of the wave profile. Singh [26] presented approximate methods in shock fitting to remedy this problem. It must be emphasized that the use of numerical techniques ignoring the existence of these shocks and not using a shock fitting technique is not a valid application of the kinematic wave theory. The problems used in this study to test different schemes and to develop and test dynamic time-steps did not result in
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
kinematic shocks. The kinematic wave model in the form presented in this paper is the basis of several commercial and research models available in the market that approximate watersheds as runoff planes. It is one of the methods available to calculate excess runoff in HEC-1, the Hydrological Simulation Program-Fortran (HSPF), the Institute of Hydrology Distributed Model (IHDM), the EPA Storm Water Management Model (SWMM) and the USDA-Water Erosion Prediction Project (WEPP) among others [25]. There has been some research to solve the kinematic wave equation for converging [28] and diverging watersheds [27]. In both cases additional terms are added to the conservation of mass equation and change it from its original form. Results that apply to the original kinematic wave equation do not necessarily apply to the modified equations. The overall objective of this paper is to improve the numerical solution of the one-dimensional kinematic wave equation for overland flow by addressing the two identified causes of oscillations: (a) the consistent mass matrix for coupling the time derivatives and (b) the central difference treatment of the advective term. The specific objectives include: 1. Evaluate the effectiveness of introducing non-symmetric upwind scheme to replace the symmetric central difference treatment of the advective term. 2. Evaluate the effectiveness of replacing the consistent scheme by a lumped matrix. 3. Develop and evaluate a dynamic time-step criteria for the most stable scheme identified from the first two objectives. 2. Numerical schemes for solving the kinematic wave equation A brief description of the kinematic wave equations for overland flow follows. Details and derivation of these equations are given by Chow et al. [6], and Henderson [13] among many others. The one-dimensional kinematic wave equation is governed by the continuity equation:
429
re ðx; tÞ is the excess rainfall rate (m/s); S0 is the element bottom slope in m/m; Sf is the friction slope in m/m; and x (m) and t (s) represent the space and time domains, respectively. Franco and Chaudhry [9] stated that interpolating h and q provides the most accurate and stable solution of this kinematic wave problem. The Galerkin finite element solution of this problem consists of expressing independent variables as: n X /ðeÞ ¼ Ni ðxÞUi ; ð4Þ i¼1
where Ui are known values of the continuous function /ðeÞ of each of the nodal points of an element (e), Ni ðxÞ are shape functions that approximate the function of /ðeÞ based on the nth nodal values [30]. The shape functions at nodes j 1 and j for one-dimensional elements are shown in Fig. 1 and are represented mathematically as: Nj1 ¼
Xj x Xj Xj1
and
Nj ¼
x Xj1 : Xj Xj1
ð5Þ
A description of the three numerical schemes considered in this paper is presented in the following section. 2.1. The consistent scheme In the consistent scheme oh=ot is assumed varying linearly between nodal values and is expressed as ( ) oUj1 oh ot ¼ Nj1 Nj : ð6Þ oUj ot ot The residual integral of the Galerkin finite element formulation is
Z Xj ðeÞ oh oq ð7Þ R ¼ ½N ðxÞT þ rðx; tÞ dx ¼ 0: ot ox Xi Integrating Eq. (7) and arranging the results yield the system of ordinary differential equations (ODE) of the form:
oh ½C þ ½Bfqg ff g ¼ 0; ð8Þ ot
oh oq þ ¼ re ðx; tÞ ð1Þ ot ox and the conservation of momentum equation for the kinematic wave assumption is reduced to: S0 ¼ Sf :
ð2Þ 2
The unit width flow (q) in m =s and the vertical flow depth (h) in m in Eq. (1) are related by the following energy equation: q ¼ ahb ;
ð3Þ 1=2 ðS0 =nÞ;
b ¼ 5=3 where, using Manning’s equation, a ¼ in which n is the Manning roughness coefficient (m1=3 =s);
Fig. 1. The weighting functions for the upwind c ¼ 1, c ¼ 0:1 and c ¼ 0.
430
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
where [C] and [B] and ff g are global matrices constructed using the direct stiffness method from the following element matrices, respectively: ðeÞ l 2 1 ðeÞ 1 1 1 C ; B ¼ ¼ 6 1 2 2 1 1
ðeÞ l rðx; tÞ and f ¼ ; ð9Þ 2 rðx; tÞ where l is the length of the element in m. 2.2. Lumped scheme The lumped scheme [23] is developed by defining the variation in oh=ot with respect to x to be constant between the midpoints of adjacent elements. The scheme uses alternative shape functions N for the time derivative defined as
Nj1 ¼ 1 H s 2l and Nj ¼ H s 2l ; ð10Þ where H ðxÞ (Heavyside function) is defined as H ðxÞ ¼ 0
if x < 0;
H ðxÞ ¼ 1
if x >¼ 0;
where l is the length of the element in m. The function that was used in this work and satisfies these conditions is (Fig. 1): F ðxÞ ¼ 3
ðx Xj Þðx Xj1 Þ ðXj Xj1 Þ2
:
ð14Þ
The residuals of the advective term in the kinematic wave equation oq=ox for one element are: Z Xj oN ðeÞ q dx Wj1 Rj1 ¼ ox Xj1 Z Xj oN ðeÞ and Rj ¼ q dx; ð15Þ Wj ox Xj1 solving these integrals for c ¼ 1:0 results in the following: ðeÞ
Rj1 ¼ 0; ðeÞ
Rj ¼ qj1 þ qj or in matrix form
Rj1 qj1 0 0 ¼ : Rj qj 1 1
ð16Þ
ð17Þ
s is the distance from node j 1. The element matrix [C] for the lumped scheme is: ðeÞ l 1 0 : ð11Þ C ¼ 2 0 1
For an upwind factor c ¼ 0:1 the results becomes
qj1 Rj1 0:45 0:55 ¼ : ð18Þ Rj qj 0:55 0:55
The [B] matrix remains the same as in the consistent method.
For c ¼ 0 the matrix reduces to the non-upwinded [B] matrix of Eq. (9).
2.3. Upwind finite element scheme 3. Methodology for dynamic time-step development An ‘‘upwind’’ finite element scheme to solve onedimensional kinematic wave equations is presented in this paper as a function of Xj1 and Xj , which are, respectively, the coordinates of the nodes of an element with nodes j 1 and j based on Heinrich et al. [12] who developed an upwind finite element scheme for convective transport equation. Element matrices for this scheme are developed in this paper. The one-dimensional upwind scheme of Heinrich et al. [12] requires that the advective term weighting function of the equation oq=ox for node j 1 be Wj1 ¼ Wj1 ðx; cÞ ¼ Nj1 þ cF ðxÞ
ð12Þ
and for node j be Wj ¼ Wj ðx; cÞ ¼ Nj cF ðxÞ;
ð13Þ
where c is the upwinding coefficient (0 for no upwinding and 1 for fully upwinded). F ðxÞ is some positive function such that it is equal to 0 at the nodes and Z Xj 1 F ðxÞ dx ¼ l; 2 Xi
3.1. Overall methodology The one-dimensional kinematic wave overland flow problem governed by Eq. (1) with specified boundary and initial conditions was solved using constant and varying rainfall and Manning’s a for a certain mesh using different time-steps. For each time-step, the average percentage error between the analytical and numerical solutions was as: " !# Pm n 1 X j¼1 qaij qnij Pm error ¼ 100; ð19Þ p i¼1 j¼1 qaij where m and p are the number of sampling points in space and in time, respectively, qa and qn are the analytical and numerical solutions for flow, respectively. Singh [26] presented a general form of the analytical solution for rainfall and Manning’s a varying in space. Eqs. (20) and (21) represent the rising and constant parts of the hydrograph while Eqs. (22) and (23) are for the receding limb of the hydrograph,
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
1=b Z x 1 hðxÞ ¼ re ðxÞ dx ; aðxÞ 0 ð1bÞ=b 1=b Z x Z 1 x 1 re ðxÞ dx dx tðxÞ ¼ b 0 aðxÞ 0 for 0 6 t 6 tr
ð20Þ
ð21Þ
and 1=b Z x
1 re ðxÞ dx ; ð22Þ aðxÞ 0 Z x
ð1bÞ=b Z x 1=b 1 1
re ðxÞ dx dx tðx; x Þ ¼ tr þ b 0 aðxÞ x
hðx; x Þ ¼
for tr 6 t 6 tf ;
ð23Þ
where re is the input rainfall rate (or rainfall excess for infiltrating slopes), b ¼ 5=3, tr is the storm duration in hours, tf is the total simulation time, and x is the intersection of the characteristic curve passing through the origin with the t ¼ tr line. Using Eqs. (20) and (21), for a fixed value of x, a hydrograph may be developed to calculate depth and discharge of water as a function of time. Similarly the water profile of the slope for a given time can also be calculated from Eqs. (20) and (21). From Eq. (23), an expression for x can be generated in function of t and x and used in Eq. (22) to solve for the receding part of the hydrograph. This part requires an iterative solution because of the non-linearity of the resulting equation. A relationship between time-step and error is generated for a series of solutions using several time-steps. As the value of the time-step increased, the error is expected to grow, provided that the spatial discretization error is small. There exists a time-step value that integrates the problem within a specified error. Using a smaller timestep results in an unnecessarily longer computational time, while a larger time-step results in large errors and leads to spurious oscillations that could destroy the solution. The specified problem is solved using several mesh sizes generating a relationship between dimensionless mesh size and optimal time-step. The results are summarized in a regression equation that is used to define a time-step estimate. The time of concentration, which depends on the storm, the slope and the roughness factor of the watershed, is introduced in the equation as a problem specific parameter so the criteria can be used for other scenarios. The time-step estimate equations are tested using a different case to ensure the validity of the new dynamic time-steps. Jaber and Mohtar [16] developed and tested dynamic time-step equations for the finite element consistent scheme of the kinematic wave equation for overland flow. The development and the testing were done for different combinations of rainfall and surface characteristics constant and varying in time and space. The
431
dynamic time-step estimates Dt, for the consistent scheme for case 1: re and a constant, case 2: re varying and a constant, case 3: re constant and a varying and case 4: re and a varying, are shown in Table 1. In [16] L0 was represented as 1=N , where N is the number of elements. The same procedure is used in this paper. Dynamic time-steps will be developed for the most stable and accurate scheme and the results are compared to the dynamic time-steps of the conventional consistent Galerkin scheme of Table 1. 3.2. Characteristics of the scenarios used in the dynamic time-step development The scenarios used for the dynamic time-step development in this research are the same as the ones used by Jaber and Mohtar [16]. In order to illustrate the solution accuracy requirement in terms of time-steps and element sizes, nodal values of water depths were computed for 5, 10, 20, 25 50, and 75 elements for each case. Various time-steps were used to integrate the system of ordinary differential equations of Eq. (8). Each scenario was solved using the following conditions: contributing area of unit width L equal to 152.4 m (500 ft), average rainfall excess intensity equal to 2.74 cm/h (2.5E ) 5 ft/s), average bottom slope equal to 0.05 m/m, and average Manning’s roughness coefficient n equal to 0:035 m1=3 =s. For cases where re and a varied in space the following relations were used: re ðxÞ ¼ 0:036x þ 5:49;
ð24Þ
aðxÞ ¼ 6:015 e0:00984x ;
where re is in cm/h and a is in s=m1=3 . The time of concentration of each simulation was calculated in seconds according to the following equation [26]: ð1bÞ=b 1=b Z x Z 1 L 1 tc ¼ re ðxÞ dx dx; ð25Þ b 0 aðxÞ 0 where b ¼ 5=3. The storm duration is set to exceed the time to concentration so that pre-equilibrium, equilibrium and post-equilibrium conditions could be evaluated. The storm length during which rainfall is continuously falling is set to be 0.4 h. The total simulation time is 1 h and the time of concentration was determined to be 0.21 h. Table 1 Consistent scheme dynamic time-step equations for the four cases Constant excess rainfall Varying excess rainfall
Constant a
Varying a
Dt ¼ tc L1:16 0 Dt ¼ 0:81tc L1:08 0
Dt ¼ 0:96tc L1:27 0 Dt ¼ 1:66tc L1:43 0
L0 is the dimensionless mesh size l=L, where l is the length of element and L is the hillslope length, tc is the time of concentration; a is the Manning’s coefficient defined in Eq. (3).
432
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
3.3. Dynamic time-step evaluation scenarios
4.2. Fourier analysis
The evaluation problems used in this study were defined as: hillslope length L ¼ 144 m (472.5 ft), S0 ¼ 0:02 m/m, roughness n ¼ 0:009 m1=3 =s. The rainfall event simulated is re ¼ 19 cm/h (1.731E ) 4 ft/s) applied for 0.133 h. The above values of re and a were used in the cases where re and a were constant. For the spatially varying re and a the following equations were used:
A Fourier analysis is performed on the consistent, lumped and upwind schemes to analyze the stability of these schemes. The Fourier or the Von Neuman stability analysis is applied to Eq. (27). Singh [26] showed that in terms of errors, Eq. (27) reduces to
re ¼ 0:2736x þ 39:51;
ð26Þ
a ¼ 8:0903 e0:0079x :
This a accounts for a variation in Manning’s roughness coefficient n from a value of 0.02 to 0:04 m1=3 =s. It also allows for a variation in slope from 0.04 to 0.13 m/m. The time of concentration ðtc Þ was found to be 0.054 h for case 1 and 0.07 h for cases 2–4. The timestep was calculated using the generated equations for each scenario as appropriate.
4. Stability analysis Stability analysis is performed for the consistent, lumped and upwind schemes using eigenvalue and Fourier analysis. It is meant to provide insights into the accuracy and stability of the schemes. Also it is an attempt to provide analytical time-step criteria for the three schemes. Eq. (1) can be linearized as follows: oh oh þ l ¼ re ðx; tÞ; ot ox
ð27Þ
where l ¼ abhb1 a is defined in Eq. (3) and b ¼ 5=3. The ODE system of equations (Eq. (8)) can be rewritten as ½Afhnew g ¼ ½P fhold g þ ff g;
ð28Þ
where ½A ¼ ½C þ 0:5Dtl½B and ½P ¼ ½C 0:5Dtl½B. Multiplying by ½A1 , we obtain 1
1
fhnew g ¼ ½A ½P fhold g þ ½A ff g:
ð29Þ
4.1. Eigenvalue analysis Some numerical oscillations that may arise in the values of h from the finite element formulation are related to the eigenvalues of the matrix product ½A1 ½P in Eq. (29) [23]. The eigenvalues of the resulting product matrix were found to be positive for all x and Dt. This is a sufficient condition to conclude that no oscillations result from the spatial integration of Eq. (1) which leads us to believe that the oscillations that arise are from the finite difference integration in time.
oe oe þ l ¼ 0; ot ox
ð30Þ
where e is the error term. Using the consistent Galerkin finite element method Eq. (30) can be written as follows: Dx 2Dx Dx l l e_j1 þ e_j þ e_jþ1 ej1 þ ejþ1 ¼ 0; 6 3 6 2 2
ð31Þ
where Dx is the element length, e_ is the time derivative of flow depth error, j 1, j, and j þ 1 represent the previous, the present and the next nodes, respectively. Integrating in time using the Crank–Nicholson finite difference scheme Eq. (31) reduces to nþ1 n nþ1 nþ1 n 2 e e ej1 ej1 ejþ1 enjþ1 j j þ þ 3 6 6 h i Dt l 1 n nþ1 n enþ1 þ e e þ e : ð32Þ ¼ j1 j1 jþ1 jþ1 Dx 2 2 Let the Courant number ðCN Þ ¼ ðlDt=DxÞ and let enj ¼ eat eikm x , where a and km are wave numbers in time and space, respectively; n 1, n and n þ 1 represent the previous, present and next time-steps, respectively. Dividing Eq. (32) by enj and using the relations ia ia ia ia e þpeffiffiffiffiffiffi ffi ¼ 2 cos a and e e ¼ 2i sin a, where i ¼ 1, results in the following amplification factor for the consistent scheme: ¼ eaDt c
1 þ ð1=2Þ cos w CN i sin w ; 1 þ ð1=2Þ cos w þ CN i sin w
ð33Þ
where W ¼ km Dx. Anderson [2] defines a stable numerical scheme as a scheme where the error is not permitted to grow in the sequence of numerical procedures as the calculation proceeds from one marching step to the next. Furthermore, since ejnþ1 ¼ eaDt enj , it is clear that jeaDt j need to be less than or equal to 1 in order for the scheme to be stable. Similar calculations for the lumped scheme result in an amplification factor equal to: eaDt L ¼
1 CN i sin w : 1 þ CN i sin w
ð34Þ
The fully upwind lumped scheme amplification factor was found to be:
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438 aDt eU ¼
CN CN CN cos w þ i sin w 1þ 2 2 2
CN CN CN 1 þ cos w i sin w : 2 2 2
ð35Þ
The amplification factors can be written in the following form: eaDt ¼ Dðcos x i sin xÞ;
ð36Þ
where x is the phase angle and D is the algorithmic damping, which is equal to the absolute value of the complex amplification factor. The algorithmic damping is a real number and needs to be smaller than 1 for the scheme to be unconditionally stable. Algorithmic damping for the consistent and the lumped method are equal to unity, while for the upwind method it is 8( 2 ) < CN 2 CN 2 2 DU ¼ þ cos w ðCN sin wÞ þ 1 : 2 2 99 ,8" 2 2 #2 ==1=2 < CN CN CN cos w þ sin w : 1þ : ;; 2 2 2
433
when there is no damping in the amplification factor [17], while a strong sharp selective damping at low wavelengths can eliminate these oscillations. Fig. 2 shows that the upwind scheme results in strong sharp damping capable of reducing the oscillations for CNs equal to one but it is accompanied by deterioration in the accuracy of the solution. Figs. 5–7 show the hydrographs for the 0.1 upwind, the fully upwind and the no upwind schemes, respectively, using the time-steps that caused the first reported oscillations. It is clear from these figures that the oscillations in overland flow problems start at the time of concentration when the flow and the flow depth are at their peak. The CNs were calculated for the peak depth when the time-step results in oscillations for 75, 50, 25 and 20 elements systems for case 1 described in Section 3.3. The CNs were found to be equal to 3 for all cases. Fig. 2 shows that for CN ¼ 3 the damping for the upwind scheme is much weaker than the damping at CN ¼ 1:0. This reduction in damping effect and sharpness makes the upwind method inadequate to reduce oscillations.
Katapodes [17] showed that the behavior of the algorithmic damping of a certain scheme in terms of the dimensionless wavelength 2P=w can describe the dissipation level per wavelength. The dimensionless wavelength can be interpreted as the number of nodes per wavelength. Fig. 2 shows the behavior of algorithmic damping in terms of the dimensionless wavelength. The results are plotted for different CNs for the consistent lumped and fully upwind lumped schemes. It is important to mention that in the case of overland flow CN is a dynamic entity that is function of the depth h. Fig. 2 shows that the three methods are unconditionally stable since D is always smaller than one. The consistent and the lumped schemes have a constant algorithmic damping of one. These are known as no damping systems. In the presence of discontinuities most discrete methods generate spurious oscillations behind numerical jumps,
5. Assessment of the numerical schemes
Fig. 2. Algorithmic damping in function of the dimensionless wavelength for a Courant number equal to 1 and Courant number equal to 3.
Fig. 3. Error variation with time-step for the consistent – C, the lumped – L, the fully upwind consistent – CU, and the fully upwind lumped – LU schemes.
In order to assess the performance of the upwind and the lumped schemes performance, several preliminary runs were performed. Fractions and multiples of the dynamic time-steps developed by Jaber and Mohtar [16] for the consistent no upwind scheme, which are the upper limits for a non-oscillatory solution with an accuracy of 95%, were used for this evaluation. 5.1. Upwind schemes Upwind schemes using factors c ¼ 0:1 and c ¼ 1:0 were used to evaluate upwinding effect on the lumped and consistent finite element solutions. The same problem used to evaluate the dynamic time-steps in Section 3.3 is used to assess the upwind finite element schemes for the overland flow problem.
434
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
As shown in Fig. 3, the full upwind scheme (c ¼ 1:0) resulted in higher error than the consistent and lumped schemes for the same time-step. The solution is nonoscillatory for time-steps equal to the dynamic time-step but it became and inaccurate and oscillations were large enough to destroy the solution when the time-step used is 1.5 times the size of the dynamic time-steps for large mesh sizes. The upwind scheme did not reduce the oscillations of this problem. The c ¼ 0:1 upwind scheme had accurate results for time-steps equal or smaller than the dynamic time-step but did not result in any reduction in the oscillations of the solution for time-steps equal to 1.5 times the dynamic time-steps as shown in Fig. 4. Figs. 5–7 show that using the time-step that
Fig. 4. Error variation with time-step for the consistent – C, the lumped – L, the upwind 0.1 consistent – CU, and the upwind 0.1 lumped – LU schemes.
Fig. 5. Upwind c ¼ 0:1 hydrograph for case 1.
Fig. 7. Consistent no upwind c ¼ 0 hydrograph for case 1.
causes oscillations in the consistent scheme results in oscillations in the upwind formulations with c ¼ 0:1 and c ¼ 1:0. It also shows that there is an improvement in the oscillatory movement but it is not enough to provide an accurate solution. This confirms the results from the Fourier analysis. The damping for a CN ¼ 3 improved the general shape of the hydrograph of the upwind scheme (Fig. 6) compared to the non-upwind scheme (Fig. 7), but it was not strong enough to provide an accurate solution. These results indicate that the upwind scheme does not deal with oscillations that result in the one-dimensional kinematic wave problems for overland flow and is not recommended for such problems. 5.2. Lumped scheme Preliminary runs of the evaluation problems using the lumped scheme instead of the consistent scheme showed significant improvement in the results. The same level of accuracy compared to the consistent scheme was reached using larger time-steps. Figs. 3 and 4 show that the lumped schemes remain accurate within 5% of the analytical solution for time-steps much larger than the consistent scheme. Since the analytical methods described above do not provide any indication of how to select a time-step that will insure an accurate nonoscillating solution, the need to develop new dynamic time-steps for the lumped scheme empirically becomes evident.
6. Results of the dynamic time-step for the lumped scheme 6.1. Dynamic time-step criteria
Fig. 6. Fully upwind c ¼ 1:0 hydrograph for case 1.
Fig. 8 shows the selected time-steps that insure accuracy within 95% and stability for different mesh sizes for the four cases, respectively. The figure also shows the regression of these time-steps against the dimensionless mesh size. The time-step used in the graphs is a dimensionless time-step. It is the actual time-step divided by the time of concentration of the system. The resulting regression equations are shown in Table 2.
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
Fig. 8. Regression for developing dynamic time-steps for all cases of the lumped scheme.
435
Fig. 10. Evaluation of dynamic time-step for case 2 lumped scheme.
Table 2 Dynamic time-step equations for the four cases of the lumped scheme Constant a Constant excess rainfall Varying excess rainfal
Dt ¼ Dt ¼
0:89tc L0:96 0 0:79tc L0:86 0
Varying a Dt ¼ 1:27tc L1:16 0 Dt ¼ 0:79tc L0
L0 is the dimensionless mesh size l=L, where l is the length of element and L is the hillslope length, tc is the time of concentration; a is the Manning’s coefficient defined in Eq. (3).
6.2. Dynamic time-step evaluation Fig. 11. Evaluation of dynamic time-step for case 3 lumped scheme.
6.2.1. Steady rainfall problems The above-developed dynamic time-steps were evaluated using different problems than the development problems and proved to be adequate for solving kinematic wave problems for the lumped scheme. The use of half the size of the dynamic time-step did not result in a significant reduction in error, while the use of 1.5 times the size of the dynamic time-step resulted in highly oscillating solutions for the four cases. The results are shown in Figs. 9–12. Fig. 13 shows the change of the dynamic time-steps with the number of elements for both the consistent and lumped schemes for the four cases, respectively. It also shows that on the average the lumped scheme dynamic time-steps are double the size of the consistent scheme. Additionally, the lumped scheme solution remains stable for time-steps double the size of the ones used in the
Fig. 9. Evaluation of dynamic time-step for case 1 lumped scheme.
Fig. 12. Evaluation of dynamic time-step for case 4 lumped scheme.
Fig. 13. Comparison between dynamic time-steps of lumped and consistent schemes for all cases.
436
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
consistent scheme solution with the error remaining under 2%. The lumped scheme reduced the oscillations of the solution without affecting the accuracy. This result agrees with Yue [31] in his analysis of the threedimensional Navier–Stokes equation for open channel flow. Yue [31] explained the better performance of the lumped scheme by the fact that the lumped matrix is a sparse matrix with diagonals playing a major role thus eliminating numerical noise contributed by non-diagonal elements and thus providing a more stable solution. 6.2.2. Unsteady rainfall problems Jaber and Mohtar [16] have shown that the dynamic time-step for the consistent scheme can be applied to unsteady rainfall using the time of concentration equation presented by Parlange et al. [22] for unsteady rainfall kinematic wave problems. The dynamic timesteps were shown to be effective for an accurate nonoscillating solution. Another type of problems that occurs often is the rainfall event that is intermittent with bursts of rain followed by periods of no rain. Using the evaluation scenario of Section 3, a rainfall event consisting of 19 cm/h was simulated from 0 to 0.1 h, followed by no rain for 0.05 h, then another burst of rain equal to 9.5 cm/h was simulated from time equal to 0.15–0.3 h. The total simulation time was 0.5 h. Since the dynamic time-step is inversely proportional to the excess rainfall rate, the larger the rainfall the smaller the time-step required to provide stable and accurate solution. This problem was simulated using the time-step that will guarantee convergence and accuracy for the most critical storm, which is the first burst of rainfall of 19 cm/h (Fig. 14). Half of the dynamic time-step, the dynamic time-step and 1.5 times the time-step were used to simulate this problem. Fig. 15 shows that the error was less than 5% for half the dynamic time-step and the dynamic time-step while 1.5 times the dynamic time-step resulted in very high error and oscillations. Fig. 16 shows the hydrographs for the numerical and analytical solution of this problem for half the dynamic time-step and the dynamic time-step, respectively. The analytical solution matches the numerical solution in both cases except for a small error in
Fig. 14. Rainfall distribution for unsteady rain problem.
Fig. 15. Error percentage in the lumped formulation solution when using fractions and multiples of the dynamic time-step for unsteady rainfall.
Fig. 16. The analytical and numerical hydrographs solutions for the lumped formulation using half the dynamic time-step and the dynamic time-step for unsteady rainfall.
the rising part of the second storm. This shows that the use of the dynamic time-step will give approximately the same result as a time-step half its size. No hydrographs are shown for a time-step equal to one and a half times the dynamic time-step because the oscillations were large enough to destroy the solution.
7. Application and implementation The dynamic time-step developed in this study can be used as the upper limit to time-step size in solving the kinematic wave problem. An adaptive time-step scheme such as the one developed by Franca and Haghighi [8] can be integrated where the dynamic time-step is used for the equilibrium phase and larger time-steps can be used for the rising and the recessing limbs of the hydrograph. The dynamic time-step equations developed could also be integrated within any hydrologic computer programs utilizing the lumped scheme of the one-dimensional overland flow kinematic wave solution. These time-steps are functions of the grid size and the time to concentration, the latter being a function of rainfall, slope, hillslope length and Manning’s roughness. All these parameters are readily available and needed to run the hydrologic processes. The time-step can then be calculated automatically depending on the user-defined mesh size.
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
8. Conclusions A non-symmetric upwind scheme to replace the symmetric central difference treatment of the advective term of the kinematic wave equation was developed and evaluated. The upwind scheme did not eliminate the oscillations of both the lumped and the consistent scheme. It required smaller or equal time-steps compared to a non-upwinded scheme for the same accuracy level. Replacing the consistent scheme with a lumped matrix was also evaluated. The lumped scheme was found to significantly improve the solution without any reduction in solution accuracy. Stability analysis shows that the oscillations in overland flow problem starts at high Courant numbers when the upwind scheme damping is not effective. The improvement in the lumped scheme results is attributed to the sparse diagonal matrix that eliminates numerical noise from off-diagonal terms. The study provided a new dynamic time-step criteria for the solution of the one-dimensional overland flow problem for the lumped scheme. A different set of problems including an intermittent rain case was used to test and evaluate the dynamic time-steps. The dynamic time-steps were found to be good criteria to determine the size of the time-step that generate accurate and stable solution. The dynamic time-step equations could be integrated within hydrologic computer programs utilizing the lumped scheme of the one-dimensional overland flow kinematic wave solution. The time-step can be calculated automatically depending on userdefined mesh size. The dynamic time-step can be the upper limit for time-step size of an adaptive time-step finite element scheme for overland flow problems. The lumped scheme of the Galerkin finite element method for transient problems is the most suitable when solving the one-dimensional overland flow equation.
Acknowledgements The authors would like to thank Purdue University Agricultural Research Program office (ARP), the Environmental Sciences and Engineering Institute (ESEI) and the Purdue Research Foundation (PRF) for their financial support of this research.
References [1] Abbott MB, Basco DR. Computational fluid dynamics: an introduction for engineers. Essex, UK: Longman Scientific and Technical; 1989. [2] Anderson DA, Tannenhill JC, Pletcher RH. Computational fluid mechanics and heat transfer. New York: McGraw-Hill; 1984.
437
[3] Bajracharya K, Barry DA. Accuracy criteria for linearized diffusion wave routing. J Hydrol 1997;197:200–17. [4] Chan EKC, Williamson S. Factors influencing the need for upwinding in two-dimensional field calculation. IEEE Trans Magn 1992;28:1611–4. [5] Christie I, Griffiths DF, Mitchell AR, Zienkiewicz OC. Finite element methods for second order differential equations with significant first derivatives. Int J Numer Meth Engrg 1976;10: 1389–96. [6] Chow VT, Maidment DR, Mays LW. Applied hydrology. New York: McGraw-Hill; 1988. [7] Fiedler FR, Ramirez JA. A numerical method for simulating discontinuous shallow flow over an infiltrating surface. Int J Numer Meth Fluids 2000;32:219–40. [8] Franca AS, Haghighi K. Adaptive finite element analysis of transient thermal problems. Numer Heat Transfer, Part B 1994; 26:273–92. [9] Franco JL, Chaudhry FH. A comparison of various finite element procedures for watershed routing. In: Computational methods in water resources XII. Southampton, UK, Boston, USA: Computational Mechanics Publications; 1998. p. 537–44. [10] Gresho PM, Lee RL. Don’t suppress the wiggles – they’re telling you something!. Comput Fluids 1981;9(2):223–55. [11] Heatwole CD, Shanholtz VO, Ross BB. Finite element model to describe overland flow on an infiltrating watershed. Trans ASAE 1982;24(1):630–7. [12] Heinrich JC, Huyakom PS, Zienkiewicz OC, Mitchell AR. An upwind finite element scheme for two-dimensional convective transport equation. Int J Numer Meth Eng 1977;11:131–43. [13] Henderson FM. Open-channel flow. New York: Macmillan; 1966. [14] Hromadka TV, DeVries JJ. Kinematic wave routing and computational error. J Hydraul Eng, ASCE 1988;114(2):207–17. [15] Hughes TJR. A simple scheme for developing upwind finite elements. Int J Numer Meth Eng 1978;13:1359–65. [16] Jaber FH, Mohtar RH. Dynamic time step for the one-dimensional overland flow kinematic wave solution. J Hydrol Eng, ASCE 2002;7(1):3–11. [17] Katapodes ND. Fourier analysis of dissipative FEM channel flow model. J Hydraul Eng, ASCE 1984;110(7):927–44. [18] Leonard BP. A survey of finite differences of opinion on numerical muddling of the incomprehensible defective confusion equation. In: Hughes TJR, editor. Finite element methods for convection dominated flows, vol. AMD 39. New York: American Society of Mechanical Engineers; 1979. [19] Mohtar RH, Segerlind LJ. Dynamic time step estimates for two dimensional transient field problems using square elements. Int J Numer Meth Eng 1998;40:1–14. [20] Mohtar RH, Segerlind LJ. Dynamic time step estimates for onedimensional linear transient field problems. Trans ASAE 1999; 42(5):1477–84. [21] Mohtar RH, Segerlind LJ. Dynamic time step and stability criteria comparison. Int J Thermal Sci 1999;38:475–80. [22] Parlange JY, Rose CW, Sander G. Kinematic flow approximation of runoff on a plane: an exact analytical solution. J Hydrol 1981;52:171–6. [23] Segerlind LJ. Applied finite element analysis. 2nd ed. New York: Wiley; 1984. [24] Sharda VN, Singh SR. A finite element model for simulating runoff and soil erosion from mechanically treated agricultural lands. 1. Governing equations and solutions. Water Resour Res 1994;30(7):2287–98. [25] Singh VP, editor. Computer models in watershed hydrology. Colorado: Water Resources Publications; 1995. [26] Singh VP. Kinematic wave modeling in water resources: surfacewater hydrology. New York: Wiley; 1996. [27] Singh VP, Agiralioglu N. Diverging overland flow. Adv Water Resour 1981;4(3):117–24.
438
F.H. Jaber, R.H. Mohtar / Advances in Water Resources 25 (2002) 427–438
[28] Singh VP, Woolhiser DA. A nonlinear kinematic wave model for watershed surface runoff. J Hydrol 1976;31:221–43. [29] Tisdale TS, Scarlatos PD, Harick JM. Streamline upwind finite element method for overland flow. J Hydraul Eng, ASCE 1998;124:350–7. [30] Vieux BE, Segerlind LJ, Mohtar RH. Surface/subsurface flow equations interactions: identifying sources of groundwater con-
tamination. Report no. 89-G1569-02, Institute of Water Research, Michigan State University, East Lansing, MI; 1990. 45 p. [31] Yue J. Selective lumping effects on depth-integrated finite element model of channel flow. Adv Water Resour 1989;12(2): 74–8. [32] Zhang W, Cundy TW. Modeling of two-dimensional overland flow. Water Resour Res 1989;25(9):2019–35.