SIC^2: Simulation and Integration of Control for Canals

Siphons

Pierre-Olivier Malaterre — 2014-10-22T16:58:49Z

It is possible to model siphons in the SIC software through the so-called "Preissmann slot " [1] [2] [3]. This method is applicable both for steady state calculations (Fluvia) or transient calculations (Sirene). The program for calculating the geometry (Talweg) automatically detects the sections which are closed at their upper end. This is possible for sections of width-elevation, abscissa-elevation, culvert (always closed), circular (if bank level high enough zbank >= zbed + 2*r) and trapezoidal (with negative slope and adhoc data) type. Anyway, internally, all these sections are converted into width-elevation format, and then the algorithms are the same for all these sections. The technique known as the Preissmann slot is to create an artificial slot at the top of the section from the ceiling to transform a closed pipe into an open free surface section. But in this artificial slot friction is not taken into account. To be more specific, the wetted perimeter P, section S, velocity V, hydraulic radius R_h=S/P used in the calculation of the pressure drop $J=\frac{Q^2}{(K S R_{h}^{2/3})^2}$ is limited to the actual section of the siphon. But elevation Z is calculated in the full section including the artificial slot.

In the calculation of the head H, the velocity V is taken in the actual section, and Z in the slot. The width-elevation tables are calculated by Talweg once for all and written in the xml file. When there is a Preissmann slot we can verify that P does not change in the slot. One can also see that S does not change either in the slot in the tables.

The advantage of this is that all calculations are the same (except for the trick on P, S, V, R_h and H) with or without slot in permanent or transitory flow conditions. The disadvantages are that the volume is stored in this slot, but it must be limited if the slot is thin.

We can show that the equations of Saint Venant pose problem with a width L=0, which one would have without the Preissmann slot from where precisely the idea of this trick.

On the other hand beware: one can have the impression that this is not the case, because when one looks at the wet section for example in the results, one does not find exactly what one would expect. For example, if we have a 9 m wide by 4 m high culvert we would expect to have a surface of 36 m². In fact we will have 36.0045 m² (36 + 9/2*0.001). This is due to the fact that the ceiling of the section is slightly modified during the treatment by Talweg, putting a section starting the Preissmann slot to 1mm more than the initial ceiling. Instead of having a flat horizontal ceiling we have a light roof hanging 1mm high, so with a small additional volume. This is to improve the numerical stability of the schema and to have a better defined $\frac{dL}{dZ}$ derivative. To better visualize this one has increased the writing format of S, P and L on the log file of F8.2 to F9.3 (cf Talweg.ans, messages 53, 54, 55 and 56). The wet perimeter P is calculated without this small modification because it only intervenes in the friction term.

The slot width must be chosen small (L=0.01 m in SIC by default, but it can be changed for each section in the section description or using the automatic modification tool) in order not to store too large volumes of water. In sewage systems, where there are manholes, some software uses a trick to reduce the section of those manholes to reduce the volume artificially stored in the slot to get a better volume conservation of the total volume of water into the network.

One can also choose this slot width L (using the automatic modification tool) having a wave speed $c=\sqrt{\frac{gS}{L}}$ close to that of sound waves in water which is about 1500 ms^-1.

This modification using the wave speed can be done also using the automatic modification tool. This can lead for example to a slot width in the order of 10^-4 m to 10^-6 m. With explicit numerical schemes this was a problem because of the very small time step required for the calculation in order to have a Courant number Cr< 1 with $C_r=\frac{(V+c)DT}{DX} <1$, where V is the velocity of the water, c the celerity of the waves, DT the computation time step and DX the space step between 2 calculation sections. But the implicit Preissmann scheme is not subject to the same constraints. It is unconditionally stable regardless of the Courant number. But a large numerical Cr (or also a very small one, which is not our case) introduces numerical bias (diffusion and phase shift, cf thesis Cunge). By default in SIC we choose a slot width of 0.01 m. For example in a culvert of 1 m², for a flow of 1 m³s^-1, we get a celerity c of the order of 30 ms^-1 with a slot of 0.01 m. For c of the order of 1500 ms^-1 we should take a slot width of approximately 0.000005 m. The advantage of a small slot is that the stored volume in the slot is limited. The drawback is that it could be less stable numerically. This is a compromise to find in case of difficulties.

