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

Double balayage - Cas général

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

Considérons un bief ayant n sections de calcul, le système à résoudre est donc :

équations de SAINT VENANT sur chaque intervalle entre 2 sections de calcul à tout instant t,
conditions aux limites amont et aval.

La discrétisation transforme ce bief en une suite de n sections de calculs reliées entre elles par les deux équations linéaires [34] et [35]. On a donc 2(n-1) équations linéaires en DQ, DZ. Les deux équations manquantes pour la résolution du système sont fournies par les conditions amont et aval que l'on linéarise à chaque pas de temps.

Condition à la limite amont :

R1.DQ1 + S1.DZ1 = T1

Dans le cas d'une loi Q(t) on a :

R1 = 1, S1 = 0 et T1 = Q(t+Dt) - Q(t)

Condition à la limite aval :

R'n.DQn + S'n.DZn = T'n, avec :
= a, R'n = 1, S'n = -a, T'n = 0

On a donc à résoudre un système linéaire de 2 n équations. Au lieu d'inverser la matrice du système on utilise la méthode du double balayage :

A11.DQi + A12.DZi = B11.DQj + B12.DZj + B13
A21.DQi + A22.DZi = B21.DQj + B22.DZj + B23

Posons DEN = A11 A22 - A21 A12

DEN est différent de zéro car les deux équations sont linéairement indépendantes.

DEN.DQi = (B11.A22 - B21.A12).DQj
+ (B12.A22 - B22.A12).DZj
+ B13.A22 - B23.A12
DEN.DZi = (B21.A11 - B11.A21).DQj
+ (B22.A11 - B12.A21).DZj
+ B23.A11 - B13.A21
[36]
avec :
A =
B =
C =
D =
E =
F =

Regardons la forme du système pour trois sections de calcul :
R1S1DQ1T11-D1-E1DZ1F11-A1-B1 .DQ2 =C11-D2-E2DZ2F21-A2-B2DQ3C2R'3S'3DZ3T'3

1er balayage :

Le premier balayage amont - aval donne :
R1S1DQ1T11-D1-E1DZ1F1R2S2 .DQ2 =T21-D2-E2DZ2F2R3S3DQ3T3R'3S'3DZ3T'3

Soit pour deux sections consécutives i et j :
Ri.DQi + Si.DZi = Ti
Ri.(Ai.DQj + Bi.DZj + Ci)
+ Si.(Di.DQj + Ei.DZj + Fi) = Ti
R*j.DQj + S*j.DZj = T*j

avec :

R*j = Ri.Ai + Si.Di
S*j = Ri.Bi + Si.Ei
T*j = Ti - Ri.Ci - Si.Fi

On cadre les coefficients R'j, S'j, T'j en normant, pour éviter la propagation d'erreurs numériques.

Rj = etc.

On a alors :

Rj.DQj + Sj.DZj = Tj

qui est la relation d'impédance amont.

Arrivé à la condition à la limite aval, on a :

qui nous donne DQn et DZn :
DQn =
DZn =
Puis on remonte le deuxième balayage.

2ème balayage

Le balayage permet de calculer les variables DQ et DZ dans chaque section au moyen des équations [37] :

[37]

Equation dynamique

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

Discrétisons chaque terme de l'équation [30] :

* =
Calcul de :

= [ - ] + [D( ) - D( )]
= [ - ] + [( - DAj) - ( - DAi)]
= [ - ] + DQj - DQi + BiDZi - BjDZj)
* = a1.DQi + a2.DQj + a3.DZi + a4.DZj + a5

avec :
a5 = [ - ]
a4 = - Bj
a3 = Bi
a2 =
a1 = -

Calcul de gA :

gA = g [ + Q].[]
= [(Ai + Aj)(Zj - Zi) + Q(DAi + DAj)(Zj - Zi) + Q(DZj - DZi)(Ai + Aj)]
= [(Ai + Aj)(Zj - Zi) + Q(BiDZi + BjDZj)(Zj - Zi) + Q(DZj - DZi)(Ai + Aj)]
* gA = [a6DZi + a7DZj + a8]

avec :

a6 = Q[Bi (Zj - Zi) - (Ai + Aj)]
a7 = Q[Bj (Zj - Zi) + (Ai + Aj)]
a8 = (Ai + Aj).(Zj - Zi)

Calcul de gASf :

gASf = g [ + [D(AiSfi) + D(AjSfj)]]
= g [ + [SfiDAi + AiDSfi + SfjDAj + AjDSfj]]

Calculons DSf :

Sf = = sgn(Q)

avec :

sgn(Q) = 1 si Q > 0
sgn(Q) = -1 si Q < 0

=> DSf = 2 sgn(Q).Q - 2 sgn (Q) .DDe
ou encore : DSf = 2 |Q| - 2 .DDe
De =

=> DDe = + R-1/3 DR
avec :

