UDK/UDC: 519.6:556.166:626/627 Prejeto/Received: 8. 1. 2015 Izvirni znanstveni članek - Original scientific paper Sprejeto/Accepted: 18. 4. 2015 A COMPARISON OF METHOD OF CHARACTERISTICS AND PREISSMANN SCHEME FOR FLOOD PROPAGATION MODELING WITH 1D SAINT-VENANT EQUATIONS PRIMERJAVA METODE KARAKTERISTIK IN PREISSMANNOVE SHEME POPLAVNIH VALOV Z 1D SAINT-VENANTOVIMI ENAČBAMI Nino Krvavica1, Vanja Travaš1^ 1 Faculty of Civil Engineering, University of Rijeka, Radmile Matejčić 3, Rijeka, Croatia Abstract The Saint-Venant equations were integrated by the explicit Method of Characteristics (MOC) and by the implicit Preissmann scheme to comparatively analyze and quantify the differences in prediction of flood wave propagation in open channels. For this purpose a hypothetical scenario was considered by defining a flood wave with a stage hydrograph at the inflow boundary of a prismatic channel. Downstream boundary was defined by a zero-gradient condition. The results are presented as stage hydrographs at equidistant sections along the channel. Comparative analysis revealed some differences between the compared methods. Discrepancies in peak wave heights are more evident in the upstream sections, while the downstream sections are more sensitive to differences in arrival times of flood wave peaks. The explicit MOC predicts lower wave heights and longer arrival times than the implicit Preissmann scheme. Keywords: Saint Venant equations, method of characteristics, Preissmann scheme, flood wave propagation. Izvleček Saint-Venantove enačbe smo povezali z eksplicitno Metodo karakteristik in implicitno Preissmannovo shemo za potrebe primerjalne analize in ocene razlik v napovedovanju širjenja poplavnih valov v odprtih kanalih. V ta namen smo opazovali hipotetičen scenarij, pri čemer smo poplavni val opredelili s faznim hidrogramom na dotočnem robu prizmatičnega kanala. Dolvodni rob je bil določen s pogojem ničnega gradienta. Rezultate študije predstavljajo fazni hidrogrami na ekvidistantnih odsekih vzdolž kanala. Primerjalna analiza odkriva nekatere razlike med obema metodama. Razlike v konicah poplavnih valov so bolj očitne v višje ležečih odsekih, medtem ko so spodnji odseki bolj izpostavljeni razlikam v času pojava konic poplavnih valov. Eksplicitna Metoda karakteristik napoveduje nižje konice poplavnih valov in kasnejše nastope konic kot implicitna Preissmannova shema. Ključne besede: Saint Venantove enačbe, metoda karakteristik, Preissmannova shema, širenje poplavnih valov. * Stik / Correspondence: vanja.travas@uniri.hr © Krvavica N., Travaš V.; Vsebina tega članka se sme uporabljati v skladu s pogoji licence Creative Commons Priznanje avtorstva -Nekomercialno - Deljenje pod enakimi pogoji 4.0. © Krvavica N., Travaš V.; This is an open-access article distributed under the terms of the Creative Commons Attribution - Non Commercial - Share Alike 4.0 Licence. Krvavica N., Travaš V.: A comparison of method of characteristics and Preissmann scheme for flood propagation modeling with 1D Saint-Venant equations - Primerjava metode karakteristik in Preissmannove sheme poplavnih valov z 1D Saint-Venantovimi enačbami Acta hydrotechnica 27/46 (2014), 1-12, Ljubljana 1. Introduction The importance of developing and studying numerical procedures for predicting geometrical, kinematical and dynamical properties of flood waves is obvious and arises from the need to predict their eventual destructive and invasive interaction with the environment. Accurate prediction becomes even more evident if a flood wave travels near or through an unprotected urban area. Numerical simulations can be used to anticipate the influence of flood wave propagation on water depths in open channels, which is crucial for planning and executing essential preventive actions. Several results of such analyses are given in terms of stage hydrographs i.e. function h(x,t) that defines the water depth h at position x at time t, and flow hydrographs i.e. function Q(x,t) that defines the spatial and temporal variations of the flow rate Q. However, note that the depended variables are usually assumed to be h(x,t) and v(x,t), since the function Q(x,t) can be expressed in terms of average flow velocity v(x,t) and depth h(x,t). The functions h(x,t) and v(x,t) must satisfy the given boundary conditions and the Saint-Venant system of differential equations, derived for the considered process under a series of usually valid assumptions (Chaudhry, 2008). The system follows from the principles of conservation of linear momentum and mass. If there are no lateral inflows along the channel and a prismatic channel with rectangular cross section area is considered, the principle of mass conservation can be expressed in differential form as (Chaudhry, 2008): dv dv dh , . Л + v~x + 3TX = 3(0-SE ), (2) h h — + v— + h— = 0. (1) On the other hand, the linear momentum equation is usually defined by considering the gravity force component in the direction of the flow, the pressure gradient and friction force. Although some other force contributions can also be included, this gives: where g is gravity acceleration (m/s2), S0 denotes the channel bed slope and SE denotes the slope of the energy gradient line. There are a variety of constitutive models that can be used to compute the energy loss S0 per unit length of channel and unit weight of fluid (Subramanya, 1997). However, in the present study the Manning equation is used (Chanson, 2004; Mohan Das, 2008): Su = n2v2R-4/3 (3) where n denotes the Manning roughness coefficient (s/m13), R denotes the hydraulic radius (m) defined as A/O, where O is the cross section perimeter (m) and A is the cross section area (m2). The latter is linearly depended on water depth h and the constant of proportionality is defined by the channel width B (m). For given initial conditions v0(x,t=0) and h0(x,t=0) and boundary conditions, equations (1), (2) and (3) can be used to compute the flow properties over spatial domain defined between points x=0 and x=L, where L denotes the channel length, and over temporal domain from t=0 up to a time of interest. However, since for a general case, the solution cannot be obtained in a closed form, different numerical techniques were developed to obtain an approximated solution (Crossley, 1999). In this paper, we focused on the Method of Characteristics (MOC) (Delphi, 2012) and the Preissmann scheme (Preissmann, 1961). A comparative analysis between the two methods was conducted and is presented hereafter. Since MOC is mainly viewed as an academic procedure (since it is difficult to adequately incorporate spatial variability) and the Preissmann scheme is widely used in practice and in several commercial software packages, the main motivation for the comparative analysis was the intention to quantify related differences, especially due to the absence of such analyses in the available literature (Shamaa, 2002; Akbari and Firoozi, 2010). Krvavica N., Travaš V.: A comparison of method of characteristics and Preissmann scheme for flood propagation modeling with 1D Saint-Venant equations - Primerjava metode karakteristik in Preissmannove sheme poplavnih valov z 1D Saint-Venantovimi enačbami Acta hydrotechnica 27/46 (2014), 1-12, Ljubljana To obtain a numerical approximation of the system defined by equations (1) and (2), the spatial and temporal domains are discretized into a finite number of points. The spatial domain is discretized by a finite number nx of computational cells with length Ax. The temporal domain is discretized by a finite number nt of time steps At. To make a discrete reference on the spatial and temporal position of relevant geometrical and mechanical quantities, the subscript i is introduced to indicate their spatial position and the superscript n is introduced to denote their temporal position. The same nomenclature is hereafter used for both methods. 2. Explicit time integration by MOC From the computational point of view, MOC is a very adequate numerical procedure that results in an efficient explicit algorithm which can be very easily implemented in any programming language. From the theoretical point of view, MOC can be used to transform a system of partial differential equations into a set of ordinary differential equations. However, note that this process is only valid for systems of hyperbolic partial differential equations such as the one that occurs in dam break problem. The essence of this method is to introduce a restriction which will specify a solution over predefined characteristics which can be geometrically interpreted as curves in the spacetime domain (usually approximated by a series of straight lines). The restriction is introduced by multiplying equation (1) with an unknown quantity (Lagrange multiplier), and adding it to equation (2). After solving those equations for the unknown multiplier, the system defined by (1) and (2) can be rewritten in a compact form for the positive C+ and negative C- characteristic as: 21 ± i.— = g(S0-SE), Dt~ с Dt ö v 0 bJ' (4) where c(x,t) denotes the wave celerity, D»/Dt is the material-time derivative and dx/dt denotes the slopes of the characteristic curve at some position in a space-time coordinate system in which the solutions of (4) are defined (Chanson, 2004). The resulting equation is based on the assumption that the celerity can be computed from the known depth by the gravity wave formula c(x,t)=(gh(x,t))05. Moreover, this equality can be used to further simplify (4) by expressing h(x,t) in terms of c(x,t). Since in this case h=c2/g, the derivative Dh/Dc is 2c/g and (4) can be rewritten as: D(v±2c) Dt = g(S0-SE ). (5) The material-time derivative D(v±2c)/Dt can be approximated by replacing the continuous spacetime domain by a grid generated with a finite number of points (Figure 1). The procedure is justified from the physical point of view if the change of the quantity (v±2c) over time interval At can be assumed to be linear. According to the previously introduced notations for spatial and temporal associations, the positive characteristic C+ illustrated in Figure 1 between points W and P can be computed as: dt At (6) and the negative characteristic C- between points E and P can be computed from: dx_ х?+1-хГ. dt At (7) The unknown values v and c at point P can be computed from the known quantities at points W and E from the previous time level n. However, note that the slopes defined by equations (6) and (7) are only valid if the flow velocity v and celerity c can be assumed as constants inside the time interval At, i.e. if the characteristics can be geometrically interpreted as straight lines (Figure 1). On the other hand, for unsteady flow conditions the quantity v±c from (5) will change in time (and space) and the start point of the characteristic curves that crosses point P will not coincide with the point W and E from the previous time step. Therefore, to compute v and c at point P, two simplifications should be introduced. The first is to X 1-1 Krvavica N., Travaš V.: A comparison of method of characteristics and Preissmann scheme for flood propagation modeling with 1D Saint-Venant equations - Primerjava metode karakteristik in Preissmannove sheme poplavnih valov z 1D Saint-Venantovimi enačbami Acta hydrotechnica 27/46 (2014), 1-12, Ljubljana assume that the adopted time step At is small enough so that the characteristics for unsteady flow can still be geometrically interpreted as straight lines. If this assumption is justified (which it usually is), the second simplification is introduced so that the slope of the characteristics can be computed. Namely, if the positive characteristic from Figure 1 is considered as an example, note that the slope at point L is unknown since the velocity vL and the celerity cL at the same position are unknown. On the other hand, the slope of the same characteristic when it passes through point P is also unknown. The problem is usually resolved by computing the slopes using the velocity vO and celerity cO, i.e. the relevant quantities at the same spatial location i as point P but from the previous time step. Accordantly, the positive characteristic is defined by AxL/At=vO+cO and the negative characteristic as AxR/At=vO-cO. For the considered case illustrated in Figure 1, the positive characteristic C+ from equation (5) can now be rewritten in discrete form as (Chanson, 2004): (vp+2cp)-(vL+2cL) At = g(S0-SE )b (8) and the negative characteristic C- can be defined ® initial conditions О boundary conditions О known • unknown © approximated t = t time marching procedure x = 0 spatial domain x = L (i -1, n) / (U n) у (i +1, n) -o- -O- W L O R as: Figure 1: Staggered space-time grid used in MOC. Slika 1: Premaknjena prostorsko-časovna mreža v metodi karakteristik. (vp-2cp)-(vR-2cR) At = — S E )r . (9) Note that the source term should also be interpolated. However, in order to retain simplicity, it can also be computed (i.e. approximated) from the known values at point L for the positive characteristic C+ (8) and at point R for the negative characteristic C- (9). Solving equations (8) and (9) simultaneously gives the unknown velocity vP and celerity cP at point P. This procedure produces an explicit algorithm that requires an additional approximation of the velocity and celerity at point L and R. However, since the relevant quantities at points W and O are known and the distance AxL is also known, a simple linear interpolation can be used for this purpose. The same procedure is adapted for negative characteristic C-. 2.1 Numerical stability Since the numerical algorithm is explicit in time, it is only conditionally stable. The essential requirement to obtain numerical stability is prescribed by the Courant condition. In the case for flood wave propagation, it can be formulated as (Huang and Song, 1985): At < Ax V0 + C0 (10) where S denotes the stabilization parameter defined between 0 and 1. Besides equation (10), to secure the numerical stability, it is also necessary to meet the almost forgotten Koren condition imposed by the inequality (Huang and Song, 1985): Krvavica N., Travaš V.: A comparison of method of characteristics and Preissmann scheme for flood propagation modeling with 1D Saint-Venant equations - Primerjava metode karakteristik in Preissmannove sheme poplavnih valov z 1D Saint-Venantovimi enačbami Acta hydrotechnica 27/46 (2014), 1-12, Ljubljana At < (/+4! - о (11) For a coordinate plane with axes Ax and At, the inequalities (10) and (11) define a spatial region in which the pairs (Ax,At) will ensure a stable numerical analysis. 2.2 Boundary conditions A stage hydrograph h(x=0,t) is specified at the upstream boundary condition at x=0 as a Dirichlet boundary condition. Numerical treatment of the flow rate Q(x=0,t) at the same cross section depends on the local Froude number Fr. Namely, for subcritical flow Fr<1 the flow rate Q(x=0,t) can be computed from the negative characteristic (9), for supercritical flow Fr>1 the flow rate Q(x=0,t) should also be specified in advance. At the downstream section the Neumann type of boundary condition is usually specified, in this case a zero gradient condition for both the water depth h and flow velocity v. 3. Implicit time integration by Preissmann scheme It is commonly accepted that amongst the more popular implicit finite difference schemes (e.g. Abbot-Ionescu scheme, Vasiliev scheme, Preissmann scheme etc.) the four point implicit scheme - or box scheme - is the most robust (Cunge et al., 1980; Szymkiewicz, 2010). Because of this it is implemented in several well-known commercial software packages such as HEC-RAS (US Army Corps of Engineers, 2010) or Bentley CivilStorm (Haestad Methods Solution Center, 2011). Preissmann scheme corresponds with the box scheme when the integration point P is in the middle of the spatial increment Ax (Szymkiewicz, 2010). With respect to the other implicit but also explicit integration methods, the main advantage of the Preissmann scheme is the ability to use variable spatial Axi and temporal increments Atn (Szymkiewicz, 2010). The spatial increment Ax denotes the distance between points i and i+1 and the temporal increment Atn denotes the time interval between the states at n and n+1. Preissmann scheme uses a non-staggered grid of points in the space-time coordinate system (Figure 2). Accordantly, a particular flow segment along the channel can be discretized more accurately (e.g. gradual change in cross section area) and the time discretization can be progressively adjusted with respect to the required temporal resolution of the results. Before considering the numerical procedure in more detail, it is opportune to rewrite the governing differential equations in terms of a flow rate Q and water stage H which is measured from a predefined datum level. Accordantly, the continuity equation, previously defined in equation (1), can be rewritten in terms of discharge as the depended variable (conservation form) as (Szymkiewicz, 2010): + t- =o. dt ox The momentum equation (2) takes the form: (12) dQ , д (ßq2\, . дн -T- + - ( — 1 + qA— = -qn dt dx\ A ) ö dx 0 2 Q\Q\ AR4/3 (13) where ß denotes the Boussinesq coefficient introduced to express the actual momentum in terms of average velocity v which is now contained in the depended variable Q (Szymkiewicz, 2010). The numerical procedure begins by specifying a computational grid (Figure 2). The unknown values H and Q at time n+1 are computed with respect to the time position of point P inside the time increment At (Figure 2). All the derivatives, functions or algebraic expressions are computed with respect to that point. Congruently, the value of an arbitrary function fP(x,t) at point P is approximated as (Preissmann, 1961): fP(i,n) ъ1^(f?+1 + fin) + \{f№ + ГГ1), Krvavica N., Travaš V.: A comparison of method of characteristics and Preissmann scheme for flood propagation modeling with 1D Saint-Venant equations - Primerjava metode karakteristik in Preissmannove sheme poplavnih valov z 1D Saint-Venantovimi enačbami Acta hydrotechnica 27/46 (2014), 1-12, Ljubljana where 9 is the weighting parameter ranging from 0 to 1 and determines how much 'weight' is attached to the values at time level n+1and how much to those at time level n. Note that for 0=1 the computation is fully implicit in time and for 0=0 it is fully explicit. Furthermore, to obtain a numerical approximation of the system defined by equations (12) and (13), time derivatives are approximated as: ŽL-lfU dt~ 2 V ■n+l.fU fU 1 + 1 +Ji+1 + li_ Atn Atn (15) and spatial derivatives as: д f / ?n+1 . fn+1\ 1Т~9(П+\1 ) дх V Axi J + (1 (16) Geometrical interpretation of the weighting parameter 0 is illustrated in Figure 2. With introduced approximations (15) and (16), the continuity equation (12) can be discretized in both time and space as: ,n+1 j+1 +"t+1 Atr, + ЧГ1+НЦ Atn ) + 0 (Qr++11+Qf+1) + (1 -Ö) = 0 (17) where the index P is introduced to denote that the width B is computed according to the formula given in (14). Note that for 0=0 the width BP can be explicitly computed and interpreted as the arithmetic mean of the width of neighbouring points in the computational grid. Similarly, the momentum equation defined in (13) can be discretized to obtain: 2 ( Atn Atn ) + ±. ((ßQ2)n+1 - (ßQ2) AXi\( A )i+1 ( A )i + i-9((ßQ2)n (ßQ2) AXi (( A )i+i ( A )i (18) +дАр(1 -0) (Я+1&) == У AR4/3)p> where the terms AP, BP and friction source term can also be computed according to the formula given in (14). Since the cross section area A is a function of water stage H, the only unknowns are the flow rate Q and water stage H in time n+1. Therefore, there are only four unknowns in equations (17) and (18). To start the computation, the state at time n is known from the initial conditions or from the previous time step. Such a pair of equations can be written for each spatial interval Ax. For a total number of nx computational cells along the channel, the procedure defines a set of 2(nx) algebraic equations with 2(nx+1) unknowns, representing the values for water stage H(x,t) and flow rate Q(x,t) inside the flow domain. The resulting two equations obtained by applying the Preissmann scheme are nonlinear, due to the convective term and the friction source term, so an iterative solution technique, such as the Newton-Raphson method (Kreyszig, 2006) is required. If parameter 0 is greater than zero, the computation results in a time marching procedure in which the step from one time line to the next is simultaneously performed for all points along the current time line (Figure 2). The non-staggered space-time grid is an attractive feature of the Preissmann scheme since it can readily be used with unequal spatial increments, which is particularly important for natural waterways where channel characteristics are highly variable even at short distances. Similarly, the applicability of unequal time steps is another important characteristic of the scheme, particularly for the case of hydrograph routing where floodwaters would generally rise relatively quickly and recess gradually in time. Krvavica N., Travaš V.: A comparison of method of characteristics and Preissmann scheme for flood propagation modeling with 1D Saint-Venant equations - Primerjava metode karakteristik in Preissmannove sheme poplavnih valov z 1D Saint-Venantovimi enačbami Acta hydrotechnica 27/46 (2014), 1-12, Ljubljana ® initial conditions O boundary conditions О known • unknown ® approximated t = t, time marching procedure = 0 spatial domain x = L Figure 2: Non-staggered space-time grid used in the Preissmann scheme. Slika 2: Nepremaknjena prostorsko-časovna mreža v Preissmanovi shemi. 3.1 Boundary conditions Since this procedure defines two unknowns for each grid point with a total of 2(nx+1) unknowns, and the total number of algebraic equations is 2(nx), there are four additional unknowns that should be determined from the boundary conditions. So, in addition to equations (17) and (18), a unique solution of the system is obtained by defining the upstream and downstream boundary conditions. Like in the previous case, only the water stage H(x=0,t) needs to be specified at the upstream cross section if Fr<1 and in the case of supercritical flow (i.e. if Fr>1), both water stage and flow rate should be specified at the same boundary. Similarly, the downstream boundary conditions are defined as zero gradient conditions for both water depth h and velocity v. 3.2 Numerical stability The Preissmann scheme is second-order accurate in both time and space if 6=0.5 and first-order accurate otherwise. Moreover, if applied to Saint-Venant system of equations, Lyn and Goodwin (1987) showed that it is unconditionally stable for 6>0.5 and neutrally stable if 6=0.5. For practical analysis it is common practice to set 6=0.55-0.65. However, it must be mentioned that in this case the scheme gives an approximation of first-order and thus generates a numerical diffusion that can affect the numerical solution. The magnitude of diffusion depends on the values of spatial and temporal increments and the weighting parameter 6 (Szymkiewicz, 2010). In other words, for 6>0.5 the scheme is computationally-diffusive. 4. Comparative analysis Major advantage of the Preissmann scheme, but also other implicit schemes, in comparison with MOC and other explicit algorithms is the fact that the numerical stability is ensured for 6>0.5 (Lyn and Goodwin, 1987). In other words, the Preissmann scheme is unconditionally stable without the need to satisfy the Courant condition (10) which defines the maximum allowable time step. In spite of the fact that implicit schemes require iterations at every time step, the requirement of explicit schemes to satisfy Courant conditions often makes them less attractive in terms of computational efficiency. However, if only computational efficiency is considered, it should also be mentioned that by increasing the number of grid points nx the situation can be reversed. Namely, by increasing nx the explicit time integration becomes more attractive, since the computational efficiency of explicit algorithm reduces linearly while efficiency of implicit algorithms reduces exponentially (since it involves a solution of a system of non-linear algebraic equations). t Krvavica N., Travaš V.: A comparison of method of characteristics and Preissmann scheme for flood propagation modeling with 1D Saint-Venant equations - Primerjava metode karakteristik in Preissmannove sheme poplavnih valov z 1D Saint-Venantovimi enačbami Acta hydrotechnica 27/46 (2014), 1-12, Ljubljana 4.1 Numerical examples To perform a comparative analysis, the considered numerical procedures were implemented in MATLAB (2011). Those numerical algorithms were used to simulate a hypothetical scenario in which a flood wave enters from the upstream boundary in a prismatic open channel with rectangular cross section. The considered channel is L=20 km long with constant width of B=30 m. Two different channel slopes where considered, moderately steep slope S0=0.05 for Test Cases 1 and 2, and very mild slope S0=0.005 for Test Cases 3 and 4. The Manning roughness coefficient was constant along the channel for all test cases and was defined as n=0.015 s/m13. The initial conditions were given as a velocity v0(x,t=0) and depth h0(x,t=0) along the flow domain (0