SIC^2 : Logiciel de Simulation Intégrée des Canaux et de leur Contrôle

Siphons

Pierre-Olivier Malaterre — 2014-10-22T09:59:56Z

Il est possible de modéliser des siphons dans le logiciel SIC grâce à la technique dite de la "fente de Preissmann" [1] [2] [3]. Cette méthode s'applique que ce soit pour les calculs en régime permanent (Fluvia) ou en régime transitoire (Sirène). Le programme de calcul de la géométrie (Talweg) détecte automatiquement les sections qui sont fermées à leur extrémité supérieure. Cela est possible pour les sections de type largeur-cote, abscisse-cote, dalot (toujours), circulaire (si la cote de berge est supérieure ou égale à la cote de fond + 2 fois le rayon) et trapézoïdales (avec un fruit négatif et données adhoc si la hauteur est suffisante pour que le plafond se referme). De toute façon, de manière interne, toutes ces sections sont transformées en format largeur-cote, et les algorithmes sont ensuite les mêmes pour toutes ces sections.

La technique dite de la fente de Preissmann consiste à créer une fente artificielle en haut de la section à partir de son plafond, pour transformer une section fermée en charge en une section ouverte à surface libre. Mais dans cette fente artificielle on ne tient pas compte de frottement. Pour être plus précis, le périmètre mouillé P, la surface S, la vitesse V et/ou le rayon hydraulique R_h=S/P utilisés dans le calcul de la perte de charge J sont limités à la section réelle du siphon (sans la fente) avec $J=\frac{Q^2}{(K S R_{h}^{2/3})^2}$

Mais la cote Z est calculée dans la section complète comprenant la fente artificielle.

Dans le calcul de la charge H on prend la vitesse V limitée à la section réelle, et le Z dans la fente. Les tables largeur-cote sont calculées par Talweg une fois pour toute et écrites dans le fichier xml. Ces tableaux sont également indiqués dans le fichier log généré par Talweg. Quand il y a une fente de Preissmann on peut vérifier que P ne change pas dans la fente. On peut aussi voir que S ne change pas non plus dans la fente dans ces tables.