=> DDe = B.DZ + R-1/3 a DZ
= ae DZ
DSf = 2 |Q| - 2 DDe
= 2 |Q| - 2 ae DZ
=> gASf = g [ + [(SfiBiDZi + SfjBjDZj)
+ Ai (2 |Qi| - 2 aei DZi)
+ Aj (2 |Qj| - 2 aej DZj)]]
* gASf = a9.DZi + a10.DZj + a11.DQi + a12.DQj + a13

avec :

a9 = [SfiBi - Ai 2 aei]
a10 = [SfjBj - Aj 2 aej]
a11 = gQ
a12 = gQ
a13 = g []

Calcul de kq :

kq = [qi + qj + Q (D(qi ) + D(qj ))]
= [qi + qj + Q (qi - qi + qj - qj )]

* kq = a14.DZi + a15.DZj + a16.DQi + a17.DQj + a18

avec :
a14 = - Q qi Bi
a15 = - Q qj Bj
a16 = Q
a17 = Q
a18 = [qi + qj ]

L'équation dynamique [30] peut donc se mettre sous la forme :

A11.DQi + A12.DZi = B11.DQj + B12.DZj + B13 [35]

avec :
A11 = + a1 + a11 - a16
A12 = a3 + a6 + a9 - a14
B11 = - ( + a2 + a12 - a17)
B12 = - (a4 + a7 + a10 - a15)
B13 = - (a5 + a8 + a13 - a18)

Equation de continuité

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

Discrétisons chaque terme de l'équation [28] :
= DAi » Bi.DZi
=> =
= + Q
q = qij (débit latéral entre les sections i et j)
L'équation [28] devient donc :
+ + Q =
+ + DQj - DQi = (qi + qj)
DQi - Bi.DZi = DQj + Bj.DZj + - (qi + qj)
qui est de la forme :
A21.DQi + A22.DZi = B21.DQj + B22.DZj + B23 [34]
avec :
A21 = 1 B21 = 1
A22 = - Bi B22 = Bj
B23 = - (...)

- Discrétisation semi-implicite

Equations générales du schéma

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

Posons :
Dfi = fA' - fA fA = fi
Dfj = fB' - fB fB = fj
L'expression d'une fonction en M s'écrit :
fM = (1-Q) + Q
fM = + (Dfi + Dfj) [31]

L'expression de la dérivée () en M s'écrit :
()M = (1-Q) + Q
=> ()M = + Q [32]
L'expression de la dérivée () en M s'écrit :
()M = ( + )
=> ()M = [33]

Nous allons utiliser les équations [31], [32] et [33] pour discrétiser les équations de ST VENANT.

Discrétisation semi-implicite

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

Les équations de SAINT VENANT n'ayant pas de solution analytique connue en géométrie réelle, on résout numériquement en discrétisant les équations, c'est-à-dire que l'on remplace les dérivées partielles par des différences finies. La discrétisation retenue est le schéma semi-implicite de PREISSMANN.
Ce schéma est implicite car il fait intervenir dans l'expression des dérivées partielles d'espace les valeurs des variables au pas de temps (...)

- Discrétisation semi-implicite

Conditions initiales et aux limites

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

Les équations aux dérivées partielles doivent être complétées par des conditions aux limites et des conditions initiales pour pouvoir les résoudre. Les conditions aux limites sont les hydrogrammes dans les noeuds amont de bief et les courbes de tarage dans les noeuds aval de modèle (car on est en écoulement fluvial). La condition initiale est la ligne d'eau fournie par le calcul en permanent (Unité (...)

- Les équations de Saint-Venant

Equation dynamique

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

Nous appliquons le théorème de quantité de mouvement à un volume d'eau compris entre les sections d'abscisses x et x + Dx.
Nous nous intéressons à la projection sur l'axe du canal :
= S Fext/x [29]

Pendant le temps Dt le volume d'eau V compris entre les sections 1 et 2 se déforme et se déplace dans le volume V' compris entre les sections 1' et 2'. Il nous faut estimer la variation de quantité de mouvement en projection sur l'axe des x terme de gauche de l'équation [29].

Estimation de la quantité de mouvement

La quantité de mouvement perdue correspond au volume V1 soit :
rV(x,t).Dt.A(x,t).V(x,t)
La quantité de mouvement gagnée correspond au volume V2 :
rV(x+Dx,t).Dt.A(x+Dx,t).V(x+Dx,t)
La variation de quantité de mouvement dans la partie commune est de :
rA(x,t+Dt).Dx.V(x,t+Dt) - rA(x,t).Dx.V(x,t)
La variation de quantité de mouvement due aux apports latéraux :

rq(x,t).Dx.Dt.V(x,t)
si débit de fuite (vitesse égale à la vitesse de l'écoulement)
+ rq(x,t).Dx.Dt.0
si débit d'apport (vitesse nulle en projection sur l'axe des x)
que l'on peut écrire de façon globale sous la forme :
krq(x,t).Dx.Dt.V(x,t)

