Abstract
Time-dependent, unsteady flow has been studied in prismatic open channels with symmetric trapezoidal and triangular cross sections and small bottom slope. The St Venant equations without lateral inflow have been discretized in explicit as well as in implicit form and solved numerically, for unsteady, subcritical flow. The inflow hydrograph used can be applied for different flood events by adjusting its parameters accordingly. The results presented are derived from the explicit schemes Lax-Diffusive, MacCormack, Lambda as well as the implicit Preissmann scheme, and are compared to those from the Muskingum-Cunge method and the widely used commercial software HEC-RAS. The peak flow computed by the Lax-Diffusive scheme was reduced at the downstream end of the channel and the arrival time of the peak increased if compared to the other methods. The Muskingun-Cunge method forecasted the shortest peak flow arrival time at the downstream end cross section. Mass conservation computed from inflow and outflow hydrographs has been confirmed, since the maximum error did not exceed 2.60%. All codes were implemented in house using Matlab®.
INTRODUCTION
The above equations are nonlinear without analytical solutions, except for special cases where the linearized form of them can be applied. Earlier studies have dealt with the numerical solution of the equations either in the direct or in the characteristic form. The direct form of equations has been solved by finite difference and finite element methods, in explicit form (Liggett & Woolhiser 1967; Strelkoff 1970; Retsinis et al. 2018), or implicit form (Preissmann 1961; Amein 1968; Baltzer & Lai 1968; Amein & Fang 1970; Fread 1973; Amein & Chu 1975; Liggett & Cunge 1975; Delis et al. 2000; Krishnappan & Altinakar 2005; Retsinis et al. 2018).
The Muskingum-Cunge method can be classified as a hydraulic method based on the kinematic wave model with added diffusion terms, representing the real behavior of the flood wave propagation, i.e., attenuation of the peak flowrate, temporal translation of the peak as well as diffusion of the flood wave (Akan 2006). The rooting behavior of the Muskingum method is based on the kinematic wave but the attenuation of flood wave is due to the numerical diffusion. The analysis showed that the applied approximation modifies the kinematic wave model to the form of advection-diffusion equation that is not consistent (Szymkiewicz 2010). The resulting equation for estimating the hydrograph downstream is the same as the one produced by the Muskingum method (Akan 2006).
Equations (1) and (2) model the unsteady flow conditions in prismatic open channels where either added rainfall or lateral inflow can be neglected. The aim of this work is the numerical solution of the complete dynamic wave model under subcritical flow conditions in a long prismatic open channel of symmetric trapezoid and triangular cross sections. Three explicit numerical schemes will be used, namely, the Lax-Diffusive, the MacCormack, and the Lambda as well as the implicit four-point Preissmann scheme. The numerical results will be compared to those obtained by the Muskingum-Cunge method and the HEC-RAS (Warner et al. 2010) river modeling software.
NUMERICAL METHODS
The governing equations may be written either in conservative or non-conservative form. The conservative form should be preferred since it better conserves the transported flow variables and simulates more accurately the flood wave celerity (Chaudhry 2008).
is the depth of the center of gravity of the wetted cross section measured from free surface. The St Venant equations are of hyperbolic type, having two real characteristic directions where the positive characteristic direction is given by λ+ = V+c, and the negative one by λ− = V−c, where c=c(x,t) is the celerity of the propagating wave. In hyperbolic problems, the number of characteristic curves inserted in the computational domain determines the required number of initial and boundary conditions for a mathematically well posed problem (Chaudhry 2008). In the present work, where subcritical flow is examined, the characteristic directions are λ+ > 0 and λ− < 0, so that at the horizontal line crossing the initial time (t= 0), we need two initial conditions for the flow depth y and the cross-sectional average velocity V, and at the vertical lines passing through each end of the channel (x= 0, x=L, L is the total length of the channel), we only need one boundary condition since only one characteristic curve enters the computational domain. The above system of partial differential equations with the appropriate auxiliary conditions has to be solved numerically using the method of finite differences.A schematic representation of the channel's geometry with both the trapezoid and triangular cross section is shown in Figure 1. The first step is the discretization of the physical domain. To do this, we divided the total length of the channel L, by a number of nx nodes (including the end nodes) with index i denoting the space location, and the total simulation time T by nt nodes with index j referring to the time level, all equally spaced for simplicity, resulting in a uniform grid with Δx and Δt being the space and temporal intervals, respectively, as shown in Figure 2. Next, the approximation of the partial derivatives and the variables met in Equation (3) was done with suitable finite difference expressions showing consistency, stability, and convergence for better overall accuracy of the method. Finally, we either computed directly all the unknown variables, or in the implicit schemes we solved the resulting algebraic system of finite difference equations with direct or iterative techniques, thus evaluating the flow variables. The different numerical schemes implemented for this work are presented next.
Schematic channel's geometry.
Computational domain.
Lax-Diffusive scheme
From Equation (4) one can compute the values of A, Q=VA, at the j + 1 time level from the flow variables at time j level known either from the initial conditions or the previous time step calculations.
MacCormack scheme
This is a two-step predictor-corrector scheme that is second-order accurate in space and time. In this algorithm, backward finite differences are used for the spatial derivatives during the predictor part, and forward finite differences for the corrector part. This sequence of differentiation may be reversed in modifications of the method. Thus, we discretize Equation (3) as follows (Chaudhry 2008).
The computed value of 
gives the variables 
and 
from which we determine the values of 
, 
. These values are then used in the corrector part to compute 
and 
.
we receive (Chaudhry 2008):
at the j + 1 time level is given by:
Lambda scheme
Preissmann scheme
Equation (16) has four unknowns, the 
, 
, 
, 
, where 
is the free surface elevation measured from a horizontal datum at grid location i and time level j + 1. Writing the two equations for each reach, we end up with 2·(nx − 1) equations. However, we have 2·nx unknowns. The remaining two equations are provided from the boundary conditions, and then the system of nonlinear algebraic equations is solved simultaneously, providing the solution of unknown flow variables at j + 1 time stage.
In the present work, the iterative Newton-Raphson method was selected for the solution of the system (Akan 2006; Retsinis et al. 2018). We guess the values of the unknown flow variables 
, 
, 
, 
at each grid point and calculate the Jacobian matrix which has as elements the partial derivatives of Equation (16) and the boundary conditions with respect to 
, 
, 
, 
. The iterations stop once the corrections 
and 
for all grid points are within the limits we have set.
All the numerical codes discussed earlier have been implemented in house, using Matlab® computational environment.
Muskingum-Cunge method
, 
, 
, Qo is a reference discharge, So is the longitudinal bottom slope of the channel, L1 is the length of the channel reach, while all variables with subscript ‘o’ correspond to those calculated from the reference discharge and m= 5/3 (Akan 2006). The peak of the inflow hydrograph was selected as the reference discharge showing small sensitivity in the final results. To improve accuracy it is suggested that Δt must be smaller than one-fifth of the time interval from the beginning until the peak of the inflow hydrograph, and the length of the channel section for Muskingum-Cunge computations must be limited according to the following relationship (Akan 2006):
RESULTS
In the Muskingum-Cunge method, the selected reference discharge was 100 m3/s yielding a maximum reach length of 1,081 m calculated by Equation (19) for the trapezoid cross section, while Lmax = 1,082 m for the triangular cross section. The spatial and temporal discretization was chosen as L1 = 1,000 m and Δt= 200 s, respectively.
For the trapezoid cross section, the temporal evolution of flow rate, flow depth, and mean velocity is presented for all methods at downstream locations x= 0, x=L/4, x=L/2, x= 3L/4, and x=L, measured from the upstream end of the channel, in Figures 3–16. Also, Figures 17–20 show the spatial variation of flow depth with distance, measured from the upstream end of the channel for the Lax-Diffusive, MacCormack, Lambda, and Preissmann methods, respectively, while Figure 21 displays the variation of the maximum flow depth for all times computed by all methods except for the Muskingum-Cunge one, along the length of the channel.
Flow depth vs time at x= 0 m, trapezoid cross section.
Average velocity vs time at x= 0 m, trapezoid cross section.
Flow rate vs time at x= 2,000 m, trapezoid cross section.
Flow depth vs time at x= 2,000 m, trapezoid cross section.
Average velocity vs time at x= 2,000 m, trapezoid cross section.
Flow rate vs time at x= 4,000 m, trapezoid cross section.
Flow depth vs time at x= 4,000 m, trapezoid cross section.
Average velocity vs time at x= 4,000 m, trapezoid cross section.
Flow rate vs time at x= 6,000 m, trapezoid cross section.
Flow depth vs time at x= 6,000 m, trapezoid cross section.
Average velocity vs time at x= 6,000 m, trapezoid cross section.
Flow rate vs time at x= 8,000 m, trapezoid cross section.
Flow depth vs time at x= 8,000 m, trapezoid cross section.
Average velocity vs time at x= 8,000 m, trapezoid cross section.
Flow depth vs distance in Lax-Diffusive, trapezoid cross section.
Flow depth vs distance in MacCormack, trapezoid cross section.
Flow depth vs distance in Lambda, trapezoid cross section.
Flow depth vs distance in Preissmann, trapezoid cross section.
Maximum flow depth-distance, trapezoid cross section.
Regarding the triangular cross section, similar results for the flow rate, depth, and average cross-sectional velocity are presented at the downstream locations x=L/2, x= 3L/4, x=L, measured from the upstream end of the channel in Figures 22–30. Figure 31 depicts the maximum value of flow depth for all times computed by all methods except for the Muskingum-Cunge one along the length of the channel.
Flow rate-time at x= 4,000 m, triangular cross section.
Flow depth vs time at x= 4,000 m, triangular cross section.
Average velocity vs time at x= 4,000 m, triangular cross section.
Flow rate vs time at x= 6,000 m, triangular cross section.
Flow depth vs time at x= 6,000 m, triangular cross section.
Average velocity vs time at x= 6,000 m, triangular cross section.
Flow rate vs time at x= 8,000 m, triangular cross section.
Flow depth vs time at x= 8,000 m, triangular cross section.
Average velocity vs time at x= 8,000 m, triangular cross section.
Maximum flow depth vs distance, triangular cross section.
For the trapezoid cross section, the flow rates computed from the MacCormack, Lambda, Preissmann, and HEC-RAS methods are almost identical, while Lax-Diffusive scheme underestimates the peak discharge. The Muskingum-Cunge method overestimates the peak discharge at all downstream locations. Both Lax-Diffusive and Muskingum-Cunge showed steeper increase in the discharge compared to the other methods. All methods except for the Lax-Diffusive one compute a small reduction of the peak discharge at the downstream end of the channel. The Lax-Diffusive computed peak is 103.2 m3/s due to friction losses, while the arrival time of the peak discharge at the downstream end for the Lax-Diffusive, MacCormack, Lambda, Preissmann, Muskingum-Cunge, and HEC-RAS method was 3.84 hr, 3.37 hr, 3.39 hr, 3.40 hr, 3.11 hr, and 3.39 hr, respectively. It is obvious that the Lax-Diffusive scheme predicts the latest arrival time of maximum flow rate, while the Muskingun-Cunge method forecasts the earliest time. The results from MacCormack, Lambda, and Preissmann schemes as well as HEC-RAS compare well regarding the flow depth, while the results from Lax-Diffusive deviate slightly. Similar conclusions can be drawn for the average cross-sectional velocity where the results from Lax-Diffusive are different from those of all other methods, while this difference is reduced further downstream. The mass conservation equation using numerical integration with trapezoid rule for the computed hydrographs downstream was confirmed. More specifically, the maximum percentage error relative to the water mass entering the upstream cross-section for Lax-Diffusive, MacCormack, Lambda, and Preissmann schemes was 2.59, 1.41, 1.64, and 1.58, respectively.
Regarding the triangular cross section, similar conclusions are drawn for the flow rate, flow depth, and the average cross-sectional velocity. The arrival times of the peak discharge at the downstream end for the Lax-Diffusive, MacCormack, Lambda, Preissmann, Muskingum-Cunge, and HEC-RAS methods are 3.85 hr, 3.37 hr, 3.38 hr, 3.39 hr, 3.11 hr, and 3.40 hr, respectively. Again, the water mass conservation was confirmed while the maximum percentage error in the case of Lax-Diffusive, MacCormack, Lambda, and Preissmann scheme was 2.61, 1.42, 1.65, and 1.44, respectively, if compared to the water mass inserted in the upstream end of the channel.
CONCLUSIONS
The flood wave propagation in long prismatic open channels of symmetric trapezoid and triangular cross sections under subcritical flow conditions has been studied. The results from four numerical methods, three explicit ones, namely, the Lax-Diffusive, MacCormack, and Lambda methods, as well as the implicit Preissmann scheme, are compared to those from the hydraulic Muskingum-Cunge method and the HEC-RAS commercial software program. All codes were implemented in Matlab environment. Due to the explicit character of Lax-Diffusive, MacCormack, and Lambda schemes, the maximum time interval was restricted according to the Courant criterion, whereas in contrast, the Preissmann scheme was unconditionally stable. The weighting factor selected for accuracy reasons was 0.5, the convergence criteria regarded the unknown discharge allowed error of 0.001 m3/s, and the linear system of algebraic equations was solved by a direct method. The structure of the implicit code was improved resulting in a low number of iterations per time level, which led to accelerated computational times. The numerical results have been compared to those obtained from the commercial HEC-RAS software and the hydraulic/hydrologic Muskingum-Cunge method.
The largest error in mass conservation equation that appeared in the Lax-Diffusive method is probably due to the fact that it is a one-step algorithm, if compared to the two-step schemes. All the results from the numerical methods, the Muskingum-Cunge method and HEC-RAS software compare well.
The major innovation for this work was the comparison of different explicit or implicit numerical schemes with the Muskingum-Cunge hydraulic/hydrologic method and the HEC-RAS software. Furthermore, great effort was made for the coding of all numerical schemes in suitable programming language, revealing practical aspects and difficulties of code implementation useful for future works. Finally it proved that the explicit schemes produced satisfactory results if compared with the results of the more demanding implicit schemes.
To conclude, the numerical results compare well to those from use of the Muskingum-Cunge method and one-dimensional unsteady flow calculations software HEC-RAS, requiring more effort than the simpler Muskingum-Cunge method which can be used for the preliminary design of long prismatic open channels under unsteady flow conditions.