L'avantage de cela c'est que tous les calculs sont les mêmes (hormis l'astuce sur P, S, V, R_h et H) avec ou sans fente, en permanent ou en transitoire. Les inconvénients sont que du volume est stocké dans cette fente, mais cela doit être limité si la fente est mince.

On peut montrer que les équations de Saint Venant posent problème avec une largeur L=0, ce que l'on aurait sans la fente de Preissmann d'où justement l'idée de cette astuce.

Par contre attention, on peut avoir l'impression que cela n'est pas le cas, car quand on regarde la section mouillée par exemple dans les résultats, on ne trouve pas exactement ce à quoi on s'attendrait. Par exemple si on a un dalot de 9 m de large par 4 m de haut on s'attendrait à avoir une surface de 36 m². En fait on aura 36.0045 m² (36 + 9/2*0.001). Cela est dû au fait que le plafond de la section est légèrement modifié lors du traitement par talweg, en mettant une section débutant la fente de Preissmann à 1 mm plus faut que le plafond initial. Au lieu d'avoir un plafond plat horizontal on a donc un léger toit pendu de hauteur 1mm, avec donc un petit volume en plus. Cela a pour but d'améliorer la stabilité numérique du schéma et d'avoir une dérivée $\frac{dL}{dZ}$ mieux définie. Pour mieux visualiser cela on a augmenté le format d'écriture de S, P et L sur le fichier log de F8.2 à F9.3 (cf Talweg.ans, messages 53, 54, 55 et 56). Le périmètre mouillé P est quand à lui calculé sans cette petite modification car il n'intervient que dans le terme de frottement.

La fente doit être choisie de largeur petite (L=0.01 m par défaut dans SIC, mais on peut la changer section par section dans les options de la section, ou globalement pour une zone sélectionnée avec l'option "baguette magique" aussi appelée "construction automatique") afin de ne pas y stocker des volumes trop importants. Dans les réseaux d'assainissement, où il y a des regards, certains logiciels utilisent une astuce pour réduire la section de ces regards pour compenser le volume stocké artificiellement dans la fente, pour avoir une meilleure conservation du volume total dans le réseau.

On peut également choisir cette largeur de fente L pour avoir une célérité des ondes proche de celle des ondes sonores dans l'eau, soit environ 1500 ms^-1 avec $c=\sqrt{\frac{gS}{L}}$ .

Cette modification peut se faire également automatiquement avec l'option "baguette magique". Cela peut conduire par exemple à des fentes de l'ordre de 10^-4 m à 10^-6 m. Dans les schémas numériques explicites cela posait un problème car impliquait des pas de temps de calcul très petits afin d'avoir un nombre de Courant Cr< 1 avec $C_r=\frac{(V+c)DT}{DX}$, où V est la vitesse de l'eau, c la célérité des ondes, DT le pas de temps de calcul et DX le pas d'espace entre 2 sections de calcul.

Mais le schéma de Preissmann implicite n'est pas soumis aux mêmes contraintes. Il est inconditionnellement stable quel que soit le nombre de Courant. Mais numériquement un Cr très grand (ou aussi très petit ce qui n'est pas notre cas ici) introduit des biais numériques (diffusion et déphasage, cf thèse de Cunge). Par défaut dans SIC nous choisissons une largeur de fente de 0.01 m. Par exemple dans un dalot de 1 m², pour un débit de 1 m³s^-1, on a une célérité c de l'ordre de 30 ms^-1 avec une fente de 0.01 m. Pour avoir c de l'ordre de 1500 ms^-1 il faudrait une largeur de fente de l'ordre de 0.000005 m. L'avantage d'une fente aussi petite c'est que du coup les volumes stockés dans la fente sont faibles. L'inconvénient c'est que numériquement cela est un peu moins bien. C'est un compromis à trouver en cas de soucis.

Références bibliographiques :

[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

Résolution du système linéaire des équations de correction

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

On a trois variables inconnues pour chaque bief n°i :

en cote : $\Delta Z^i_1$ à l'amont
en cote : $\Delta Z^i_n$ à l'aval
en débit : $\Delta Q^i$ sur tout le bief

Soit N le nombre de biefs de maille. On pourrait résoudre le système des 3 N équations. Mais, de toute façon, il sera nécessaire de refaire, après correction des débits, un calcul général de ligne d'eau pour obtenir la cote dans chacune des sections des biefs.

On voit donc que les seules inconnues qui nous intéressent directement sont les $\Delta Q^i$.

Si l'on écrit les équations dans un ordre convenable on voit apparaître des blocs carrés d'ordre N :

$\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}$

Matrice . Vecteur = Constante

Soit : $M.V = K$ [8]

D, F sont des matrices diagonales
I est la matrice unité
H est une matrice triangulaire supérieure avec une diagonale principale nulle

Dans la première ligne de blocs du système, on range les équations de continuité des débits aux noeuds qui ne sont pas noeud aval de maille (remplissage du bloc A).

On complète ensuite par les relations d'égalité entre les cotes amont aux diffluences (remplissage des bloc B et C1).

Dans la deuxième ligne on a les équations de condensation.

Dans la troisième ligne on trouve les lois aval aux biefs aval de maille, ainsi que les relations d'égalité entre les cotes amont et les cotes aval des biefs aux nœuds qui ne sont pas nœud aval de maille.

En multipliant [8] par une matrice de passage P convenable :
$(P.M).V = P.K$

On obtient un système de N équations en $\Delta Q$ :

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

avec :
$\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{ Matrice unite}\cr}}$

La méthode n'est valable, que si la matrice $\chi_1$ est inversible. Or celle-ci présente la particularité de toujours pouvoir se mettre sous la forme d'une matrice triangulaire supérieure avec la diagonale composée uniquement de 1. On peut alors inverser $\chi_1$ par la méthode de Gauss-Jordan. Le système [9] est résolu par la méthode Gauss avec un pivot partiel.

Équations de correction

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

La ligne d'eau corrigée sera solution de l'équation :

$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]

En développant [2], en série de Taylor limitée aux termes du premier ordre :

[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]

avec :
$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}}$

On sait que $a_{j}^{i}$ n'est pas nul puisque la résolution de [1] est faite par la méthode de Newton où intervient le rapport $\frac{f}{a^{i}_{j}$.

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

Pour un bief i de n sections, les (n-1) équations [3] deviennent :

$\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} $

En condensant ces équations relatives à deux sections, on obtient une relation caractéristique du bief :

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

avec :

i = indice du bief
1 = indice de la section amont du bief
n = indice de la section aval du bief
$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}$ ( avec : $d^i_n = 1$ et $e^i_n = 0$ )

Les autres relations dont on dispose sont :

Pour tous les nœuds de la maille qui ne sont pas nœuds aval de maille :

Continuité des débits :

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

avec :

k = nombre de biefs de maille reliés au nœud
$\varepsilon ^i = 1$ pour un bief amont
$\varepsilon ^i = -1$ pour un bief aval

Égalité des cotes :

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

soit k - 1 relations, avec :

u, v = 1 si le bief part du nœud et n si le bief arrive au nœud
i, l = biefs reliés au nœud

Pour tous les nœuds aval de maille :

Si la condition aval est une cote donnée : $\Delta Z^i_n = 0$ [7]
Si la condition aval est une courbe de tarage aval (dans ce cas un seul bief par nœud aval) : $\Delta Z^i_n = g(\Delta Q^i) $

Classement des biefs

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