avec :
k = 1 pour un débit de fuite et
k = 0 pour un débit d'apport.

On a donc :
D(mVx) = (rV.Dt.A.V) (x+Dx,t) - (rV.Dt.A.V)(x,t) + (rA.Dx.V) (x,t+Dt)

(rA.Dx.V)(x,t) - k(rq.Dx.Dt.V) (x,t)
que l'on peut écrire sous la forme :
= rDx [ + - kqV]

Il faut évaluer à présent le terme de droite de l'équation [29] qui représente la résultante des forces extérieures en projection sur l'axe des x. On considère uniquement l'effet des forces de gravité, de pression et de frottement.

Estimation des forces extérieures

La résultante des forces de gravité est :
Fgx = r.g.A.Dx.sin(S0)

comme on est en écoulement quasi horizontal sin(S0) » S0 (pente du canal), donc :
Fgx = rg.A.Dx.S0

La résultante des forces de pression est obtenue grâce à l'hypothèse de répartition hydrostatique des pressions. La résultante des forces de pression s'exerçant sur la masse d'eau comprise entre x et x+Dx est la même qu'en statique, c'est-à-dire quand la surface libre est horizontale nous pouvons alors écrire :

Fg = - rg.A.Dx = -Fp (résultante des forces nulle)
Donc : |Fp| = rg.A.Dx et Fp est perpendiculaire à la surface libre.

La normale à la surface libre a pour composantes :
n
car l'équation de la surface libre dans des axes liés au fond est h - y (x,t) = 0
On suppose que |n| = 1 car la surface libre est presque horizontale. La résultante des forces de pression s'écrit donc :
Fp
Donc :
Fpx = -rg.A.Dx
Les forces de frottement sont estimées à partir de la formule de MANNING-STRICKLER :
Ffx = - rg.A.Dx.Sf
avec : Sf =
La résultante des forces extérieures en projection sur l'axe du canal est donc :
S Fext/x = rg.A.Dx.S0 - rg.A.Dx. - rg.A.Dx.Sf
= rg.A.Dx [S0 - Sf - ]
L'équation [29] se met donc sous la forme :
rDx [ + - kqV] = rg.A.Dx [S0 - Sf - ]
Soit :
+ + g.A = g.A [S0 - Sf] + kqV
que nous écrivons sous la forme :
+ + g.A = - g.A Sf + kqV [30]

Equation de continuité

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

Cette équation traduit la conservation de la masse d'eau. On considère la variation du volume d'eau compris entre deux sections à l'abscisse x et x + Dx pendant le temps Dt.

Masse d'eau entrante :
r.Q(x,t).Dt + r.q.Dx.Dt
Masse d'eau sortante :
r.Q(x+Dx,t).Dt
Variation du stock :
r Vt+Dt - r Vt = r(A.Dx)t+Dt - r(A.Dx)t
Le bilan exprimant la conservation de la masse d'eau s'écrit :
r(A.Dx)t+Dt - r(A.Dx)t = r(Q.Dt)x + rq.Dx.Dt - r(Q.Dt)x+Dx
Par passage à la limite, on obtient l'équation [28] :

+ = q [28]

Les équations de Saint-Venant

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

Le canal étant découpé en zones homogènes, les biefs, nous examinons comment se présente l'étude de la ligne d'eau transitoire sur un bief. Nous verrons ensuite comment résoudre le problème pour l'ensemble du réseau hydraulique.
On prend les mêmes hypothèses que pour l'Unité 2. De plus, on ne considère que des phénomènes transitoires lisses donc la propagation de fronts raides ne peut être (...)

- Les équations de Saint-Venant

Liste des symboles

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

A = surface mouillée
B = largeur au miroir
C_{F =} coefficient de débit pour une vanne de fond dénoyée
C_{G =} coefficient de débit de référence pour une vanne
C_{S =} coefficient de débit pour une vanne de fond noyée
g = accélération de la gravité
H = charge totale
h = hauteur de l'eau au dessus du fond
h_{s =} hauteur de la vanne (pour calculer les débordements)
K = coefficient de Strickler
k = 0 pour un débit entrant, 1 pour un débit sortant
k_{F =} coefficient de réduction pour une vanne noyée
L = longueur du seuil
n = coefficient de Manning
Q = débit
q = débit latéral
R = rayon hydraulique
S_{o =} pente du fond
S_{f =} pente de frottement
V = vitesse moyenne
W = ouverture d'une vanne de fond ou d'une prise
x = abscisse dans le sens de l'écoulement
y = tirant d'eau
Z = cote
Z_c= cote critique
e = faible incrément
m = coefficient de débit pour un écoulement en charge dénoyé
m' = coefficient de débit pour un écoulement en charge noyé
m_F= coefficient de débit pour un écoulement à surface libre dénoyé
m_S= coefficient de débit pour un écoulement à surface libre noyé
S = symbole de sommation