Bibliographical References:

[1] Cunge J.A., Wegner M., 1964. "Intégration numérique des équations d'écoulement de Barré de Saint Venant par un schéma implicite de différences finies. Application au cas d'une galerie tantôt en charge, tantôt à surface libre". La Houille Blanche, n°1-1964, pp 33-39

[2] Vasconcelos J.G., Wright S.J., 2004. "Numerical modeling of the transition between free surface and pressurized flow in storm sewers". Innovative Modeling of Urban Water Systems, Monograph 12, Chap 10, W. James, Ed., pp 189-214

[3] Ukon T., Shigeta N., Watanabe M., Shiraishi H., Uotani M., 2008. "Correction methods for dropping of simulated water level utilising Preissmann and MOUSE slot models". 11th International Conference on Urban Drainage, Edinburgh, Scotland, UK, pp 1-9

Resolution of the linear system of the correcting equations

David Dorchies — 2012-02-01T16:45:00Z

Three unknown variables are attached to each reach n°i:

in terms of water elevation: $\Delta Z^i_1$ at the upstream node
in terms of water elevation: $\Delta Z^i_n$ at the downstream node
in terms of discharge: $\Delta Q^i$ anywhere in the reach

Let N be the number of reaches in the loop. It is possible to solve the linear system composed of the 3N equations. But an entire backwater curve computation will be necessary, after correction of the discharges, in order to get the water elevation in each cross section of the reaches.

Therefore, the directly interesting unknowns are the $\Delta Q^i$ variables.

If the equations are written in a judicious order, some (N.N) blocs appear:

$\pmatrix{A&B&0\cr D&I&F\cr G&H&I\cr}.\pmatrix{\Delta Q\cr \Delta Z_1\cr \Delta Z_n\cr} = \pmatrix{C1\cr C2\cr C3\cr}$

Matrix . Vector = Constant

That can be written: $M.V = K$ [8]

D, F are diagonal matrices
I is the unity matrix
H is an upper triangular matrix with a nil principal diagonal

The discharge continuity equations at the nodes which are not downstream nodes of a loop, are placed in the first block line of the linear system (filling the A block). We then complete the block line with the equality relations between the water elevations at the distributory nodes (filling the B and C1 blocks).

The condensation equations are placed in the second block line.

The downstream boundary conditions at the downstream reaches of the loop, and the equality relations between the upstream and the downstream water elevations of the reaches at the nodes that are not downstream nodes of a loop, are placed in the third block line.

By multiplying [8] with a proper transition matrix P:

$(P.M).V = P.K$

we get a N equations system in the $\Delta Q$ variable:

$A_0.\Delta Q = C_0$ [9]