Les noms des nœuds permettent de décrire l'enchaînement des biefs. La position d'un bief dans le réseau est entièrement définie par la connaissance des noms des nœuds amont et aval. On fixe ainsi en même temps le sens de l'écoulement. On décrit ainsi simplement la topologie du réseau qui est donc représenté sous la forme d'un graphe orienté. Les nœuds sont les sommets et les biefs sont les arcs de ce graphe. Ils sont numérotés automatiquement par le programme en fonction de l'ordre dans lequel ils sont introduits dans le fichier de données.

Un écoulement fluvial étant régi par l'aval, le calcul d'une ligne d'eau s'effectue en remontant l'écoulement, de l'aval vers l'amont. Afin de démarrer le calcul, il faut donc disposer d'une relation reliant la cote et le débit à la condition limite aval.

Un réseau hydraulique peut être représenté par un graphe orienté connexe sans circuit, sur lequel on effectue un parcours à l'aide de l'algorithme suivant :

On cherche un sommet sur lequel n'arrive aucun arc et d'où il ne part qu'un seul arc. On sélectionne ce dernier dans une liste puis on le supprime du graphe. On applique récursivement le procédé sur le sous graphe obtenu tant qu'il existe des sommets ainsi définis. Ensuite soit le problème est terminé, soit il reste des arcs non récupérés. Dans ce cas, on cherche dans le sous-graphe restant un sommet sur lequel n'arrive aucun arc et d'où il part plus d'un arc. On sélectionne ces derniers, puis on les supprime du graphe et on applique récursivement le procédé sur le sous-graphe obtenu tant qu'il existe des sommets ainsi définis. Ensuite soit le problème est terminé, soit il reste des arcs non récupérés. Dans ce cas, on recommence sur le sous-graphe restant l'algorithme depuis le début.

En lisant à l'envers la liste des arcs ainsi obtenue, on trouve une relation d'ordre partiel satisfaisante (Remarque : seuls les nœuds amont de biefs sont concernés par cet algorithme). De la liste, on tire également le repérage des biefs de maille. Ces biefs sont groupés en commençant par la première diffluence (les deux premiers biefs successifs ayant le même nœud amont). Si le modèle a plusieurs conditions aval en cote, alors tous les biefs suivant la première diffluence de la liste appartiennent à la maille. Si, en revanche, la maille se referme, elle s'arrête à la dernière confluence (le dernier bief ayant le même nœud aval qu'un autre bief situé dans la liste et localisé après le début de la maille).

Appliquons le 1er critère :
Le premier sommet sur lequel n'arrive aucun arc et d'où il ne part qu'un seul arc est le sommet DAM. On retient le bief [1] dans la liste et on le supprime du graphe.
On a alors :

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

On ne peut plus trouver de nœud étant nœud amont d'un unique bief. On cherche alors un nœud n'ayant jamais été nœud aval, mais étant nœud amont de plusieurs biefs. Le nœud DC1 vérifie ces conditions. Les biefs 2 et 4 sont alors ajoutés à la liste :

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

Le même algorithme est appliqué à nouveau. On obtient ainsi le nœud DC2 et les biefs 3 et 6 sont ajoutés à la liste.

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

Il n'existe plus de nœud ayant deux biefs aval. On applique alors le premier algorithme. On obtient ainsi le nœud FC1, et le bief 5 est ajouté à la liste.

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

La maille commence à partir du premier nœud de diffluence (biefs 2 et 4, nœud DC1). Tous les biefs suivants appartiennent à la maille, car cette maille n'est pas fermée (il existe plusieurs conditions limite aval).

Classement des biefs de maille pour le calcul en régime permanent - Solution initiale

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

On fixe un débit initial dans chaque bief de maille en respectant la continuité des débits aux nœuds, en prenant par exemple une répartition au prorata de la section moyenne du bief (rapport du volume du bief aux cotes de débordement sur sa longueur).

