SIC^2: Simulation and Integration of Control for Canals

Numerical methods

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

In all the cases analyzed, we obtained an equation of the form f(Z)=0 or f(W)=0 to be solved.
In this section we examine in detail how those equations are solved numerically.

- Numerical methods

Recapitulation of the Newton's method

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

- Numerical methods

Resolution of the water surface profile equation

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

The water surface profile equation [11] is solved by Newton's method. The initial value is taken equal to:

Zio =

In our case, the function f(Z) does not have a continuous first derivative, because the geometry of the computational section is only known at certain points. In order to avoid any adverse effects on the computation convergence, it is limited by a bisection algorithm.

Starting bisection algorithm values are:

Zmin = Zci
Zmax = Hj + DH(Zci)

In a Newton iteration l two cases can occur:

1) The calculated value Zil is in the interval (Zmin, Zmax).

This value is retained
if f(Zil) > 0 then Zmin = Zil
if f(Zil) < 0 then Zmax = Zil

since f is a decreasing function.

2) The calculated value Zil is not in the interval (Zmin, Zmax). Replace Newton's iteration by a bisection algorithm:

Zil =

Similarly, if one detects a breakdown of the Newton algorithm due to the non continuity of the derivative of f(Zi), one also adopts the bisection algorithm.

Computation of the critical depth

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

We previously observed that it was necessary to calculate the critical depth in each computational section in order to determine if a sub-critical solution exists. The critical depth in a section is defined by:

= 1

This equation can be transformed as follows:

f(Zci) = Log ( ) = 0
which is a decreasing function.

which can be solved by Newton's method.

The initial value is determined by taking the critical water depth in the rectangle equivalent to the section corresponding to the overtopping water elevation (ZDi).

Zcio = ( )1/3 + ZFi

with: X =

Starting values for the bisection algorithm are:

Zmin = ZFi
Zmax = ZDM + 10

with: ZDM = elevation of the last point of the section.

These values are modified in the following way:

If f(Zcil) > 0 : Zmin = Zcil

If f(Zcil) < 0 : Zmax = Zcil

since f is a decreasing function.

Computation at a singular section

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

At a singular section, one can have as many weirs (five maximum) or gates (five maximum) as needed.

One must then solve equation [23] where the fk(Zi, Zj) terms are the discharges going through each device number k. Zi is the water surface elevation upstream of the device and Zj is the water surface elevation downstream of it. This equation is solved by Newton's method.

The initial value is determined as follows: Consider only one weir with length XL equal to the sum of the lengths of all the devices, and located at the elevation Zmin of the lowest one. If the downstream water elevation Zj is less than Zmin, one initializes with the free-flow weir formula:

Zio = Zmin + ( )2/3

If Zj is greater than Zmin, one considers:

Zio = Zj +

Starting values for the bisection algorithm are:

Zmin = elevation of the lowest device, or Zj if the latter value is greater.

Zmax = elevation of the highest device + 10 meters.

These values are modified as follows:

if f(Zil) > 0 : Zmax = Zil
if f(Zil) < 0 : Zmin = Zil

where f = f_k(Z_i, Z_j) - Q, which is a increasing function.

Let us examine the case of a regulator gate in a singular section. One has to solve the equation [24] which is of the type f(W) = 0, W being the regulator gate opening allowing to achieve the target upstream water elevation Zi. One starts calculating the discharges going through all the fixed gates, and one obtains:

fr(Zi, Zj, W) = Q - fk(Zi, Zj)

This equation is solved by bisection algorithm by adopting as the initial opening, half of the maximum authorized opening.

The calculation is stopped when the distance between the limits of the bisection algorithm is less than 1 mm.

Computation of the offtake opening

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

Equation [27] which is of the type:
fp(Z1, Z2, W) - Qp = 0
has to be solved. Z2 is calculated in the case when the offtake downstream condition varies:
Z2 = ZD + ( )2/3
or Z2 = ZD + (Zo-ZD) ( )1/n
One looks for the solution by bisection algorithm.
The initial value of W is taken equal to half of the authorized maximum opening. Calculation is stopped when the difference between the offtake discharge and the targeted discharge is less than (...)

- Numerical methods