with:
$\left\{ {\matrix{A_0 = A - B.D + \chi_0.(H.D - G)\cr C_0 = C_1 - \chi_0.C_3\cr \chi_0 = B.F.{\chi_1}^{-1}\cr \chi_1 = H.F - I\cr I = \mbox{ Unity matrix}\cr}}$

This method can only be applied if the $\chi_1$ matrix can be inverted. Fortunately this matrix can always be transformed into a triangular superior matrix with the diagonal only composed of 1. Therefore $\chi_1$ can easily be inverted through the Gauss - Jordan method. The system [9] is solved through the Gauss method with a partial Gauss elimination.

Correcting equations

David Dorchies — 2012-02-01T16:40:00Z

The corrected backwater curve will be solution of the equation:

$f(Z^{i}_{j} + \Delta Z^{i}_{j}, Z^{i}_{j+1} + \Delta Z^{i}_{j+1},Q^{i}_{j} + \Delta Q^{i}) = 0$ [2]

After a first order Taylor development of [2], we get:

[2] - [1] => $a^{i}_{j}.\Delta Z^{i}_{j} + b^{i}_{j}.\Delta Z^{i}_{j+1} + c^{i}_{j}.\Delta Q^{i} = 0 $ [3]

with:

$a_{j}^{i} = \frac{\partial f}{\partial Z^{i}_{j}}$
$b_{j}^{i} = \frac{\partial f}{\partial Z^{i}_{j+1}}$
$c_{j}^{i} = \frac{\partial f}{\partial Q^{i}}$

We know that aij cannot be nil, since [1] is solved through the Newton method using the ratio $\frac{f}{a^{i}_{j}$.

Let:

$d^{i}_{j} = \frac{b^{i}_{j}}{a^{i}_{j}}$ et $e^{i}_{j} = \frac{c^{i}_{j}}{a^{i}_{j}}$

For a reach i with n cross sections, the (n-1) equations [3] give:

$\Delta Z^{i}_{1} = d^{i}_{1}.\Delta Z^{i}_{2} + e^{i}_{1}.\Delta Q^{i} $

...

$\Delta Z^{i}_{j} = d^{i}_{1}.\Delta Z^{i}_{j+1} + e^{i}_{j}.\Delta Q^{i} $

...

$\Delta Z^{i}_{n-1} = d^{i}_{n-1}.\Delta Z^{i}_{n} + e^{i}_{n-1}.\Delta Q^{i} $

After condensation of these equations, we can get one relation representative of the total reach:

$\Delta Z^{i}_{1} = D^{i}_{n}.\Delta Z^{i}_{n} + E^{i}.\Delta Q^{i} $ [4]

with:

i = index of the reach
1 = index of the upstream section of the reach
n = index of the downstream section of the reach
$D^i = \prod_{j=1}^{n} d^i_j$
$E^i = D^i \prod_{k=1}^{n} \frac{e^i_k}{\prod_{j=k}^{n} d^i_j}$ ( with: $d^i_n = 1$ and $e^i_n = 0$ )

The other available relations are:

For all the nodes of the loop that are not downstream nodes of the loop:

Discharge continuity

$\sum_{i=1}^{k} \varepsilon ^i \Delta Q^i =0$ [5]

with:

k = number of reaches of the loop connected to the node
$\varepsilon ^i = 1$ upstream reach
$\varepsilon ^i = -1$ downstream reach

Water elevations equality

$Z^i_u + \Delta Z^i_u = Z^l_v + \Delta Z^l_v$ [6]

giving k - 1 equations, with:

u, v = 1 if the reach leaves the node and n if the reach joins the node
i, l = reaches connected to the node

For all the downstream nodes of the loop:

$\Delta Z^i_n = 0$ [7]

if the downstream boundary condition is a given elevation, or

$\Delta Z^i_n = g(\Delta Q^i) $

if the downstream boundary condition is a rating curve (in this case only one reach can join the corresponding downstream node).

Sort algorithm for the reaches of the loop before the steady flow computation - Initial solution

David Dorchies — 2012-02-01T16:33:00Z

An initial discharge value is chosen in each reach of the loop, compatible with the discharge continuity at the nodes, using for example a discharge division proportional to the mean cross area of the reach (reach volume at the overtopping bank elevation divided by its length).

Then, the backwater curve is computed in each reach, in sequence and following the above described order (Cf. Classification of reaches). The water elevation at a distributory node (boundary condition for the upstream reach) is taken as the average computed elevation of the downstream reaches starting from this node.

For a given portion of the reach n°i, defined by two sections, the Zⁱ_j elevation (in the upstream cross section n°j) is to be computed, knowing the Zⁱ_j+1 elevation (in the downstream cross section n°j+1) and the Qⁱ_j discharge flowing through the cross section n°j.

This is done by solving (numerically) the backwater equation that can be written as:

$f(Z^{i}_{j}, Z^{i}_{j+1},Q^{i}_{j}) = 0$ [1]

Classification of reaches

David Dorchies — 2012-02-01T16:33:00Z

Reaches are identified by their nodes. The position of a reach in the network is entirely defined by the names of its upstream and downstream nodes. The direction of flow is defined at the same time. The network topology is simply described as an oriented graph.

The reaches constitute the arcs of that graph, delineated by the nodes, upstream and downstream. They are automatically numbered by the program according to the order in which they are input in the data file.

Subcritical flow being controlled by the downstream conditions, the calculation of a water surface profile proceeds upwards, commencing at the downstream end.

Therefore, a relationship between water surface elevation and discharge is needed as a downstream boundary condition to start the calculation.

A hydraulic network can be represented by a connected oriented graph, without closed circuits, along which one can move with the help of the following algorithm:

One must look for a node where no arc arrives, and from where only one arc leaves.

The latter arc is added to a list, and is then removed from the graph.

The process is recursively applied on the resulting sub-graph, as long as nodes (as defined above) exist.

Then, either the problem is over, or arcs not added to the list still remain.

In the latter case, one searches in the remaining sub-graph a node where no arcs arrive, and from where more than one arc leaves.

Those are selected, then they are removed from the graph, and the process is recursively applied to the resulting sub-graph, as long as such nodes exist. Then, either the problem is over, or non-listed arcs still remain.

In the latter case, the algorithm is applied again, from the beginning, on the remaining sub-graph.

Reading backwards, the list of arcs thus obtained, one can find a satisfactory partial relationship (Remark: only the upstream nodes of the reaches are processed by this algorithm). From the list, one can also establish the identity of reaches of a given loop. These reaches are gathered starting from the first distributory point (the two first successive reaches with the same upstream node). If the model has several downstream boundary conditions in terms of elevations, then all the reaches following the first distributory point of the list belong to the loop. On the opposite, if the loop is being closed, it stops with the last junction (the last reach with the same downstream node of another reach belonging to the list and located after the beginning of the loop). We will illustrate the use of the algorithm on a simple example, as shown in figure 3.

Convention for the reaches of the nodes :

REACH > 0 if it starts from the node
REACH < 0 if it arrives at the node

We have:

DAM [1]
DC1 [2] [4] [-1]
DC2 [3] [6] [-2]
FC1 [5] [-4]
FC2 [-5]
FC3 [-6]

Apply the first criterion:

The first node where no arc arrives and from where only one arc starts is the upstream node.

Add reach [1] to the list, and remove it from the graph. Then one obtains:

NODES

DAM ~~[1]~~
DC1 [2] [4] ~~[-1]~~
DC2 [3] [6] [-2]
FC1 [5] [-4]
FC2 [-5]
FC3 [-6]
Liste : [1]

We cannot find any node being the upstream node of one single reach. We look for one node never being a downstream node but being the upstream node of several reaches. Node DC1 is fitting these conditions. Then the reaches number 2 and 4 are added to the list:

DC1 ~~[2]~~ ~~[4]~~ ~~[-1]~~
DC2 [3] [6] ~~[-2]~~
FC1 [5] ~~[-4]~~
Liste : [1], [2], [4]

The same algorithm is applied once again. We get the node named DC2 and the reaches number 3 and 6 are added to the list:

DC2 ~~[3]~~ ~~[6]~~ ~~[-2]~~
FC1 [5] ~~[-4]~~
Liste : [1], [2], [4], [3], [6]

We cannot find any other node with two downstream reaches. We then use the first algorithm. We get the node named FC1 and the reach number 5 is added to the list.

Liste : [1], [2], [4], [3], [6], [5]

The loop starts from the first distributory point (reaches 2 and 4, node DC1). All the following reaches of the list belong to the loop, since this loop is not closed (several downstream boundary conditions exist).

Interpolation of singular sections

David Dorchies — 2009-10-22T15:33:50Z

In the case of a singular section it is necessary to have two computational sections at the same abscissa. If the model user entered two data sections at the same abscissa, both these sections are retained as computational sections. If only one section was entered, the downstream computational section is interpolated using the singular data section and the data section immediately downstream with a 1 m step. The interpolated section is then placed at the same abscissa as the singular section (see figure 14).

Figure 14

Section k is interpolated between sections i and j with a 1m step and is then placed at the section i abscissa. Therefore, one has to take care to enter two sections (i and j) at the same abscissa, especially if the bed elevation upstream and downstream of the device is different or if the section dimensions are different.

Figure 15

In figure 15, one will enter the two cross sections (i and j), as shown, but j will be given the same abscissa as i, so that there will not be any interpolation.

Interpolation of computational sections

David Dorchies — 2009-10-22T15:33:50Z

If the distance between two data sections is more than one computational (space) step, the program interpolates computational sections between these two data sections, in accordance with the step decided by the user.

In reality, the space step is adjusted in order to give a whole number of equal computational intervals between the two sections considered. The interpolation is performed at constant water depth.

In the figure below, the interpolated widths in each data section that enable to obtain all the widths of the computational section are shown in broken lines.

A computational section could thus have up to double the number of points of the data sections. One must therefore adopt some criteria to eliminate some points and to avoid storing too many points.

Any point of the computational section which does not modify the section area by more than 1 % is eliminated. Similarly, all points which do not modify the wetted perimeter by more than 5 % are excluded. Therefore, whatever the water depth may be, we can expect a precision of 1 % regarding the estimation of the area and 5 % regarding the wetted perimeter.

The computational sections are completed vertically by a fictitious point located 500 m above the canal bank elevation in order to allow calculation even if overtopping occurs.

The computational sections are numbered within each reach, so that if the computational space step is modified in any given reach, the numbering of the sections within other reaches will remain unchanged.

Transformation of a data section into width-elevation format

David Dorchies — 2009-10-22T15:33:50Z

At any cross section the program looks for all the high and low points. It determines right and left banks and eliminates points outside the bed. Using the high and low points it divides the cross section into channels and for each channel it effects the transformation into width-elevation (see figure 12) format. Then for each elevation it adds the widths of all the channels.

Figure 12

Computational sections

David Dorchies — 2009-10-22T15:33:50Z

Data sections may be unequally distributed along the canal. In fact, the model user should select sections which best represent the canal dimensions, the changes of slopes and so on. Depending on the regularity of the canal, the spacing of the data sections may then be small or large.

For the hydraulic calculation, the spacing between computational sections should be such that a reliable estimation of the water surface profile is possible. This spacing is chosen by the model user depending on his knowledge of the canal hydraulic behavior (when the data sections are too far away, the model interpolates supplementary computational sections in order to allow a better simulation of the water surface profile).

All the entered data sections are retained as computational sections. Irrespective of the manner in which a data section was defined, the program transforms it into width-elevation data for storage and interpolation. Only groups of 4 characteristic values: elevation, width, wetted perimeter (if a section had been entered in terms of abscissa-elevation, the wetted perimeter would take into account the section asymmetry) and area are finally retained.

Parametric form

David Dorchies — 2009-10-22T15:33:50Z

Sections of special geometrical shape can be input in parametric form.

Circle: defined by its radius R, the bottom elevation ZF, and the bank elevation ZB. If ZB is not entered, ZB = ZF + 2R.

Figure 6

Culvert: defined by the bed width L, the side slope m (which may be negative) and by the bed elevation ZF and bank elevation ZB.

Figure 7

Power Relationship: defined by XL0, a, the bank elevation and the bed elevation.

L = XL0 * (Y/Y0)a

with:

Y0 = ZB - ZF
Y = Z - ZF water depth

The parabolic profile corresponds to a = 0.5

Figure 8

Rectangle: defined by its width, bed and bank elevations.

Figure 9

Trapezium or triangle: Defined by its bed width L (nil for the triangle) and the canal bed and bank elevations. Side slope may be negative in the trapezium case. In case of negative slope, if the section gets closed at a level below ZB, the program makes the necessary adjustment to the bank elevation.

Figure 10