Puis, on calcule la ligne d'eau dans chaque bief successif, en suivant l'ordre établi (Cf. Classement des biefs) et en prenant pour condition aval (dans le cas d'un nœud défluent) la moyenne des cotes obtenues lors du calcul des biefs aval partant de ce nœud.

Soit une partie du bief n°i, définie par deux sections, on cherche la cote Zⁱ_j (dans la section amont n°j) connaissant la cote Zⁱ_j+1 (dans la section aval n°j+1) et le débit Qⁱ_j s'écoulant dans la section n°j.
Cela revient à résoudre numériquement l'équation de la courbe de remous qui peut se mettre sous la forme :

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

Sections paramétrées

David Dorchies — 2009-10-22T15:06:46Z

On peut introduire les sections de forme géométrique particulière sous une forme paramétrée.
Cercle : défini par le rayon R, la cote du fond ZF et la cote de berge ZB. Si ZB n'est pas précisé l'on prend ZB = ZF + 2R.

Dalot : défini par la largeur au plafond L et le fruit m (qui peut être négatif), ainsi que par la cote du fond ZF et de couverture ZB.

Puissance : défini par XLo, a et la cote de berge et du fond.
L = XLo * (Y/Yo)a
avec :
Yo = ZB - ZF
Y = Z - ZF tirant d'eau
Le profil parabolique correspond à a = 0.5

Rectangle : défini par la largeur, les cotes du fond et de berge.

Trapèze ou triangle : défini par la largeur au plafond L (nulle pour le triangle) et les cotes du fond et de berge. Le fruit peut être négatif dans le cas du trapèze. En cas de fruit négatif si avant ZB la section se referme, le programme déplace la cote de berge.

Sections singulières

David Dorchies — 2009-10-22T15:06:46Z

Les sections comportant des ouvrages en travers sont appelées singulières. Dans ces sections, les lois générales de calcul de la ligne d'eau ne s'appliquent pas. Ces lois sont remplacées par les formules de débit des ouvrages.

Au niveau de la géométrie, on ne décrit pas les ouvrages hydrauliques, mais on indique les sections où ils sont situés. Il faut donc, à ce stade, introduire le gabarit du canal et non celui des ouvrages.


Profil structure en travers

Une singularité est considérée comme ponctuelle par rapport à l'ensemble du bief. Il faut donc connaître 2 sections, que l'on supposera à la même abscisse, correspondant à l'amont et à l'aval de la singularité.

La section amont caractérise le gabarit du canal à l'amont des ouvrages, la section aval caractérise le gabarit du canal à l'aval des ouvrages. Si l'on n'introduit qu'une section pour décrire une section singulière, le programme génère automatiquement une section supplémentaire (voir § ci-dessus).

Remarque : attention, il ne faut pas (pour un calcul en régime transitoire) que les premières sections de tous les biefs partant d'un noeud amont de modèle soient singulières, ni que les dernières sections de tous les biefs arrivant à un noeud aval de modèle soient singulières. Sinon dans ce cas on aura un message d'erreur "Erreur numérique (pivot système réduit)".

Sections de calcul

David Dorchies — 2009-10-22T15:06:46Z

Les sections de données sont réparties inégalement le long du canal. En effet, le projeteur choisit les sections qui sont le mieux représentatives du gabarit du canal, des changements de pentes, etc. Donc suivant que le canal est régulier ou pas les sections de données sont peu ou très espacées.

Le calcul hydraulique nécessite un pas d'espace permettant une bonne estimation de la ligne d'eau. Celui-ci est fixé par le projeteur en fonction de la connaissance qu'il a du comportement hydraulique du canal (quand les sections de données sont trop éloignées, le modèle interpole des sections de calcul supplémentaires pour permettre une meilleure simulation de la ligne d'eau).

Toutes les sections de données sont gardées comme section de calcul. Quel que soit le type retenu pour définir une section de donnée, le programme la transforme en largeur cote pour le stockage et l'interpolation. On ne conserve que des quadruplés : cotes, largeurs, périmètres mouillés, surfaces. Si la section a été introduite en abscisse-cote, le périmètre tient compte de la dissymétrie de la section.

Transformation d'une section de données en largeur-cote

David Dorchies — 2009-10-22T15:06:46Z

Le programme recherche sur une section en travers tous les points hauts et tous les points bas. Il détermine donc les berges droites et gauches et élimine les points hors du lit. Ensuite à partir des points hauts et bas, il découpe la section en travers en chenal et pour chaque chenal il fait la transformation en largeur-cote. Puis il somme les largeurs de tous les chenaux pour chaque cote.

Interpolation des sections de calcul

David Dorchies — 2009-10-22T15:06:46Z

Si deux sections de données sont distantes de plus d'un pas d'espace, le programme interpole à partir de ces deux sections de données des sections de calcul, au pas fixé par l'utilisateur.

En réalité, le pas est ajusté pour donner un nombre entier d'intervalles de calcul égaux entre les deux sections considérées. Cette interpolation est faite à tirant d'eau constant.

La section de calcul peut donc avoir jusqu'à 2 fois plus de points que la section de donnée. Il faut donc prévoir des critères pour éliminer certains points afin de ne pas alourdir le nombre de points à stocker.

On élimine tout point de la section de calcul qui ne modifie pas la surface de la section de plus de 1 %, de même pour les points qui ne modifient pas le périmètre mouillé de plus de 5 %. Ainsi quel que soit le tirant d'eau, on est sûr d'avoir une estimation de la surface au moins à 1 % près et du périmètre à 5 %.

Les sections de calcul sont complétées verticalement par un point à 500 m de la cote de débordement pour permettre le calcul après débordement.

Les sections de calcul sont numérotées par bief, ainsi si le pas de calcul est modifié sur un bief la numérotation des autres biefs reste inchangée.