Modélisation 1D des Réseaux d'Assainissement
Modélisation 1D des Réseaux d'Assainissement
THESE
Marc Buyer
Membres du jury
2 PLAN DE LA THESE.................................................................................................................................... 4
1 INTRODUCTION .......................................................................................................................................... 6
1 INTRODUCTION ........................................................................................................................................ 51
ANNEXES
LISTE DES FIGURES
Figure 1 : Section en travers d’un canal........................................................................................7
Figure 2 : Volume de contrôle.......................................................................................................8
Figure 3 : Grandeurs caractéristiques du ressaut hydraulique.....................................................10
Figure 4 : Surface d’application des forces de pressions ............................................................13
Figure 5 : Support d’une fonction φ quelconque et domaine d’intégration ................................14
Figure 6 : Boule centrée sur P et courbe régulière Γ...................................................................15
Figure 7 : Transitions d’un écoulement à surface libre à un écoulement en charge ...................20
Figure 8 : Distributions de pression et de vitesse dans le cas d’un écoulement en charge (a) et
d’un écoulement à surface libre (b).....................................................................................21
Figure 9 : Lignes de charges dans le cas d’un écoulement en charge (a) et d’un écoulement à
surface libre (b) ...................................................................................................................22
Figure 10 : Représentation d’un collecteur avec l’artifice numérique de la fente de Preismann24
Figure 11 : Jonction de collecteurs..............................................................................................27
Figure 12 : Schéma de principe d’une confluence ......................................................................33
Figure 13 : Lignes d’eau possibles dans un déversoir latéral prismatique..................................34
Figure 14 : Types de profils de surface pour des nombres de froude F<1 et F>1 dans les cas
prismatique (θ>0) et non prismatique (θ<0) .......................................................................35
Figure 15 : Exemple d’abaques : canaux rectangulaires avec entonnement. ..............................38
Figure 16 : Vitesses longitudinale et latérale dans le déversoir ..................................................38
Figure 17 : Comparaison des profils d’écoulement sur un déversoir frontal et latéral ...............40
Figure 18 : Effet de la vitesse latérale .........................................................................................41
Figure 19 : Application du théorème de Bernoulli......................................................................41
Figure 20 : Effet de l’entonnement .............................................................................................42
Figure 21 : Exemple d’un déversoir de Sélestat..........................................................................44
Figure 22 : Ecoulement à travers un élargissement brusque.......................................................47
Figure 23 : Surface autour du point P(xP,tP) projetée sur un plan (x,t) (grille irrégulière)..........55
Figure 24 : Croisement des caractéristiques dans le cas du ressaut ............................................57
Figure 25 : Ondes issues d’une rupture de barrage .....................................................................58
Figure 26 : Quadrillage du plan (x,t)...........................................................................................60
Figure 27 : Schéma de type explicite ..........................................................................................61
Figure 28 : Schéma de type implicite..........................................................................................62
Figure 29 : Schéma de discrétisation dans le plan (x,t)...............................................................63
Figure 30 : Représentation du problème de Riemann .................................................................68
Figure 31 : Formation d’une onde de choc. ................................................................................69
Figure 32 : Formation d’une onde de choc de raréfaction ..........................................................70
Figure 33 : Modification des caractéristiques dans le cas de l’onde de choc de raréfaction ......71
Figure 34 : Faisceau de caractéristiques dans le cas du système d’équations.............................73
Figure 35 : Représentation d’un volume de contrôle ..................................................................75
Figure 36 : Décomposition du domaine d’intérêt en cellules de calculs et flux aux interfaces. .76
Figure 37 : Représentation de la valeur d’une variable u dans une cellule de calcul..................86
Figure 38 : Variables discontinues au niveau des interfaces.......................................................86
Figure 39 : Forme des caractéristiques........................................................................................87
Figure 40 : Représentation de la fonction associée à un limiteur de flux ...................................88
Figure 41 : Structure générale des logiciels réalisés ...................................................................98
Figure 42 : Schématisation des cellules de calcul du système de collecteurs réels et des tronçons
fictifs..................................................................................................................................100
Figure 43 : Conditions aux limites applicables aux extrémités des tronçons............................101
Figure 44 : Choix de la hauteur normale pour condition à la limite amont ..............................102
Figure 45 : Choix de la hauteur critique pour condition à la limite amont ...............................102
Figure 46 : Choix d’une hauteur imposée pour condition à la limite amont.............................103
Figure 47 : Choix de la hauteur normale pour condition à la limite aval..................................104
Figure 48 : Choix de la hauteur critique pour condition à la limite aval...................................104
Figure 49 : Choix d’une hauteur imposée pour condition à la limite aval ................................105
Figure 50 : Vue du canal expérimental d’Obernai ....................................................................107
Figure 51 : Lignes d’eau obtenues numériquement et expérimentalement pour la courbe S1 .109
Figure 52 : Lignes d’eau obtenues numériquement et expérimentalement pour la courbe S3 .111
Figure 53 : Lignes d’eau obtenues numériquement et expérimentalement pour la courbe M1 112
Figure 54 : Lignes d’eau obtenues numériquement et expérimentalement pour la courbe M2 113
Figure 55 : Lignes d’eau obtenues numériquement et expérimentalement pour la courbe M3 115
Figure 56 : Lignes d’eau obtenues numériquement à plusieurs instants dans le cas 1 .............117
Figure 57 : Lignes d’eau obtenues numériquement à plusieurs instants dans le cas 2 .............119
Figure 58 : Lignes d’eau obtenues numériquement à plusieurs instants dans le cas 3 .............121
Figure 59 : Lignes d’eau obtenues numériquement à plusieurs instants dans le cas 4 .............123
Figure 60 : Vue du dispositif de mesure de débit......................................................................126
Figure 61 : Vue de dessus du Venturi .......................................................................................126
Figure 62 : Ligne d’eau dans le Venturi pour un débit de 26 m3h-1..........................................127
Figure 63 : Ligne d’eau dans le Venturi pour un débit de 44.8 m3h-1.......................................128
Figure 64 : Ligne d’eau dans le Venturi pour un débit de 57.3 m3h-1.......................................128
Figure 65 : Ondes croisées apparaissant dans la rectangulaire du Venturi ...............................129
Figure 66 : Schéma de principe du banc d’essai physique de déversoir. ..................................130
Figure 67 : Vue du déversoir du banc d’essai ...........................................................................131
Figure 68 :Chaîne d’acquisition des valeurs de débit ...............................................................132
Figure 69 : Méthode manuelle de mesure des hauteurs d’eau ..................................................133
Figure 70 : Vues du dispositif 3D d’acquisition des hauteurs d’eau.........................................134
Figure 71 : Visualisation 3D des hauteurs d’eau dans le déversoir après interpolation............134
Figure 72 : Comparaison hauteurs d’eau et débits expérimentaux et numériques dans un cas
fluvial pour 1 crête déversante et w=7 cm ........................................................................139
Figure 73 : Comparaison hauteurs d’eau et débits expérimentaux et numériques dans un cas
torrentiel pour 1 crête déversante et w=7 cm ....................................................................139
Figure 74 : Comparaison hauteurs d’eau et débits expérimentaux et numériques dans un cas de
ressaut pour 1 crête déversante et w=7 cm........................................................................140
Figure 75 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰,
Iav=8‰, Q=88.7 m3h-1, 1 crête déversante et w=5 cm ......................................................143
Figure 76 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰,
Iav=8‰, Q=65.2 m3h-1, 1 crête déversante et w=5 cm ......................................................143
Figure 77 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰,
Iav=0.5‰, Q=73.3 m3h-1, 1 crête déversante et w=12.5 cm ..............................................145
Figure 78 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰,
Iav=8‰, Q=117.7 m3h-1, 1 crête déversante et w=12.5 cm ...............................................146
Figure 79 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰,
Iav=0.5‰, Q=73.3 m3h-1, 2 crêtes déversante et w=5 cm..................................................148
Figure 80 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰,
Iav=8‰, Q=73.9 m3h-1, 2 crêtes déversante et w=5 cm.....................................................148
Figure 81 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰,
Iav=8‰, Q=84 m3h-1, 2 crêtes déversante et w=7.5 cm.....................................................150
Figure 82 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰,
Iav=0.5‰, Q=89 m3h-1, 2 crêtes déversante et w=7.5 cm..................................................151
Figure 83 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰, .153
Figure 84 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=0.5‰,
Iav=8‰, Q=108 m3h-1, 1 crête déversante et w=12.5 cm ..................................................154
Figure 85 : Comparaison des hauteurs d’eau expérimentales et numériques pour Iam=6‰,
Iav=1.33‰, Q=103 m3h-1, 1 crête déversante et w=9.4 cm ...............................................156
Figure 86 :Définition des variables utilisées dans le logiciel DO. ............................................172
Figure 87 : Détermination de la hauteur d’eau admissible à l’amont (Tram) ...........................173
Figure 88 : Forme de la ligne d’eau au débit décennal .............................................................174
Figure 89 :Déversoir à crête basse : Forme de la ligne d’eau au débit décennal ......................176
LISTE DES TABLEAUX
Tableau 1 : Conditions d’apparition des six lignes d’eau décrites et détermination empirique de
la vitesse latérale .................................................................................................................35
Tableau 2 : Formules empiriques et conditions d’application ....................................................36
Tableau 3 : Divers types de schémas numériques.......................................................................61
Tableau 4 : Configuration utilisée pour la courbe S1................................................................109
Tableau 5 : Configuration utilisée pour la courbe S3................................................................110
Tableau 6 : Configuration utilisée pour la courbe M1 ..............................................................111
Tableau 7 : Configuration utilisée pour la courbe M2 ..............................................................113
Tableau 8 : Configuration utilisée pour la courbe M3 ..............................................................114
Tableau 9 : Configuration utilisée dans le cas 1 pour la validation en transitoire ....................116
Tableau 10 : Temps de propagation du ressaut obtenus dans le cas 1 ......................................116
Tableau 11 : Configuration utilisée dans le cas 2 pour la validation en transitoire ..................118
Tableau 12 : Temps de propagation du ressaut obtenus dans le cas 2 ......................................118
Tableau 13 : Configuration utilisée dans le cas 3 pour la validation en transitoire ..................120
Tableau 14 : Temps de propagation du ressaut obtenus dans le cas 3 ......................................120
Tableau 15 : Configuration utilisée dans le cas 4 pour la validation en transitoire ..................122
Tableau 16 : Temps de propagation du ressaut obtenus dans le cas 4 ......................................122
Tableau 17 : Synthèse des résultats concernant la validation pour les collecteurs ...................124
Tableau 18 : Configurations de pentes étudiées pour 1 crête déversante et w=7 cm................135
Tableau 19 : Comparaison débits déversés simulés - débits déversés mesurés pour 1 crête
déversante et w=7 cm........................................................................................................137
Tableau 20 : Configurations de pentes étudiées pour 1 crête déversante et w=5 cm................141
Tableau 21 : Comparaison débits déversés simulés - débits déversés mesurés pour 1 crête
déversante et w=5 cm........................................................................................................142
Tableau 22 : Configurations de pentes étudiées pour 1 crête déversante et w=12.5 cm...........144
Tableau 23 : Configurations de pentes étudiées pour 2 crêtes déversante et w=5 cm ..............147
Tableau 24 : Comparaison débits déversés simulés - débits déversés mesurés pour 2 crêtes
déversante et w=5 cm........................................................................................................147
Tableau 25 : Configurations de pentes étudiées pour 2 crêtes déversante et w=7.5 cm ...........149
Tableau 26 : Comparaison débits déversés simulés - débits déversés mesurés pour 2 crêtes
déversante et w=7.5 cm.....................................................................................................150
Tableau 27 : Configurations de pentes étudiées pour 2 crêtes déversante et w=12.5 cm .........151
Tableau 28 : Comparaison débits déversés simulés - débits déversés mesurés pour 2 crêtes
déversante et w=12.5 cm...................................................................................................152
Tableau 29 : Configurations de pentes étudiées pour le déversoir à crête haute ......................155
Tableau 30 : Comparaison débits déversés simulés - débits déversés mesurés pour le déversoir à
crête haute .........................................................................................................................155
LISTE DES SYMBOLES
a: Valeur propre
a% : Valeur propre de la matrice jacobienne approchée
A: Section mouillée.
b: Valeur de Bd pour x=0
B: Largeur au miroir
Bd : Largeur du canal dans le déversoir
c: Célérité des ondes
C: Coefficient de Chézy.
Cd : Coefficient de débit dans le cas du déversoir frontal.
Ccfl : Nombre de Courant.
d: Elément du vecteur correctif
D: Vecteur correctif
Dam : Diamètre de la conduite amont
Dav : Diamètre de la conduite aval
Dh : Profondeur hydraulique
e: Vecteur propre associé à la matrice jacobienne
e% : Vecteur propre associé à la matrice jacobienne approchée
F: Vecteur flux
Fi : Flux calculé dans la cellule de calcul i.
Fi+1/2 : Flux à l’interface des cellules i et i+1.
Fpression : Force de pression
Fr : Nombre de Froude
G: Vecteur source
g: Accélération de la pesanteur
h: Tirant d’eau
ham : Hauteur d’eau sur le seuil à l’amont du déversoir (m)
hav : Hauteur d’eau sur le seuil à l’aval du déversoir (m)
hm : Hauteur d’eau moyenne sur le seuil
H: Charge totale
Ham : Energie spécifique à l’amont du seuil (m)
Hav : Energie spécifique à l’aval du seuil (m)
Hs : Energie spécifique.
Iam : Pente de la conduite amont (m)
Iav : Pente de la conduite aval (m/m)
I1 : Pression hydrostatique
I2 : Pression latérale
Imot : Perte de charge dans la conduite aval
J: Matrice jacobienne
J% : Matrice jacobienne approchée
k=n*c : Coefficient associé à la géométrie du déversoir. Avec c coefficient de forme du
déversoir (c=1 dans le cas ou le déversoir a une paroi mince) et n* nombre de
parois déversantes (n*=1 ou 2).
L: Limiteur de flux pour les méthodes TVD amont, symétrique, MUSCL et Mac-
Cormack.
Ldo : Longueur du déversoir (m)
n: Coefficient de Manning
pam : Hauteur du seuil à l’amont (m)
pav : Hauteur du seuil à l’aval (m)
P: Périmètre mouillé
Pentedo : Pente du déversoir (m/m)
p: Pression
q: Rapport entre deux débits
Q: Débit
Qdev : Variation de débit du au déversement
Q10 : Débit décennal (m3/s)
Qts : Débit de temps sec (m3/s)
Qcr : Débit critique (m3/s)
R: Matrice des vecteurs propres
Rh : Rayon hydraulique.
S0 : Pente du canal
Sf : Terme de frottement
tam : Tirant d’eau dans la conduite amont (m)
Tram : Hauteur de remplissage dans la conduite amont (≠ tam si régime torrentiel) (m)
TV : Variation totale
u: Vitesse moyenne dans la direction principale de l’écoulement
u% : Vitesse dans la direction principale de l’écoulement obtenue par la moyenne de
Roe
U: Vecteur écoulement
v: Vitesse dans la direction latérale
n
Vmax : Vitesse de l’onde la plus rapide dans le domaine de calcul
w: Hauteur de crête du déversoir.
W: Variable adimensionnelle de hauteur de crête.
X: Variable adimensionnelle de distance suivant la longueur du déversoir.
y: Variable adimensionnelle de hauteur d’eau.
ξ: Coefficient de perte d’énergie
ρ: Masse volumique
ν: Coefficient défini en fonction du rapport entre charge et diamètre aval
Θ: Variable d’entonnement
θ: Angle associé au changement longitudinal de largeur
Φ: Angle entre vitesse latérale et vitesse longitudinale
φ 1 : Fonction de limitation de flux
i+
2
PARTIE 1
INTRODUCTION
GENERALE
1
Partie 1 : Introduction
2
Partie 1 : Introduction
• Les phénomènes présents évoluent très rapidement (des variations de débit importantes
peuvent être observées à l’échelle de quelques minutes).
• Ils contiennent des confluences et des défluences.
• Le sens de l'écoulement peut s'inverser.
• Le régime d'écoulement est en permanence perturbé par des changements de pentes ou
de forme de conduites, des apports latéraux, etc.
• Les ouvrages hydrauliques (seuils, siphons, singularités qui entraînent des pertes de
charges, bassins et déversoirs d’orage) nombreux dans les réseaux demandent chacun
une modélisation particulière.
• L’écoulement dans une conduite peut se produire à surface libre ou en charge et passer
d’un type d’écoulement à l’autre très rapidement.
A l’heure actuelle, les équations de la mécanique des fluides (Navier –Stokes, Reynolds,
Barré de Saint-Venant, Euler, etc.…) ont donné naissance à de nombreux logiciels qui visent à
les résoudre. C’est le domaine de la C.F.D. [Versteeg and Malalasekera-1995]. Ce domaine
d’étude est très vaste. En effet, celui-ci concerne aussi bien les fluides visqueux que non
visqueux, compressibles aussi bien qu’incompressibles ou encore monophasique aussi bien que
multiphasique. De plus, ces logiciels et bien sûr les équations qui sont à leur base peuvent être à
caractère 1D, 2D ou 3D. En ce domaine, on peut aller du plus globalisant au plus fin. Tout
dépend du phénomène à décrire, de sa complexité physique mais également géométrique et
encore de la finesse des résultats attendus.
Dans le domaine de l’assainissement, les méthodes numériques issues de la CFD ne sont
pas utilisées. On en trouve des applications dans le domaine de l’hydraulique fluviale avec des
logiciels tels que RUBAR [Paquier-1995] ou TELEMAC qui pour le premier permettent de
modéliser le comportement d’une rivière et pour le second de modéliser le comportement de
grandes étendues d’eau soumise par exemple aux effets du vent et des marées. Ces logiciels
résolvent les équations de Barré de Saint-Venant 1D ou 2D. Pour l’étude plus précise de
phénomènes localisés on trouve des logiciels tels que FLUENT ou CFX qui résolvent les
équations de Reynolds. Ce sont des logiciels 3D qui tiennent compte des phénomènes de
turbulence et qui sont capables de prendre en charge des écoulements multiphasiques.
Pour l’assainissement, on trouve plusieurs logiciels dont, les plus employés sont SWMM (le
plus ancien et le plus employé car du domaine public), HYDROWORKS, MOUSE et CANOE.
Au niveau de l’hydraulique, ils sont tous bâtis sur la résolution des équations de Barré de Saint-
Venant 1D par des méthodes numériques classiques du type différences finies. Ces logiciels
contiennent des modules supplémentaires qui permettent de représenter certains phénomènes
spécifiques (mises en charges, singularités, ouvrages mobiles, déversoirs, etc..). Ces logiciels
montrent des limites dans leurs aptitudes à gérer les changements de régimes d’écoulements
car, dans ce cas, on peut voir apparaître de fortes discontinuités sur la ligne d’eau (ressaut
hydraulique), mais également dans leurs capacités à reproduire de manière correcte le
comportement de certains ouvrages tels que le déversoir d’orage. C’est dans le but de pallier
certaines de ces difficultés que nous avons entrepris ce travail.
3
Partie 1 : Introduction
2 PLAN DE LA THESE
Dans cette étude, nous nous intéressons aux équations de Barré de Saint-Venant en
régime transitoire et à la mise au point d’un modèle associé bien adapté à la problématique de
l’assainissement. Dans cette optique, le modèle doit être capable de :
Dans la première partie de ce travail, nous commençons par décrire les équations de
Barré de Saint-Venant 1D sous forme conservative dont la formulation inclue le ressaut
hydraulique et qui sont au cœur des modèles développés. Puis, nous faisons le bilan les divers
modèles simplifiés de Barré de Saint-Venant en discutant leurs domaines d’application et de
validité. Pour finir, nous décrivons quelques phénomènes et ouvrages qui peuvent apparaître
dans le réseau d’assainissement et que nous avons pris en compte dans nos modèles.
Dans la seconde partie, nous décrivons les méthodes de résolution classiques des
équations de Barré de Saint-Venant en mettant en évidence leurs limites dans le cadre d’une
application au domaine de l’assainissement. Puis, nous présentons, les méthodes numériques à
capture de chocs de type Godunov, possédant le caractère TVD (Total Variation Diminishing)
c’est à dire non oscillatoire, qui sont capables de résoudre les équations de Barré de Saint-
Venant en tenant compte des fortes discontinuités qui peuvent apparaître dans les variables
étudiées. En fin de partie, nous décrivons les équations et méthodes numériques que nous avons
mises en œuvre dans les modèles que nous avons développés.
Enfin, dans la troisième partie qui est la plus importante, nous présentons dans un
premier temps, les modèles de collecteurs et déversoirs que nous avons réalisés en précisant les
conditions initiales et aux limites applicables. Puis, nous décrivons les phases expérimentales
qui nous ont permis les valider.
4
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
PARTIE 2
ECOULEMENTS
DANS
LES
COLLECTEURS
ET
OUVRAGES
DU
RESEAU
D’ASSAINISSEMENT
5
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
1 INTRODUCTION
Dans cette partie, nous décrivons dans un premier temps le système d’équations de Barré de
Saint-Venant qui permet de représenter un écoulement à surface libre dans un collecteur ou
canal ainsi que les équations régissant le ressaut hydraulique par l’intermédiaire des conditions
de Rankine-Hugoniot.
Dans le but d’unifier dans un seul système d’équations les deux systèmes précédents nous
définissons la notion de solution faible d’un système d’équations et de condition d’entropie.
Ceci nous permettra de définir le système d’équations englobant à la fois l’écoulement à
surface libre en collecteur et le ressaut hydraulique.
Dans un deuxième temps, nous précisons le phénomène du passage en charge.
Pour finir, on détaille le comportement et les méthodes de calculs appliquées à divers ouvrages
tels que confluence, déversoir, seuil ou vanne.
Dans ce paragraphe, nous allons établir les équations de Barré de Saint-Venant qui
gouvernent les phénomènes d’écoulement (Avalanches, ruptures de barrages, écoulements en
rivières ou en canaux artificiels, mais encore comportement des grandes étendues d’eau).
Les équations de Barré de Saint-Venant établies en 1871 sont les équations les plus
utilisées pour modéliser les écoulements non stationnaires graduellement varié à surface libre.
Ces équations sont non linéaires et de type hyperbolique. Elles constituent en fait une
simplification des équations générales de la mécanique des fluides c’est à dire les équations de
Navier-Stokes. [Paquier-1995]. Ces équations peuvent être 1D, 2D ou 3D. Ceci dépend de la
complexité du phénomène que l’on veut décrire mais également de son échelle. Ainsi, dans le
domaine de l’assainissement, compte tenu de la complexité géométrique et hydraulique des
ouvrages tels que les bassins de stockage, déversoirs d’orages ou confluences, il serait
nécessaire de tenir compte des phénomènes de turbulence par l’intermédiaire d’une
modélisation 3D. Cette démarche est envisageable pour un ouvrage mais pas pour l’ensemble
du réseau. Le temps de calcul, la capacité des ordinateurs mais surtout la difficulté de
convergence des équations de Reynolds limite considérablement cette approche. C’est
pourquoi, dans le domaine de l’hydraulique urbaine, les modèles utilisés sont en général 1D.
Dans ce qui suit, nous n’allons pas déduire les équations de Saint-Venant des équations
de Navier-Stokes mais uniquement montrer comment elles peuvent être établies dans leur
forme 1D en distinguant équations continues et équations du ressaut hydraulique puis écrire
celles-ci sous forme conservative afin d’unifier les deux systèmes. L’intérêt de cette démarche
est d’avoir dans un seul système d’équations ces deux phénomènes hydrauliques afin d’être le
plus efficace dans la résolution des équations.
Les notations adoptées dans le cadre de l’établissement des équations de Barré de Saint-
Venant sont les suivantes :
6
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
• La distribution des pressions est considérée comme hydrostatique dans une section. Celle-ci
est obtenue grâce à la relation suivante :
h
Fpression hydrostatique = ρgI1 =ρg ∫ ( h-η ) Bdη (2.1.)
0
• Liquide incompressible.
• On effectue comme fait dans le cas des équations de Navier-Stokes, deux bilans [Vila-
1986] :
∂A ∂Q
+ =0 (2.2.)
∂t ∂x
7
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
On utilise un volume de contrôle Ω situé entre deux sections du canal ou collecteur situées aux
abscisses x et x+dx.
x
x+dx
Fond du canal
dx
d r r rr ∂ r r
∫ ρudω = ∫ ρu ( u.n ) dσ + ∫ ρudω = ∑ Fext (2.3.)
dt Ω ∂Ω Ω ∂t
Fext=Fpression+Fgravité+Ffrottement (2.4.)
• La force engendrée par l’accélération de la pesanteur (Fgravité) agit sur Ω et est donnée par :
avec :
S0 ( x ) sin ( α )
• La force de frottement (Ffrottement) est due à une contrainte de cisaillement et agit sur le
périmètre mouillé.
uu uu uu
Ffrottement = ∫ -τ 0 dσ = ∫ -g dσ = -g Pdx = -Ag dx = -gASf dx (2.6.)
Σ Σ C C CR h
Exprimé sous la forme précédente Sf est donné par une relation de type Chézy. On
rencontre plus communément ce terme exprimé sous la forme de la relation de Manning-
Strickler :
8
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Q2
Sf = 4
(2.7.)
K S2 A 2 R h 3
Cette expression suppose une hypothèse fondamentale qui consiste à agréger dans un même
mode de calcul les pertes d’énergie dues aux frottements contre les parois et celles dues à la
viscosité du fluide. Lorsque le flux se partage en plusieurs branches, ces deux paramètres
jouent des rôles différents : la viscosité contrôle la cohérence de l’écoulement et la rugosité
l’adhérence aux parois
Or : p ≡ p(η) = ρg (h − η)
∂h
Fpression = −ρg Adx (2.9.)
∂x
r rr ∂ (ρQu )
∫ ρu (u.n )dσ = ∫ div((ρu ) ⊗ u )dω = dx (2.10.)
∂Ω Ω ∂x
et
∂ ∂ (ρQ )
∫ Ω ∂t
ρudω =
∂t
dx (2.11.)
On arrive à :
∂Q ∂ (Qu ) ∂h
+ + gA = gA(S 0 (x ) − S f ) (2.12.)
∂t ∂x ∂x
∂A ∂Q
+ =0
∂t ∂x
(2.13.)
∂Q + ∂ (Qu ) + gA ∂h = gA(S − S )
∂t ∂x ∂x
0 f
C’est le système d’équations hyperboliques le plus classiquement utilisé dans la plupart des
logiciels de simulation dans le cadre de l’hydrologie urbaine (SWMM, CANOE, MOUSE,
HYDROWORKS).
9
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Le ressaut hydraulique est une surélévation brusque de la surface libre d’un écoulement
permanent qui se produit lors du passage du régime torrentiel au régime fluvial. Ce passage
peut être dû à une rupture de pente ou à la présence d’un ouvrage comme le déversoir d’orage
sur le réseau. Il est accompagné d’une agitation marquée avec entraînement d’air et de grandes
pertes d’énergie.
=Perte de charge
Charge spécifique=
Pour Fr compris en 1.7 et 2.5, on constate le même phénomène mais plus accentué.
Dans ce cas, il se produit des petits tourbillons superficiels.
10
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Pour Fr entre 2.5 et 4.5 l’écoulement est pulsatoire. La plus grande turbulence se vérifie
soit près du fond soit à la surface. Chaque pulsation produit une onde de période irrégulière.
Cette onde peut se propager sur une très grande distance.
11
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Qapp = 0 (2.14.)
Or :
Qapp = ρA ( s - u )
⇒ ρA ( s - u ) = 0
D’où :
[Q] = s [ A ] (2.15.)
où :
[-] : désigne saut d’une quantité notée « - » de part et d’autre de la discontinuité.
u : correspond à la vitesse moyenne de l’écoulement.
s : correspond à la vitesse de propagation du ressaut hydraulique.
uQapp = Fpression
⇒ uA ( s-u ) = Fpression
⇒ s [ Q ] = Fpression + uQ (2.16.)
Il est nécessaire de préciser les forces de pression qui agissent sur les frontières du
volume de contrôle étudié :
rr
Fpression = − ∫ [Link]σ = Fp1 + Fp2 (2.17.)
Σ = ∂Ω
où Fp1 agit en x et x+dx sur la section mouillée et Fp2 agit sur le périmètre mouillé. La
force de pression est décomposée en deux termes dont le premier (Fp1) correspond au terme de
pression hydrostatique. Celle-ci est due à la hauteur de colonne d’eau au-dessus de la position
verticale η. Le second terme (Fp2) exprime qu’une modification de la section utile aura pour
conséquence de gêner ou de faciliter l’écoulement. Cette variation est prise en compte par une
variation de largeur et donc du périmètre mouillé du canal. Nous tennons compte de la
composante de cette force uniquement dans la direction de l’écoulement.
∂ ∂ h
Fp1 = − ∫ ρg (h − η)dσ dx = − ∫ ρg(h − η)B(x, η)dη dx (2.18.)
∂x A ∂x 0
h ∂
Fp2 = + ∫ ρg(h − η) B(x, η)dη dx (2.19.)
0 ∂x
12
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
En ajoutant Fp1 et Fp2 on retrouve l’expression de Fpression donnée dans le cas des
équations continues. Nous illustrons les surfaces d’actions de ces deux forces sur la figure
suivante :
Fp2 Fp1
Dans le cas du ressaut, le terme [-] sur Fp2 n’apporte pas de contribution. Puisque le ressaut est
supposé localisé dans une section donc les largeurs à gauche et à droite du ressaut sont
supposées constantes. Il ne reste que [Fp1] et donc :
avec :
h
I1 (A, x ) = ∫ (h − η)B(x, η)dη (2.21.)
0
On considère un problème aux valeurs initiales dans le cas scalaire pour t>0.
u t +f ( u ) x = 0, u ( x,0 ) = u 0 ( x ) (2.22.)
Nous allons généraliser la notion de solution pour les équations de la forme précédente.
Supposons pour le moment que u est une solution classique de l’équation (2.22.). Soit C10 la
classe de fonctions C1 φ qui disparaissent en dehors d’un ensemble compact en t ≥ 0 c’est à dire
(spt φ)∩(t ≥ 0) ⊆ D) où D est le rectangle 0 ≤ t ≤ T, a ≤ x ≤ b et choisies telles que φ=0 en
dehors de D et sur les lignes t=T, x=a et x=b.
13
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Support de φ
x=a x=b x
b T
∫∫ ( u +f ) φ dxdt = ∫∫ ( u +f ) φ dxdt = ∫ ∫ ( u +f ) φ dxdt = 0
t>0
t x
D
t x
a 0
t x
∫∫ ( uφ +f ( u ) φ ) dxdt+∫
t ≥0
t x
t=0
u 0φdx=0 (2.23.)
Nous avons ainsi montré que si u est une solution classique de (2.22.) alors (2.23.) est vérifiée
pour toute fonction φ ∈ C10 .
Définition 10 :
Une fonction mesurable et bornée u (x , t ) est appelée solution faible du problème aux valeurs
initiales (2.22.) avec des valeurs initiales bornées et mesurables u 0 si la relation (2.23.) est
vérifiée pour toutes fonctions test φ ∈ C10 .
Nous pouvons également montrer que tous les types de discontinuités ne sont pas
admissibles [Smoller-1984]. En fait, la condition (2.23.) apporte de sévères contraintes sur les
courbes de discontinuité. Soit Γ une courbe régulière au passage de laquelle u a une
discontinuité de saut c’est à dire que u a des limites bien définies sur les deux cotés de Γ et que
u est régulière ailleurs que sur Γ. Soit P un point quelconque sur Γ et soit D une boule centrée
sur P (cf. figure suivante). Nous supposons que dans D, Γ est donné par x = x (t ) . Soient D1 et
D2 les composantes de D qui sont séparées par Γ. Soit φ ∈ C10 et si nous supposons cette
fonction nulle en t=0 ; alors de (2.23.) nous tirons :
∫∫ (uφ
D
t + fφx )dxdt = ∫∫ (uφ t + fφx )dxdt + ∫∫ (uφ t + fφx )dxdt =0
D1 D2
14
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
D Γ
D1 P
Q2
Q1
D2
En utilisant le fait que u est de classe C1 dans Di, le théorème de la divergence donne :
Tant que φ=0 sur ∂D ces intégrales curvilignes sont nulles sur Γ. Donc si u l = u (x (t ) − 0, t ) et
u r = u (x (t ) + 0, t ) alors on a :
Ainsi,
0 = ∫ φ(− [u ]dx + [f (u )]dt )
Γ
s[u ] = [f (u )] (2.24.)
15
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Si nous nous intéressons aux solutions faibles d’un système d’équations hyperbolique mis sous
forme matricielle :
∂U ∂F(U, x )
+ = G (U, x ) avec U ∈ R m (2.25.)
∂t ∂x
Une solution discontinue de part et d’autre d’une ligne se déplaçant à la vitesse s dans le plan
(x,t) est caractérisée par la condition d’entropie qui est une solution faible physiquement
admissible :
Les équations de Saint-Venant peuvent-elles se mettre sous cette forme (de manière à ce que
l’on retrouve les conditions de Rankine-Hugoniot à partir des équations continues) ? Il s’agit de
mettre les équations (2.13.) sous la forme (2.25.), afin que (2.15.)et (2.20.) donnent (2.26.).
∂Q ∂ Q 2
+ + gI1 = G (U, x ) (2.27.)
∂t ∂x A
Nous calculons G de manière à vérifier (2.13.) :
Il vient :
∂I1 ∂h
L(U, x ) = g − gA (2.29.)
∂x ∂x
Nous introduisons en fait une variation de I1 par rapport à x qui en fait est due à la variation de
la géométrie du canal et en particulier de la largeur. En fait cette variation de pression due à la
variation longitudinale de la géométrie et donc de la section mouillée représente I2. Ceci est
montré dans ce qui suit.
Or :
∂I1 ∂h
h
∂ 2 A ( x,η )
g = gA + ∫ g ( h − η ) dη (2.30.)
∂x ∂x 0 ∂x∂η
Remarque :
La relation (2.30.) est établie en utilisant le théorème de Leibnitz qui permet la dérivation sous
le signe d’intégration. Soit :
16
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
u2 (α)
φ(α) = ∫ f ( x,α ) dx
u1 ( α )
u (α)
dφ 2 ∂f du du
= ∫ dx+f ( u 2 ,α ) 2 -f ( u1 ,α ) 2
dα u1 ( α ) ∂α dα dα
Il vient ainsi :
h
∂ 2 A ( x,η ) h
∂B ( x,η )
L ( U,x ) = ∫ g ( h − η ) dη = ∫ g ( h − η ) dη = gI 2 (2.31.)
0
∂x∂η 0
∂x
∂A ∂Q
∂t + ∂x = 0
(2.33.)
∂Q + ∂ Q + gI1 = g(S 0 − S f ) + gI 2
2
∂t ∂x A
Cette forme est compatible avec les équations du ressaut mobile (2.15.) et (2.20.).
Par commodité d’écriture, ce système peut être mis sous la forme matricielle suivante :
U t + F(U )x = G (U ) ⇔ U t + J (U )U x = G (U )
avec :
A Q 0
U = ,
F(U ) = Q 2 , G (U) =
Q
A
+ gI1
gI 2 + gA ( S0 − Sf )
U est appelé vecteur écoulement, F(U) vecteur flux et G(U) vecteur source.
17
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Dans le domaine de l’assainissement, les modèles qui ont été mis au point peuvent être
regroupés en deux grandes familles :
¾ Onde quasi-permanente,
¾ Onde,
¾ Onde cinématique,
¾ Onde simple ou onde de gravité.
Pour plus de détails concernant ces divers types de modèles on pourra consulter [Buyer-1999].
2.6 Choix d’un modèle pour une modélisation fine des réseaux
Le choix d’un modèle mécaniste se fait en général en fonction de la finesse attendue sur
les résultats et d’autre part de l’ordre de grandeur des différents termes des équations de BSV.
Il faut comparer les ordres de grandeurs des termes d’inertie et de frottement. Par exemple, un
modèle d’onde cinématique ou diffusante pourra être choisis dans le cas d’écoulement à
frottements prépondérants [Moussa–1996b]. A l’inverse, on pourra prendre un modèle d’onde
simple lorsque la pente du bief est faible c’est à dire lorsque les termes de frottements sont
négligeables.[Kovacs–1988].
Il existe des méthodes quantitatives d’aide au choix du degré de simplification des
équations de Barré de Saint-Venant à utiliser pour modéliser les écoulements en rivière. L’une
d’entre elles est basée sur l’analyse en petites perturbations de la propagation d’onde et
introduit un critère de choix faisant intervenir deux paramètres sans dimension qui sont le
nombre de Froude de l’écoulement moyen et la période du signal d’entrée.[Moussa–1996b].
Les modèles simples ont l’avantage de nécessiter moins de paramètres et par conséquent
moins de mesures que les modèles plus complexes d’où une préférence pour les modèles
d’onde diffusante, cinématique ou linéaire pour de nombreux problèmes. Mais ces modèles ne
sont pas valables dans le cas d’écoulements qui font intervenir des ondes se propageant vers
l’amont ou lorsque la profondeur et la vitesse de l’écoulement changent rapidement.
De plus, si l’on cherche à obtenir une vision globale du comportement du réseau
l’utilisation d’un modèle conceptuel est suffisante.
18
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Les écoulements dans les systèmes de collectes des eaux usées ou de pluie montrent les
particularités suivantes :
• L’écoulement en réseau d’assainissement est de type « transcritique » c’est à dire que l’on
peut observer dans un même tronçon les régimes d’écoulements fluvial et torrentiel. Ces
transitions peuvent être dues à des changements de pentes entre les collecteurs, à la
présence d’ouvrages ou encore être provoquées par un apport d’eau subit dans le réseau lors
d’un événement pluvieux. Le passage du régime fluvial au régime torrentiel ou inversement
peut être détecté par évaluation du nombre de Froude. Cette particularité doit être prise en
compte. En effet, de celle-ci dépend la localisation des conditions aux limites qui seront
utilisées lors du calcul des débits et hauteurs d’eau dans les collecteurs.
• Le passage du régime torrentiel au régime fluvial peut engendrer un ressaut hydraulique qui
doit être localisé et pris en compte lors de la résolution numérique des équations de Barré
de Saint Venant afin de déterminer avec précision les débits et hauteur d’eau en amont et en
aval de celui-ci.
Compte tenu de ces particularités et en vue d’obtenir une solution fine pour le
comportement de l’eau dans ce système nous avons décidé d’utiliser le modèle de Barré de
Saint-Venant pris sous sa forme conservative. On rappelle ci-dessous la forme matricielle
de ces équations :
U t + F(U )x = G (U )
avec :
A Q 0
U = ,
F(U ) = Q 2 , G (U) =
Q
A
+ gI1
gI 2 + gA ( S0 − Sf )
19
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Hewitt et Hall-Taylor (1970) ont distingué huit régimes possibles pour des écoulements
mêlant gaz et liquide :
• Ecoulement stratifié dans lequel la phase liquide est en-dessous de la phase gazeuse (a).
• Ecoulement ondulé qui possède une interface ondulée entre phase liquide et phase
gazeuse (b).
• Ecoulement en bouchon avec une surface de nature très ondulée qui atteint le fond du
tuyau et qui sépare ainsi la phase gazeuse en cellules indépendantes ( c )
• Ecoulements en bulles avec des bulles et des poches de gaz qui sont toutes distribuées
sur la partie supérieure de la conduite (d).
• Ecoulements en gouttes avec une distribution quasi uniforme de bulles de gaz dans le
phase liquide (e)
• Ecoulement annulaire avec une large portion de gaz qui pousse le liquide (f)
• Ecoulement en spray (g)
• Ecoulement en mousse (h)
Dans ces deux derniers cas (g et h) l’écoulement est de nature diphasique. Le mélange eau-air
est homogène avec dans le premier cas (g) de fines bulles d’air clairement séparée de l’eau et
dans le second cas (h) une émulsion formée par les deux phases dans laquelle l’eau et l’air ne
peuvent plus clairement être distinguées.
Ces divers types de transition d’un écoulement à surface libre à un écoulement en charge sont
représentées sur la figure suivante :
20
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Dans ce qui suit, je présente quelques éléments bibliographiques dus à Hager [Hager-
1999] qui permettent de prendre conscience de la similitude qui existe entre un écoulement en
charge et un écoulement à surface libre.
Figure 8 : Distributions de pression et de vitesse dans le cas d’un écoulement en charge (a)
et d’un écoulement à surface libre (b)
En principe, les deux types d’écoulements peuvent être traités de façon similaire
excepté que :
• Pour une conduite sous pression, la ligne piézométrique ne coïncide pas avec la
génératrice supérieure de la conduite.
• Pour un écoulement en canal ouvert, la surface libre n’est pas connue à l’avance.
Les processus mécaniques vont de pair avec des phénomènes de dissipation. Par
conséquent, il n’est pas possible de maintenir l’énergie initiale plus d’un certain temps et au-
delà d’une certaine distance sans dépenser d’énergie additionnelle. Ces pertes d’énergie ou
encore « pertes de charges » sont dues aux phénomènes de dissipation. Ces phénomènes sont
principalement dus :
• Aux frottements des couches de liquide sur les parois des collecteurs. Ce phénomène est
étroitement lié à l’état de surface du collecteur et à la viscosité du fluide. Ces pertes par
frottement sont dues à des contraintes de cisaillement entre les couches de liquide.
• A la présence de singularités coudes, vannes, seuils, entonnements. Ce phénomène est dû à
une modification de la géométrie dans la direction de l’écoulement et donc de la section
mouillée. A chaque fois que l’écoulement est accéléré ou décéléré des contraintes de
cisaillement internes apparaissent et l’énergie mécanique de l’écoulement est diminuée.
Ces pertes d’énergie se manifestent par des remous ou des mouvements de type vortex et ne
peuvent être évaluées qu’indirectement
Dans le but de répondre à cette question on écrit dans un premier temps les équations
suivantes en conduite fermée et pour le cas d’un canal à surface libre.
21
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Figure 9 : Lignes de charges dans le cas d’un écoulement en charge (a) et d’un écoulement
à surface libre (b)
u2 dH
H = hp + , = So − Sf (2.34.)
2g dx
u2 dH
H=h+ , = So − Sf (2.35.)
2g dx
La section mouillée A pour une conduite sous pression est donnée par A=A(x) et dans le
cas à surface libre la section mouillée en la position x=x* varie avec la hauteur d’eau et
A=A(x,h).
Q
Si u = en différentiant les deux équations précédentes, on obtient :
A
Pour le cas en charge :
dh p Q 2 dA
− = S0 − Sf (2.36.)
dx gA 3 dx
dh Q 2 dA
− = S0 − Sf (2.37.)
dx gA 3 dx
Pour une géométrie de conduite donnée c’est à dire connaissant une fonction A(x), une
inclinaison de la conduite S0(x), le débit Q et la pente énergétique Sf l’équation pour
l’écoulement en charge peut être résolue pour la distribution de pression locale hp. Au contraire,
la relation concernant la surface libre ne peut pas être résolue car A dépend également de la
hauteur d’eau.
dA ∂A ∂A dh
= + (2.38.)
dx ∂x ∂h dx
22
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
dh
dx
[ ]
1 − Fr 2 −
Q 2 ∂A
gA 3 ∂x
= S0 − Sf (2.39.)
La similitude entre les deux équations encadrées montre que l’écoulement en charge et
l’écoulement à surface libre sont identique à la condition que :
Q 2 ∂A U
Fr =
2
= =0 (2.40.)
gA 3 ∂h c
Ceci signifie que le nombre de Froude doit être nul. Effectivement, dans une conduite
en charge, les informations hydrauliques transitent « instantanément » vers l’amont et l’aval de
la conduite. Ceci n’est pas le cas dans en hydraulique à surface libre.
Ceci montre bien que l’écoulement en charge est un cas particulier d’un écoulement à
surface libre et que les coefficients de perte de charge (obtenus pour les écoulements en charge)
peuvent être appliqués au cas correspondant d’un écoulements à surface libre tant que le
nombre de Froude reste petit. La limite supérieure est 1 c’est à dire lorsque le terme de pression
dans l’équation à surface libre disparaît. En pratique le concept de transposition du coefficient
de perte de charge ne doit pas être appliqué au-delà de Fr=0.7 [Hager-1999].
23
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
ΖΟΝΕ 6
ΖΟΝΕ 5
ΖΟΝΕ 4
ΖΟΝΕ 3
ΖΟΝΕ 2
ΖΟΝΕ 1
Zone 1 : Il s’agit d’une fente inférieure qui possède une faible épaisseur. Cette fente permet au
logiciel de prendre en compte l’altimétrie des ouvrages et de modifier la représentation de la
ligne d’eau. Elle permet également de traiter le problème du fond sec.
Zone 3 : Elle correspond à la section normale du tronçon ou la relation A=f(h) est régulière.
Dans cette zone, l’écoulement fonctionne à surface libre sans anomalies susceptibles
d’influencer les grandeurs hydrauliques.
Zone 5 : La fente supérieure de la conduite qui permet au modèle d’assimiler les écoulements
en charge à des écoulements à surface libre. Avec cet artifice de calcul, les écoulements en
charge peuvent être traités par le même modèle que les écoulements à surface libre.
Cette zone peut être utilisée pour tenir compte des volumes de stockage qui existent
effectivement dans les systèmes d’assainissement entre le toit de la conduite et le niveau du sol
(cheminées, branchements, avaloirs).
Zone 6 : Cette zone est créée artificiellement pour représenter le phénomène de débordement et
la quantité d’eau débordée.
Une étude expérimentale permettant la comparaison des résultats obtenus sur un banc
d’essai physique avec ceux obtenus grâce à un logiciel de modélisation basé sur la résolution
des équations de barré de Saint-Venant 1D sous forme conservative, en introduisant l’artifice
numérique de la fente de Preismann et utilisant un schéma numérique à capture de choc comme
ceux décrits dans la partie 3 de cette thèse a été menée par B. Trajkovic et al.[ Trajkovic-1999].
Les pentes étudiées étaient de 1.4 et 2.7 %. La comparaison des résultats expérimentaux et
numériques semble fournir des résultats proches. Des problèmes numériques apparaissent
lorsque la vanne permettant le passage de l’écoulement en charge à l’écoulement à surface libre
est manipulée trop vite.
24
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Compte tenu des instabilités numériques que peut générer la fente de Preismann, nous nous
sommes inspirés des similitudes entre les équations des écoulements en charge et à surface
libre. Pour cela, dans un premier temps, nous avons choisi de limiter la section A à la pleine
section. Dans un deuxième temps, nous prenons en compte le fait que dans une conduite en
charge les informations hydrauliques transitent « instantanément » vers l’amont et l’aval.
A
La difficulté réside dans le calcul de la célérité c = gD h = g qui devrait tendre vers
B
l’infini.
Pour rendre cette approche valide, il faut que c prenne la valeur la plus grande possible
( c U ). Pour se faire, on choisit la valeur de la largeur au miroir B de telle sorte que c soit la
plus grande du système modélisé au point de mise en charge.
4 LA CONFLUENCE
4.1 Introduction
1. L’écoulement est entièrement fluvial dans les tronçons respectivement amont, latéral et aval
(FLUV-FLUV-FLUV).
2. L’écoulement est fluvial dans les branches amont et latérale et devient torrentiel à l’aval
(FLUV-FLUV-TOR).
3. L’écoulement est torrentiel dans l’une des branches amont ou latérale et fluvial dans le
tronçon aval. (FLUV-TOR-FLUV) ou (TOR-FLUV-FLUV). La formation d’un ressaut
hydraulique est observée dans l’une des branche latérale ou amont.
4. L’écoulement est fluvial dans l’une des branches amont ou latérale mais devient torrentiel à
l’aval. (FLUV-TOR-TOR) ou (TOR-FLUV- TOR)
5. L’écoulement est torrentiel dans les branches amont et latérales et devient fluvial dans la
branche aval. (TOR-TOR-FLUV)
6. L’écoulement est torrentiel dans les trois branches (TOR-TOR-TOR).
25
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Dans la littérature, on trouve presque exclusivement des études pour le cas tout fluvial ou le cas
tout torrentiel. Les cas transcritiques sont très peu étudiés. Pour ce type d’étude, la phase
expérimentale afin de formuler, puis de valider les relations théoriques est indispensable.
De cette manière, il paraît difficile de mettre au point un traitement unique de la
confluence pour toutes les configurations et régimes d’écoulements possibles puisqu’elles
nécessitent une anticipation du comportement hydraulique ce qui peut s’avérer hasardeux.
Le cas le plus souvent rencontré dans les cours d’eau naturels est caractérisé par un
régime fluvial dans toutes les branches. Les premières études concernant les confluences ont
d’ailleurs été menées pour l’étude de cette configuration.
Dans presque toutes les études réalisées seul le canal latéral forme un angle par rapport
à l’axe formé par les deux canaux amont et aval.
L’une des premières étude répertoriée concernant la confluence est due à Taylor qui en
1944, supposa que :
1. Les largeurs des trois branches étaient égales
2. Une pente de radier nulle
3. Des frottements négligeables
4. Une répartition hydrostatique des pressions et une distribution des vitesses uniforme
5. Des lignes de courants parallèles aux parois latérales autour du volume de contrôle
6. Des hauteurs d’eau égales dans les conduites amont et latérale
7. Des forces de pression latérales négligeables.
26
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
configurations ont été étudiées pour des nombres de Froude inférieurs à 0.6 et un rapport entre
le débit dans le canal amont et celui du canal aval compris entre 0 et 1.
P*
En particulier, Hager, dans son livre [Hager-1999] mène une étude qui permet de
déterminer les hauteurs d’eau dans les tronçons connaissant les débits amont et latéral et une
hauteur d’eau. Celle-ci est basée sur la conservation de l’énergie, le principe fondamental de la
dynamique et la conservation des débits.
Ici, on se contentera de citer les résultats obtenus par Hager. Les trois canaux qui
constituent la confluence ont la même largeur b. µ correspond à un coefficient de contraction
du à la présence de la conduite latérale.
Q 2u Ql2 Qc2
pu + Qu + pl + Ql = pc + Qc (2.41.)
2gb 2 2gb 2 2gµ 2 b 2
Les indices u, l, c et d sont respectivement associés aux canaux amont, latéral, contracté
et aval.
Q 2u Q 2 cos δ Q2
pu b + + p l b cos δ + l = p d b + d + p∗ b cos δ (2.42.)
gb gb gµb
27
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
p l + npc
p* = avec 0≤n<∞ (2.43.)
1+ n
p* ( n = 0 ) = p l et p* ( n → ∞ ) = p c
Les distributions de pressions dans les canaux amont et latéraux doivent être identiques
pour assurer la continuité des pressions à travers le point de confluence.(pu=pl). Ce résultat a été
montré par Visher en 1958.
Qc = Qd = Qu + Ql (2.44.)
Ainsi, on peut éliminer les termes pu, pl et pc pour obtenir le coefficient de contraction µ
Q
comme une fonction du rapport des débits q = u , de l’angle de la jonction δ et du coefficient
Qd
n. On obtient :
n cos δ 2 n cos δ
1
−1
n cos δ 2 2
−1
µ = 1 ± 1 − 1 +
1+ n
2q + 2 (1 − q )
2
cos δ − (1 − 3q + 3q ) 1 + n 1 + n
1 + 1 +
D’autre part, les conditions suivantes doivent être respectées par la relation précédente :
Pour q=1 c’est à dire pas de débit latéral, le coefficient de contraction vaut 1,
Pour δ=0 c’est à dire des branches parallèles, le coefficient de contraction est
µ(δ = 0 ) = 1 .
La première condition est toujours respectée alors que la seconde conduit à prendre une valeur
n=1/2 pour le coefficient. En insérant cette valeur dans l’expression précédente on obtient pour
le coefficient de contraction :
1
2 −1
2 1 2 1 2 1
µ −1
= 1 + (1 − q )(2 − q )1 − cos δ − cos δ + cos δ 1 + cos δ (2.45.)
3 3 9 3
Les coefficients de perte de charges pour les branches amont et latérales sont données
par :
ξ u = (µ −1 − 1) − 1 + 3q − 2q 2
2
(2.46.)
ξ l = (µ −1 − 1) + q − 2q 2
2
28
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Hager a également donné une expression des coefficients de perte d’énergie dans le cas
d’un confluence de forme courbe.
hu − hd =
(1 + ξu ) Qd2 − Q2u
2gSd2
(2.47.)
hl − hd =
(1 + ξl ) Qd2 − Ql2
2gSd2
Remarque : Si la largeur de la branche aval est trop faible par rapport à la somme des deux
largeurs en amont de la jonction, une transition du fluvial au torrentiel apparaît à l’entrée de la
branche aval.
C’est Bowers (1950) qui le premier, analysa un écoulement torrentiel dans une jonction
de canaux. En fonction de la géométrie de la jonction et de la valeur du nombre de Froude dans
le canal amont un ressaut hydraulique apparaissait ou non dans l’un des confluents. Dans le cas
ou aucun ressaut n’était présent, des ondes croisées similaires à celles présentes dans les
rétrécissements apparaissaient. Dans tous les cas, Bowers conseillait pour le dimensionnement
d’adopter des hauteurs de canaux bien plus grandes que les hauteurs usuelles pour
dimensionner ce genre de confluences.
Une autre étude associée aux écoulements torrentiels est due à Schnitter et al. (1955). Il
étudia une configuration dans laquelle l’angle de confluence était pratiquement nul mais pour
laquelle le canal latéral avait une pente très forte. Il constata la présence d’ondes croisées très
fortes dans la branche aval qui étaient encore amplifiées dans le cas où les deux branches amont
avaient de nombre de Froude différents. Néanmoins, l’installation d’un mur de séparation dont
l’une des extrémités se situait au niveau du point de confluence et de hauteur décroissante dans
la branche aval permettait de faire quasiment disparaître les ondes croisées.
Behlke et Pritchett (1966) menèrent une analyse des écoulements torrentiels dans les
jonctions de canaux trapézoïdaux et rectangulaires pour des angles de confluences de 15°, 30°
et 45° et des nombres de Froude compris entre 2 et 7. Ils observèrent deux sauts obliques dont
les origines se situaient au point de confluence. La hauteur de ces deux sauts au voisinage du
29
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
point de confluence est la même. Ils supposèrent que l’effet de la branche latérale sur
l’écoulement dans la branche principale est comparable à la perturbation causée par la présence
d’un mur oblique. Un effet similaire est attribué à l’écoulement dans la branche latérale. Ainsi,
ils déterminèrent les angles que font les sauts par rapport à la direction principale mais ils ne
tentèrent pas de déterminer la hauteur maximale de ces sauts.
Un modèle pour la prédiction de l’angle de l’onde principale dans une jonction simple
est du à Greated (1968) qui étudia un écoulement dans une jonction à 60° horizontale pour des
nombres de Froude compris entre 6 et 11. Il n’essaya pas de calculer la hauteur des ondes
suggérant qu’une étude expérimentale beaucoup plus approfondie était nécessaire en raison de
l’influence prépondérante de la forme des canaux. Sa formulation a été reprise par Hager
[Hager-1989b] qui s’intéressa non seulement à la localisation de ces ondes mais également à la
hauteur maximale qu’elles pouvaient atteindre. Celle-ci peut atteindre à l’opposé du canal
latéral 10 à 20 fois la hauteur d’eau dans le canal amont.
Cette forte augmentation de la hauteur d’eau peut engendrer de gros problèmes dans les
conduites fermées en raison de l’apparition de phénomènes tels que pulsations, entraînement
d’air, transition vers une mise en charge ou apparition de ressauts hydrauliques.
30
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
∂ (uh ) ∂ (vh )
∂x + ∂y = 0
u ∂u v ∂u ∂h
+ + = S 0 x − S fx (2.48.)
g ∂x g ∂y ∂x
u ∂v v ∂v ∂h
+ + = S 0 y − S fy
g ∂x g ∂y ∂y
Les pertes par frottement sont exprimées par l’intermédiaire de l’équation de Manning-
Strickler :
u (u 2 + v 2 ) v(u 2 + v 2 )
1 1
2 2
S fx = , S fx = (2.49.)
(1 n )2 h (1 n )2 h
4 4
3 3
∂h ∂ (hu ) ∂ (hv )
∂t + ∂x + ∂y = 0
∂ (hu ) ∂ 2 1 ∂ (uvh )
+ u h + gh 2 + =0 (2.50.)
∂ t ∂x 2 ∂y
∂ (hv ) ∂ (uvh ) ∂ 2 1
+ + v h + gh 2 = 0
∂t ∂x ∂y 2
U t + F(U )x + G (U )y = 0
(2.51.)
hu
h hv
1 2
avec : U = hu , F (U ) = hu + gh et G (U ) =
2
huv
2 1 2
hv v 2
h + gh
huv 2
Pour l’étude des changements de sections, les pentes de fond et pertes par frottement sont prises
en compte [Rahman-1997] en ajoutant un vecteur source au second membre. Le système de
Barré de Saint-Venant devient :
31
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
U t + F ( U ) x + G ( U ) y = S(U)
hu
h
1 2
avec: U = hu , F ( U ) = hu + gh
2
, (2.52.)
hv 2
huv
hv 0
G ( U) = huv et S(U) = gh(S0x − Sfx )
2 1 2 gh(S0y − Sfy )
v h + gh
2
u (u 2 + v 2 ) v(u 2 + v 2 )
1 1
2 2
S fx = , S fx = (2.53.)
(1 n ) (1 n )
2 4 2 4
h 3
h 3
Les matrices jacobiennes obtenues dans le cas du système sans second membre sont les
suivantes :
0 1 0
∂F(U ) 2
A(U ) = = − u + gh 2u 0
∂U − uv
v u
(2.54.)
0 0 1
∂G (U )
B(U ) = = − uv v u
∂U − v 2 + gh 0 2v
λ1 = u et λ 2 = λ 3 = u ± gh (2.55.)
λ1 = v et λ 2 = λ 3 = v ± gh (2.56.)
32
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Tronçon amont 1
α1 Tronçon aval
α2
Tronçon amont 2
Q2 Q 22 Q32
cos α1 1 + gI1 ( h ) + cos α 2 + gI 2 ( h ) = + gI3 ( h ) (2.57.)
S1 ( h ) S2 ( h ) S2 ( h )
Bi h 2
Ii ( h ) = pour i ∈ {1, 2,3} (2.59.)
2
33
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
5 LES DEVERSOIRS
On s’intéresse à présent aux différentes approches existantes du fonctionnement
hydraulique d’un déversoir d’orage latéral. Dans ce paragraphe, on expose dans un premier
temps le fonctionnement hydraulique de l’ouvrage. Puis, on présente quelques relations
empiriques classiquement utilisées pour évaluer les débits déversés. Ensuite, on évoque le
principe du logiciel DO qui est un outil d’aide au dimensionnement des déversoirs assez
répandu dans les DDAF. Les deux paragraphes suivants détaillent les méthodes développées,
pour le premier, sur le concept de conservation de l’énergie spécifique constante, et pour le
second, sur l’équation de la quantité de mouvement. Ces deux méthodes permettent l’obtention
des lignes d’eau ainsi que les débits déversés. Pour finir, on détaille les expressions des
variations débits déversés associées au deux formulations précédentes.
hn1 hn1
34
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Dans le cas du déversoir non prismatique le nombre de forme de lignes d’eau susceptibles
d’apparaître dans le déversoir est encore plus important. En effet, la figure ci-dessous montre
que dans le cas d’un canal non prismatique (θ<0) en fluvial la ligne d’eau peut non seulement
monter comme dans le cas prismatique mais également baisser. Dans le cas torrentiel, elle peut
baisser mais également monter si le canal n’est plus prismatique (θ<0).
Figure 14 : Types de profils de surface pour des nombres de froude F<1 et F>1 dans les
cas prismatique (θ>0) et non prismatique (θ<0)
Initialement, les débits déversés par l’intermédiaire des déversoirs d’orages ont été
évalués à travers l’utilisation de relations empiriques. Ces équations sont toutes bâties à partir
de résultats expérimentaux. On trouve par exemple les formules de Engels (1917), de Coleman
et Smith (1923), de Balmaceda et Gonzales (1930) ou encore de Dominguez (1945) qui
permettent le calcul du débit déversé en fonction des valeurs de hauteur d’eau à l’amont et/ou à
l’aval du déversoir [El Khashab-1975]. Ces relations ne sont applicables que pour certains
types d’écoulement et uniquement pour certaines géométries de déversoir.
On trouve beaucoup de formules empiriques pour les seuils frontaux. Pour les seuils
latéraux à crête basse, ces formules sont beaucoup plus rares. Ils en existent cependant
quelques-unes unes, valables moyennant un certain nombre d’hypothèses :
35
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Cependant, dès que la longueur du seuil devient importante, il est nécessaire de tenir
compte des variations de la ligne d’eau. De plus, un régime fluvial à l’amont et à l’aval du
déversoir ne signifie pas qu’il le reste dans le déversoir. Enfin, toutes ces formules contiennent
un coefficient déterminé expérimentalement ; il permet de tenir compte des approximations
choisies mais fait toujours référence à une géométrie et à des conditions d’écoulement précises.
Or, il est impossible de trouver une formule pour chaque cas de figure ; c’est pourquoi on
cherche à modéliser le fonctionnement par des équations plus fondamentales pour disposer
d’outils plus généraux.
5.3 Approche du dimensionnement de déversoir par Giersch
Le logiciel DO a été crée par la DDAF du Bas-Rhin en 1985 [Giersch-1985] pour aider
les ingénieurs et les techniciens au dimensionnement et au diagnostic des déversoirs d’orage et
il est encore aujourd’hui un des seuls outils à disposition.
Il permet le calcul des seuils latéraux à crête haute et basse. Malheureusement, cet outil
est affecté par plusieurs imprécisions dans la démarche. Nous citerons qu’il n’est pas capable
de traiter les configurations faisant apparaître un ressaut hydraulique dans l’ouvrage et impose
un régime d’écoulement fluvial à l’extrémité amont du déversoir. D’autre part, la relation
permettant le calcul du débit déversé est proche de la formule de Poléni qui n’est valable que
pour le déversoir frontal.
36
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
De plus, pour le déversoir à crête basse le calcul est basé sur l’hypothèse que la ligne
d’eau diminue sur le seuil pour arriver à une hauteur d’eau nulle à l’aval de la crête ce qui à
priori est faux.
Nous avons regroupé dans l’annexe 1 de ce travail quelques détails et commentaires
concernant le principe de calcul à la base de cet outil.
5.4 Raisonnement à énergie spécifique constante
Une approche plus physique initiée par Ackers en 1957 basée sur un raisonnement à
énergie constante a permis de progresser dans la connaissance du comportement hydraulique
du déversoir. En particulier, cette approche a permis de s’intéresser non seulement à
l’évaluation du débit déversé mais également à la forme de la ligne d’eau sur la crête du
déversoir. Malheureusement, comme le montre l’étude de El Khashab, cette méthode tombe en
défaut dans certains cas car les équations s’avèrent inadaptées. En particulier dans le cas de
l’apparition d’un ressaut hydraulique dans l’ouvrage.
Cette approche est plus théorique. Elle se base sur des équations phénoménologiques
avec une hypothèse importante qui est celle d’une énergie restant constante le long du seuil. En
fait, l’énergie varie très lentement car ses variations sont dues essentiellement aux pertes de
charge linéaires, sauf dans le cas d’un ressaut hydraulique où la dissipation d’énergie devient
importante.
De plus, il propose une formule pour le débit déversé qui est basée sur la formule de
Poléni mais corrigée, qui permet de prendre en compte l’effet de la vitesse latérale et celui de
l’entonnement dans le cas d’un seuil oblique.
Tous les calculs sont menés pour aboutir à une solution qui puisse être appliquée
facilement. On dispose de formules pour le dimensionnement et le diagnostic des déversoirs
latéraux à crête haute et à crête basse, pour des conduites rectangulaires et circulaires.
L’ensemble des résultats a été mis sous forme d’abaques faisant intervenir les variables
adimensionnelles suivantes :
kx h θ w
X= , y= , Θ= , W=
b Hs k Hs
37
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Ils ont été vérifiés expérimentalement et l’erreur commise reste acceptable. Elle est
d’environ 5 % si on excepte la formation de ressaut, les calculs sont en bonne concordance
avec les modèles réduits. Dans le cas de l’apparition, dans l’ouvrage, d’un ressaut hydraulique,
lieu de dissipation d’énergie, l’approche à énergie constante ne peut plus être appliquée.
El Khashab et Smith [El Khashab-1976] ont proposé un modèle basé sur l’équation de
la quantité de mouvement. On considère que l’écoulement est unidirectionnel selon l’axe
principal de l’écoulement (Ox) ; on suppose, de plus, que dans une section droite la répartition
des pressions est hydrostatique et que la vitesse est uniforme.
u
(Ox)
v
x x+dx
38
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Qdev vucosϕ − 2v 2 Q 2 ∂A
S0 − Sf + + .
dh Q g gA 3 ∂x
h ′(x) = = (2.60.)
dx Q² ∂A
1 − 3 .
gA ∂h
Cette expression permet de déterminer les lignes d’eau apparaissant dans le déversoir
dans le cas uniforme permanent et uniquement dans le cas où il n’y a pas de ressaut dans le
déversoir. Dans le cas ou un ressaut apparaît dans le déversoir, il faut rajouter les équations de
Rankine-Hugoniot au système d’équations précédent afin de connaître la ligne d’eau [Carleton,
1985].
La connaissance préalable des conditions aux limites est nécessaire pour obtenir la
convergence du calcul. Tous les cas possibles de ligne d’eau doivent être inventoriés. Le calcul
se fait en émettant plusieurs hypothèses de départ (régime d’écoulement : fluvial, torrentiel)
concernant le comportement hydraulique du déversoir. En effet, on choisit un cas parmi les 6
exposés à la figure 12. Après calcul, on vérifie la validité de ces hypothèses. La difficulté réside
dans le choix des hypothèses de départ. En effet, certaines d’entre elles ne sont pas acceptables
et peuvent faire diverger le calcul rapidement.
L’établissement du débit latéral est basé sur la constatation que les profils d’écoulement
sur un déversoir frontal et sur un déversoir latéral sont semblables. Souvent, on utilise la
relation de Poléni valable pour le déversoir frontal. Celle-ci est donnée par :
dQ
= C d 2g (h − w ) 2
3
(2.61.)
dx
Hager a travaillé sur les déversoirs latéraux et a établi une loi de déversement basée sur
la loi de Poléni [Hager-1987]. Celui-ci l’a corrigé pour tenir compte du déversement latéral.
Hager a adapté cette relation au cas du déversoir latéral en affectant à la loi de Poléni une série
de coefficients qui permettent de tenir compte des effets de :
39
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
u
u Φ
v
Crête Crête
v
u
h
u w
On utilise donc une formule de type Poléni et on ajoute des coefficients permettant de
prendre en compte l’effet de la vitesse latérale et de sa direction (ωv), ainsi que l’effet de
l’entonnement du déversoir (ωΦ).
dQ
Qdev = − = −0.6n *.c w gH 3 (y − W)3/ 2 × ωu × ωφ (2.62.)
dx
Dans le cas d’un seuil frontal, la vitesse principale u est nulle ; si on exprime la vitesse v
à partir de la définition de la charge spécifique H, on obtient, en supposant que la hauteur d’eau
sur le seuil est faible :
v = 2g(H − w) ; Or cette vitesse est faible, on peut donc considérer que H est proche de h. On
obtient :
v = 2g(h − w) (2.63.)
u²
Dans le cas d’un seuil latéral, la charge spécifique s’exprime par H = +h.
2g
u²
Si << h , on se rapproche du cas du déversoir frontal, si ce n’est pas le cas, cette vitesse a
2g
une influence sur le débit déversé. On a alors, en supposant toujours que la hauteur d’eau est
proche de celle du seuil, u = 2g(H − w) .
On définit alors ωu comme le rapport de la vitesse v du cas du seuil latéral et de la vitesse u du
cas du seuil frontal :
u H−w
ωu = = (2.64.)
v h−w
40
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
u
dx
Or :
u = vcosΦ
u2 u2
⇒ cos 2 Φ = ⇒ sinΦ = 1 − (2.67.)
v2 v2
1
2
P1 u1 ² P v ²
z1 + + = z2 + 2 + 2
ρg 2g ρg 2g
u1 ² v 2 ²
h+ = + w + ε(h − w)
2g 2g
41
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
2
ε=
3
u1 ² (h − w)(ε − 1) (y − W)(ε − 1)
−1 = = (2.68.)
v2 ² H − w − ε(h − w) 1 − ρy − W(1 − ε)
1/2
y−W
ω Φ = sinΦ = (2.69.)
3 − 2y − W
L’effet de l’entonnement doit être pris en compte quand la crête fait un angle θ (θ<0)
avec la direction principale du canal (Ox).
θ u
Φ
v
θ
Π
dQ = − vcos − θ − Φ × (h − w).dx (2.70.)
2
Or :
Π Π
cos − θ − Φ = cos − (θ + Φ) = sin(θ + Φ) (2.71.)
2 2
42
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
1/2
y−W
sin(θ + Φ) = = F = sinθinθ.c + sinΦinΦ.c (2.72.)
3 − 2y − W
On considère des valeurs de θ <<1, donc on peut remplacer les valeurs de sin θ et cos θ par
leurs équivalents : θ et (1- θ2/2). En développant chaque terme et en négligeant les termes en θ
d’ordre supérieur à 2, on trouve :
1
1/2
ωΦ = sinΦ = A. 1 + θ 2 − 1 (2.73.)
A
dQ 1− W y−W 1
1/2
= −0.6n*c w gH (y − W)
3 3/2
× × 1 − θ 2 − 1
dx y−W 3 − 2y − W A
où :
dQ 1− W 3(1 − y)
1/2
= −0.6n * c w gH (y − W) ×
3 3/2
× 1 − θ (2.74.)
dx 3 − 2y − W y − W
Par rapport à une formule classique, on trouve des valeurs du débit déversé plus fortes à
cause de la prise en compte des effets de la vitesse latérale et de l’entonnement. C’est cette
relation donnant le débit déversé qui a été intégrée à notre modèle.
5.7 Modélisation tridimensionnelle
43
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Crête déversante
Compte tenu des difficultés et limites associées aux diverses méthodes d’obtention de la
ligne d’eau et/ou du débit déversé décrites dans les paragraphes précédents (Méthodes
empiriques, énergie spécifique constante,…), nous avons choisi d’ajouter au système
d’équation de Barré de Saint-Venant sous forme conservative un terme de déversement pour
tenir compte du débit déversé par l’intermédiaire du déversoir d’orage. Ainsi, la forme des
équations telles que nous les avons utilisées est la suivante :
∂A ∂Q
∂t + ∂x = Q dev
(2.75.)
∂Q + ∂ Q + gI1 = gI 2 + gA(S 0 − S f ) + QQ dev
2
∂t ∂x A
A
Par commodité d’écriture, ce système peut être mis sous la forme matricielle suivante :
U t + F(U )x = G (U ) ⇔ U t + J (U )U x = G (U )
avec :
44
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
A Q Q dev
U = ,
F(U ) = Q 2 ,
G (U ) = QQ dev
Q + gI1 gI 2 + gA(S 0 − S f ) +
A A
où :
dQ
Qdev = − : Terme d’apport qui pour nous correspond au débit déversé par
dx
l’intermédiaire du déversoir d’orage. Ce terme de déversement est explicité par la relation de
Hager (cf. Partie 2 § 5.6.4).
Q2
Sf = 4
(2.76.)
2 2
K A Rh
S
3
Les forces de pression sont calculées en tenant compte des propriétés géométriques du canal. I1
correspond au terme de pression hydrostatique et I2 représente la force de pression due à la
variation de la largeur du canal.
Ces forces s’expriment de la manière suivante :
h
I1 = ∫ (h − η)Bdη (2.77.)
0
h
∂B
I 2 = ∫ (h − η) dη (2.78.)
0
∂x
où h=hauteur d’eau ; η=variable d’intégration indiquant la distance depuis le fond du canal ;
B(x,η) largeur du canal à la distance η du lit du canal qui s’exprime :
∂A(x, η)
B(x, η) = (2.79.)
∂η
45
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
6 CHANGEMENT DE SECTION
6.1 Elargissement
∆H12
ξe = (2.80.)
V12 / 2g
2
∆H A
ξ e 90° = 2 12 = 1 − 1 (2.81.)
V1 / 2g A 2
Par analogie avec ce cas particulier, l’influence de l’angle d’élargissement est pris en
compte par la relation :
δ
φ e (δ ) = + sin (2δ ) 0 ≤ δ ≤ 30
90
(2.83.)
5 δ
φ e (δ ) = − 30 ≤ δ ≤ 90
4 360
6.2 Contraction
46
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
L’écoulement contracté est essentiellement influencé par l’angle de contraction et le rapport des
A
sections mouillées φ = 2 < 1 .
A1
Peu de résultats expérimentaux sont disponibles concernant cet ouvrage. On notera
l’article de Gardel [Gardel-1962]. La relation obtenue pour le coefficient de perte de charge
est :
1.83(1− φ )0.4
∆H 1 δ
ξ ν = 2 12 = (1 − φ) (2.84.)
V2 / 2g 2 90
Il est à noter que le coefficient de perte de charge est toujours inférieur à celui de
l’élargissement correspondant et ξ ν augmente de manière conséquente avec l’angle de
contraction
Des investigations expérimentales ont été réalisées pour obtenir les hauteurs d’eau maximales
atteintes [Hager-1999]. Mais pour obtenir les hauteurs d’eau dans toute la structure il faut
passer par des logiciels de modélisation. C’est ce qu’ont réalisé par exemple Rahman et
Chaudhry en 1997 [Rahmann-1997] en résolvant le système d’équations de Barré de Saint-
Venant 2D par une méthode numérique à capture de choc de type Mac-Cormack pour les
élargissements.
47
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
Dans les modèles numériques unidimensionnels que nous avons réalisés, la présence du
changement de section est prise en charge par un terme de pression supplémentaire qui
correspond à la pression exercée par l’eau sur les parois de l’ouvrage. C’est en fait un terme de
pression latérale [Delis-2000]. Celui-ci est exprimé par le terme I2 que l’on retrouve au second
membre du système d’équations de Barré de Saint-Venant. On rappelle son expression ci-
dessous :
h
∂B
I 2 = ∫ (h − η) dη (2.85.)
0
∂x
∂A(x, η)
B(x, η) = (2.86.)
∂η
48
Partie 2 : Ecoulements dans les collecteurs et déversoirs du réseau d’assainissement
49
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
PARTIE 3
METHODES
MISES EN ŒUVRE
DANS LES
OUTILS NUMERIQUES
REALISES
50
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
1 INTRODUCTION
Dans cette partie de la thèse, on expose dans un premier temps la méthode des
caractéristiques associée à la théorie des ondes qui a mené par la suite au développement des
méthodes à capture de choc. Les limites de cette méthode sont illustrées grâce à quelques
exemples représentatifs des situations qui peuvent être rencontrées dans le cadre de l’étude du
fonctionnement hydraulique d’un système. Ensuite, on présente les méthodes aux différences
finies classiques les plus couramment employées dans le cadre de l’hydraulique urbaine et plus
particulièrement dans le domaine de l’assainissement. Pour finir, on expose les méthodes
numériques de type Godunov à capture de choc qui servent de base aux logiciels de simulation
qui ont été construits. Les méthodes numériques effectivement utilisées dans nos logiciels sont
présentées en fin de partie.
• Le second est celui du traitement des crues (onde de crue diffusante ou encore
onde cinématique ou de continuité pour un modèle encore plus simplifié).
Dans le cas des crues, on cherche dans le cadre de la protection des populations contre
les inondations, à connaître à partir de mesures réalisées en amont l’évolution de la crue. Dans
ce cas, on met en œuvre des modèles d’onde de crue diffusante ou d’onde cinématique.
Pour l’onde cinématique, Seddon étudia les crues du Mississipi et établi en 1899 que les
ondes se propagent à une vitesse supérieure à celle de l’écoulement. La nature de ces ondes est
différente de celles produites dans le cas à inertie prépondérante. Ceci permit de donner un
ordre de grandeur plus réaliste de la vitesse de propagation des crues.
Pour l’onde de crue diffusante, la méthode consiste à déterminer l’hydrogramme de
sortie connaissant l’hydrogramme amont en tenant compte de l’atténuation de l’onde au cours
du temps [Viollet-1998]. Le décalage temporel dû au temps mis par le signal pour traverser le
collecteur ou le canal est directement lié à la vitesse de l’eau qui ne doit pas être confondue
avec la vitesse de propagation de l’onde [Mottiee–1996].
L’approche la plus précise, à utiliser si l’on veut par exemple dimensionner un ouvrage
et évaluer son impact en terme d’inondation, passe par la modélisation numérique des équations
de BSV. Dans un premier temps, on remplace le système d’équations aux dérivées partielles
par un système d’équations algébrique ; c’est la phase de discrétisation des équations. Puis, on
applique à ce système une méthode numérique de résolution qui permet de calculer les valeurs
de hauteurs d’eau et de vitesse.(Différences finies, volumes finis, éléments finis)[Viollet-1998].
Actuellement, pour l’étude des crues, on utilise de plus en plus souvent des modèles
bidimensionnel.
51
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Le cas de la présence d’un ressaut hydraulique ne peut pas être résolu par la méthode
des caractéristiques. En effet, la transition d’un régime d’écoulement torrentiel à un régime
fluvial se traduit par le fait que l’une des courbes caractéristique devient localement verticale.
Ceci signifie qu’une partie de l’information ne se propage plus ni vers l’amont, ni vers l’aval. Il
y a donc une discontinuité physique entre deux points très proches qui interdit la résolution
locale du système de Barré de Saint-Venant. Une solution passe alors par l’intervention des
bilans de masse et de quantité de mouvement à travers la discontinuité afin d’obtenir la vitesse
de propagation du ressaut ainsi que les hauteurs conjuguées. Ceci rend beaucoup plus lourde
une résolution par cette méthode. L’application d’une méthode de résolution du type de celle
des caractéristiques nécessite le découplage des équations du ressaut et des équations continues.
De manière plus efficace, les méthodes de résolutions numériques actuelles permettent
d’assurer le couplage des équations du ressaut avec les équations de BSV par l’intermédiaire de
conditions de saut dites de Rankine-Hugoniot. Ainsi, le ressaut est en fait contenu dans la
formulation des équations de BSV et une méthode numérique capable de prendre en compte les
discontinuités permet de déterminer une solution sans plus avoir à traiter les équations du
ressaut à part.
Dans ce cas, on peut obtenir des solutions simples dites de Ritter [Viollet-1998] et une
interprétation en terme de caractéristiques. Mais la simulation correcte des ondes de
submersion et de crues associées à une rupture de barrage passe par l’utilisation d’un schéma
numérique capable de prendre en compte de fortes discontinuités et en particulier par la
résolution d’un problème de Riemann.
∂h ∂u ∂h
∂t + h ∂x + u ∂x = 0
1 ∂u u ∂u ∂h (3.87.)
+ + = S0 − Sf
g ∂t g ∂x ∂x
Avec :
52
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
∂h ∂h
dh = dx + dt
∂x ∂t
(3.88.)
du = ∂u dx + ∂u dt
∂x ∂t
Où les deux équations du système (3.88) sont les différentielles totales des hauteur d’eau, dh, et
de la vitesse moyenne, du, respectivement.
Le système d’équation constitué des équations (3.87) et (3.88) donne un jeu de quatre équations
∂h ∂h ∂u ∂u
aux dérivées partielles avec quatre inconnue : , , , .
∂x ∂t ∂x ∂t
Si l’on exprime la profondeur d’eau, h, par la célérité de l’onde, c, on écrira dans le cas
rectangulaire :
c 2 = gh (3.89.)
d ( c 2 ) = 2cdc = d ( gh )
∂c ∂u ∂c
2
∂t + c + 2u =0
∂x ∂x
(3.90.)
∂u + u ∂u + 2c ∂c = g ( S − S )
∂t ∂x ∂x
0 f
∂ ∂
∂t + ( u + c ) ∂x ( u + 2c ) = g ( S0 − Sf )
(3.91.)
∂ ∂
∂t + ( u − c ) ∂x ( u − 2c ) = g ( S0 − Sf )
Les arguments de gauche peuvent représenter des dérivées totale ; on peut donc écrire :
d dx
dt ( u + 2c ) = g ( S0 − Sf ) avec ( u + c) =
dt
(3.92.)
d ( u − 2c ) = g ( S − S ) avec ( u − c) =
dx
dt 0 f
dt
dx
= cw = ( u ± c ) (3.93.)
dt
53
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
qui représente la vitesse par rapport au sol. Il y a deux ondes dont l’une se propage
obligatoirement vers l’aval et l’autre peut se propager soit vers l’amont, soit vers l’aval.
Les quatre équations du système (3.92) forme le système d’équations aux dérivées
totales des caractéristiques ; elles remplacent le système d’équations aux dérivées partielles de
Barré de Saint-Venant (2.13).
dx
A condition qu’un observateur suive les courbes définies par = ( u ± c ) , en admettant
dt
que le canal soit sans frottement et horizontal, alors ( S0 − Sf ) = 0 , prennent les formes simples
dx
U + 2c = Cte sur les caractéristiques positives C+ définies par = ( u + c)
dt
dx
U − 2c = Cte sur les caractéristiques négatives C- définies par = ( u − c)
dt
Les courbes caractéristiques, C+ et C-, peuvent être tracées dans un plan, x et t et dans ce
cas précis sont des droites puisque le canal est supposé sans frottement et horizontal. Ceci
signifie qu’une perturbation de la surface libre est simplement translatée au cours du temps sans
aucune atténuation. Un écoulement tel qu’une famille de caractéristiques est rectiligne est
appelé une onde simple.
dx
On a ainsi défini les courbes caractéristiques = ( u ± c ) dans le cas du système de
dt
Barré de Saint-Venant ainsi que les invariants de Riemann du système u ± 2c = Cte .
On recherche dans le plan (x,t) des courbes caractéristiques le long desquelles l’équation
aux dérivées partielles devient une différentielle totale. Par exemple, lorsque la section du
collecteur est rectangulaire les courbes caractéristiques sont définies dans le plan (x,t) par les
deux équations différentielles [Graff-1996]:
dx
dt = u + c courbe notée C+
(3.94.)
dx = u − c courbe notée C
dt -
On note c la célérité de l’onde dynamique définie par c 2 = gh . La hauteur d’eau est exprimée
par la célérité.
Ces courbes définissent les célérités de propagation des ondes dites dynamiques.
Lorsque le nombre de Froude Fr<1 la vitesse absolue de l’onde (courbe C- ) est dirigée vers
l’amont et la courbe C+ descend le courant. On est dans le cas de l’écoulement fluvial. Lorsque
le nombre de Froude Fr>1 les courbes C- et C+ descendent le courant. On est dans le cas de
l’écoulement torrentiel.
54
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
g g
du + dh + g(Sf − S0 )dt = 0 du − dh + g(Sf − S0 )dt = 0
h h (3.95.)
dx = (u + gh )dt dx = (u − gh )dt
Une solution des équations précédentes peut être établie de façon numérique ; elle
donne la description d’un mouvement non permanent d’une position à l’autre dans le canal.
t
C-
C+
P(xP,tP)
L R 104
101
11 12 13
14
1 2 3 4 5 x
Figure 23 : Surface autour du point P(xP,tP) projetée sur un plan (x,t) (grille irrégulière)
De plus, pour compléter la méthode, il est nécessaire de connaître des conditions aux
limites. En effet, pour xP proche de l’une des deux extrémités de la canalisation, l’une au moins
des caractéristiques ne coupe pas la droite t=0. Les conditions aux limites s’expriment par une
relation du type f (u , h)=0 qui peut être, par exemple, la donnée de u + 2 gh sur la droite x=0.
Si les caractéristiques ont des pentes de signe différent (ce qui est le cas lorsque l’écoulement
est fluvial), il est nécessaire de connaître une condition à la limite amont et une condition à la
limite aval. Si elles ont des pentes de même signe (ce qui est le cas lorsque l’écoulement est
55
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
torrentiel), il faut connaître deux relations entre u et h à l’amont, aucune n’étant indispensable à
l’aval [Kovacs-1988].
Les schémas résolvant la grille des caractéristiques demandent une grande quantité
d’interpolations et exigent une puissance de calcul relativement importante. De plus, ils sont en
général non conservatifs puisque les caractéristiques sont assimilées à des segments de droite
entre les points de discrétisation. Ils sont adaptés essentiellement aux canaux de forme simple
et peu variable, ce qui signifie qu’ils ne sont pas adaptés à l’assainissement. Les réseaux
contiennent en général de nombreuses singularités et sont géométriquement complexes.
La méthode des caractéristiques n’est pas valable pour les problèmes où apparaît une
discontinuité. Ainsi, les problèmes du ressaut hydraulique et de la rupture de barrage
(discontinuité dans les conditions initiales) ne peuvent pas être résolus par cette méthode.
Néanmoins, comme nous allons le voir, on peut en avoir une interprétation en terme de
caractéristiques.
Une équation hyperbolique non linéaire peut admettre des chocs ne figurant pas dans les
données par focalisation des caractéristiques [Euvrard-1990]. C’est le cas du ressaut
hydraulique. Dans les données initiales du problème la discontinuité sur la hauteur d’eau n’est
pas présente. Au cours du temps un front d’onde se forme et la courbure des lignes
d’écoulement devient très importante. Ainsi, les courbes caractéristiques de même nature se
resserrent pour finalement se croiser. Ceci n’est pas acceptable car, sur chaque caractéristique,
la solution est constante et égale à la pente de la droite. En conséquence, au point
d’intersection, le problème aurait une infinité de solutions. Cette situation est représentée sur la
figure suivante :
56
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
h t
?
Ainsi, il faut admettre que les solutions peuvent être de nature discontinues et la
méthode des caractéristiques tombe en défaut. Le ressaut hydraulique se propage à une certaine
vitesse différente de celle de l’écoulement et de celle des ondes de surface. C’est pourquoi, la
détermination de la vitesse de propagation du ressaut ainsi que de sa position nécessite
l’utilisation d’une technique de calcul plus complexe : mise en œuvre des équations de Barré de
Saint-Venant sous forme conservative et utilisation d’une méthode numérique supportant les
discontinuités.
Dans le cas de la rupture de barrage, le problème est différent car la discontinuité est
présente au temps initial. C’est cette discontinuité initiale qui ne peut pas être gérée par la
méthode des caractéristiques car on ne peut pas définir les caractéristiques au point de
discontinuité.
Pour contourner le problème on peut initialiser la solution à t=0+ en utilisant une
solution analytique du type solution de Ritter [Viollet-1998].
57
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
x
Vitesse à t+∆t
x
t
t*
Raréfaction Choc
x
• Apparition d’une onde choc dans la zone de faible tirant d’eau qui fait remonter la ligne
d’eau de manière abrupte (Ressaut).
• Apparition d’une onde de détente dans la zone à fort tirant d’eau qui fait baisser la hauteur
de la surface libre.
Dans le travail qu’il a mené, Benayada a construit une solution analytique du système
de Barré de Saint-Venant complet par le calcul opérationnel à savoir en utilisant la transformée
de Laplace dans le but d’étudier les écoulements dans les rivières et les canaux avec fond
mobile [Benayada–1994]. Dans un premier temps, celui-ci réarrange le système de Barré de
Saint-Venant en le linéarisant autour du régime permanent par introduction de nouvelles
58
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
variables qui représentent les écarts aux valeurs moyennes de la hauteur et de la vitesse et en
utilisant des variables sans dimensions permettant de généraliser les résultats.
Puis, il introduit dans ce nouveau système les transformées de Laplace l [ u ] et l [ y ] des
variables sans dimensions de hauteur y et de vitesse u dont les expressions sont ensuite
recherchées sous la forme l [ u ] = Ae rX et l [ y ] = Be rX où X correspond à l’abscisse
adimensionnelle le long de l’écoulement. Après avoir déterminé les expressions de U et Y,
celui-ci retourne à l’original en utilisant le théorème d’inversion de Mellin-Fourier :
1
f (x, t) = ∫
2πi C
e pt F(X, p)dp (3.96.)
Les expressions précédentes ont été établies pour une excitation échelon de la vitesse
mais Benayada les a généralisées à une excitation quelconque afin de pouvoir simuler des
hydrogrammes de crue.
Les expressions de u et de y comportent deux parties. La première a une forme hyperbolique et
la seconde est de forme trigonométrique. Ceci signifie que le mouvement peut être périodique.
59
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
On commence par quadriller le plan (x, t) afin d’obtenir des mailles de taille (∆x, ∆t) où
∆x est le pas d’espace et ∆t le pas de temps. Le but est de calculer les valeurs de débit et
hauteur d’eau en chaque point du maillage. Les méthodes aux différences finies sont toutes
basées sur les développements de Taylor des fonctions continues et dérivables. Plus les pas de
temps et d’espace sont petits plus les développements limités sont proches des valeurs exactes.
A un instant donné, on peut écrire :
∂f ∂ 2f ∆x 2 ∂ 3f ∆x 3
f i +1 = f i + ∆x + 2 + + ε(∆x 4 )
∂x ∂x 2! ∂x 3 3!
(3.98.)
∂f ∂ 2f ∆x 2 ∂ 3f ∆x 3
f i −1 = f i − ∆x + 2 − 3 + ε(∆x 4 )
∂x ∂x 2! ∂x 3!
Les dérivées partielles par rapport à x sont comptées à l’abscisse i∆x. Les équations sont
valables quel que soit le temps t. On exprime alors à l’aide des développements limités les
∂f ∂ 2 f
termes , ,….. en fonction des valeurs de la fonction f aux nœuds du maillage avec les
∂x ∂x 2
notations suivantes :
f in = f ( i∆x, n∆t ) 0 ≤ i ≤ If ; 0 ≤ n ≤ Nf
(n+1)∆t
n∆t
(n-1)∆t
60
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Pondération Schéma
α3=0 explicite
α3=0.5 implicite centré dans le temps
α3=1 totalement implicite
α1=1 α2=1 Progressif
α1=0.5 α2=0.5 Centré
α1=0 α2=0 Régressif
α1=1 α2=0 Mixte décentré
Tableau 3 : Divers types de schémas numériques
Si l’équation différentielle à résoudre contient d’autre variables que f nous pouvons les
exprimer également sous forme d’un développement de Taylor. Une fois que toutes les dérivées
partielles sont approchées par les valeurs de f aux différents nœuds du maillage, on substitue
leurs expressions dans l’équation différentielle à résoudre.
- Schéma explicite :
Seule la dérivée de f par rapport au temps s’exprime en fonction des valeurs de f au pas
de temps n+1.
La différentielle par rapport à x s’exprime en fonction des valeurs de f au pas de temps
précédent de calcul (pas n). On calcule f in +1 connaissant les valeurs de f calculées au temps n∆t
pour, en général, les pas d’espaces i-1, i et i+1 (schéma à trois points). Un schéma explicite est
caractérisé par le fait que l’on puisse exprimer explicitement une valeur inconnue en fonction
de valeurs connues.
t
(n+1)∆t
n∆t
(n-1)∆t
Les schémas explicites sont les plus simples. Cependant, il est nécessaire de prendre un
pas de temps très court, du fait de la condition de stabilité numérique qui impose :
61
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
∆x
∆t ≤ (3.101.)
u +c
- Schéma implicite :
(n+1)∆t
n∆t
(n-1)∆t
Pour les schémas explicites, on peut calculer les valeurs de f de proche en proche en
progressant par balayage sur tous les pas d’espace i de pas de temps en pas de temps
(Progression horizontale).
Pour les schémas implicites, les valeurs de f ne peuvent être calculées au pas de temps n+1
qu’en résolvant le système formé des If ou If–1 équations (selon la définition des conditions aux
limites) linéaires algébriques du schéma.
(PM).(f)n+1 = (SM)
où (PM) est une matrice bande de taille (If –1, If –1) qui contient autant de bandes qu’il y a de
valeurs inconnues de f (valeurs de f au pas de temps n+1) par équation.
Le système est représenté par une matrice bande uniquement dans le cas d’un réseau linéaire ou
arborescent (si les équations sont correctement ordonnées). Si le réseau est maillé les matrices
sont beaucoup moins bien conditionnées.
(SM) est la matrice colonne second menbre qui contient tous les termes connus.
En général, le système formé des If ou If –1 équations est non linéaire puisque les coefficients
des dérivées partielles de f dans l’équation à résoudre sont souvent fonction des valeurs de f.
62
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Pour être résolu, un tel système matriciel doit être linéarisé (les coefficients des matrices (PM)
et (SM) doivent être choisis constants), ce qui suppose des itérations successives.
On retrouve la définition du schéma explicite lorsque la matrice (PM) est diagonale. Lorsque la
matrice (PM) est triangulaire inférieure, ce qui est le cas lorsque l’on a un schéma explicite
centré dans le temps et l’espace, on peut à l’aide de la condition aux limites amont f 0n +1 ,
calculer les valeurs de f aux différents pas d’espace i (au pas de temps n+1) de proche en
proche sans être obligé de résoudre le système en inversant la matrice (PM). On parle alors de
schéma semi-explicite.
t ν∆x (1−ν)∆x
(n+1)∆t
(1−α3)∆t
∆t
α3∆t
n∆t
∆x
i∆x (i+1)∆x x
Figure 29 : Schéma de discrétisation dans le plan (x,t)
représentent, quand ils sont présents, la partie principale de l’erreur commise par le schéma
puisque ce sont les termes d’ordre le plus bas. Cette partie principale de l’erreur est connue
sous le nom de diffusion numérique des schémas aux différences finies [Kovacs–1988].
Le schéma est consistant si lorsque l’on fait tendre ∆x et ∆t vers zéro alors, on retrouve
les équations aux dérivées partielles de départ.
Le schéma de Preismann est du premier ordre en x (c’est à dire consistant en 0(∆x)) et du
second ordre en t (c’est à dire consistant en 0(∆t2)) pour α3 ≠ 0,5. Si α3 = 0,5 le schéma est du
second ordre en x et en t. [Malaterre–1994].
L’étude théorique de la stabilité d’un schéma numérique peut être menée en utilisant la
technique de l’analyse de Fourrier. Cette technique n’est valable que pour les systèmes
linéaires, il faut donc linéariser les équation discrétisées de BSV..
On remplace dans les équations discrétisées la vitesse ou le débit et la hauteur d’eau par
des séries de Fourrier puis on met le système sous forme matricielle. On détermine ensuite les
valeurs propres de la matrice d’amplification.
Le schéma est stable si ces valeurs propres sont inférieures à 1 [Szymkiewicz–1996].
Pour que le schéma soit stable, il faut que l’erreur soit bornée si le temps tend vers l’infini.
Nous venons de voir que la méthode des caractéristiques n’est pas apte à résoudre les
problèmes faisant apparaître des solutions discontinues ou pour lesquels des conditions initiales
discontinues doivent être considérées.
Pour la méthode des différences finies, ses inconvénients sont, d’une part le présupposé
de la dérivabilité de la solution et d’autre part la difficulté de traiter des géométries réalistes
donc complexes [Monthe-1992]. Cependant, cette méthode est très performante dans le cas
d’écoulements lents. En particulier, le nombre de Froude doit rester inférieur à 1 voire inférieur
à 0.5 [Paquier-1995] et peut également résoudre les situations dans lesquelles le régime
d’écoulement est purement torrentiel.
64
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Dans le cadre de l’assainissement, la méthode aux différences finies la plus utilisée reste
celle de Preismann dans laquelle l’algorithme du double balayage (les débits sont d’abord
calculés de l’amont vers l’aval puis les hauteurs d’eau de l’aval vers l’amont) permet
d’effectuer l’inversion de matrices de manière simple. Une étude récente [Meselhe et al.–1997]
montre que le schéma de Preismann devient instable lorsque le régime d’écoulement est proche
du régime critique ou lorsque le régime devient « transcritique » c’est à dire que dans le même
bief ou tronçon coexistent régimes fluviaux et torrentiels. Ainsi, cette méthode n’est pas
utilisable lorsqu’un ressaut hydraulique est présent ou lorsqu’il y a assèchement ou remplissage
d’un bief.
Les logiciels les plus courant utilisent alors des approximations plus ou moins adaptées
pour résoudre ces types de problèmes [Meselhe et al.-1997].
Une autre méthode numérique performante est celle des éléments finis. Dans le cas de
son application à la résolution des équations de Barré de Saint-Venant, la démarche consiste à
utiliser un schéma aux différences finies sur le temps, puis à linéariser et enfin à résoudre à
chaque pas de temps le problème elliptique ainsi obtenu par la méthode des éléments finis.
Nous n’avons pas considéré celle-ci dans cette étude car nous n’avons pas estimé nécessaire de
s’y intéresser dans le cadre d’une résolution d’équations 1D au service de l’hydraulique urbaine
et de l’assainissement. Elle a déjà été testée pour la résolution des équations de Barré de Saint-
Venant par l’utilisation d’éléments de type Galerkine [Szymkiewicz-1991]. Sa mise en œuvre
dans le cadre de la résolution des équations de Barré de Saint-Venant par l’intermédiaire
d’éléments finis discontinus vient d’être lancée au laboratoire S.H.U. dans le cadre d’une thèse
de doctorat réalisée par Maher Abdallah. Elle a l’avantage de permettre une prise en compte
très fine de la géométrie du support et pourrait être bien adaptée à des problèmes du type
jonction de collecteurs.
65
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
66
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
En fait, c’est Lax et Wendroff en 1960 [Euleterio-1997] qui ont montré que ce type de
méthodes si elles sont convergentes mènent à des solutions correctes des équations dites
solutions faibles. Comme les chocs se propagent à une vitesse différente de celle de la matière,
ils induisent une augmentation d’entropie [Hanich-1996]. Comme la solution n’est pas unique,
il a fallu construire un critère d’entropie dans le schéma numérique qui permet d’assurer la
convergence automatique vers la solution physiquement admissible du problème.
Dans sa première méthode numérique Godunov utilisa une solution exacte du problème
de Riemann. Ensuite, sa méthode a été rendue plus attractive par une extension à un second
ordre de précision par Van Leer [Van Leer-1979] puis par le développement de solveurs de
Riemann approchés comme celui de Roe ou celui de Osher et Salomon.[Delis et Skeels–2000].
L’extension au second ordre de la méthode rappela avec pertinence un théorème que
Godunov avait établi :
Tout schéma numérique, d’un ordre de précision supérieur à 1, produira des oscillations
au niveau des discontinuités.
La difficulté associée au fait que les méthodes du premier ordre sont trop imprécises et
que les méthodes d’ordre plus élevé produisent des oscillations resta un obstacle au
développement de « bonnes » méthodes numériques pour les équations de conservations
hyperbolique.
Une voie possible pour rendre le schéma non oscillatoire serais d’utiliser des méthodes
monotones. Mais celles-ci sont au plus précises au premier ordre et sont donc d’utilisation
limitée. Une meilleure voie pour résoudre cette contradiction est celle de la construction de
méthodes numériques TVD (Total Variation Diminishing). Elles sont apparues dans les années
80. Les travaux fondateurs sont dus à Harten [Harten-1983] [Harten-1984], Sweby [Sweby-
1984] et de nombreux autres. Ces travaux ont mené à la mise au point d’une classe de méthodes
numériques qui produit des résultats non oscillatoires au voisinage des discontinuités et précise
au moins au second ordre. Il existe deux méthodes de construction de méthodes numériques
TVD :
C’est ce type de méthodes numériques que nous avons utilisées pour développer notre
modèle dans l’optique de l’appliquer au domaine de l’hydraulique urbaine et plus
particulièrement à l’assainissement.
3.1.1 Définition
u t + au x = 0
uL si x < 0 (3.102.)
u ( x, 0 ) = u 0 ( x ) =
u R si x > 0
u0(x)
uL
uR
x=0 x
On considère dans ce qui suit le cas de l’une des équations de type hyperbolique les plus
simple pouvant mener à des solutions discontinues. C’est l’équation de Bürger associé à des
valeurs initiales discontinues. On donne une solution de cette équation dans deux cas. Le
premier donne naissance à une onde de compression et le second à une onde de dépression.
On considère donc le problème aux valeurs initiales pour les équations de Bürger :
1 2
u t + f ( u )x = 0 , f ( u ) = 2 u
(3.103.)
u si x < 0
u ( x, 0 ) = u 0 (x) = L
u R si x > 0
68
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Comme le flux est convexe c’est à dire que λ′ ( u ) = f ′′ ( u ) > 0 les vitesses des
caractéristiques à gauche sont plus grandes que celles de droite λ L > λ R . Dans ce cas, les
données sont compressives. Les caractéristiques vont se croiser immédiatement et on observe la
formation d’un choc qui se déplacera à une vitesse S.
t
uL
uR
Onde de choc de
vitesse S
u si x-St < 0
u ( x, t ) = L (3.104.)
u R si x-St > 0
∆f 1
S= = (uL + uR ) (3.105.)
∆u 2
Cette solution discontinue est une onde de choc et est de nature compressive. Elle
satisfait la condition :
69
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
t
uR
uL
Onde de détente
uL si x ≤ x L
u − uL
u 0 ( x ) = u L + R ( x − xL ) si x L < x < x R (3.107.)
xR − xL
uR si x ≥ x R
t
uR
uL
Queue de l’onde
Tête de l’onde
Le front droit ou tête de l’onde est donné par la caractéristique qui émane de xR :
x = xR + λ (uR ) t (3.108.)
Le front gauche ou queue de l’onde est donné par la caractéristique qui émane de xL :
x = xL + λ (uL ) t (3.109.)
Ce phénomène de dispersion des ondes n’est pas rencontré dans l’étude des systèmes
d’équations hyperboliques linéaires à coefficients constants. La solution entière est :
x-x L
u ( x,t ) = u L si
t
≤ λL
x − xL x-x L
u0 ( x ) = λ ( u ) = si λ L < < λR (3.110.)
t t
x-x R
u ( x,t ) = u R si
t
≥ λR
71
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
x
u ( x,t ) = u L si ≤ λL
t
x x
λ (u) = si λ L < < λ R (3.111.)
t t
x
u ( x,t ) = u R si ≥ λ R
t
u si x-St < 0
u ( x, t ) = L
u R si x-St > 0 si u > u
L R
1
S = (uL + uR )
2
x (3.112.)
uL si ≤ uL
t
x x
u ( x, t ) = si u L < < u R si u L ≤ u R
t t
x
uR si ≥ u R
t
La solution exacte est une onde simple qui émane de l’origine cette onde est également
une onde de choc lorsque uL>uR ou une onde de raréfaction si uL≤uR.
On vient donc d’étendre le concept de solution pour pouvoir tenir compte des
discontinuités. Le problème est de distinguer une solution physiquement admissible d’une
solution non admissible. Comme on l’a vu dans le développement précédent qu’une
discontinuité physiquement admissible respectant la relation de Rankine-Hugoniot doit aussi
obligatoirement satisfaire la condition d’entropie.
Les notions que l’on vient d’étudier peuvent être directement étendues au cas du
système hyperbolique de lois de conservation.
t
ai
a2 am-1
am
a1
m m
U L = ∑ αi K ( ) , U R = ∑ βi K ( )
i i
(3.115.)
i =1 i =1
I m
U ( x, t ) = ∑ αi K ( ) + ∑ β K( )
i i
i (3.116.)
i =1 i = I +1
D’après l’équation (3.115) il est facile de voir que le saut de U à travers la globalité de
la structure d’onde dans la résolution du problème de Riemann est :
∆U = U R − U L = ( β1 − α1 ) K ( ) + K + ( β m − α m ) K (
1 m)
(3.117.)
Ceci est une combinaison linéaire des vecteurs propres avec des coefficients qui
pondèrent les ondes présentes dans le problème de Riemann. Le coefficient de pondération
pour l’onde i est βi − αi et le saut en U au passage de l’onde i noté ( ∆U )i est :
( ∆U )i = ( βi − αi ) K (i) (3.118.)
73
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
U t + F ( U )x = 0
U L si x < 0 (3.119.)
U ( x, 0 ) = U 0 ( x ) = U si x > 0
R
Ce système peut par exemple représenter le système d’équations de Barré de Saint-Venant sans
second membre.
(
U in +1 = U in − λ Fin+1 2 − Fin−1 2 ) (3.120.)
∆t
où λ = et Fi +1 2 = F ( U nj−l1 ,K , U nj+ l2 ) où l1 et l2 sont deux entiers non négatifs et cette
∆x
fonction est appelée flux numérique qui est une approximation du flux physique de (3.119).
Définition 9 : Consistance
Fi +1 2 ( v,K , v ) = F ( v ) (3.121.)
Ceci signifie que si tous les arguments de Fi +1 2 sont égaux à v alors le flux numérique est le
même que le flux physique en U=v.
74
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
U t + F(U) x = 0 (3.122.)
∫
volume
(Udx − F(U)dt) = 0 (3.123.)
t2
Volume de contrôle
t1
xL xR
D’où :
xR xR t2 t2
∫ U(x,t
xL
2 )dx = ∫ U(x,t1 )dx + ∫ F(U(x L ,t))dt − ∫ F(U(x R ,t))dt
xL t1 t1
(3.124.)
En posant :
x
1 i+1/2
∆x xi∫−1/2
U in = U(x,0)dx
x
(3.125.)
1 i+1/2
∆x xi∫−1/2
U in +1 = U(x,∆t)dx
75
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
∆t
1
∆t ∫0
Fi −1/2 = F(U(x L ,t))dt
∆t
(3.126.)
1
∆t ∫0
Fi +1/2 = F(U(x R ,t))dt
x t
Avec ∆x le pas d’espace, ∆t le pas de temps, i = et n = entiers.
∆x ∆t
i représente le nombre de cellules de calcul.
On nomme i + 1/2 ( respectivement i − 1/2 ) l’interface des cellules i et i+1 (respectivement i-1
et i). Fi +1/2 représente donc le flux entre les cellules i et i+1.
Cellules de calcul
Fi −1/2 Fi +1/2
i-1 i i+1
i-1/2 i+1/2 x
Figure 36 : Décomposition du domaine d’intérêt en cellules de calculs et flux aux
interfaces.
La figure ci-dessus représente la manière dont est décomposé le domaine de calcul. Les cellules
de calcul notées i, dans lesquelles les valeurs des variables recherchées sont constantes par
morceau, ont pour limites les interfaces i-1/2 et i+1/2. On calcule les flux numériques traversant
ces interfaces et permettant d’évaluer l’influence des valeurs contenues dans i sur ses voisines.
De nombreux schémas numériques conventionnels peuvent être mis sous une forme
conservative. Le problème est, pour chacun d’entre eux, de déterminer la fonction flux
adéquate. Pour tous ces schémas numériques de nature explicite se pose le problème du choix
du pas de temps à utiliser pour la résolution. De ce choix dépend la stabilité du schéma
numérique considéré.
76
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
n
En appelant Vmax la vitesse de l’onde la plus rapide dans le domaine de calcul au pas de temps
n, on définit le nombre de Courant maximum :
∆tVmax
n
Ccfl = (3.128.)
∆x
Ccfl ∆x
∆t = n
(3.129.)
Vmax
Dans le cas de l’étude du système de Barré de Saint-Venant, le pas de temps sera donné par :
Ccfl ∆x
∆t = (3.130.)
u +c
Le plus utilisé des solveurs de Riemann approché est sans doute celui introduit par Roe
en 1981 [Roe-1981]. Depuis, cette méthode a été raffinée et appliquée à un large panel de
problèmes physiques. Les raffinements de la méthode de Roe ont été introduit par Roe et Pike
en 1984 [Roe and Pike-1984]. Grâce à ces raffinements, il n’était plus nécessaire d’exprimer
explicitement la matrice jacobienne approchée de Roe. Cette méthodologie est plus simple et
utile pour la résolution de problème de Riemann pour des problèmes physiques complexes tel
que celui de Barré de Saint-Venant [Delis-1998].
77
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
U t + F ( U )x = 0
U L si x < 0 (3.131.)
U ( x, 0 ) = U 0 ( x ) = U si x > 0
R
∂F ( U )
J (U) = (3.132.)
∂U
Ut + J ( U) Ux = 0 (3.133.)
Dans l’approche de Roe, on remplace la matrice jacobienne exacte par une matrice jacobienne
constante :
J% = J% ( U L , U R ) (3.134.)
où les indices (L,R) représentent les numéros des cellules de calcul à l’interface desquelles le
problème de Riemann est à résoudre. Cette matrice est une fonction des états constants UL et
UR.
U t +JU
% =0
x
UL si x<0 (3.135.)
U ( x,0 ) = U si x>0
R
Ce problème de Riemann est résolu de manière exacte. Le problème approché est obtenu en le
système initial de lois de conservation non linéaire par un système linéaire à coefficients
constants dans lequel intervient les valeurs initiales.
Pour un système d’équations quelconque, la matrice jacobienne de Roe doit satisfaire les
propriétés suivantes :
a% ( ) ≤ a% ( ) ≤ K ≤ a% ( )
1 2 m
(3.136.)
e% ( ) , e% ( ) , K , e% ( )
1 2 m
(3.137.)
78
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
J% ( U, U ) = J ( U ) (3.138.)
F ( U R ) − F ( U L ) = J% ( U R − U L ) (3.139.)
La propriété (1) est évidente car le problème approché doit préserver les propriétés
mathématiques du problème non linéaire originel.
La propriété (2) assure la consistance avec les lois de conservation.
La propriété (3) assure le caractère conservatif de la méthode. Il garantit également la
reconnaissance exacte de discontinuités isolées c’est à dire que si les données UL et UR sont
connectées par une onde simple isolée, alors le solveur de Riemann approché reconnaîtra cette
onde de manière exacte.
La construction de matrices qui satisfont ces trois propriétés pour un problème hyperbolique
général peut être très complexe et peu intéressante en terme de programmation. C’est pourquoi,
Roe et Pike ont construit une méthode qui rend inutile la construction explicite de cette matrice
jacobienne approchée.
Dans la méthode qu’ils ont proposée, la différence des valeurs constantes par morceau est
projetée sur la base des vecteurs propres :
m
∆U = U R − U L = ∑ α( ) e% ( )
k k
(3.140.)
k=1
D’autre part, le flux à l’interfaces des cellules de calcul i et i+1 est exprimé sous la forme :
et la différences des valeurs constantes par morceau des flux numériques est donnée par :
m
∆F = FR − FL = ∑ α(i+1) 2 a% ( ) e% ( )
k k k
(3.142.)
k=1
Pour mettre en œuvre la méthode de Roe il faut calculer les valeurs de e% ( ) , a% ( ) et α% ( ) . Puis, on
k k k
utilise ces valeurs pour appliquer un schéma numérique parmi ceux étudiés dans le paragraphe
4.3.
C’est la méthode de Roe que nous avons appliqué au système de Barré de Saint-Venant. On
présente cette application dans le paragraphe suivant.
79
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Les valeurs et vecteurs propres de la matrice jacobienne approchée sont exprimée sous la
forme:
a% 1,2
i+ 1
= u% 1,2
i+ 1
± c%1,2
i+ 1
2 2 2
1 (3.143.)
e%1,2 = % 1,2
i+ 1
2 a i+ 1
2
Les valeurs des coefficients de pondération introduits par la projection du vecteur écoulement
A
U = dans la base des vecteurs propres sont donnés par :
Q
α 1,2
= 2
(
∆ i + 1 Q + −u% i + 1 ± c% i + 1 ∆ i + 1 A
2 2
) 2
(3.144.)
i+ 1
2 ±2c% i + 1
2
Qi A i + Qi +1 A i +1
u% i + 1 = (3.145.)
2 A i +1 + A i
et :
( I1 )i +1 − ( I1 )i
g si A i +1 − A i ≠ 0
c% 2
i+ 1
= A i +1 − A i (3.146.)
2
c% 2 = c% 2 si A i +1 − A i = 0
i i +1
Ces valeurs permettent de calculer la fonction flux aux interfaces Fi+1 2 pour la mise en œuvre
du schéma numérique.
Dans ce qui suit, on expose le principe d’autres types de solveurs qui peuvent être
appliqués à la résolution des équations de Barré de Saint-Venant
C’est le solveur de Riemann le plus simple. La forme de la fonction flux numérique est
la suivante [Delis-2000] :
80
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
1 1
F% LF1 = FR + FL − ( U R − U L ) (3.147.)
i+
2 2 λ
Récemment, le solveur de Lax Friedrich a été amélioré pour devenir le solveur local de
Lax Friedrich (LLF) [Nujic-1995]. Cette forme dépend uniquement de la célérité de l’onde
locale la plus rapide et sa forme est donnée par :
1
F% LLF
1 = FR + FL − a max ( U R − U L ) (3.148.)
i+
2 2
(
a max = max a1i , a i2 , a1i +1 , a i2+1 ) (3.149.)
Cette forme est moins diffusive que la forme classique (3.147) de la méthode de Lax-
Friedrich et a toujours l’avantage de satisfaire la condition d’entropie.
Ce solveur a été mis au point par Harten, Lax et Van Leer [Delis-2000]. Il consiste à approcher
la solution de tout problème de Riemann par l’utilisation des trois états constants suivants :
n
U i pour x < bl 1 t
i+
n
2
U x 1 , t = U 1 pour b l
1 t<x<b r
1 t (3.150.)
i+ 2 i+ 2 i+
2
i+
2
U in+1 pour br 1 t < x
i+
2
br 1 = max u i +1 + ci +1 , a% 1 1 (3.151.)
i+ i+
2 2
bl 1 = min u i − ci , a% 2 1 (3.152.)
i+ i+
2 2
L’état moyen U n 1 est défini de telle sorte que le solveur de Riemann soit consistant avec la
i+
2
forme intégrale des lois de conservations, c’est à dire :
81
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
b r 1 U in+1 − bl 1 U in
i+ i+ Fi +1 − Fi
Un 1 = 2 2
− (3.153.)
i+
2 b r
1 −b l
1 br 1 − bl 1
i+ i+ i+ i+
2 2 2 2
et b% r,l sont les approximations numériques pour les vitesses des caractéristiques maximum et
i +1 2
b + 1 Fi − b − 1 Fi +1 b+ 1 b− 1
i+ i+ i+ i+
FHLLE
1 = 2 2
+ 2 2
( Ui+1 − Ui ) (3.154.)
i+
2 b+ 1 − b −
1 b+ 1 − b− 1
i+ i+ i+ i+
2 2 2 2
avec :
b + 1 = max b r 1 , 0 (3.155.)
i+
2 i+ 2
b − 1 = max b l 1 , 0 (3.156.)
i+
2 i+ 2
1 2
k
F HLLE = i +1 i ∑ q i + 1 αi+ 1 e% i + 1
+ − k k
1 F F (3.157.)
i+
2 2 k =1 2 2 2
avec :
b+ 1 + b− 1 b+ 1 b− 1
i+ i+ i+ i+
q1,21 = 2 2
a% 1,21 − 2 2 2
(3.158.)
i+
2 b+ 1 − b− 1 i+
2 b+ 1 − b− 1
i+ i+ i+ i+
2 2 2 2
Les bases de cette approche peuvent être trouvées dans les articles de Osher et Osher et
Salomon. Comme pour le solveur de Roe, le solveur de Osher est basé sur l’utilisation de la
matrice jacobienne du système d’équations de conservation à résoudre [Euleterio-1997]. En
particulier, la matrice jacobienne doit être diagonalisable.
Le principe de cette méthode consiste à décomposer la matrice jacobienne. Cette
décomposition est réalisée en distinguant les valeurs propres de valeurs positives et celles de
valeurs négatives. Ainsi, la matrice jacobienne est scindée sous la forme :
82
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
avec :
∂F + ∂F −
J+ (U) = ; J− (U) = (3.160.)
∂U ∂U
J + ( U ) possède des valeurs propres positives ou nulles et J − ( U ) des valeurs propres négatives
ou nulles. Cette décomposition peut s’appliquer indifféremment à des systèmes linéaires ou non
linéaires.
L’approche de Osher consiste ensuite à supposer qu’il existe des fonctions vectorielles flux
F + ( U ) et F − ( U ) qui vérifient :
F ( U ) = F+ ( U ) + F− ( U ) (3.161.)
Les valeurs initiales U L et U R du problème de Riemann pour les lois de conservation sont
notées :
U 0 ≡ U L = U in ; U1 ≡ U R = U in+1 (3.162.)
Dans ce cas, le flux numérique à utiliser pour la mise en œuvre d’un schéma numérique de type
Godunov est obtenu par l’utilisation des relations intégrales :
U1
∫ J ( U ) dU = F ( U ) − F ( U )
− − −
1 0 (3.163.)
U0
et
U1
∫ J ( U ) dU = F ( U ) − F ( U )
+ + +
1 0 (3.164.)
U0
U
1 1 1
F 1 = F ( U 0 ) + F ( U1 ) − ∫ J ( U ) dU (3.165.)
i+
2 2 2 U0
L’intégration par rapport à U est réalisée dans l’espace des phases où les composantes
de U sont des nombres réels et l’obtention d’un résultat dépend du domaine d’intégration
choisi.
Cette méthode a encore été peu appliquée au domaine de l’hydraulique. Pour sa mise en
œuvre dans le cas des ondes simples c’est à dire de compression et dépression, on renvoie à
[Euleterio-1997].
83
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Un grand nombre de limiteurs ont été mis au points (MINBEE, SUPERBEE, SUPRA, …)
[Roe-1983], [Van Leer-1973]. Nous ne les décrirons pas ici.
• Variation totale
1 +∞
TV ( u ) = lim sup u ( x + δ ) − u ( x ) dx
δ ∫−∞
(3.166.)
δ→0
TV ( u n ) = ∫
+∞
u ′ ( x ) dx (3.167.)
−∞
En fait, même pour les fonctions discontinues, la relation (3.167.) reste valable au sens
de la théorie des distributions [Smoller-1984]. Si u = u ( x, t ) , on peut généraliser les définitions
précédentes. En fait si l’on s’interesse à la convergence, il suffit de définir la variation totale à
un temps t fixé t = t n
On définit la variation totale de u(x,t) à un temps t fixé que l’on note TV ( u ( t ) ) par :
84
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
+∞
TV ( u n ) = ∑ u in+1 − u in (3.168.)
i =−∞i
• Schéma TVD
TV ( u n +1 ) ≤ TV ( u n ) (3.169.)
TV ( u n ) ≤ TV ( u n −1 ) ≤ K ≤ TV ( u 0 ) (3.170.)
La notion TVD a été introduite pour la première fois par Harten [Harten-1983] [Harten-
1984]. Une caractéristique essentielle d’un schéma TVD est qu’il définit des solutions sans
oscillations, en particulier à la traversée d'une région à forte variation de gradient. Un tel
schéma ne crée ni n’amplifie les extrémas, puisqu’il préserve la monotonie, c’est une propriété
que ne partagent pas l’ensemble des schémas numériques.
Le critère TVD d’un schéma est lié à l’existence d’une quantité minimale de viscosité
numérique nécessaire dans ce schéma. Les fondements théoriques des méthodes TVD ne sont
disponibles que pour les problèmes scalaires unidimensionnels. En pratique, l’expérience
accumulée dans de nombreuses applications à des problèmes non linéaires et/ou
multidimensionnels, a montré que la théorie unidimensionnelle scalaire sert bien de ligne de
conduite pour étendre ces idées.
85
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Ces méthodes numériques sont basées sur une approche géométrique. Les valeurs des
variables contenues dans les cellules de calculs ne sont plus considérées comme constantes par
morceaux mais comme des fonctions linéaires par morceaux comme représenté sur la figure ci-
dessous :
Les valeurs constantes par morceaux u in sont remplacées par une fonction linéaire u i ( x ) :
u i ( x ) = u in +
( x − xi ) ∆ (3.171.)
i
∆x
∆i
où est la pente de u i ( x ) dans la cellule de calcul i. Les valeurs de u i ( x ) au niveau des
∆x
extrémités de la cellule de calcul jouent un rôle important. Elles sont données par :
1 1
u iL = u i ( 0 ) = u in − ∆ i ; u iR = u i ( ∆x ) = u in + ∆ i (3.172.)
2 2
Ces valeurs sont appelées valeurs aux limites extrapolées. On note que l’intégrale de u i ( x )
dans la cellule Ii est identique à celle de u in ; en conséquence, la propriété de conservation est
préservée. La figure ci-dessous représente la linéarisation par morceaux dans trois cellules de
calculs successives.
86
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
u t + f ( u )x = 0
u i ( x ) x<0 (3.173.)
u ( x, 0 ) = u x
i +1 ( ) x>0
qui permet de calculer le flux de type Godunov au niveau des interfaces et d’obtenir une
solution du problème de Riemann généralisé. La solution ne contient plus de région uniforme
comme dans le cas du problème de Riemann classique, mais les caractéristiques sont
maintenant courbées comme le montre la figure ci-dessous :
On peut ainsi obtenir des schémas numériques d’un ordre de précision très élevé. Ces
méthodes numériques sont dites méthodes de type MUSCL (Monotone Upstream-centered
Scheme for Conservation Law).
Comme pour toutes les méthodes d’ordre élevé, se pose le problème de l’oscillation au
voisinage des zones à forts gradients. La technique classiquement utilisée pour supprimer ces
oscillations consiste à remplacer la pente ∆ i par une pente limitée ∆i par certaines contraintes
et permet de rendre ces méthodes TVD.
Théorème :
1
∆i = sign ∆ i − 1 + sign ∆ i + 1 × min βi − 1 . ∆ i − 1 ; βi + 1 . ∆ i + 1
2 2 2 2 2 2 2
(3.174.)
2 2
β 1 = , β 1 =
i−
2 1+ c i+
2 1− c
L’analyse TVD précédente est suffisante pour construire des schémas TVD amont et centré à
limitation de pente. Une autre interprétation du théorème précédent est très utile. On introduit le
limiteur de pente ξi tel que :
∆ i = ξi ∆ i (3.175.)
87
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
1 1
∆i = (1 + ω) ∆u i− 1 + (1 − ω) ∆u i+ 1 (3.176.)
2 2 2 2
où :
2β 1 r
i−
ξL ( r ) = 2
1 − ω + (1 + ω ) r
2β 1 r
i+
ξR ( r ) = 2
(3.178.)
1 − ω + (1 + ω ) r
∆ 1
i−
r= 2
∆ 1
i+
2
0, r≤0
ξ (r) = (3.179.)
min {ξ L ( r ) ,ξ R ( r )} , r>0
Ce limiteur de pente est analogue au limiteur de flux ULTRABEE. Il existe d’autre type de
limiteurs de pente. On citera les limiteurs de type Van-Leer, Albada, MINBEE ou encore
SUPERBEE., …[Euleterio-1997]
88
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Une des méthodes permettant de construire des schémas numériques qui présentent le
caractère TVD est l’approche à limitation de flux. Elle nécessite un flux numérique associé à un
schéma numérique d’ordre supérieur à un f HI1 et un flux numérique associé à un schéma
i+
2
monotone c’est à dire d’ordre un f LO1 . On peut définir un flux d’ordre élevé TVD par :
i+
2
HI LO
1 =f 1 +φ
f TVD fi+ 1 − fi + 1
LO
1 (3.180.)
i+ i+ i+
2 2 2 2 2
Méthode de résolution :
U t + F ( U ) x = G ( U )
(3.181.)
U ( x, t ) = U
n n
On veut faire évoluer la solution Un du temps t=tn à la solution Un+1 au temps t=tn+1 pour
le pas de temps ∆t = t n +1 − t n . Le schéma de base que l’on peut utiliser consiste à résoudre la
succession de problèmes aux valeurs initiales suivant :
89
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Les conditions initiales pour le problème homogène sont les conditions initiales du
problème complet et la solution après un temps ∆t est U(adv) qui peut être vue comme une
solution envisagée. Le second système peut ensuite être résolu pour un temps ∆t et avec la
condition initiale donnée par la résolution du premier système c’est à dire U(adv). La solution
ainsi obtenue peut être vue comme la solution approchée du problème complet.
∆t
A : U (i
adv )
= U in − Fi + 1 − Fi − 1
∆x 2 2 (3.183.)
S: U in +1 = U (i
adv )
(
+ ∆tG U (i
adv )
)
le schéma complet s’écrit finalement :
∆t
U in +1 = U in − ( adv )
Fi + 1 − Fi − 1 + ∆tG U i
∆x 2 2
( ) (3.184.)
90
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
∆t % % %n
U in +1 = U in − Fi + 1 − Fi − 1 + ∆tG i (3.185.)
∆x 2 2
Les flux aux interfaces sont calculés en déterminant les valeurs des quantités présentées
dans la partie 3 § [Link] grâce à la technique du solveur de Roe et l’association du terme
source au schéma numérique est réalisée par la méthode décrite dans la partie 3 § 3.6.
Dans ce qui suit, on précise la méthode utilisée pour calculer le terme source et on
détaille les schémas numériques utilisés pour la résolution du système de Barré de Saint-Venant
rappelé ci-dessous :
U t + F(U )x = G (U ) ⇔ U t + J (U )U x = G (U ) (186.)
avec :
A Q
U = ,
F(U ) = Q 2
Q + gI1
A
On expose dans ce qui suit la manière dont le terme source, dont on rappelle l’expression,
est intégré au schéma numérique :
Qdev
G= (3.187.)
gI 2 + gA ( S0 − Sf ) + QQdev
A
Avec :
Q2 n 2
Sf = 4
: Relation de Manning-Strickler ;
A2R h 3
h
I1 = ∫ (h − η)Bdη : Terme de pression hydrostatique ;
0
h
∂B ∂I
I2 = ∫ ( h − η) dη = 1 : Terme de pression latérale ;
0
∂x ∂x
91
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
dQ 1− W 3(1 − y)
1/2
Qdev = = −0.6n*c w gH (y − W) ×
3 3/2
× 1 − θ variation de débit du
dx 3 − 2y − W y−W
au déversement.
La valeur du terme source dans une cellule de calcul G % n est approchée par une
i
moyenne arithmétique des valeurs obtenues en utilisant les solutions des deux problèmes de
% n et G
Riemann (i-1, i) et (i, i+1) c’est à dire. G % n . Ces valeurs sont calculées en faisant
i+ 1 i− 1
2 2
et en supposant que le canal soit localement rectangulaire le terme source est approché par :
%
Q %
( %
dev A i ± 1 , Bi ± 1 , u i ± 1
2 2
%
2
)
%n
G = %
( ) (3.188.)
))
% %
(
QQ
( ) (
dev A i ± 1 , Bi ± 1 , u i ± 1
i± 1 %
2
% % % %
gI2 A i ± 1 , Bi ± 1 + gA i ± 1 S0 − Sf A i ± 1 , Bi ± 1 +
% % % 2 2 2
2 2 2 2 2 A %
i± 1
2
où :
%
A = A i A i +1 (3.189.)
i+ 1
2
B% n = Bi + Bi +1 (3.190.)
i+ 1
2 2
Q A i + Qi +1 A i +1
u% i + 1 = i (3.191.)
2 A i +1 + A i
Les flux aux interfaces obtenus par application de la méthode du solveur de Roe (Partie 3 §
[Link]) sont modifiés par l’ajout d’un terme correctif introduit dans le vecteur D. L’expression
du flux est donc :
1 % %
F% 1 = Fi + Fi ±1 + R i ± 1 Di ± 1 (3.192.)
i±
2 2 2 2
92
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
La fonction ψ est la correction d’entropie pour les valeurs propres a ik+1/2 et de la forme :
a si a ≥ δ
ψ (a ) = 2 2 (3.194.)
( a + δ ) /2δ si a <δ
où δ est un petit nombre positif que l’on prend ainsi [Harten et Hyman-1983] :
La fonction limiteur de flux Lki±1/ 2 contrôle les termes de second ordre de façon à ce qu’un
résultat lisse et non oscillant soit garanti même en présence de discontinuités :
où :
d ik+1/ 2 = σ ( a% ik+1/ 2 )( Lki +1 + Lki ) − ψ ( a% ik+1/ 2 + γ ik+1/ 2 ) α ik+1/ 2 pour k={1, 2} (3.198.)
Avec :
93
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
Nous écrirons le schéma sous la forme proposée par Garcia-Navarro [Garcia-Navarro and
Alcrudo-1992] :
1 λ
U in +1 = (U i(1) + U i(2) ) + [R i +1/ 2 Di +1/ 2 − R i −1/ 2 Di −1/ 2 ] (3.203.)
2 2
(d ik+1/ 2 ) PC = ψ(a% ik+1/ 2 )[1 − λ a% ik+1/ 2 ][1 − Lki +1/ 2 ]α ik+1/ 2 pour k={1, 2} (3.204.)
où :
α ik+12 −s
rik+1/ 2 = et s = sgn(a ik+1/ 2 ) (3.206.)
α ik+1/ 2
R i +1/ 2 et Di +1/ 2 peuvent être évalués en U n ou en U (2) , nous avons choisi d’utiliser les valeurs
de U (2) .
La formule d’interpolation MUSCL, présentée par Van Leer, fournit une famille de
schémas paramétrés du second ordre et un du troisième ordre. Plutôt qu’un limiteur de flux,
l’approche plus géométrique suggère un limiteur de pente qui limite les fluctuations des
variables dépendantes.
Dans les cellules de calcul, les valeurs ne sont plus constantes. On introduit une pente
sur la valeur du vecteur écoulement. On remplace les arguments U i +1 et U i du flux numérique
par U iR+1 et U iL+1 où les arguments U R et U L sont calculés ainsi :
94
Partie 3 : Méthodes mises en œuvre dans les outils numériques réalisés
1
U iL+1/ 2 = U i + [(1 − m)∆ i+−1/ 2 + (1 + m)∆ i−+1/ 2 ] (3.207.)
4
1
U iR+1/ 2 = U i +1 − [(1 − m)∆ i−−3/ 2 + (1 + m)∆ i++1/ 2 ] (3.208.)
4
où :
3− m
1≤ β ≤ avec m ≠1 (3.210.)
1− m
1
La précision en espace est déterminée par la valeur de m, ici on pose m = et on obtient une
3
précision du troisième ordre.
Ici â ik+1 , R̂ i +1/ 2 et αˆ ik+1/ 2 sont évalués comme avant, mais les arguments U i +1 et U i ont été
remplacés respectivement par U iR+1/ 2 et U iL+1/ 2 .
Ce sont ces méthodes numériques et en particulier la méthode amont que nous avons
utilisé pour la construction de nos modèles d’écoulements en collecteurs et déversoirs.
5 CONCLUSION
Dans cette partie de la thèse nous avons mis en évidence les limites des méthodes
numériques classiques comme celles des caractéristiques et des différences finies dans le but de
résoudre les équations de Barré de Saint-Venant appliquées au domaine de l’assainissement.
Ensuite, nous avons décrit une classe de méthodes numériques (méthodes à capture de chocs de
type TVD) qui permettent de résoudre des systèmes d’équations hyperboliques pouvant faire
apparaître des discontinuités comme c’est le cas des équations de Barré de Saint-Venant. Pour
finir, nous avons décrit, les méthodes effectivement mises en œuvre dans les modèles
numériques réalisés.
95
Partie 4 : Validation des outils numériques réalisés
PARTIE 4
VALIDATION
DES
OUTILS DE CALCULS
REALISES
96
Partie 4 : Validation des outils numériques réalisés
1 INTRODUCTION
Dans cette partie de la thèse on a pour objectif, dans un premier temps de présenter les
modèles d’écoulements que nous avons construit mais également de montrer l’aptitude des
ces logiciels à reproduire les débits et hauteurs d’eau mesurés sur plusieurs dispositifs
expérimentaux. Ces dispositifs sont un canal de travaux pratiques, un banc d’essai physique
de déversoir et pour finir un Venturi.
Les hauteurs d’eau ont été mesurées par un dispositif caméra permettant d’obtenir une
vision 3D de la surface libre et les débits par l’association d’un capteur ultrason et d’un
Venturi. Dans d’autres cas, les hauteurs d’eau ont été mesurées par limnimétrie.
Dans un deuxième temps, on expose dans ce chapitre les résultats expérimentaux et
numériques obtenus dans le cadre de trois campagnes de mesures qui ont été menées :
• Validation des modèles de collecteurs programmés à l’aide de 4 méthodes numériques
différente grâce à un canal de TP. Les résultats nous ont permis de choisir l’un de ces
modèles pour la suite du développement.
• Validation de la prise en compte du changement de section par les modèles grâce à un
Venturi.
• Validation du modèle de déversoir grâce au banc d’essai physique de déversoir.
97
Partie 4 : Validation des outils numériques réalisés
De manière simplifiée, les valeurs de la hauteur d’eau et du débit qui transitent dans un
système de canaux et collecteurs sont obtenues par l’intermédiaire des relations décrites dans
les chapitres précédents. Tous ces modèles peuvent être scindés en deux strates :
De manière générale, ces outils de calcul que nous avons réalisé ont tous une structure
qui peut être représentée de la manière suivante :
TRONCON
CALC. TRON.
GRAPHIQUE
Le module HYDRO a pour but de fixer l’hydrogramme d’entrée dans le système. Ce débit
peut être constant (temps sec) ou fortement variable (temps de pluie). De manière générale,
dans le cadre d’une application aux réseaux d’assainissement, les hydrogrammes mettent en
évidence des pics de débit qui peuvent avoir de graves conséquences (débordements du
réseaux, inondations
98
Partie 4 : Validation des outils numériques réalisés
Le module COND. INI. permet d’appliquer une hauteur d’eau dans le canal au temps initial.
La valeur de la hauteur d’eau dans un tronçon peut être imposée ou fixée à la hauteur
normale.
Le module COND. LIM. permet de fixer les conditions aux limites du système. C’est le point
le plus délicat de la mise en œuvre du modèle. Celles-ci sont décrites dans le paragraphe 2.3
de ce chapitre.
Le module AF. COND. LIM. permet d’affecter les valeurs des conditions aux limites aux
cellules des deux tronçons fictifs. En effet, ces valeurs sont recalculées à chaque pas de temps.
Dans le cas d’une hauteur normale ou critique, on applique cette valeur calculée au pas de
temps suivant. Dans le cas d’une hauteur imposée on analyse la possibilité, en fonction du
point de contrôle, de l’apparition d’un ressaut hydraulique.
Dans le module CALC. CEL, sont déterminées les valeurs constantes par morceau dans les
cellules de calcul. Il s’agit de la pression hydrostatique, de la largeur au miroir, des
composantes du vecteur flux et des valeurs propres de la matrice jacobienne. (cf. partie 2). La
définition de ces valeurs doit tenir compte du problème étudié. La géométrie du problème, en
particulier la pente et la nature du terrain, définissent ces valeurs et permettent la
modélisation du problème.
Dans le module CALC INTERF, les valeurs de toutes les quantités nécessaires au calcul des
flux aux interfaces ainsi que le vecteur source sont calculées.
Dans le module CALC TRON, sont évaluées les valeurs du vecteur écoulement au pas de
temps t+1 puis représentées graphiquement par le module GRAPHIQUE.
Les valeurs de hauteurs normales et critiques mais également les relations liant section
mouillée et hauteur d’eau, rayon hydraulique et hauteur d’eau, largeur au miroir et hauteur
d’eau et leurs inverses ne peuvent pas dans la plupart des cas être obtenues de manière directe.
Par exemple, dans le cas de la hauteur normale, celle-ci est obtenue grâce à la relation de
Manning-Strickler. Pour résoudre cette relation afin d’obtenir la hauteur normale il faut
utiliser un solveur. Or l’utilisation de solveurs coûte cher en temps de calcul. D’autres
exemples sont ceux des canaux circulaires ou ovoïdes pour lesquels on ne dispose pas de
relation directe liant rayon hydraulique et hauteur d’eau ou encore largeur au miroir et hauteur
d’eau.
C’est pourquoi, nous avons calé des polynômes permettant d’obtenir directement ces
valeurs : hn=f(Q) ou Rh=f(h) ou h=f(S), …
Les polynômes qui ont été utilisés dans les logiciels que nous avons réalisés sont
regroupés dans l’annexe 4.
Les conditions initiales correspondent à une hauteur d’eau imposée dans toutes les
cellules de calculs de tous les tronçons du système modélisé. On peut choisir entre la hauteur
normale et une hauteur d’eau imposée dans chaque tronçon.
99
Partie 4 : Validation des outils numériques réalisés
[Link] Introduction
Les conditions aux limites sont affectées à deux tronçons dits « fictifs » constitués de
deux cellules de calculs qui sont situés à l’amont et à l’aval du système de collecteurs et
déversoirs « réels » connectés en série. Ces tronçons fictifs sont nécessaires au logiciel car le
schéma numérique mis en œuvre est un schéma à 5 points. Il nécessite donc la donnée des
valeurs contenues dans deux cellules de calcul à l’amont et à l’aval de la cellule dans laquelle
on veut obtenir la hauteur d’eau et le débit. Ainsi, pour le calcul des variables dans la
première cellule du premier tronçon réel on aura besoin des valeurs contenues dans les deux
cellules du tronçon fictif amont. De même, pour le calcul des variables dans la dernière cellule
du dernier tronçon réel, on aura besoin des valeurs contenues dans le tronçon fictif aval. Ceci
est illustré par la figure suivante.
Tronçon Tronçon
Tronçon réel 1 Tronçon réel n
fictif amont fictif aval
On donne le choix à l’utilisateur de fixer les conditions aux limites amont et aval du
système de collecteurs. Les valeurs applicables sont les suivantes :
Dans le cas de l’application d’une hauteur imposée, celle-ci peut être de plusieurs types.
Elle peut être constante et représentera une hauteur d’eau invariable comme dans le cas d’un
bassin ou d’un plan d’eau. Elle peut également représenter une hauteur de seuil ou la hauteur
d’ouverture d’une vanne. Dans ce cas, la présence de ces ouvrages a été intégrée au logiciel
grâce à des lois de seuil et de vannes. Ces lois de seuil et de vannes ont été utilisées dans la
phase de validation du logiciel sur un canal expérimental. Pour faire apparaître les divers
types de courbes de remous et les ressauts hydrauliques, nous avons parfois mis en place une
vanne rectangulaire à l’amont du canal et un seuil à l’aval.
100
Partie 4 : Validation des outils numériques réalisés
Q, géom.
hn hc himp
Nous avons imaginé une méthode de gestion des conditions aux limites adaptées à la
physique du problème considéré. Nous n’avons pas utilisé les conditions aux limites de type
Neuman ou Dirichlet de manière classique. Les conditions utilisées peuvent plutôt être
qualifiées de « mixtes ». Nous avons mis en œuvre une méthode d’analyse de l’applicabilité
ou non d’une condition à la limite amont ou aval qui est fonction des conditions hydrauliques
dans les canaux et collecteurs.
En effet, avant d’affecter les valeurs données par l’utilisateur aux tronçons fictifs amont
et aval, on teste si ces conditions sont applicables en fonction des conditions hydrauliques de
l’écoulement. Par conditions hydrauliques, on entend les valeurs des hauteurs normales et
critiques et ainsi de la présence ou non d’un point de contrôle au niveau de la limite
considérée. De plus, dans le cas d’une hauteur d’eau imposée invariable, on teste si celle-ci
est de nature à donner naissance à un ressaut hydraulique.
Si cette analyse permet d’affirmer que la condition à la limite est applicable, dans ce cas
la valeur est affectée aux cellules du tronçon fictif. C’est dans ce cas une condition de type
Dirichlet.
Si cette condition n’est pas applicable c’est à dire qu’au niveau hydraulique le point de
contrôle ne se situe pas en l’extrémité considérée du domaine de calcul alors une progression
linéaire entre les cellules du tronçon fictif et la valeur contenue dans la cellule du tronçon réel
la plus proche est calculée. Cette progression linéaire qui constitue une condition à la limite
de type Neuman est calculée à partir de la valeur contenue dans la cellule de calcul du tronçon
réel.
Cette méthode de gestion des conditions aux limites est nécessaire car on travaille en
régime transitoire. Ainsi, une condition à la limite qui était applicable à une extrémité du
domaine à un instant donné peut ne plus l’être au pas de temps suivant.
101
Partie 4 : Validation des outils numériques réalisés
hn
Linéarisation hn
Figure 44 : Choix de la hauteur normale pour condition à la limite amont
Si l’on veut appliquer la hauteur normale comme condition à la limite amont, on teste
si le canal est à pente forte ou à pente faible c’est à dire si la pente critique est supérieure ou
inférieure à la pente du canal ou encore si la hauteur normale est respectivement inférieure ou
supérieure à la hauteur critique.
Si le canal est caractérisé par une pente faible, le régime d’écoulement sera fluvial s’il
n’est pas perturbé. En conséquence, le point de contrôle est à l’aval et la hauteur normale ne
peut être appliquée à l’amont et on effectue une linéarisation entre les deux cellules du
tronçon fictif amont et la première cellule du premier tronçon réel.
Si le canal est à pente forte, le régime d’écoulement est à priori torrentiel et le point de
contrôle est à l’amont. En conséquence, la hauteur normale peut être appliquée en ce point.
hc
Oui h>hc
Non
Linéarisation hc
Figure 45 : Choix de la hauteur critique pour condition à la limite amont
102
Partie 4 : Validation des outils numériques réalisés
himp
Oui Non
h>hc
103
Partie 4 : Validation des outils numériques réalisés
hn
Linéarisation hn
Figure 47 : Choix de la hauteur normale pour condition à la limite aval
Si l’on veut appliquer la hauteur normale comme condition à la limite aval, on teste si
le canal est à pente forte ou à pente faible.
Si le canal est caractérisé par une pente faible, le régime d’écoulement sera fluvial s’il
n’est pas perturbé. En conséquence, le point de contrôle est à l’aval et la hauteur normale peut
être appliquée en ce point.
Si le canal est à pente forte, le régime d’écoulement est a priori torrentiel et le point de
contrôle est à l’amont. En conséquence, la hauteur normale ne peut pas être appliquée en ce
point et on effectue une linéarisation entre les deux cellules du tronçon fictif aval et la
dernière cellule du premier tronçon réel.
hc
Linéarisation hc
Figure 48 : Choix de la hauteur critique pour condition à la limite aval
104
Partie 4 : Validation des outils numériques réalisés
himp
Oui Non
h<hc
On peut illustrer l’application de ces conditions aux limites grâce au logiciel nommé
hsl qui permet de reproduire les courbes de remous susceptibles d’apparaître dans un unique
canal rectangulaire en écoulement à surface libre. Le logiciel est présenté dans l’annexe 5
avec quelques exemples.
105
Partie 4 : Validation des outils numériques réalisés
106
Partie 4 : Validation des outils numériques réalisés
Dans le but de valider les schémas numériques programmés nous avons tenté de
reproduire différentes courbes de remous sur un canal de travaux pratiques.
En faisant varier débit, pente, en posant ou ne posant pas des vannes à l’entrée et des
vannes ou des seuils à la sortie, nous avons obtenu différentes courbes de remous appelées S1,
S3, M1, M2, M3 (décrites plus haut ), ainsi que différents ressauts. Seule la courbe de remous
de type S2 n’a pas pu être reproduites en raison de la proximité des valeurs des hauteurs
normale et critique dans le cas du canal utilisé.
De plus, nous avons placé selon les cas, pour se mettre dans des conditions propices à
l’apparition d’une courbe de remous donnée, une vanne à 25cm de l’amont et/ou un seuil à
l’aval. Au niveau numérique, la présence du seuil et/ou de la vanne a été traduite par
l’application de conditions aux limites obtenues grâce à une loi de vanne à l’amont et une loi
de seuil à l’aval.
Nous avons comparé, après stabilisation, les hauteurs réelles avec celles calculées par
les différentes méthodes, et la position de l’éventuel ressaut hydraulique. Nous avons
également comparé les temps de propagation des ressauts réels et calculés.
107
Partie 4 : Validation des outils numériques réalisés
4.3 Erreur
- L’aiguille du débitmètre oscillant beaucoup, nous avons estimé le débit fourni par la
pompe à l’aide d’un seau et d’un chronomètre ; des erreurs de l’ordre de 10% de la valeur
mesurée sont donc possibles.
- Les pentes considérées étant parfois très faibles, et leur obtention quelque peu délicate,
l’erreur sur la pente mesurée n’est peut-être pas négligeable.
- Les seuils et les vannes n’étant pas tout à fait étanches, des différences entre la distance
théorique du ressaut de l’entrée du canal, et celle effectivement observée sont constatées.
- Une erreur d’estimation des coefficients intervenant dans les lois de seuil et de vannes est
possible. Concernant la loi de seuil, le coefficient de débit a été calé par rapport aux
résultats expérimentaux.
- Des oscillations verticales de la surface libre dues à la présence du ressaut.
Ces sources d’erreurs que nous avons limités le plus possible nous on conduit à
considérer qu’un écart de ± 2mm par rapport à la valeur lue était représentatif de la réalité.
4.4.1 Validation des schémas numériques pour le régime non uniforme permanent
Les tableaux regroupant les valeurs de hauteur d’eau obtenues par les quatre méthodes
numériques et expérimentalement pour les cas envisagés dans ce paragraphe sont réunis dans
l’annexe 6. Les quatre logiciels utilisés pour cette phase de validation étaient du même type
que celui présenté en annexe 5 programmé à l’aide de la méthode upwind mais trois autres
programmés à l’aide des méthodes symmetric, MUSCL et Mac-Cormack. Pour obtenir les
résultats exposés dans les paragraphes suivants, les nombres de Courant que nous avons
utilisé étaient de 1 pour les méthodes upwind, symmetric et MUSCL et de 0.1 pour la
méthode Mac-Cormack en raison de sa tendance à l’oscillation.
[Link] Courbe S1
Cette courbe apparaît dans le cas d’un canal à pente forte lorsque l’on franchit la
hauteur critique pour passer en régime fluvial. En particulier, on voit apparaître une telle
courbe à l’aval d’un ressaut hydraulique présent dans un canal à pente forte. C’est ce que nous
avons fait expérimentalement en essayant de repousser le ressaut vers l’amont en réglant la
vanne et le seuil. Les conditions expérimentales sont regroupées dans le tableau suivant :
108
Partie 4 : Validation des outils numériques réalisés
Courbe S1
Longueur (m) Largeur (m) Strickler Débit (m3/h) Pente (m/m)
Remarques
On impose une condition à la limite aval grâce à un seuil de 5 cm et une condition à la limite amont grâce
à une vanne de 3,5 cm placée à 0,25 m.
[Link].2 Résultats
La figure suivante représente les valeurs lues avec et sans les erreurs ainsi que les
valeurs calculées par les 4 schémas numériques.
Courbe de remous
+ Valeur expérimentale.
0.12
+ Valeur expérimentale avec
erreur ajoutée ou retranchée.
0.1
0.08
Hauteur Méthode Mac Cormack
d’eau (m)
0.06
Méthode Symmetric
Méthode upwind
0.04
Méthode MUSCL
0.02
0
0 1 2 3 4 5 6
Position (m)
109
Partie 4 : Validation des outils numériques réalisés
- Upwind : 2,25m
- MUSCL : 2 ,25m
- Symmetric : 1,5m
- McCormack : 1,75m
[Link] Courbe S3
Cette courbe apparaît dans le cas d’un canal à pente forte lorsque l’on se situe en-
dessous de la hauteur normale. Au niveau expérimental, nous avons forcé la hauteur d’eau à
se situer en dessous de la hauteur normale grâce à une vanne placée à l’amont du canal. Les
données expérimentales sont regroupées dans le tableau suivant :
Courbe S3
3
Longueur (m) Largeur (m) Strickler Débit (m /h) Pente (m/m)
Remarques
On impose une condition à la limite amont grâce à une vanne de 1 cm placée à 0,25 m.
[Link].2 Résultats
La figure suivante représente les valeurs lues avec et sans les erreurs ainsi que les
valeurs calculées par les 4 schémas numériques.
110
Partie 4 : Validation des outils numériques réalisés
Courbe de remous
0.02 + Valeur expérimentale.
0.014
0.012
Hauteur Méthode Mac Cormack
0.01
d’eau (m)
Méthode Symmetric
0.008
Méthode upwind
0.006 Méthode MUSCL
0.004
0.002
0
0 1 2 3 4 5 6
Position (m)
Figure 52 : Lignes d’eau obtenues numériquement et expérimentalement pour la courbe S3
Les quatre schémas numériques produisent des résultats satisfaisants pour la courbe
S3. Néanmoins, on note une nouvelle fois que le schéma Mac-Cormack montre une tendance
à l’oscillation.
[Link] Courbe M1
Cette courbe apparaît dans le cas d’un canal à pente faible lorsque l’on se situe au-
dessus de la hauteur normale. Au niveau expérimental, nous avons forcé la hauteur d’eau à se
situer au-dessus de la hauteur normale grâce à un seuil placé à l’aval du canal. Les données
expérimentales sont regroupées dans le tableau suivant :
Courbe M1
3
Longueur (m) Largeur (m) Strickler Débit (m /h) Pente (m/m)
Remarques
On impose une condition à la limite aval grâce à un seuil de 11 cm placée en bout de canal.
111
Partie 4 : Validation des outils numériques réalisés
[Link].2 Résultats
La figure suivante représente les valeurs lues avec et sans les erreurs ainsi que les valeurs
calculées par les 4 schémas numériques.
Courbe de remous
0.16
+ Valeur expérimentale.
0.14
Hauteur
d’eau(m) 0.13
Méthode Mac Cormack
0.12 Méthode Symmetic
Méthode upwind
0.11 Méthode MUSCL
0.1
0 1 2 3 4 5 6
Position (m)
[Link] Courbe M2
Cette courbe apparaît dans le cas d’un canal à pente faible lorsque l’on se situe entre
les hauteurs normale et critique. En particulier, elle apparaît à l’aval d’un ressaut hydraulique
avec en bout de canal une chute brusque. Ce sont ces conditions que l’on a reproduit au
niveau expérimental.
112
Partie 4 : Validation des outils numériques réalisés
Courbe M2
Longueur (m) Largeur (m) Strickler Débit (m3/h) Pente (m/m)
Remarques
On impose une condition à la limite amont grâce à une vanne de 2 cm à 0,25 m de l'amont du canal.
[Link].2 Résultats
La figure suivante représente les valeurs lues avec et sans les erreurs ainsi que les
valeurs calculées par les 4 schémas numériques :
Courbe de remous
0.06
+ Valeur expérimentale.
+ Valeur expérimentale
0.05 avec erreur ajoutée ou
retranchée.
0.04
Hauteur
d’eau (m) 0.03
Méthode upwind
0.01
Méthode MUSCL
0
0 1 2 3 4 5 6
Position (m)
113
Partie 4 : Validation des outils numériques réalisés
[Link] Courbe M3
Cette courbe apparaît dans le cas d’un canal à pente faible lorsque l’on se situe en-
dessous de la hauteur critique. En particulier, elle apparaît à l’amont d’un ressaut hydraulique.
Au niveau expérimental, on a fait apparaître un ressaut en imposant à l’amont une hauteur
inférieure à la hauteur critique par l’intermédiaire d’une vanne que l’on a réglé de telle sorte
que le ressaut se situe le plus à l’aval possible du canal.
Courbe M3
Longueur (m) Largeur (m) Strickler Débit (m3/h) Pente (m/m)
Remarques
On impose une condition à la limite amont grâce à une vanne de 3 cm à 0,25 m de l'amont du canal.
[Link].2 Résultats
La figure suivante représente les valeurs lues avec et sans les erreurs ainsi que les
valeurs calculées par les 4 schémas numériques :
114
Partie 4 : Validation des outils numériques réalisés
Courbe de remous
0.06
+ Valeur expérimentale.
+ Valeur expérimentale
0.05 avec erreur ajoutée ou
retranchée.
0.04
Hauteur
d’eau (m) 0.03
Méthode Mac Cormack
Méthode upwind
0
0 1 2 3 4 5 6
Position (m)
Figure 55 : Lignes d’eau obtenues numériquement et expérimentalement pour la courbe
M3
Comme pour la courbe M2, les méthodes upwind et symmetric produisent les valeurs
les plus proches des valeurs expérimentales tandis que les résultats obtenus grâce aux
méthodes MUSCL et Mac-Cormack sont plus éloignés des valeurs attendues.
propagations des ressauts irréalistes. Ceux-ci étaient beaucoup plus courts (Au moins divisé
par deux).
Dans ce qui suit, nous présentons les cas de quatre ressauts. Pour deux d’entre eux, le
canal a été réglé en pente faible et pour les deux autres en pente forte.
Le tableau suivant regroupe les données expérimentales appliquées au canal :
Cas 1
3
Longueur (m) Largeur (m) Strickler Débit (m /h) Pente (m/m)
Remarques
On impose une condition à la limite amont grâce à une vanne de 2,9 cm à 0,25 m de l'amont du canal.
Expérimental 65 s
Méthode Symetric 66 s
Méthode MUSCL 60 s
Méthode Upwind 62 s
Les trois figures suivantes présentent l’évolution de la ligne d’eau au cours du temps.
Celles-ci ont été obtenues avec les schémas numériques symetric, upwind et MUSCL.
Les courbes représentent la ligne d’eau après 5s, 10s, 20s, 30s, 40s et pour finir à la
stabilisation du ressaut.
Après 5s :
Après 10s :
Après 20s :
Après 30s :
Après 40s :
Stabilisation :
116
Partie 4 : Validation des outils numériques réalisés
Méthode muscl
0.07
0.06
0.05
0.04
Hauteur
0.03
d’eau (m)
0.02
0.01
-0.01
0 1 2 3 4 5 6
Position (m)
Méthode upwind
0.07
0.06
0.05
0.04
Hauteur
d’eau (m) 0.03
0.02
0.01
-0.01
0 1 2 3 4 5 6
Position (m)
M éthode sym etric
0.07
0.06
0.05
0.04
H auteur
0.03
d’eau (m )
0.02
0.01
-0.01
0 1 2 3 4 5 6
Position (m )
117
Partie 4 : Validation des outils numériques réalisés
Cas 2
Longueur (m) Largeur (m) Strickler Débit (m3/h) Pente (m/m)
Remarques
On impose une condition à la limite amont grâce à une vanne de 2,2 cm à 0,25 m de l'amont du canal.
Temps de propagation
Expérimental 56 s
Méthode Symetric 54 s
Méthode MUSCL 57 s
Méthode Upwind 56 s
Les trois figures suivantes présentent l’évolution de la ligne d’eau au cours du temps.
Celles-ci ont été obtenues avec les schémas numériques symetric, upwind et MUSCL.
Les courbes représentent la ligne d’eau après 3s, 5s, 10s, 20s, 30s et pour finir à la
stabilisation.
Après 5s :
Après 10s :
Après 15s :
Après 20s :
Après 25s :
Stabilisation :
118
Partie 4 : Validation des outils numériques réalisés
0.06
0.05
0.04
Hauteur
0.03
d’eau (m)
0.02
0.01
-0.01
0 1 2 3 4 5 6
Position (m)
Méthode muscl
0.07
0.06
0.05
0.04
Hauteur 0.03
d’eau (m)
0.02
0.01
-0.01
0 1 2 3 4 5 6
Position (m)
Méthode upwind
0.07
0.06
0.05
0.04
Hauteur 0.03
d’eau (m)
0.02
0.01
-0.01
0 1 2 3 4 5 6
Position (m)
119
Partie 4 : Validation des outils numériques réalisés
Cas 3
3
Longueur (m) Largeur (m) Strickler Débit (m /h) Pente (m/m)
Remarques
On impose une condition à la limite aval grâce à un seuil de 7,9 cm placé à l'extrêmité du canal.
Temps de propagation
Expérimental 33 s
Méthode Symetric 35 s
Méthode MUSCL 35 s
Méthode Upwind 35 s
Les trois figures suivantes présentent l’évolution de la ligne d’eau au cours du temps.
Celles-ci ont été obtenues avec les schémas numériques symetric, upwind et MUSCL.
Les courbes représentent la ligne d’eau après 5s, 10s, 15 s, 20s, 25 s et pour finir à la
stabilisation.
Après 5s :
Après 10s :
Après 15s :
Après 20s :
Après 25s :
Stabilisation :
120
Partie 4 : Validation des outils numériques réalisés
Méthode upwind
0.1
0.05
Hauteur 0
d’eau (m)
-0.05
-0.1
0 1 2 3 4 5 6
Position (m)
Méthode symetric
0.1
0.05
Hauteur
0
d’eau (m)
-0.05
-0.1
0 1 2 3 4 5 6
Position (m)
Méthode muscl
0.1
0.05
Hauteur 0
d’eau (m)
-0.05
-0.1
0 1 2 3 4 5 6
Position (m)
Cas 4
3
Longueur (m) Largeur (m) Strickler Débit (m /h) Pente (m/m)
Remarques
On impose une condition à la limite aval grâce à un seuil de 6 cm placé à l'extrêmité du canal.
Temps de propagation
Expérimental 30 s
Méthode Symetric 32 s
Méthode MUSCL 28 s
Méthode Upwind 29 s
Les trois figures suivantes présentent l’évolution de la ligne d’eau au cours du temps. Celles-
ci ont été obtenues avec les schémas numériques symetric upwind et MUSCL.
Les courbes représentent la ligne d’eau après 5s, 10s, 15s, 20s et pour finir à la stabilisation.
Après 5s :
Après 10s :
Après 15s :
Après 20s :
Stabilisation :
122
Partie 4 : Validation des outils numériques réalisés
Méthode symetric
0.1
0.05
Hauteur
0
d’eau (m)
-0.05
-0.1
0 1 2 3 4 5 6
Position (m)
Méthode upwind
0.1
0.05
Hauteur 0
d’eau (m)
-0.05
-0.1
0 1 2 3 4 5 6
Position (m)
Méthode muscl
0.1
0.05
Hauteur 0
d’eau (m)
-0.05
-0.1
0 1 2 3 4 5 6
Position (m)
123
Partie 4 : Validation des outils numériques réalisés
On présente à l’aide du tableau suivant une synthèse des résultats obtenus concernant
les divers types de courbes de remous obtenues en régime non uniforme permanent ainsi que
ceux obtenus en régime transitoire pour les 4 cas de ressauts hydrauliques que nous avons
envisagé.
Signalétique
La méthode (Mc Cormack) choisit d’intégrer le terme source dans les calculs de flux à
chaque itération, mais la précision n’est pas très bonne et on observe des oscillations dans les
courbes représentant les hauteurs et les débits ; de plus, la convergence n’est obtenue qu’avec
des C.F.L. très petits (0,1) ce qui engendre un temps de calcul important.
La méthode MUSCL converge très rapidement pour des C.F.L. proches de 1 ; elle semble
donc adaptée dans le cadre de calculs rapides de moindre précision.
124
Partie 4 : Validation des outils numériques réalisés
- La méthode Mac Cormack produit des oscillations dans les courbes de débit et de hauteurs
d’eau, les autres non. Cependant, ces oscillations ont une amplitude de plus en plus faible
au cours des itérations.
- En règle générale, la méthode Upwind est la plus précise.
La comparaison des résultats numériques avec l’expérience donne des erreurs de
moins de 5% (par rapport à la longueur du canal) sur la position du ressaut et sur les hauteurs
d’eau, ce qui est satisfaisant. Les allures des courbes de hauteur d’eau obtenues sont conforme
à notre attente
4.6 Conclusion
Pour la validation des schémas numériques en régimes transitoires les trois méthodes
numériques donnent des résultats en bon accord avec la réalité.
Les résultats sont donc satisfaisants en ce qui concerne les écoulements à surface libre
dans des canaux rectangulaires. Les conclusions précédentes nous ont conduit à ne conserver
que le schéma numérique upwind pour le développement des logiciels de modélisation des
déversoirs et de la confluence.
125
Partie 4 : Validation des outils numériques réalisés
Capteur
ultrasons
Venturi
40cm
Sens de l’écoulement 20cm
venturi
On mesure manuellement les hauteurs d’eau au sein du venturi à l’aide d’une règle,
pour différents débits.
On mesure la hauteur totale du venturi à intervalles réguliers sur son axe central à
partir du sommet de l’ouvrage qui sert de point de départ à la mesure. On prend ensuite la
126
Partie 4 : Validation des outils numériques réalisés
distance entre ce point de référence et la surface libre puis par différence on obtient la hauteur
d’eau dans le canal.
5.2 Résultats obtenus
On a comparé la hauteur d’eau mesurée et celle calculée pour trois débits différents. On
présente ci-dessous les résultats obtenus. Les valeurs de hauteurs d’eau associées à ces trois
courbes sont regroupées dans l’annexe 7.
0.09
h mesurée +3 mm
0.08 h calculée modèle 1D
h mesurée
0.07 h mesurée -3 mm
0.06
hauteur d'eau (m)
0.04
partie rectangulaire
0.03
0.02
0.01
0
1.5 2 2.5 3 3.5
position (m)
127
Partie 4 : Validation des outils numériques réalisés
0.12
h obernai
h obernai erreur+3mm
h obernai erreur-3mm
0.1
h calculée modèle 1D
0.08
hauteur d'eau (m)
fin du venturi
début du venturi
0.06
partie rectangulaire
0.04
0.02
0
1.5 2 2.5 3 3.5
position (m)
0.14
h obernai
h obernai erreur+3mm
0.12 h obernai erreur-3mm
h calculée modèle 1D
0.1
hauteur d'eau (m)
0.06
partie rectangulaire
0.04
0.02
0
1.5 2 2.5 3 3.5
position (m)
5.3 Interprétation
pour la ligne d’eau expérimentale. Un premier avant l’entrée dans le venturi, un autre dans la
partie rectangulaire et un dernier à l’aval. C’est la forme du venturi qui impose un écoulement
fluvial à l’amont et torrentiel à l’aval.
Les résultats du logiciel font bien apparaître ces paliers comme le montre la figure
suivante. On a globalement une ligne d’eau qui recoupe celle obtenue expérimentalement.
Cependant en comparant plus finement les résultats, il apparaît des écarts entre le logiciel et
l’observation.
En effet, le logiciel produit une ligne d’eau correcte en amont et en aval de la partie
rectangulaire du venturi. Ainsi, la courbe simulée se trouve dans la plage d’incertitude de la
mesure manuelle. Cependant pour la partie rectangulaire, des écarts plus grands sont observés.
Ceux-ci sont de l’ordre de 10 à 15 % d’erreur suivant les débits simulés. Le logiciel a
tendance à surévaluer les hauteur d’eau.
Une première explication de ces écarts pourrait être la différence de rugosité entre les
canaux amont et aval et le canal à Venturi. Mais ces différences sont a priori faibles et le canal
à Venturi est d’une longueur trop petite pour que des erreurs sur les coefficients de perte de
charge puisse engendrer des écarts aussi significatifs.
L’explication de ces écarts est plus certainement liée au fait que l’entonnement est un
ouvrage qui engendre des phénomènes de nature bidimensionnelle. En effet, le passage de
l’écoulement dans la zone contractée induit des ondes qui interfèrent et créent des courants
dans d’autres directions que celle de l’écoulement. Les lignes de courant ne sont plus
unidirectionnelles dans la partie rétrécie.
Apparitition
de vagues
Toutefois, les écarts moyens observés sont acceptables par rapport à la précision désirée. Pour
résoudre ce problème, on peut envisager d’ajouter une perte de charge au modèle pour
diminuer la hauteur d’eau dans le palier. En général, une perte de charge a tendance à
engendrer une augmentation de la ligne d’eau à l’amont et de la diminuer à l’aval.
En conclusion, la comparaison des surfaces libres rend compte de l’efficacité du logiciel
1D concernant le calcul de la surface libre en cas d’entonnement.
129
Partie 4 : Validation des outils numériques réalisés
Dans ce qui suit sont présentés le banc d’essai ainsi que les dispositifs de mesures
associés puis la comparaison des résultats expérimentaux et logiciels obtenus.
Dans ce qui suit, on présente les résultats de la phase de validation du modèle de Barré
de Saint-Venant concernant le partage des débits au niveau du déversoir et les valeurs de
hauteur d’eau dans l’ouvrage.
Celle-ci a été réalisée à l’aide d’un banc d’essai physique de déversoir représenté sur
les figures ci-dessous :(fig.49 et 50).
Réservoir
Venturi
Vanne papillon
11 m
DO
16 m
130
Partie 4 : Validation des outils numériques réalisés
Le banc d’essai est constitué d’une réserve d’eau enterrée dans laquelle est placée une
pompe immergée capable de débiter 250 m3/h. Celle-ci alimente un second bac d’une capacité
d’environ 1 m3 placé en hauteur et dans lequel le niveau de l’eau reste constant. On garantit
ainsi une alimentation du banc d’essai avec un débit constant. L’arrivée de l’eau dans le
collecteur circulaire amont ainsi que son retour vers le réservoir enterré après passage dans le
déversoir et les canaux circulaires conservé et déversé sont assurés par des canaux
rectangulaires d’une largeur de 40 cm. Le diamètre des canaux circulaires est de 20 cm. Leurs
pentes ainsi que celle du déversoir sont réglables grâce à des systèmes de supports montés sur
tiges filetées.
Les débits sont mesurés dans les branches rectangulaires amont, aval et déversée grâce
à l’association d’un capteur à ultrason et d’un Venturi.
En fait, chaque canal est équipé d'un Venturi couplé à un capteur à ultrasons situé à
environ 1 m en amont du Venturi. Les capteurs à ultrasons mesurent la hauteur d'eau dans le
canal puis le débit est calculé grâce à la loi du Venturi (de type Q = [Link]) fournie par le
constructeur. Les capteurs sont reliés à un ordinateur muni du logiciel HP Vee Lab qui permet
le traitement des données, le calcul des débits et le tracé des courbes d’évolution des débits en
fonction du temps.
Une étape supplémentaire fait la somme des débits conservés et déversés et la compare
au débit d'entrée en admettant une fourchette d'erreur de 5%. Ainsi, on dispose d’une mesure
pour le débit d’entrée, le débit conservé et le débit déversé ce qui permet de vérifier la
conservation des débits à tout moment.
131
Partie 4 : Validation des outils numériques réalisés
Transmetteurs
PC
Capteur ultrasons HP VEE LAB
Venturi
Dans un premier temps, les mesures de hauteurs d’eau sur la crête du déversoir ont été
réalisées manuellement grâce à un limnimètre. C’est cette méthode qui a été utilisée dans le
cas du déversoir à crête basse de hauteur de crête 7 cm.
132
Partie 4 : Validation des outils numériques réalisés
Puis le laboratoire s’est doté d’un instrument de mesure capable d’obtenir une image
en 3D de la nappe de surface dans le déversoir [HOLO3-2000] et ainsi par traitement
numérique d’obtenir la hauteur d’eau en tout point du déversoir et en particulier, pour ce qui
nous intéresse, sur une ligne parallèle à la crête joignant l’amont et l’aval du déversoir.
Ces mesures sont réalisées en lumière structurée et utilisent le principe de triangulation qui
consiste à projeter l’image d’un réseau de franges parallèles sur l’objet à mesurer, puis à
observer grâce à une caméra numérique cet objet depuis un point faisant un angle non nul
avec la direction de projection [Lipeme Kouyi-2001]. Un angle non nul impose une
déformation de l’image. En effet, la position du capteur entraîne une déformation des pixels
qui n’ont donc pas toujours la même taille. Cette déformation du pixel dans les plans (x,y) et
(y,z) a été quantifiée et corrigée numériquement afin de connaître la position d’un point en
unité de longueur (mm) depuis les bords du déversoir [Anstett-2001].
Le système de mesure est composé d’un projecteur stroboscopique muni d’un réseau
de franges d’égale épaisseur envoyant un flash de 100 µs toutes les secondes et d’une caméra
numérique de résolution 755x567 pixels, de temps d’acquisition de 30 µs. Le stroboscope et
la caméra numérique sont synchronisés.
Avant de procéder aux mesures, le dispositif est calibré par l’observation de plans de
référence horizontaux dont les hauteurs sont connues. Le calcul de la surface libre c’est à dire
des hauteurs d’eau est ensuite fait à partir de ces plans de références.
Le réseau de franges est déformé par le relief de l’objet et l’analyse de cette
déformation nous permet de remonter à la forme initiale de l’objet. La caméra transmet
l’image à un ordinateur qui restitue, après traitement numérique basé sur l’analyse de Fourier
[Chardon-2001], une mesure en trois dimensions de la surface libre sur la crête à partir de
l’image du réseau. Le logiciel utilisé est donné par le fournisseur (HOLO3) et est programmé
sous MATLAB.
Sur la figure ci-dessous, on présente le dispositif de mesure utilisé :
133
Partie 4 : Validation des outils numériques réalisés
Au final, après interpolation sur plusieurs observations afin d’obtenir une valeur
moyenne des hauteurs d’eau observées en tout point de la surface libre c’est à dire qui soit
représentative de la nappe liquide observée, on obtient des graphiques du type de celui
représenté ci-dessous :
134
Partie 4 : Validation des outils numériques réalisés
Il s’agit, du tracé 3D réalisé à partir des coordonnées des points obtenues grâce au
dispositif caméra en vue de confronter les hauteurs d’eau observées expérimentalement aux
valeurs obtenues par notre logiciel de simulation 1D. De ce graphique on extrait les trois
lignes d’eau qui sont mises en évidence en noir. Elles représentent dans le cas d’un déversoir
à une crête déversante les lignes d’eau successivement proches de la crête, au milieu du
déversoir et au fond du déversoir. Dans le cas du déversoir à crête double, la ligne d’eau au
fond est remplacée par la ligne d’eau près de la seconde crête déversante.
Elles permettent de visualiser la variabilité de la hauteur d’eau dans le déversoir selon la
direction transversale.
Ces trois lignes sont transposées, en plus de celle obtenue par le logiciel de simulation,
sur des graphiques 2D représentant la hauteur d’eau en fonction de la position dans le
déversoir. Toutes les configurations ont été étudiées par cette méthode sauf celle du déversoir
à crête basse simple de hauteur de crête 7 cm car nous ne disposions pas encore du dispositif
de mesure 3D des hauteurs d’eau.
6.3 Déversoir à crête basse simple prismatique
Pente amont Pente aval Pente amont Pente aval Pente amont Pente aval
(‰) (‰) (‰) (‰) (‰) (‰)
1 1 1
2 3 3
4 5 5 9 5
1
6 7 7
8 9 9
10
D’autre part, dans chaque configuration de pentes, nous avons fait varier le débit entre
deux valeurs extrêmes qui correspondent pour la première à une valeur légèrement supérieure
au débit critique (Débit pour lequel le déversoir commence à déverser) et pour la seconde la
valeur pour laquelle la conduite possédant la pente la plus faible passe en charge. En effet,
lorsque nous avons réalisé cette campagne de mesure, le logiciel n’était pas encore capable de
simuler un écoulement en charge. Ce débit maximum se situe pour tous les cas autour d’une
valeur de 100 m3/h.
135
Partie 4 : Validation des outils numériques réalisés
De plus, les coefficients de Strickler ont été évalué sur le banc d’essai et varient selon
la pente et les débits injectés. Les mesures effectuées ont mis en évidence une plage pour le
coefficient de Strickler comprise entre 95 et 115. Ainsi, nous avons utilisé une valeur de 95
dans le cas d’une pente supérieure à 5 ‰ et de 115 dans le cas d’une pente inférieure ou égale
à 5 ‰.
D’autre part, la vanne papillon située à l’aval du canal conservé perturbe l’écoulement
dans le cas d’un régime fluvial. Cette influence aval peut remonter dans le déversoir et influer
sur le débit déversé. C’est pourquoi, au niveau de la simulation, nous avons appliqué en
condition à la limite aval une loi de vanne établie expérimentalement. Celle-ci permet
d’associer la hauteur d’eau mesurée à 1 m à l’amont de la vanne au débit transitant dans le
canal aval. Les valeurs mesurées nous ont permis d’établir la relation suivante :
Nous avons sur le site expérimental mesuré les débits injecté (Q I ) et déversé (Q DE )
pour une certaine configuration de pentes. Ensuite, nous avons, au laboratoire, effectué une
simulation en considérant le même débit d’entré et les mêmes pentes pour obtenir le débit
déversé simulé (Q DS ) . Pour finir nous avons calculé pour chaque cas les pourcentages de
débit déversé simulé (PS ) et déversé expérimental (PE ) par les relations :
Expérimentation :
Q
PE = 1 − DE
QI
Simulation :
Q
PS = 1 − DS
QI
Pour finir, nous avons déterminé l’erreur commise sur le pourcentage de déversement
par :
E = PE − PS
Les tableaux ci-dessous regroupent les valeurs issues des mesures ainsi que des calculs
évoqués ci-dessus :
136
Partie 4 : Validation des outils numériques réalisés
Pente amont de 1 ‰ (I am = 1 ‰ ) : Pente amont de 3 ‰ (I am = 3 ‰ ) :
Iav QI QDE QDS PE PS E Emoy
3 3 3
(‰) (m /h) (m /h) (m /h) (%) (%) (%) (%)
Iav QI QDE QDS PE PS E Emoy 30.1 6.0 7.1 20.1 23.6 3.5
3 3 3
(‰) (m /h) (m /h) (m /h) (%) (%) (%) (%) 1 50.3 15.0 16.1 29.8 32.0 2.2 3.2
24.3 5.1 6.3 21.0 25.9 4.9 89.4 35.7 32.2 40.0 36.0 3.9
37.6 10.0 11.0 26.7 29.1 2.4 39.1 6.9 9.9 17.6 25.3 7.7
1 3.4
59.2 20.2 20.3 34.1 34.4 0.3 2 52.3 13.0 15.2 24.8 29.1 4.2 5.3
79.6 31.3 29.2 39.3 36.7 2.6 91.0 33.8 30.2 37.1 33.2 3.9
38.7 5.2 8.1 13.3 21.0 7.7 40.2 4.8 6.4 11.9 15.8 4.0
3 62.9 16.2 17.2 25.8 27.3 1.6 4.1 4 52.0 10.3 7.4 19.8 14.1 5.6 5.1
80.2 26.5 24.0 33.0 29.9 3.1 84.3 28.6 23.8 33.9 28.2 5.6
48.4 8.5 8.1 17.7 16.6 1.0 53.3 10.7 9.4 20.1 17.6 2.5
5 60.1 14.2 14.0 23.7 23.3 0.3 3.0 6 76.8 24.0 19.9 31.3 25.9 5.3 4.2
79.2 25.8 19.8 32.6 25.0 7.6 97.5 36.1 31.3 37.0 32.1 4.9
52.6 10.5 8.6 20.0 16.4 3.5 48.2 8.3 6.8 17.2 14.1 3.1
7 62.7 15.6 12.6 24.9 20.2 4.8 5.2 8 65.6 17.5 13.4 26.6 20.4 6.2 4.9
77.8 24.6 18.8 31.6 24.2 7.4 81.9 26.9 22.5 32.8 27.5 5.4
53.8 10.6 8.8 19.8 16.4 3.4 40.1 4.3 4.1 10.7 10.2 0.4
10 71.2 20.4 15.4 28.7 21.6 7.1 5.8 10 65.0 16.3 13.9 25.1 21.4 3.7 3.0
93.7 34.3 27.8 36.6 29.7 6.9 92.3 32.8 28.2 35.5 30.6 5.0
Pente amont de 9‰ (I am = 9 ‰ ) :
Pour le calcul de l’erreur, nous avons utilisé une valeur absolue. En effet, selon le cas
envisagé la valeur expérimentale peut être supérieure ou inférieure à la valeur simulée mais
nous n’avons pas pu établir de règle stricte concernant la tendance du modèle à sur évaluer ou
sous évaluer les débits déversés selon que le débit injecté est plus ou moins important. Il
semble néanmoins que dans la majorité des cas traités le modèle ait tendance à sous évaluer le
débit déversé.
On constate que les erreurs sur les pourcentages de déversement n’excèdent jamais 8.5
%. De plus, il faut tenir compte du fait que les valeurs de débits obtenues expérimentalement
sont fiable à ± 5 % . Cette erreur est due à la précision des sondes à ultrasons utilisées en
association avec les Venturi pour évaluer les débits.
En conclusion, on peut considérer que les valeurs obtenues par le modèle sont en bon
accord avec les valeurs attendues et que les erreurs qui affectent les valeurs sont tout à fait
admissibles dans le cadre d’une application au domaine de l’assainissement.
Voici pour terminer les lignes d’eau obtenues dans plusieurs configurations qui
mettent en évidence quelques types de courbes de remous que l’on peut voir apparaître à
l’intérieur d’un déversoir prismatique. Ces trois cas ont été choisis parmi les configurations
envisagées précédemment. Les simulations sont réalisées en régime transitoire. Ainsi, on peut
observer l’évolution de la ligne d’eau au cours du temps jusqu’à la stabilisation qui
correspond au régime permanent. Nous ne présentons pas ici l’évolution de la ligne d’eau
mais uniquement les lignes d’eau obtenues après stabilisation.
Pour se rapprocher le plus possible de la réalité du banc d’essai physique et pour éviter
d’avoir à appliquer des conditions aux limites au niveau des extrémités du déversoir nous
avons modélisé le trio, collecteur amont d’une longueur de 6 m, déversoir d’une longueur de
1.5 m et collecteur aval de 5 m. La hauteur de crête considérée est de 7.15 cm. Ainsi les
conditions aux limites sont appliquées aux extrémités des collecteurs. A l’amont, la condition
a la limite appliquée est la hauteur normale et à l’aval la loi de vanne décrite dans le
paragraphe précédent.
On a sélectionné un cas fluvial, un cas torrentiel et un cas faisant apparaître un ressaut
hydraulique dans le déversoir. Les valeurs de hauteurs d’eau mesurée et simulées sont
regroupées dans l’annexe 9 pour les trois cas envisagés.
Nous avons ici considéré une configuration dans laquelle la pente du collecteur amont
est de 3 ‰ et des pentes de 1 ‰ pour le déversoir et le canal aval. Le débit injecté en continu
à l’amont du système est de 24,3 m3/h.
Les figures ci-dessous montrent la forme de la ligne d’eau dans le déversoir et la
courbe d’évolution du débit dans le système déversoir plus collecteur. Les valeurs de hauteurs
d’eau expérimentales ont été obtenues en six points du déversoir (tous les 30 centimètres)
grâce à un limnimètre.
138
Partie 4 : Validation des outils numériques réalisés
Water profile in the sewer side weir Flow rate variation curve
9 25
Qupstream=24.3 m3/h
Simulated water profile 24
8.5 Measured water depth
23
22
8 Simulated flow rate
21 Experimental flow rate
7.5
20
Qdownstream=19.2 m3/h
19
7
Sewer side weir crest 18
Qdownstream=17.6 m3/h
6.5 17
6 6.5 7 7.5 5.5 6 6.5 7 7.5 8
Dans le cas d’un écoulement fluvial la simulation prévoit une remontée de la ligne
d’eau dans le déversoir en allant de l’amont vers l’aval ce qui est conforme aux observations
ainsi qu’à la théorie.
L’erreur commise sur les valeurs de hauteur d’eau est ici en moyenne de 4 mm. Ceci
est acceptable si on tient compte de la variabilité de la ligne d’eau dans le déversoir et donc de
la difficulté d’obtenir une valeur fiable pour le tirant d’eau. On estime que l’erreur de mesure
sur la hauteur d’eau est de ±5 mm.
Nous avons cette fois considéré une configuration dans laquelle la pente du collecteur
amont est de 3 ‰, celle du déversoir de 1 ‰ et celle du collecteur aval de 7 ‰. Le débit
injecté en continu à l’amont du système est de 62,7 m3/h.
Water profile in the sewer side weir Flow rate evolution curve
12 64
Qupstream=62.7 m3/h
62
11
60
Simulated water profile
Measured water depth 58 Simulated flow rate
10
56 Experimental flow rate
9 54
52
8
50
48
7
Sewer side weir crest Qdownstream=45.3 m3/h
46
Qdownstream=45.0 m3/h
6 44
6 6.5 7 7.5 5.5 6 6.5 7 7.5 8
139
Partie 4 : Validation des outils numériques réalisés
Dans le dernier cas présenté ici, nous avons considéré une configuration dans laquelle
la pente du collecteur amont est de 5 ‰ et des pentes de 1 ‰ pour le déversoir et le canal
aval. Le débit injecté en continu à l’amont du système est de 47 m3/h.
Water profile in the sewer side weir Flow rate evolution curve
12 48
Qupstream=47.0 m3/h
Simulated water profile 46
11
Measured water depth 44
10 42
8 36 Qdownstream=34.9 m3/h
34
7
Sewer side weir crest 32
Qdownstream=31.0 m3/h
6 30
6 6.5 7 7.5 5.5 6 6.5 7 7.5 8
140
Partie 4 : Validation des outils numériques réalisés
Pour l’étude de cette configuration et de celles qui suivent, le dispositif expérimental a été
amélioré. En effet, les crêtes ont été moulées par thermoformage et non plus obtenues de
manière artisanale. La méthode mise en œuvre pour la mesure des hauteurs d’eau a été
modifiée. C’est celle décrite dans le paragraphe 6.2 concernant le dispositif caméra. Les
résultats obtenus sont présentés dans ce qui suit et les valeurs de hauteurs d’eau simulées et
expérimentales sont regroupées dans l’annexe 10 :
On considère dans ce qui suit le cas d’un déversoir à crête basse dont les caractéristiques
sont les suivantes :
• 1 crête déversante
• Hauteur de crête w=5 cm
• longueur=1.50 m.
Les configurations de pentes envisagées sont présentées dans le tableau ci-dessous. Quant
à la pente du déversoir, elle a été prise comme constante et égale à 0.5 ‰.
Pente amont de 0.5 ‰ ; Pente aval de 0.5 ‰ : Pente amont de 0.5 ‰ ; Pente aval de 8 ‰ :
141
Partie 4 : Validation des outils numériques réalisés
QI QDE Q DS PE PS E Emoy
(m /h) (m /h) (m 3/h)
3 3
(%) (%) (%) (%)
20.5 1.0 1.8 4.9 8.8 3.9
36.6 5.5 8.6 15.0 23.5 8.5
54.3 15.2 18.4 28.0 33.9 5.9
65.7 22.5 25.2 34.2 38.4 4.1
3.4
78.9 30.5 32.8 38.7 41.6 2.9
91.2 39.9 40.0 43.8 43.9 0.1
108.5 50.3 50.0 46.4 46.1 0.3
121.0 58.7 56.5 48.5 46.7 1.8
142
Partie 4 : Validation des outils numériques réalisés
12.0
11.0
10.0
Hauteur d'eau (cm)
9.0 h_sim
h_mes crête
h_mes milieu
8.0 h_mes fond
7.0
6.0
5.0
0 20 40 60 80 100 120 140 160
Position (cm)
• Les résultats suivants ci-dessous concernent le même déversoir pour une pente amont de
Iam=0.5‰, une pente aval de Iav=0.5‰. Le débit considéré était de Q=65.2m3h-1
11.0
10.0
9.0
Hauteur d'eau (cm)
h_sim
h_exp crête
8.0
h_exp milieu
h_exp fond
7.0
6.0
5.0
0 20 40 60 80 100 120 140 160
Position (cm)
143
Partie 4 : Validation des outils numériques réalisés
On considère dans ce qui suit le cas d’un déversoir à crête basse dont les caractéristiques
sont les suivantes :
• 1 crête déversante
• Hauteur de crête w=12.5 cm
• longueur=1.50 m.
Les configurations de pentes envisagées sont présentées dans le tableau ci-dessous. Quant
à la pente du déversoir, elle a été prise comme constante et égale à 0.5 ‰.
Pente amont de 0.5 ‰ ; Pente aval de 0.5 ‰ : Pente amont de 0.5 ‰ ; Pente aval de 8 ‰ :
144
Partie 4 : Validation des outils numériques réalisés
• Les résultats ci-dessous concernent le déversoir précédemment décrit pour une pente
amont de Iam=0.5‰, une pente aval de Iav=0.5‰. Le débit considéré était de Q=73.3m3h-
1
:
16.0
15.0
hauteur d'eau (cm)
14.0
h_sim
h_mes crête
h_mes milieu
h_mes fond
13.0
12.0
11.0
0 20 40 60 80 100 120 140 160
Position (cm)
La ligne d’eau monte lorsque l’on progresse sur la crête du déversoir. Ceci est
caractéristique d’un écoulement de type fluvial dans l’ouvrage.
•
Les résultats suivants ci-dessous concernent le même déversoir pour une pente amont de
Iam=0.5‰, une pente aval de Iav=8‰. Le débit considéré était de Q=117.7m3h-1
145
Partie 4 : Validation des outils numériques réalisés
14.0
13.0
Hauteur d'eau (cm)
h_mes crête
h_mes milieu
h_mes fond
h_sim
12.0
11.0
0 20 40 60 80 100 120 140 160
Position (cm)
La ligne d’eau est cette fois quasi horizontale avec une tendance à la baisse en se
dirigeant vers l’aval. On a dans ce cas une ligne d’eau de type torrentielle.
La grande variabilité des courbes expérimentale n’est qu’apparente puisque
l’amplitude de l’axe des ordonnées correspondant aux hauteurs d’eau n’est que de trois
centimètres. D’autre part, on peut remarquer que si le graphique était réalisé dans un repère
orthonormé l’amplitude de ces variations prendrait une toute autre échelle.
Les résultats suivants concernent des déversoirs à doubles crêtes. Les hauteurs de crêtes
étudiées sont de w=5, 7.5 et 12.5 cm.
• 2 crêtes déversante
146
Partie 4 : Validation des outils numériques réalisés
Les configurations de pentes envisagées sont présentées dans le tableau ci-dessous. Quant
à la pente du déversoir, elle a été prise comme constante et égale à 0.5 ‰.
Pente amont de 0.5 ‰ ; Pente aval de 0.5 ‰ : Pente amont de 0.5 ‰ ; Pente aval de 8 ‰ :
La comparaison des débits simulés et mesurés révèle encore une fois une excellente
corrélation. Ceci montre l’aptitude de notre modèle à évaluer les débits déversés dans le cas
du déversoir à crête double.
•
Les résultats ci-dessous concernent le déversoir précédemment décrit pour une pente
amont de Iam=0.5‰, une pente aval de Iav=0.5‰. Le débit considéré était de Q=92.6m3h-1
147
Partie 4 : Validation des outils numériques réalisés
11.0
10.0
9.0
Hauteur d'eau (m)
8.0 h_sim
h_mes crête
h_mes milieu
7.0 h_mes fond
6.0
5.0
4.0
0 20 40 60 80 100 120 140 160
Position (cm)
10.0
9.0
8.0
Hauteur d'eau (cm)
h_sim
h_mes crête
7.0
h_mes milieu
h_mes fond
6.0
5.0
4.0
0 20 40 60 80 100 120 140 160
Position (cm)
148
Partie 4 : Validation des outils numériques réalisés
Les deux lignes présentées dans ce paragraphe sont de type torrentiel et comme dans
les cas des déversoirs précédents, les hauteurs d’eau obtenues numériquement surévaluent les
valeurs expérimentales.
• 2 crêtes déversante
• Hauteur de crête w=7.5 cm
• Longueur de crête L=1.5 m.
Les configurations de pentes envisagées sont présentées dans le tableau ci-dessous. Quant
à la pente du déversoir, elle a été prise comme constante et égale à 0.5 ‰.
Les tableaux suivants donnent les valeurs des débits déversés pour différents débits injectés :
149
Partie 4 : Validation des outils numériques réalisés
• Les résultats ci-dessous concernent le déversoir précédemment décrit pour une pente
amont de Iam=0.5‰, une pente aval de Iav=8‰. Le débit considéré était de Q=84m3h-1.
11.0
10.0
Hauteur d'eau (cm)
h_sim
h_mes crête
9.0
h_mes milieu
h_emes fond
8.0
7.0
0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0
Position (cm)
Dans ce cas, la ligne d’eau est de type torrentiel. Bien que les résultats numériques
surévaluent la ligne d’eau, on observe la même tendance.
• Les résultats suivants ci-dessous concernent le même déversoir pour une pente amont de
Iam=0.5‰, une pente aval de Iav=0.5‰. Le débit considéré était de Q=89 m3h-1.
150
Partie 4 : Validation des outils numériques réalisés
12.0
11.0
10.0
Hauteur d'eau (cm)
9.0 h_sim
h_mes crête
h_mes milieu
8.0 h_mes fond
7.0
6.0
5.0
0 20 40 60 80 100 120 140 160
Position (cm)
• 2 crêtes déversante
• Hauteur de crête w=12.5 cm
• Longueur de crête L=1.5 m.
151
Partie 4 : Validation des outils numériques réalisés
Les tableaux suivants donnent les valeurs des débits déversés pour différents débits
injectés :
• Les résultats ci-dessous concernent le déversoir précédemment décrit pour une pente
amont de Iam=0.5‰, une pente aval de Iav=0.5‰. Le débit considéré était de Q=100m3h-1..
152
Partie 4 : Validation des outils numériques réalisés
18.0
17.0
16.0
Hauteur d'eau (cm)
15.0 h_sim
h_mes crête
h_mes milieu
14.0 h_mes fond
13.0
12.0
11.0
0 20 40 60 80 100 120 140 160
Position (cm)
• Les résultats suivants ci-dessous concernent le même déversoir pour une pente amont de
Iam=0.5‰, une pente aval de Iav=8‰. Le débit considéré était de Q=108 m3h-1.
153
Partie 4 : Validation des outils numériques réalisés
14.0
13.0
Hauteur d'eau (cm)
h_sim
h_mes crête
h_mes milieu
h_mes fond
12.0
11.0
0 20 40 60 80 100 120 140 160
Position (cm)
La ligne d’eau simulée est torrentielle. Les valeurs expérimentales montrent dans la
partie amont diminution de la hauteur d’eau puis une élévation en allant vers l’aval ce qui
semblerait indiquer l’apparition d’un ressaut. Or avec une pente de 8 ‰ dans la conduite aval,
ceci semble improbable.
Dans le cas du déversoir à crête haute, la hauteur de crête est supérieure au diamètre de la
conduite aval. L’écoulement se produit donc en charge dans la conduite aval. Ainsi, c’est
l’étude de cette configuration qui nous a permis de valider la méthode de modélisation que
nous avons employée pour représenter un écoulement sous pression.
Le déversoir à double crête haute que nous avons étudié avait une longueur de 1.5 m et
une hauteur de crête de 94 mm. Les caractéristiques du système de collecteurs sont les
suivantes :
• Diamètre amont de 188 mm et diamètre aval de 67 mm.
• Pente du canal amont de 6 ‰ et pente du canal aval de 1.33 ‰.
154
Partie 4 : Validation des outils numériques réalisés
QI QDE Q DS PE PS E Emoy
3 3 3 (%)
(m /h) (m /h) (m /h) (%) (%) (%)
11.3 3.3 2.2 29.2 19.5 9.7
28.1 19.4 17.3 69.0 61.6 7.5
42.9 33.2 31.3 77.4 73.0 4.4
54.6 44.8 41.8 82.1 76.6 5.5
4.4
74.4 61.4 61.2 82.5 82.3 0.3
90.1 75.7 75.2 84.0 83.5 0.6
103.5 88.9 83.2 85.9 80.4 5.5
122.1 108.3 106 88.7 86.8 1.9
Les résultats montrent une bonne corrélation entre valeurs expérimentales et valeurs
simulées. Ceci montre l’aptitude de notre modèle à rendre compte de manière correcte du
comportement du déversoir à crête haute. Ceci valide la méthode utilisée dans le logiciel afin
de prendre en compte la mise en charge.
• Les résultats ci-dessous concernent le déversoir précédemment décrit pour une pente
amont de Iam=6‰, une pente aval de Iav=1.33‰. Le débit considéré était de Q=103m3h-1
155
Partie 4 : Validation des outils numériques réalisés
13
12
Hauteur d'eau (cm)
11
h_sim
h_mes crête
h_mes milieu
h_mes fond
10
8
0 20 40 60 80 100 120 140 160
Position (cm)
On observe dans un premier temps une forte diminution de la hauteur d’eau puis une
stabilisation pour ensuite voir apparaître un ressaut hydraulique à l’aval en raison de la
contraction de la conduite aval. On observe un décalage du même ordre de grandeur que dans
les cas à crête basse mais cette fois, la simulation sous-évalue la hauteur d’eau alors que
c’était l’inverse dans les cas à crête basse.
D’autres tests réalisés nous ont permis de vérifier que la différence entre les débits
déversés par un déversoir à crête haute simple et double pour une même hauteur de crête est
très faible. Seule les hauteurs d’eau sont plus élevées dans le cas à crête simple que dans le
cas à crête double.
156
Partie 4 : Validation des outils numériques réalisés
∂A ∂Q
+ =0 (4.213.)
∂t ∂x
∆t
U in +1 = U in + λ Fin+ 1 − Fin− 1 + ∆tG in avec λ = (4.214.)
2
2 ∆x ( j)
Si on somme ensuite sur toutes les cellules de calculs d’un tronçon j c’est à dire de 1 jusqu’à
nbdx(j), on obtient :
nbdx ( j ) nbdx ( j )
∆x ( j) ∑ [UNEW(1, i, j) − U(1, i, j)] = ∆t ∑ [FP(1, i, j) − FN(1, i, j)]
i =1 i =1
(4.216.)
157
Partie 4 : Validation des outils numériques réalisés
D’où :
nbdx ( j )
et finalement, on obtient :
nbdx ( j )
∆x ( j) ∑ [UNEW(1, i, j) − U(1, i, j)] = ∆t[FP(1, nbdx ( j), j) − FN(1,1, j)]
i =1
(4.218.)
Au niveau expérimental, nous n’avions aucun dispositif qui aurait pu nous permettre de
comparer les valeurs obtenues grâce au logiciel avec des valeurs mesurées. Nous avons
néanmoins comparé les résultats obtenus grâce à notre modèle avec ceux produits par un
logiciel de simulation 3D qui est FLUENT.
Il s’avère que dans les cas que nous avons testé, les lignes d’eau obtenues grâce à notre
logiciel ne correspondent pas du tout à celles prévues par FLUENT et ceci quelle que soit la
méthode de prise en compte de la confluence envisagée (Projection des flux numériques sur la
direction de la conduite aval ou bilan de quantité de mouvement permettant de calculer une
hauteur d’eau dans la confluence).
C’est pourquoi le laboratoire a lancé un sujet de thèse concernant cette problématique.
158
Partie 4 : Validation des outils numériques réalisés
8 CONCLUSION
Les résultats obtenus grâce à nos modèles, hormis celui concernant la confluence,
donnent des résultats satisfaisant. La confluence nous a tout de même permis de vérifier
le caractère conservatif du schéma numérique utilisé.
Les meilleurs résultats obtenus concernent les débits déversés. Les erreurs constatées
entre les données expérimentales et les résultats de simulations n’excèdent en général pas
10 % ce qui est suffisant dans le cadre d’une application au domaine de l’assainissement.
Les cas de grosses erreurs adviennent lorsque le débit injecté est proche du débit de début
de déversement de l’ouvrage. Dans ce cas, on peut être proche de la limite de sensibilité
des capteurs ultrasons.
Les erreurs constatées sur les hauteurs d’eau peuvent être imputées d’une part aux
variations transversales de la ligne d’eau et d’autre part à l’exclusion des phénomènes de
turbulence dans les équations de base du modèle. Cependant, la comparaison des courbes
de remous obtenues expérimentalement et numériquement met en évidence les mêmes
sens de variation et confirme l’aptitude du modèle à localiser le ressaut hydraulique. Ceci
montre que nos modèles sont aptes à représenter les phénomènes d’écoulement en
collecteurs avec entonnement et dans les déversoirs, même dans le cas du déversoir à
crête haute ou à double crête.
D’autres campagnes de mesures pourront encore être menées dans l’avenir concernant le
déversoir. En effet, on dispose de nombreuses crêtes moulées qui n’ont pas encore été testées.
On citera les crêtes permettant d’obtenir des déversoirs à crête basse avec entonnement ou
encore des déversoirs à crête simple ou double, haute ou basse pour d’autres hauteurs,
diamètres et longueurs de crêtes.
159
Partie 5 : Conclusion générale
PARTIE 5
CONCLUSION
GENERALE
160
Partie 5 : Conclusion générale
L’étude bibliographique que nous avons mené a permis de mettre en évidence que les
modèles utilisés classiquement n’étaient pas aptes à reproduire de manière fine le
comportement hydraulique d’un réseau d’assainissement. En effet, des discontinuités telles
que le ressaut hydraulique peuvent apparaître. Pour les prendre en compte et tenir compte de
la variabilité des hauteurs d’eau dans le système, il faut mettre en œuvre le système de Saint-
Venant sous sa forme conservative.
De plus, l’étude bibliographique a également montré que les méthodes numériques
classiques mises en œuvre pour la résolution d’un tel système d’équations tombent en défaut
dans les situations qui font apparaître de fortes discontinuités sur les débits et les hauteurs
d’eau ainsi que dans le cas d’un nombre de Froude proche de 1. Pour une résolution efficace
de ces équations des méthodes numériques à capture de chocs doivent être utilisées.
En ce qui concerne les déversoirs d’orages, nous avons mis en évidence l’intérêt de
considérer cet ouvrage comme faisant partie intégrante du système de collecteurs. De plus, les
débits déversés sont évalués par l’intermédiaire d’une relation dite de Hager qui permet de
tenir compte d’un éventuel entonnement et du changement de direction de la vitesse causé par
le déversement latéral.
Ensuite, nous avons mis en œuvre une méthode de prise en compte du phénomène de
mise en charge inspirée de la méthode de la fente de Preismann et adaptée à l’utilisation des
schémas numériques à capture de choc.
Pour finir, l’étude bibliographique nous a conduit à intégrer à notre modèle de Barré
de Saint-Venant le changement de section en introduisant un terme de pression latérale.
Les logiciels que nous avons construits ainsi que les méthodes numériques sur
lesquelles ils sont basés ont de multiples intérêts. En effet, les conditions aux limites
applicables aux extrémités du domaine de calcul sont basées sur l’hydraulique et les
paramètres de calages sont quasiment inexistants. Pour les méthodes numériques, elles
possèdent toutes le caractère TVD qui les rend non oscillatoires (Sauf pour la méthode Mac-
Cormack) et sont toutes au moins du second ordre de précision. Le caractère TVD est obtenu
soit par la méthode de limitation de pente ( méthode MUSCL) soit par la méthode à limitation
de flux (méthodes Upwind et Symmetric).
161
Partie 5 : Conclusion générale
permettant de décrire le comportement d’un tel système lorsque le régime d’écoulement est
fluvial ou torrentiel dans les trois branches mais le cas d’écoulements mixtes est peu envisagé.
Les études visant à globaliser ce type de situations même en simplifiant la géométrie de la
confluence sont rares et consistent en général à considérer les hauteurs d’eau dans les trois
branches égales au niveau de la confluence.
Dans le cas particulier de la confluence, l’écoulement pourrait être pris en compte comme
un écoulement de type 2D. Ainsi, un maillage 2D et une méthode numérique à capture de
choc pourraient être appliqués à cette structure. Ce travail fait l’objet de la thèse de M. Maher
Abdallah qui vient de débuter au laboratoire SHU.
Plus généralement, dans les ouvrages comme les déversoirs qui possèdent plusieurs
entrées et/ou sorties, les phénomènes de turbulence sont à prendre en compte d’une manière
plus fine que par la simple utilisation d’un modèle simplifié comme celui de Manning-
Strickler. Pour se faire, le passage par une étude individuelle de chaque ouvrage pour obtenir
une loi de fonctionnement et l’intégrer ensuite dans un modèle unidimensionnel ou une
modélisation de type 3D paraît indispensable tant la géométrie et les phénomènes de
turbulences sont importants en vue de l’évaluation des hauteurs d’eau et des capacités de
transport des structures.
D’autres ouvrages fréquemment rencontrés dans les réseaux d’assainissements comme
les bassins d’orages notamment à travers les quantités de matières stockées et les volumes
d’eau qui y transitent restent à étudier. En arrière plan, se situe le problème plus général du
transport solide et l’étape suivante de ce travail pourrait consister à coupler une équation de
transport au système de Barré de Saint-Venant afin d’évaluer les quantités de matières
déversées vers le milieu naturel par l’intermédiaire du déversoir d’orage.
162
Bibliographie
BIBLIOGRAPHIE
M. Anstett
Validation d’un capteur de mesure 3D de surface libre en lumière structurée.
Mémoire de stage DUT Mesure Physique
Juillet/Août 2001
K. Babaeyan-Koopaei
Choke-free flow in ovoïdal sewers with increase in bed elevations
Journal of Hydraulic Engineering, Vol. 127, N°2, February, 2001, pp.149-153.
L. Benayada
Construction d’une solution analytique et critique d’algorithmes numériques afférents
aux équations de Saint Venant en vue d’application aux écoulements transitoires sur fond
mobile.
Institut National Polytechnique de Toulouse, 1994, 167 p., Thèse de doctorat, Spécialité
mécanique des fluides.
M. Buyer
Etude bibliographique et modélisation des écoulements à surface libre en canaux et
collecteurs.
Mémoire de DEA “Mécanique et Ingénierie” option “Sciences de l’Eau”, 1999.
M. Carleton
Contribution à l’analyse et à la modélisation du fonctionnement des déversoirs d’orages.
Thèse docteur ingénieur, INSA, Lyon, Septembre 1985.
L. Chardon
Mise au point d’un algorithme de mesure de surface d’eau en lumière structurée.
Maîtrise de Mathématiques Option Ingénierie, 2001
B. Chocat (coordinateur)
Encyclopédie de l’hydrologie urbaine et de l’assainissement.
Paris, Ed. Tec&Doc Lavoisier, 1997
H. Combes
Modélisation numérique 3D des écoulements dans les déversoirs d’orage
Mémoire de fin d’étude ENGEES, 2000.
163
Bibliographie
A. El Khashab
Hydraulics of flow over side weirs.
Thèse de doctorat-Department of civil engineering of the University of Southampton
Juin 1975
P. Giersch
Les déversoirs d’orages : principe, données constructives et calcul
Direction Départementale de l’Agriculture et de la Forêt (DDAF) du Bas-Rhin, Service du
Génie Rural des Eaux et Forêts, Bureau d’études, 1985.
Y. Kovacs
Modèles de simulation d’écoulement transitoire en réseau d’assainissement.
CERGRENE, Ecole Nationale des Ponts et Chaussées, 1988, 328 p., Thèse de doctorat,
Spécialité : Sciences et Techniques de l’Environnement.
D. Euvrard
Résolution numérique des équations aux dérivées partielles de la physique, de la
mécanique et des sciences de l’ingénieur.
Différences finies, éléments finis, méthode des singularités.
Ed. Masson 1990, ISBN: 2-225-82128-3.
164
Bibliographie
A. Gardel
Perte de charge dans un étranglement conique. (Head loss in a conical contraction).
Bulletin technique de la suisse romande 88(21),1962, pp.313-320 ; 88(22), 1962, pp. 325-337.
S.K Godunov
A difference method for the numerical computation of continuous solutions of
hydrodynamic equations.
Mat. Sbornik, 47, 271-306 (1959) (Translated as JPRS by US DEPT of Commerce, 1960).
W.H. Hager
Die Hydraulik von Vereinigungsbauwerken.
Gas-Wasser-Abwasser, 62. Jahrgang 1982, N°7, pp 148-149
W.H. Hager
L’écoulement dans les déversoirs latéraux (Flow in side weirs).
Canadian Journal of Civil Engineering, 1986, 13(5):501-509.
W.H. Hager
Lateral outflow over side weirs.
Journal of Hydraulic Engineering, Vol. 113, N°4, April 1987, pp.491-504.
W.H. Hager
Transitional flow in channel junctions.
Journal of Hydraulic Engineering, Vol. 115, 1989a, N°2, pp 243-259
W.H. Hager
Supercritical flow in channel junctions.
Journal of Hydraulic Engineering, Vol. 115, 1989b, N°5, pp 595-616
W. H. Hager
Wastewater Hydraulics-Theory and Practice.
Ed. Springer, 1999
L. Hanich
Résolution des équations de la mécanique des fluides par des méthodes TVD en
coordonnées généralisées.
Université de Caen, 1996, 171 p., Thèse de doctorat, Spécialité mécanique.
165
Bibliographie
A. Harten
High-resolution schemes for hyperbolic conservation laws.
Journal of Computational Physics 49, 357-393, (1983).
A. Harten
On a class of high-resolution total variation stable finite difference schemes.
SIAM Journal Numerical Analysis, Vol. 21, 1-23 (1984).
E. Herouin
Modélisation des écoulements complexes à surface libre en milieu naturel (étude
bibliographique)
Rapport de DEA de Mécanique des fluides. Université Claude Bernard (Lyon 1), 1991
Hjelmfelt, A.T.
Flow in elliptical channels.
Water Power 19 (10), 1967, pp. 429-431.
HOLO 3
Fonctions Matlab de mesure de hauteur d’eau par méthode le lumière structurée pulsée.
Saint Louis, 2000, 10 p.
Y. Kovacs
Modèles de simulation d’écoulement transitoire en réseau d’assainissement.
CERGRENE, Ecole Nationale des Ponts et Chaussées, 1988, 328 p., Thèse de doctorat,
Sciences et Techniques de l’Environnement.
G. Lipeme Kouyi
Modélisation numérique 3D des écoulements dans les déversoirs d’orage.
Mémoire de DEA, Ecole Nationale du Génie de l’Eau et de l’Environnement (ENGEES), 2001,
80 p.
M. Louaked et L. Hanich
TVD scheme for the shallow water equations.
Journal of Hydraulic Research, Vol. 36, 1998, N°3, pp. 363-378.
166
Bibliographie
[Link]
Etudes du CEMAGREF- Equipements pour l’eau et l’environnement n°14. Modélisation,
analyse et commande optimale LQR d’un canal d’irrigation.
Edition CEMAGREF-DICOVA et CEMAGREF Montpellier, 1994, 220 p.
G. Martinet
Contribution à la modélisation numérique des avalanches de neige dense et des laves
torrentielles.
Université Joseph Fourier, Grenoble I, 1992, 218 p., Thèse de doctorat, Spécialité mécanique.
L.A. Monthe
Etude des équations aux dérivées partielles hyperboliques. Application aux équations de
Saint-Venant.
Institut National des Sciences Appliquées de Rouen, 1997, 181 p., Thèse de doctorat,
Discipline : Mathématiques Appliquées, Spécialité : Analyse numérique.
G. Moretti
The λ-scheme
Computers and fluids,vol. 7, 191-205, 1979
H. Mottiee
Un modèle numérique pour la simulation des réseaux d’assainissement pluvial fondé sur
le concept de stockage.
Institut National des Sciences Appliquées de Lyon, 1996, 260 p., Thèse de doctorat, Spécialité :
Conception en Bâtiment et Techniques Urbaines, N° d’ordre : 96 ISAL 0029.
M. Nujic
Efficient implementation of non oscillatory schemes for the computation of free surface
flows.
167
Bibliographie
A. Paquier
Modélisation et simulation de la propagation de l’onde de rupture de barrage.
Université Jean Monnet de Saint-Etienne, 1995, 192 p., Thèse de doctorat, Spécialité analyse
numérique.
C.E. Rice
Open channel junctions with Supercritical Flow.
United States Department of Agriculture. Agricultural Research Service
ARS-14, January 1985
P.L. Roe
Approximate Riemann solvers, parameter vectors, difference schemes.
Journal of Computational Physics 43, 357-372, 1981.
P.L. Roe
Some contributions to the modelling of discontinuous flows.
In proceedings of the SIAM/AMS Seminar, San Diego, 1983.
D. Rollet
Instrumentation des déversoirs d’orage : Mise en place de l’autosurveillance à Sélestat
Mémoire de fin d’étude ENGEES, 2002.
H. H. G. Savenije.
Analytical expression for tidal damping in alluvial estuaries.
Journal of Hydraulic Engineering, vol. 124, N° 6, June 1998, pp.615-618.
168
Bibliographie
J. Smoller.
Shock waves and reaction-diffusion equations. Second edition.
Ed. :Springer-Verlag-1984-632 p .
P.K. Sweby
High resolution schemes using flux limiters for hyperbolic conservation laws.
SIAM Journal of Numerical Analysis 21, 995-1011 (1984).
R. Szymkiewicz
Finite element method for the solution of the Saint-Venant equations in an open channel
network.
Journal of Hydrology–1991-Vol. 122–pp. 275-287.
C. Thirriot et L. Benayada
Solution analytique approchée des équations de l’écoulement transitoire à surface
libre en canal avec perte de charge et pente de fond.
C. R. Acad. Sci. Paris, t. 320, Série IIb, p. 325-330, 1995.
E. F. Toro
Riemann solvers and numerical methods for fluid dynamics.
A practical introduction.
Ed. Sringer-592 p-1997
E. F. Toro
Shock-capturing methods for free surface flows
Ed. John Whiley & Sons, 2000
B. Van Leer
Towards the ultimate conservation difference scheme I. The quest of monotonicity.
Lecture Notes in Physics, Vol. 18,163-168, 1973.
B. Van Leer
Towards the ultimate conservation difference scheme V. A second order sequel to
Godunov’s method.
Journal of Computational Physics 32, 101-136, 1979
169
Bibliographie
J-P Vila
Sur la théorie et l’approximation numérique de problèmes hyperbolique non linéaires.
Application aux équations de Saint Venant et à la modélisation des avalanches de neige
dense.
Université de Paris VI, 1986, 480 p., Thèse de doctorat, Spécialité Sciences Mathématiques.
170
Annexes
ANNEXES
171
Annexes
Annexe 1
Présentation de l’approche de dimensionnement du
déversoir par Giersch.
hav
ham
tam Dam
Tram pam pav
Dav
C’est une démarche de dimensionnement pour des déversoirs latéraux à crête haute :
• On impose au départ une hauteur de crête amont (pam) et une longueur de déversoir
(Ldo) ;
172
Annexes
d’écoulement telles que le ressaut se produise dans la conduite amont avant le déversoir.
Donc, on prendra pour le calcul :
¾ Tram=hauteur d’eau réelle à l’amont (tam) si le régime est fluvial dans la
conduite amont.
¾ Tram ≥ la hauteur conjuguée de tam, s’il est torrentiel.
• On fait ensuite une vérification des hypothèses de départ au débit décennal, basée
sur la conservation de l’énergie sur le seuil, en admettant uniquement une perte de
charge linéaire de 1% et en calculant la hauteur d’eau moyenne nécessaire pour
déverser Q10-Qcr par la formule de Poléni.
• Hypothèses initiales :
4 × Qcr
Longueur du déversoir (Ldo) : en 1ère approximation on peut prendre Ldo =
Dam
173
Annexes
Le calcul est basé sur le théorème de Bernoulli pour le débit critique. On est en régime
fluvial (imposé par le choix de Tram) et la ligne d’eau monte sur le seuil ; au débit critique, il
ne doit pas y avoir déversement. On suppose que la ligne d’eau arrive juste à la hauteur de la
crête à l’aval et on considère que la charge hydraulique sur la conduite étranglée est égale à la
hauteur de crête aval :
Par des considérations de pertes de charge sur la conduite aval, on peut montrer que :
Vav 2
Pav − ν − Dav
2g
Lav =
Imot − Iav
avec :
Imot : perte de charge dans la conduite aval (valeur tablée)
ν : coefficient défini en fonction du rapport entre charge et diamètre aval.
On suppose que la seule perte d’énergie sur le déversoir est due à une perte de charge linéaire
de 1 cm/m. Soit Ham l’énergie spécifique à l’amont du déversoir (par rapport au radier amont)
et Hav l’énergie spécifique à l’aval du déversoir (toujours par rapport au radier amont), on doit
avoir :
Hav ≈ Ham – 0,01×Ldo.
Ligne d’énergie
hav
ham
Ham Hav
pam
pav
Dav
On connaît le tirant d’eau amont admissible donc la hauteur sur le seuil à l’amont :
On calcule la hauteur d’eau moyenne (hm) sur le seuil pour évacuer le débit à déverser à
l’aide de la formule de Poléni :
174
Annexes
4hm − ham
On calcule la hauteur d’eau sur le seuil à l’aval du déversoir : hav = ;
3
D’où la charge spécifique à l’aval du déversoir (par rapport au niveau du radier amont) en
supposant que la vitesse est nulle :
Dans tous les cas, on doit avoir Hav < Ham – 0,01×Ldo. Si cette condition n’est pas vérifiée,
on recommence les calculs en augmentant Ldo ou Tram jusqu’à ce que les valeurs d’énergie
soient les plus proches possibles.
1.4 Remarques
• Dans le calcul de la hauteur d’eau moyenne hm : la formule utilisée pour le calcul du débit
déversé est proche de la formule de Poléni (qui s’applique normalement à des déversoirs
frontaux).
• Dans la valeur de la perte d’énergie sur le seuil de 0.01 m/m qui n’est pas justifiée.
175
Annexes
cherche à calculer ensuite la longueur du déversoir nécessaire pour avoir hav = 0 à l’aval du
déversoir au débit décennal :
hav=0
Tram
pam
pav
Cependant, pour obtenir hav=0, on est obligé de supposer que la ligne d’eau diminue
sur le seuil alors qu’on est en régime fluvial (on a imposé Tram pour cela) ; cette hypothèse est
à priori fausse. Cette ligne d’eau correspond à un régime torrentiel ; si l’on suppose qu’à
l’amont du déversoir on est en régime fluvial, on peut envisager les cas (b) ou (c) mais on ne
peut pas prendre alors comme condition à la limite la hauteur normale.
Conclusion :
Pour les seuils latéraux à crête haute, la démarche consiste à déterminer la hauteur de crête et la
longueur du déversoir à partir du calcul de perte de charge ; le principe du calcul prend en
compte une ligne d’eau réaliste mais quantifie son évolution de manière empirique :
l’utilisation de la formule de Poléni pour le calcul du débit déversé et la définition de la hauteur
d’eau moyenne sur le seuil restent discutables. Pour les seuils latéraux à crête basse, le calcul
est basé sur l’hypothèse que la ligne d’eau diminue sur le seuil pour arriver à une hauteur aval
nulle, ce qui est à priori faux.
176
Annexes
Annexe 2
Généralités sur les systèmes d’équations hyperboliques
Les lois de conservation sont des équations aux dérivées partielles issues de la physique
qui peuvent être écrites sous la forme :
u t + f ( u )x = 0
U t + F ( U )x = 0
où :
u1 f1
u f
U = 2 , F (U) = 2
M M
u f
m m
U est appelé vecteur des variables conservées et F est le vecteur flux dont chaque composante
Fi est une fonction des composants uj de U.
Les entrées de la matrice J sont les dérivées partielles des composantes fi du vecteur flux F par
rapport aux composantes uj du vecteur des variables conservées c’est à dire a ij = ∂fi ∂u j
Les équations de conservation présentées dans la définition 1 peuvent également être écrites
sous une forme quasi-linéaire lorsqu’elles sont sans second membre. Celle-ci est obtenue en
décomposant le terme de flux de la manière suivante :
∂F ( U ) ∂F ∂U
=
∂x ∂U ∂x
177
Annexes
Ut + J ( U ) Ux = 0
qui est un cas particulier d’un système d’équations sans second membre.
Les valeurs propres λi de la matrice J sont les solutions des polynômes caractéristiques :
J − λ i I = det ( J − λ i I ) = 0
où :
I est la matrice identité. Les valeurs propres de la matrice jacobienne d’un système de
lois de conservations avec second membre sont appelées les valeurs propres du système.
Un vecteur propre à droite d’une matrice J correspondant à une valeur propre λi de J est
T
un vecteur K ( ) = k1( ) , k (2 ) ,L , k (m) qui satisfait la relation JK ( ) = λ i K ( ) . De manière similaire,
i i i i i i
Un système d’équation avec second membre est dit hyperbolique en un point du plan (x,
t) si J a m valeurs propres λ1 ,L , λ m et un ensemble correspondant de vecteurs propres
K ( ) ,L , K ( ) . Le système est dit strictement hyperbolique si les valeurs propres λi sont toutes
1 m
distinctes.
Une matrice J est dite diagonalisable si elle peut être exprimée sous la forme :
J = KΛK −1 où Λ = K −1JK
en terme d’une matrice diagonale Λ et d’une matrice K. Les éléments diagonaux de Λ sont les
valeurs propres λ i de J et la colonne K(i) de K sont les vecteurs propres à droite de J
correspondant aux valeurs propres λ i , c’est à dire :
178
Annexes
λ1 L 0
0 L 0
Λ= , K = ( K (1) ,L , K (m) ) , JK (i) = λ i K (i)
M M M
0 L λm
179
Annexes
Annexe 3
Calculs concernant la mise en œuvre de la méthode de Roe
dans le cas du système de Saint-Venant
∂Q ∂Q
∂F(U) ∂A ∂Q
J= =
∂U
∂ Q2
∂ Q 2
+ gI + gI
∂A A ∂Q A
Les détails des calculs des éléments de la matrice dans le cas du canal prismatique sont
présentés ci-après :
∂Q
=0 Car Q et A sont des variables indépendantes.
∂A
∂Q
=1
∂Q
h (x ,t )
∂ Q2 Q2 ∂
+ gI = − 2 + g ∫ (h − η)b(x , η)dη
∂A A A ∂A 0
h (x ,t )
Q2 ∂ ∂h
=− + g ∫ (h − η)b(x , η)dη
A 2
∂h 0 ∂A
Q2 h ( x , t ) ∂h
=− + g ∫ b(x, η)dη
A 2
0 ∂A
Q2 ∂h
= − 2 + gA
A ∂A
−1
A ∂h ∂A
= − u + g car
2
= puisque l' on est dans le cas d' un canal
b ∂A ∂h
prismatique A n' est fonction que de h
h (x ,t )
∂ Q2 2Q ∂I ∂
+ gI = +g = 2u + g ∫ (h − η)b(x , η)dη = 2u
∂Q A A ∂Q ∂Q 0
0 1
J= A
g − u
2
2u
b
Q
avec u = .
A
180
Annexes
gA
a 1, 2 = u ± c avec c=
b
1
e 1,2 =
u ± c
∂U ~ ∂U
+J =G
∂t ∂x
~
où J correspond à la matrice jacobienne approchée.
On rappelle les conditions que doit satisfaire la matrice jacobienne approchée dans le cadre du
solveur de Roe :
~ u 1i +, 21 ± ~ci1+, 21
a i1+, 21 = ~
2 2 2
~e 1, 2 = 1, 2
1
i+ 1 ~
a i+ 1
2
2
~
u et de ~c qui satisfont
Le problème de la détermination de J est transféré en celui du calcul de ~
simultanément les première et cinquième conditions précédemment énumérées.
~
Soit R i + 1 la matrice des vecteurs propres. Alors, J peut être diagonalisée sous la forme :
2
~
J = R i + 1 diag ~
a ik+ 1 R i−+11
2 2 2
La quatrième condition que doit vérifier le jacobien approché implique que ∆ i+ 1 U peut
2
s’exprimer sous la forme d’une combinaison linéaire des vecteurs propres :
181
Annexes
2
~k
∆ i + 1 U = ∑ α ik+ 1 e i+ 1
2 k =1 2 2
On détermine grâce à la relation précédente les valeurs des coefficients α1i + 1 et α i2+ 1 :
2 2
2
∆ i + 1 U = ∑ αik+ 1 e% ki+ 1
2 k =1 2 2
A − Ai 1 1
⇒ i +1 = α i + 1 2 u
1
+ αi2+ 1
Qi +1 − Qi +c 2 u i + 1 − ci + 1
i+ 12 i+ 12 2 2
∆ i + 1 A = α1i + 1 + α i2+ 1
2 2 2
⇒
2 2
( 2 2
)
∆ i + 1 Q = α i + 1 u i + 1 + ci + 1 + α i + 1 u i + 1 − ci + 1
1 2
2 2
(
2
)
α1i + 1 = ∆ i + 1 A − αi2+ 1
2 2 2
⇒
2 2
(
2 2
2
)
2 2 2
(
2
2
) 2
(
∆ i + 1 Q = ∆ i + 1 A u i + 1 + ci + 1 − αi + 1 u/ i + 1 + ci + 1 + αi + 1 u/ i + 1 − ci + 1
2
)
α1i + 1 = ∆ i + 1 A − αi2+ 1
2 2 2
⇒
α 2 1 = 2 2
(
∆ i + 1 Q - u i + 1 + ci + 1 ∆ i + 1 A
2 2
)
i+
2 −2ci + 1
2
1
α = 2 2 2 2 2
(
/ i + 1 ∆ i + 1 A - ∆ i + 1 Q + u i + 1 + c/ i + 1 ∆ i + 1 A
−2c
2
)
i+ 12 −2ci + 1
2
⇒
αi2+ 1 = 2 2
(
∆ i + 1 Q - u i + 1 + ci + 1 ∆ i + 1 A
2 2
)
2 −2ci + 1
2
1
α = 2 2
(
−∆ i + 1 Q + u i + 1 + ci + 1 ∆ i + 1 A
2 2
)
i+ 12 −2ci + 1
2
⇒
α2 = i+ 2 2
(
∆ 1 Q - u i + 1 + ci + 1 ∆ i + 1 A
2 2
)
i + 1
2 −2ci + 1
2
182
Annexes
1
∆ i + 1 Q + − u i + 1 + ci + 1 ∆ i + 1 A ( )
α i + 1 2 =
2 2 2 2
2ci + 1
2
⇒
α 2 = i + 2
∆ 1 Q + − u i + 1 − ci + 1 ∆ i + 1 A
2 2 2
( )
i+ 12 −2ci + 1
2
⇒α 1,2
= 2
(
∆ i + 1 Q + − u i + 1 ± ci + 1 ∆ i + 1 A
2 2
) 2
i+ 1
2 ±2ci + 1
2
Démonstration :
On a :
~
• Ji + 1 ∆ i + 1 U = ∆ i + 1 F
2 2 2
2
• ~k
∆ i + 1 U = ∑ α ik+ 1 e i+ 1
2 k =1 2 2
~
• J = R i + 1 diag ~
a ik+ 1 R i−+11
2 2 2
Alors :
~
∆ i + 1 F = Ji + 1 ∆ i + 1 U
2 2 2
2
~ ~k
⇒ ∆ i + 1 F = Ji + 1
2
∑α
2 k =1
k
i+ 1
2
e i+ 1
2
~ ~ 1 + ~J α 2 e ~2
⇒ ∆ i+ 1 F = Ji + 1 α1i + 1 e i+ 1 i+ 1 i+ 1 i+ 1
2 2 2 2 2 2 2
⇒ ∆ i + 1 F = α1i + 1 R i + 1 diag ~
a k R −1 e ~ 1 + α 2 R diag ~ −1 ~ 2
a i + 1 R i + 1 e i+ 1
k
2 2 2 i + 12 i + 12 i+ 12 i+ 1
2
i+ 1
2 2 2 2
avec k = 1,2.
a b 1 1 d − b
Soit A = alors A -1 = adj(A ) =
c d A A − c a
On a :
183
Annexes
1 1 1 ~
a2 − 1
R i+ 1 = ~1 et avec la remarque précédente : R −11 = i + 12
~ ~ − ~
a i+ 1 − ~ 1
2
2 a i+ 1 a i+ 1 i+
2
2
a i1+ 1 a1
2 2 2 2 i + 12
Calcul de :
R i + 1 diag ~
a ik+ 1 R i−+11 e~1
2 2 2
i+ 1
2
~
a2 − 1 1
~ 1 i + 12
= R i + 1 diag a ~ 1
k
~ 2 ~1 − ~
2 2 a
i+ 1
i + 12 − a i + 12 a1 1 a i + 1
i + 12 2
~a i1+ 1 0 1
= R i+ 1 2
2 0 a i + 1 0
~ 2
2
1 1 a 1 ~ 1
= ~1 ~ i + 2
a i + 1 a i + 1 0
2
2 2
~ a1
i + 12 ~ 1 1 ~ 1 ~ 1
= 2 = a 1 =a e
~
a 1 i+ 1 ~
2 a i+ 1
i+ 1 i+ 1
i+ 1 2
2 2
2
D’où :
~ ~ 1 + ~J α 2 e ~2
∆ i + 1 F = Ji + 1 α1i + 1 e i+ 1 i+ 1 i+ 1 i+ 1
2 2 2 2 2 2 2
2
⇒ ∆ i + 1 F = ∑ α ik+ 1 ~ ~k
a ik+ 1 e i+ 1
2 k =1 2 2 2
184
Annexes
2
∆ i + 1 F = ∑ α ik+ 1 ~ ~k
a ik+ 1 e i+ 1
2 k =1 2 2 2
∆Q + (−u + c)∆A
∆Q = (u + c ) + ∆Q + (−u − c)∆A (u − c)
2c − 2c
⇒ 2
Q ∆Q + (− u + c)∆A ∆Q + (−u − c)∆A
∆ + g∆I = (u + c )(u + c ) + (u − c )(u − c )
A 2c − 2c
∆Q c/ 2 − u/ 2 ∆Q u/ 2 − c/ 2
∆Q = (u/ + c ) + ∆A − (u/ − c) + ∆A
2c 2c 2c 2c
⇒
Q2 ∆Q
∆ + g∆I = (u + c )2 − ∆Q (u − c )2 + ∆A (− u + c )(u + c )2 + ∆A (u + c )(u − c )2
A 2c 2c 2c 2c
∆Q = ∆Q
Q
2
∆Q
∆ + g∆I = (u + c/ + u − c/ )(u/ + c − u/ + c )
⇒ A 2c
∆A 2
+ {( ) (
u + c 2 + 2uc (− u + c ) + u 2 + c 2 − 2uc (u + c ) ) }
2c
Q2 ∆A
⇒ ∆ + g∆I = ∆Q + {
− u/ 3 + u 2 c − c/ 2 u/ + c 3 − 2u 2 c + 2/ u/ c/ 2 + u/ 3 + u 2 c + c/ 2 u/ + c 3 − 2u 2 c − 2/ u/ c/ 2 }
A 2c
Q2 ∆A
⇒ ∆
A
+ g∆I = ∆Q +
2c
{
2c 3 − 2u 2 c }
Q2
⇒ ∆ + g∆I = ∆Q − u 2 ∆A + c 2 ∆A
A
Q2
⇒ u 2 ∆A − 2u∆Q + ∆ = c 2 ∆A − g∆I
A
∆I Q2
On pose c 2 = g d' où u 2 ∆A − 2u∆Q + ∆ = 0
∆A A
Q2
∆ = 4Q 2 − 4∆A∆
A
185
Annexes
Q2
∆ = 4Q 2 − 4∆A∆
A
Q2 Q2
⇒ ∆ = 2 (Q i +1 − Q i ) − (A i +1 − A i ) i +1 − i
2
A i +1 A i
2 A i +1Q i2 A i Q i2+1 2
⇒ ∆ =2 Q 2
/ i +1 − 2Q Q (
i i +1 + Q 2
−
/ i / i +1
Q −
Ai
− )
A i +1
/ i
+Q
− 2Q i Q i +1 A i A i +1 + A i2+1Q i2 + A i2 Q i2+1
⇒ ∆ =2
A i A i +1
(A i+1Q i − A i Q i+1 )
⇒ ∆ =2
A i A i +1
A i +1 Ai
2/ (Q i +1 − Q i ) ± 2/ Qi − Q i +1
~
u i+ 1 = Ai A i +1
2 2/ (A i +1 − A i )
− A i ± A i +1 A i +1 m A i
Qi + Q i +1
Ai A i +1
⇒~
u i+ 1 =
2 ( A i +1 + A i )( A i +1 − A i )
( A i +1 − A i )Q + ( A i +1 − A i )Q
i i +1
Ai A i +1
⇒~
u i+ 1 =
2 ( A i +1 + A i )( A i +1 − A i )
Qi A i + Q i +1 A i +1
⇒~
u i+ 1 =
2 A i +1 + A i
186
Annexes
Annexe 4
Caractéristiques hydrauliques des collecteurs circulaires et
ovoïdes
Les caractéristiques hydrauliques à calculer sont les rayon hydraulique (Rh), surface mouillée
(S), périmètre mouillé (P), hauteur en fonction de la surface (h), largeur au miroir (B), pression
hydrostatique (I), hauteur normale (hn), hauteur critique (hc) et débit normal (Qn).
Nous avons considéré les conduites rectangulaires, circulaires en fer a cheval et ovoïdes. Les
conduites de formes arrondies sont caractérisés par le rapport de leur hauteur (T) sur leur
largeur (B). Les calculs ont été réalisés dans le cas des conduites ovoïdes 1.5 et 4/3, fer à cheval
0.75 et circulaire 1.
Pour connaître les caractéristiques hydrauliques des ovoïdes T=1.5 B, circulaires T=B et
en fer à cheval T=0.75 B, on s’est appuyé sur les travaux réalisés par Hager [Hager-1999]. En
effet, ce dernier a établi des relations hydrauliques adimensionnelles pour ces conduites.
Pour ce qui est des ovoïdes T=4/3 B, à partir de la géométrie de ces conduites, on a
construit les équations hydrauliques adimensionnelles sous le même modèle que celui de Hager
en s’inspirant des travaux de K. Babaeyan-Koopaei pour les expressions des grandeurs
hydrauliques. Les expressions de ces grandeurs utilisées par K. Babaeyan-Koopaei [Babaeyan-
Koopaei-2001] sont présentées sur la figure ci-dessous :
187
Annexes
Puis nous avons tracé sous excel les graphes donnant les courbes Rh, S, B, I en fonction
de h/T et finalement calé des polynômes sur ces courbes. Ce sont ces polynômes qui ont été
intégré au logiciel pour permettre le calcul des grandeurs hydrauliques.
Pour ce qui concerne les termes de pression hydrostatiques dans le cas des conduites
ovoïdes et en fer à cheval, celles-ci ont été calculées en mettant en œuvre des relations valables
pour le cas elliptique. Les conduites ovoïdes ont été approchées par des conduites elliptiques
comme représenté ci-dessous dans le cas de l’ovoïde T=1.5B et du fer à cheval T=0.75B.
I 1
(
= (2 y − 1)arccos(1 − 2 y ) − 2(1 − 2 y ) y − y 2
) 1
2
3
(
+ 2 y − y2 ) 3
2
avec y =
h
ρgBT 2
8 T
Pour les conduites circulaire et rectangulaire les intégrales permettant d’obtenir le terme de
pression ont été calculées.
188
Annexes
h 4 h h h 4 h 3 h 2
T( − 0 .8181( ) + 1.5399 ( )3 − 1.1872 ( ) 2 T( − 0 .7813 ( ) + 1.4252 ( ) − 1.1446 ( )
T T T T T T
Rh h h
+ 0 .667 ( ) − 0 .0001 ) + 0 .719 ( ) − 0 .0013 )
T T
h h h h h h
S T 2 (−0.5455( ) 3 + 0.9431( ) 2 + 0.1193( ) + 0.001) T 2 (−0.5455( ) 3 + 0.9431( ) 2 + 0.1193( ) + 0.001)
T T T T T T
h 6 h h
) + 158.54( )5 − 179.74( )4
h 6 h h
T( − 54.474( T( − 59.561( ) + 172 .59( )5 − 195 .84( ) 4
T T T T T T
B h h h
h 3 h 2 h + 109 .39( )3 − 32.359( ) 2 + 5.7314( ) + 0.0668 )
+ 99.223( ) − 28.411( ) + 4.8798( ) + 0.0402)
T T T T T T
h h h h h h
I T 3 (0.1313( )3 + 0.1435( ) 2 − 0.0.126 ( ) + 0.0004 ) T 3 ( 0 .052 ( ) 3 + 0 .2619 ( ) 2 − 0 .0335 ( ) + 0 .0013 )
T T T T T T
Qps 0 . 171 K s I 0 . 5 T 8 / 3 0 . 221 K s I 0 .5
T 8/3
Q Q
Hn 1 . 09 T (1 − (1 − 0 . 884 ) 0 .5 1 . 0159 T (1 − (1 − ) 0 .5
( 0 . 171 K s I 0 .5T 8 / 3 ) 0 .5 ( 0 . 221 K s I 0 . 5 T 8 / 3 0 .5
)
hn 2 hn 2
0 . 171 K s I 0 . 5 T 8 / 3 1 . 9 ( ) × 0 . 221 K s I 0 . 5 T 8 / 3 1 . 942 ( ) ×
T T
Qn
hn 2 hn 2
(1 − 0 . 421 ( ) ) (1 − 0 . 487 ( ) )
T T
Q Q
Hc 1 . 34 ( ) 0 .5 1 . 28 ( ) 0 .5
( 9 . 81 T ) 0 . 5 ( 9 . 81 T ) 0 . 5
h 4 h h
T( − 0.5257 ( ) + 0 .7707 ( ) 3 − 0 .7120 ( ) 2 h 3 h h
T T T T( − 0.2748 ( ) − 0 . 17 ( ) 2 + 0.7505 ( )
Rh h
T T T
+ 0 .7300 ( ) − 0.00155 ) − 0.0049 )
T
h h h h h h
S T 2 (−0.9053( ) 3 + 1.3570( ) 2 + 0.3401( ) + 0.003527 T2 (−0.9333( )3 + 1.1659( )2 + 0.8522( ) − 0.0228)
T T T T T T
h h h
T( − 74.49( ) 6 + 222.6( ) 5 − 260( ) 4 h 6 h h
T T T T( − 60.168( ) + 184 .57 ( )5 − 227 .09( ) 4
T T T
B h h h
+ 150( ) 3 + 46.71( ) 2 + 8.603( ) + 0.09642) h 3 h 2 h
+ 143 .45( ) − 50.724 ( ) + 9.8533( ) + 0.1616 )
T T T T T T
h h h h h
I T 3 (0.4527 ( ) 2 − 0.06008 ( ) + 0.001745 ) T 3 ( 0 .0688 ( ) 3 + 0 .5592 ( ) 2 − 0 .073 ( ) + 0 .0034 )
T T T T T
π
Qps 5
K s I 0 .5 T 8 / 3 0 . 457 K s I 0 . 5 T 8 / 3
4 3
189
Annexes
Q Q
Hn 0 .926 T (1 − (1 − 3 .11 ) 0 . 5 ) 0 .5 0 . 850 T (1 − (1 − 2 . 19 0 .5
) 0 .5 ) 0 .5
K sI T 0 .5 8/3
K sI T 8/3
h 2 h
0 . 457 K s I 0 . 5 T 8 / 3 2 . 8 ( ) (1 − 0 . 71 ( ) 2 )
T T
5
pour h/T<0.8
2
π π δ − sin δ cos δ
3 3 h 2 h
K s I 0 .5 T 8 / 3 0 . 457 K s I 0 .5 T 8 / 3 2 . 8 ( ) (1 − 0 . 8 ( ) 2
Q 5
δ π T T
4 3
h 6
+ 0 . 25 () )
T
pour h/T<0.93
0 .5
Q Q
Hc 0 . 787 ( ) 0 .5
( 9 . 81 T )
0 .5
( 9 . 81 T ) 0 . 5
Pour le canal rectangulaire, les seules difficultés résident dans le calcul de la hauteur normale et
de la définition d’un débit à pleine section. Pour la hauteur normale, on utilise la relation de
Manning-Strickler dans laquelle on remplace les section mouillée et rayon hydraulique par
leurs expressions en fonction de la hauteur d’eau et on calcule la hauteur normale en
déterminant le zéro de la fonction grâce à un solveur. Pour le débit à pleine section, on a
simplement limité la hauteur d’un canal rectangulaire à 4 fois sa largeur.
On a ainsi pu créer dans le logiciel des fonctions permettant le calcul de ces grandeurs
hydrauliques sans que ces calculs soient trop coûteux en temps.
190
Annexes
Annexe 5
Présentation du logiciel hsl
Un bon moyen d’illustrer l’application des conditions aux limites est le logiciel nommé
hsl qui permet de reproduire les courbes de remous susceptibles d’apparaître dans un unique
canal rectangulaire en écoulement à surface libre.
Les valeurs des conditions aux limites appliquées sont matérialisées par deux triangles
rouges à l’amont et à l’aval du canal.
3 PRESENTATION DU LOGICIEL
On présente ci-dessous la fenêtre d’accueil du logiciel hsl
Les conditions aux limites applicables à l’amont et à l’aval sont les hauteurs normales,
critiques ou imposées.
4 QUELQUES EXEMPLES
4.1 Courbe M2
191
Annexes
Dans le cas suivant, il s’agit d’un canal à pente faible dans lequel on a fait apparaître
une courbe de type M2 en imposant la hauteur normale à l’amont et la hauteur critique à l’aval.
L’application de la hauteur critique à l’aval est possible puisque l’on est dans le cas d’un
écoulement fluvial c’est à dire que le point de contrôle est à l’aval. L’application de la hauteur
normale à l’amont n’est pas prise en compte et il y a linéarisation. Sur les 200 m de longueur du
canal, la courbe M2 ne se développe pas complètement et on ne retrouve pas la hauteur
normale à l’amont. Le fait que la hauteur critique ne soit pas atteinte à l’aval est du à
l’affectation de la condition limite dans les cellules du tronçon fictif aval.
4.2 Courbe M3
192
Annexes
conjuguée du ressaut hydraulique, dans ce cas le ressaut apparaît et se propage vers l’aval. On a
ainsi à l’amont du ressaut une courbe de remous de type M3 et à l’aval on retrouve la hauteur
normale.
Dans la même configuration, si on applique une hauteur imposée de 1.3 m cette fois le
ressaut n’apparaît pas car la hauteur imposée est supérieure à la valeur de la hauteur conjuguée
du ressaut. Dans ce cas, c’est une linéarisation sur les deux cellules fictives et la cellule amont
du tronçon réel qui est réalisée. En fait, on reste à la hauteur normale.
4.3 Courbe S1
On considère un canal de pente 1%, de largeur 1.2 m, de longueur 200m, de rugosité 100 et
un débit injecté de 6 m3s-1.
193
Annexes
On considère un canal de 100 m de long de rugosité 100, de pente 1%, de largeur 1.2 m et un
débit injecté de 10 m3s-1.
194
Annexes
effectivement appliquée à l’amont mais pas à l’aval. A l’aval, on a linéarisation. Ainsi, on voit
apparaître une courbe de remous de type S2.
195
Annexes
Annexe 6
Tableaux de valeurs concernant la validation du logiciel
pour les courbes de remous.
• Courbe S1 :
COURBE S1
Abscisse Hauteurs mesurées Hauteurs simulées
Méthode Méthode Mac-
Méthode Upwind Méthode MUSCL
symétrique Cormack
0,25 2,2 2,3 2,3 2,3 2,3
0,75 2,5 2,4 2,4 2,4 2,5
1,25 2,6 2,8 2,5 2,5 2,8
1,75 2,8 7,6 2,7 2,8 4,4
2,25 6,9 8,5 5,6 5,2 7,6
2,75 7,7 9,1 7,6 7,6 8,1
3,25 8,3 9,7 8,1 8,1 8,7
3,75 8,8 10,4 8,7 8,7 9,2
4,25 9,3 10,9 9,2 9,2 9,7
4,75 9,8 11,6 9,7 9,7 10,3
5,25 10,3 12,0 10,3 10,3 10,8
5,75 10,7 12,7 10,8 10,8 11,3
• Courbe S3 :
COURBE S3
Abscisse Hauteurs mesurées Hauteurs calculées
Méthode Méthode Mac-
Méthode Upwind Méthode MUSCL
symétrique Cormack
0,25 0,7 0,8 0,8 0,8 0,8
0,75 0,8 0,9 0,9 0,9 0,9
1,25 1,0 1,1 1,1 1,1 1,1
1,75 1,1 1,2 1,2 1,2 1,2
2,25 1,3 1,3 1,3 1,3 1,3
2,75 1,4 1,5 1,3 1,3 1,5
3,25 1,4 1,5 1,5 1,5 1,5
3,75 1,5 1,6 1,5 1,5 1,5
4,25 1,4 1,6 1,5 1,5 1,5
4,75 1,5 1,6 1,6 1,6 1,6
5,25 1,5 1,6 1,6 1,6 1,6
5,75 1,5 1,6 1,6 1,6 1,6
196
Annexes
• Courbe M1 :
COURBE M1
Abscisse Hauteurs mesurées Hauteurs calculées
Méthode Méthode Mac-
Méthode Upwind Méthode MUSCL
symétrique Cormack
0,25 12,3 13,5 12,5 12,5 12,7
0,75 12,3 13,7 12,7 12,7 12,8
1,25 12,4 13,9 12,8 12,8 13,1
1,75 12,6 14,1 12,9 12,9 13,1
2,25 12,8 14,3 13,1 13,1 13,3
2,75 13,1 14,5 13,2 13,3 13,5
3,25 13,3 14,8 13,5 13,5 13,6
3,75 13,5 14,9 13,5 13,6 13,7
4,25 13,6 15,1 13,7 13,7 13,9
4,75 13,7 15,3 13,9 13,9 14,0
5,25 13,9 15,5 14,0 14,0 14,1
5,75 13,9 15,6 14,1 14,1 14,4
• Courbe M2 :
COURBE M2
Abscisse Hauteurs mesurées Hauteurs calculées
Méthode Méthode Mac-
Méthode Upwind Méthode MUSCL
symétrique Cormack
0,25 1,3 1,3 1,3 1,3 1,3
0,75 1,5 1,6 1,6 1,6 1,6
1,25 1,6 1,7 1,7 1,7 1,7
1,75 1,9 2 2 2 2,0
2,25 2,2 2,3 2,1 2,1 2,3
2,75 2,4 2,5 2,4 2,4 2,7
3,25 4,6 4,8 4,8 3,3 4,8
3,75 4,7 4,7 4,7 4,7 4,7
4,25 4,6 4,5 4,5 4,5 4,5
4,75 4,5 4,4 4,4 4,4 4,4
5,25 4,2 4,3 4,3 4,3 4,3
5,75 3,7 3,9 3,9 4,1 3,9
197
Annexes
• Courbe M3 :
COURBE M3
Abscisse Hauteurs mesurées Hauteurs calculées
Méthode Méthode Mac-
Méthode Upwind Méthode MUSCL
symétrique Cormack
0,25 1,9 1,9 1,9 1,9 1,8667
0,75 2 2,0 2,0 2,0 2,0000
1,25 2,3 2,3 2,3 2,3 2,2667
1,75 2,5 2,4 2,4 2,4 2,4000
2,25 2,7 2,7 2,7 2,7 2,6667
2,75 2,9 2,8 2,8 2,8 2,8000
3,25 3,1 3,1 3,1 3,1 3,0667
3,75 3,4 3,3 3,3 3,3 3,3333
4,25 3,6 3,6 3,6 3,6 3,6000
4,75 4,2 4,1 4,3 3,9 4,9333
5,25 4,9 4,9 4,8 5,3 4,8000
5,75 4,5 4,7 4,7 4,9 4,6667
198
Annexes
hauteur hauteur
hauteur hauteur hauteur hauteur hauteur
Position expérimentale expérimentale - Position Position Position Position
expérimentale calculée calculée calculée calculée
(m) + 3mm 3mm (m) (m) (m) (m)
(m) (m) (m) (m) (m)
(m) (m)
1.8 0.078 0.081 0.075 1.5 0.07643 2.135 0.06964 2.7275 0.05659 3.49 0.01916
2 0.077 0.08 0.074 1.525 0.07641 2.155 0.06841 2.751 0.05634 3.515 0.01928
2.04 0.072 0.075 0.069 1.55 0.07639 2.175 0.06713 2.7745 0.05599 3.54 0.01939
2.08 0.067 0.07 0.064 1.575 0.07636 2.195 0.06586 2.798 0.05551 3.565 0.01948
2.16 0.061 0.064 0.058 1.6 0.07632 2.215 0.06463 2.8215 0.05481 3.59 0.01957
2.24 0.056 0.059 0.053 1.625 0.07629 2.235 0.06352 2.845 0.0517 3.615 0.01966
2.32 0.051 0.054 0.048 1.65 0.07628 2.255 0.06257 2.87475 0.04634 3.64 0.01975
2.4 0.052 0.055 0.049 1.675 0.07627 2.275 0.06182 2.9045 0.04127 3.665 0.01984
2.48 0.052 0.055 0.049 1.7 0.07627 2.295 0.06125 2.93425 0.03792 3.69 0.01993
2.56 0.052 0.055 0.049 1.725 0.07627 2.315 0.06083 2.964 0.03532 3.715 0.02002
2.64 0.052 0.055 0.049 1.75 0.07626 2.335 0.06049 2.99375 0.03317 3.74 0.02013
2.72 0.052 0.055 0.049 1.775 0.07625 2.355 0.06016 3.0235 0.03133 3.765 0.02024
2.8 0.049 0.052 0.046 1.8 0.07623 2.375 0.05987 3.05325 0.02974 3.79 0.02035
2.88 0.045 0.048 0.042 1.825 0.07618 2.3985 0.05963 3.083 0.02835 3.815 0.02044
2.96 0.04 0.043 0.037 1.85 0.07613 2.422 0.0594 3.11275 0.02715 3.84 0.02052
3.04 0.029 0.032 0.026 1.875 0.07610 2.4455 0.05918 3.1425 0.02606 3.865 0.02059
3.12 0.019 0.022 0.016 1.9 0.07609 2.469 0.05897 3.17225 0.02487 3.89 0.02066
3.2 0.019 0.022 0.016 1.925 0.07608 2.4925 0.05874 3.202 0.02406 3.915 0.02071
3.28 0.019 0.022 0.016 1.95 0.07608 2.516 0.05851 3.23175 0.02332 3.94 0.02092
3.36 0.019 0.022 0.016 1.975 0.07579 2.5395 0.05828 3.2615 0.02257
3.44 0.022 0.025 0.019 1.995 0.07507 2.563 0.05805 3.29125 0.02186
3.52 0.018 0.021 0.015 2.015 0.07499 2.5865 0.05782 3.321 0.02125
3.6 0.016 0.019 0.013 2.035 0.07430 2.61 0.05761 3.35075 0.02066
2.055 0.07359 2.6335 0.0574 3.3805 0.02012
2.075 0.07277 2.657 0.0572 3.41025 0.01963
2.095 0.07184 2.6805 0.057 3.44 0.0192
2.115 0.07079 2.704 0.0568 3.465 0.01905
199
Annexes
• Cas Q=44.8 m3h-1
hauteur hauteur
hauteur hauteur hauteur hauteur hauteur
Position expérimentale expérimentale - Position Position Position Position
expérimentale calculée calculée calculée calculée
(m) + 3mm 3mm (m) (m) (m) (m)
(m) (m) (m) (m) (m)
(m) (m)
1.8 0.108 0.111 0.105 1.5 0.10453 2.135 0.09757 2.7275 0.079 3.49 0.02612
2 0.106 0.109 0.103 1.525 0.10453 2.155 0.09578 2.751 0.07846 3.515 0.02619
2.04 0.107 0.11 0.104 1.55 0.10453 2.175 0.09383 2.7745 0.07786 3.54 0.02627
2.08 0.105 0.108 0.102 1.575 0.10453 2.195 0.0919 2.798 0.07717 3.565 0.02634
2.16 0.099 0.102 0.096 1.6 0.10453 2.215 0.0901 2.8215 0.07628 3.59 0.02643
2.24 0.088 0.091 0.085 1.625 0.10453 2.235 0.08852 2.845 0.07294 3.615 0.02652
2.32 0.079 0.082 0.076 1.65 0.10453 2.255 0.0872 2.87475 0.06678 3.64 0.02662
2.4 0.072 0.075 0.069 1.675 0.10453 2.275 0.08617 2.9045 0.05958 3.665 0.02672
2.48 0.068 0.071 0.065 1.7 0.10453 2.295 0.0854 2.93425 0.05419 3.69 0.02682
2.56 0.068 0.071 0.065 1.725 0.10453 2.315 0.08483 2.964 0.05019 3.715 0.02691
2.64 0.068 0.071 0.065 1.75 0.10453 2.335 0.08435 2.99375 0.04696 3.74 0.027
2.72 0.066 0.069 0.063 1.775 0.10453 2.355 0.08385 3.0235 0.04426 3.765 0.02708
2.8 0.067 0.07 0.064 1.8 0.10453 2.375 0.08342 3.05325 0.04194 3.79 0.02717
2.88 0.062 0.065 0.059 1.825 0.10453 2.3985 0.08311 3.083 0.03993 3.815 0.02726
2.96 0.054 0.057 0.051 1.85 0.10453 2.422 0.08289 3.11275 0.03813 3.84 0.02735
3.04 0.045 0.048 0.042 1.875 0.10453 2.4455 0.08274 3.1425 0.03651 3.865 0.02744
3.12 0.036 0.039 0.033 1.9 0.10453 2.469 0.08259 3.17225 0.03504 3.89 0.02755
3.2 0.032 0.035 0.029 1.925 0.10453 2.4925 0.0824 3.202 0.03373 3.915 0.02765
3.28 0.027 0.03 0.024 1.95 0.10453 2.516 0.08218 3.23175 0.03256 3.94 0.02775
3.36 0.027 0.03 0.024 1.975 0.10453 2.5395 0.08193 3.2615 0.03147
3.44 0.027 0.03 0.024 1.995 0.10453 2.563 0.08167 3.29125 0.03035
3.52 0.028 0.031 0.025 2.015 0.10505 2.5865 0.08139 3.321 0.02938
3.6 0.025 0.028 0.022 2.035 0.10423 2.61 0.08109 3.35075 0.02857
2.055 0.103 2.6335 0.08076 3.3805 0.02774
2.075 0.10174 2.657 0.08038 3.41025 0.027
2.095 0.10048 2.6805 0.07995 3.44 0.02634
2.115 0.09913 2.704 0.07949 3.465 0.02605
200
Annexes
• Cas Q=57.3 m3h-1
hauteur hauteur
hauteur hauteur hauteur hauteur hauteur
Position expérimentale expérimentale - Position Position Position Position
expérimentale calculée calculée calculée calculée
(m) + 3mm 3mm (m) (m) (m) (m)
(m) (m) (m) (m) (m)
(m) (m)
1.8 0.128 0.131 0.125 1.5 0.127 2.135 0.11565 2.7275 0.09419 3.49 0.03058
2 0.13 0.133 0.127 1.525 0.12697 2.155 0.11353 2.751 0.09362 3.515 0.03067
2.04 0.127 0.13 0.124 1.55 0.12695 2.175 0.11136 2.7745 0.09295 3.54 0.03078
2.08 0.125 0.128 0.122 1.575 0.12693 2.195 0.10919 2.798 0.09216 3.565 0.0309
2.16 0.116 0.119 0.113 1.6 0.1269 2.215 0.10709 2.8215 0.09117 3.59 0.03101
2.24 0.105 0.108 0.102 1.625 0.12686 2.235 0.10514 2.845 0.08613 3.615 0.0311
2.32 0.093 0.096 0.09 1.65 0.12683 2.255 0.10345 2.87475 0.07721 3.64 0.03118
2.4 0.085 0.088 0.082 1.675 0.1268 2.275 0.10212 2.9045 0.06868 3.665 0.03126
2.48 0.083 0.086 0.08 1.7 0.12679 2.295 0.10117 2.93425 0.06301 3.69 0.03134
2.56 0.082 0.085 0.079 1.725 0.12678 2.315 0.10049 2.964 0.05863 3.715 0.03142
2.64 0.082 0.085 0.079 1.75 0.12678 2.335 0.09996 2.99375 0.05501 3.74 0.03151
2.72 0.081 0.084 0.078 1.775 0.12678 2.355 0.09945 3.0235 0.05192 3.765 0.03161
2.8 0.077 0.08 0.074 1.8 0.12678 2.375 0.09895 3.05325 0.04922 3.79 0.03172
2.88 0.074 0.077 0.071 1.825 0.12678 2.3985 0.09851 3.083 0.04685 3.815 0.03182
2.96 0.067 0.07 0.064 1.85 0.12676 2.422 0.09814 3.11275 0.04473 3.84 0.03192
3.04 0.054 0.057 0.051 1.875 0.12672 2.4455 0.09784 3.1425 0.04283 3.865 0.03201
3.12 0.044 0.047 0.041 1.9 0.1267 2.469 0.09755 3.17225 0.04114 3.89 0.03209
3.2 0.038 0.041 0.035 1.925 0.1267 2.4925 0.09727 3.202 0.03962 3.915 0.03217
3.28 0.033 0.036 0.03 1.95 0.12675 2.516 0.097 3.23175 0.03822 3.94 0.03225
3.36 0.03 0.033 0.027 1.975 0.12626 2.5395 0.09673 3.2615 0.03684
3.44 0.029 0.032 0.026 1.995 0.1251 2.563 0.09647 3.29125 0.03557
3.52 0.03 0.033 0.027 2.015 0.12486 2.5865 0.09623 3.321 0.03452
3.6 0.03 0.033 0.027 2.035 0.12381 2.61 0.096 3.35075 0.03355
2.055 0.12262 2.6335 0.09575 3.3805 0.03249
2.075 0.12117 2.657 0.09545 3.41025 0.03164
2.095 0.11951 2.6805 0.09511 3.44 0.03085
2.115 0.11765 2.704 0.09468 3.465 0.03051
201
Annexes
Annexe 8
Présentation du logiciel déversoir d’orage
Le logiciel déversoir d’orage à été conçu dans le but de répondre à la problématique du
déversement de l’eau du réseau d’assainissement vers la rivière par l’intermédiaire du déversoir
d’orage. L’objectif de cet outil numérique est le diagnostic des déversoirs d’orage. En
introduisant les caractéristiques géométriques de l’ouvrage (hauteur de crête, longueur,…) ainsi
que les caractéristiques des conduites amont et aval, le logiciel est capable de fournir le débit
déversé en fonction du débit amont, ainsi que la ligne d’eau le long de la crête déversante.
Jusqu’à présent, aucun outil de calcul n’était capable de décrire correctement le
comportement hydraulique de cet ouvrage dans toutes les configurations possibles. Comme
l’ont fait Robinson et McGhee [Robison-1993], la résolution numérique de ces équations est
toujours basée sur un algorithme qui décrit tous les cas possibles en fonction des régimes
d’écoulement et des conditions hydrauliques dans le déversoir. Il est donc nécessaire de
connaître ou d’émettre une hypothèse sur le régime d’écoulement avant calcul. Or, l’ouvrage
perturbe l’écoulement en raison de la variabilité du débit. Ainsi, faire l’hypothèse que le régime
d’écoulement sera toujours le même à l’amont ou à l’aval du déversoir peut s’avérer hasardeux.
De plus, l’une des particularités du déversoir d’orage est que le ressaut hydraulique y apparaît
fréquemment.
Bien sûr, les cas complexes de déversoirs possédant plusieurs entrées et/ou sorties, des
crêtes courbes ou encore par exemple un coude à l’amont ne peuvent pas être traités avec notre
logiciel. Les situations décrites précédemment doivent être traitées avec des outils plus
performant de nature 2D voire 3D tel que FLUENT.
203
Annexes
Déversoir latéral
Conduite
Conduite
Aval
Amont
Conduite
Déversée
Déversoir frontal
Conduite
Conduite Amont Déversée
Conduite
Aval
Conduite Aval
Conduite Amont
Angle d’entonnement
♦ haute :
Déversoir Crête déversante
204
Annexes
205
Annexes
♦ basse :
• La mise en charge des conduites amont et aval est prise en compte dans le calcul.
• Dans certains cas, les déversoirs sont équipés d’une vanne à l’amont du collecteur aval dans
le but de réguler le débit conservé. Celle-ci peut être prise en compte par le logiciel.
Vanne
Conduite Aval
Conduite Amont
Déversoir
Déversoir
Conduite Amont θ1
Conduite Aval
θ2
206
Annexes
a)
b)
2 INTERFACE UTILISATEUR
La figure ci-dessous présente la fenêtre d’accueil du logiciel déversoir d’orage :
207
Annexes
Cette étape permet de renseigner l’utilisateur sur les caractéristiques hydrauliques des
conduites amont et aval ; en particulier, on obtient les informations concernant le débit à pleine
section en régime permanent et uniforme et le régime d’écoulement.
De plus, des renseignements concernant les caractéristiques géométriques du déversoir
sont donnés ; comme le type de crête (haute ou basse) ou la longueur de crête (courte/longue).
208
Annexes
hn/hc
Fluvial (pente faible)
Q/Qmax
2.2 Calculs
209
Annexes
Fenêtre du logiciel déversoir permettant le choix du débit injecté pour le cas d’un unique
calcul
Voici, représenté sur la figure le type de résultats que l’on observe après convergence.
Sur la figure du haut sont représentés la ligne d’eau en rouge, le radier en noir et les
génératrices supérieures des collecteurs amont et aval et la crête pour le déversoir. La figure du
bas représente la courbe d’évolution du débit
210
Annexes
Il faut définir le débit maximal ainsi que le nombre de pas de calcul entre 0 et le débit maximal.
Le nombre de pas de calculs définit le nombre points à partir desquels sera construite la courbe
de fonctionnement.
Fenêtre du logiciel déversoir permettant le choix du débit injecté pour le cas d’une courbe
de fonctionnement
On représente ci-dessous le type de courbe de fonctionnement susceptibles d’être obtenues:
Qconservé
Début du Qinjecté
déversement
Courbe de fonctionnement obtenue grâce au logiciel déversoir
211
Annexes
L’ensemble des calculs effectués précédemment est consultable dans cette rubrique. On
peut choisir soit les résultats d’un calcul pour un débit déversé en fonction du débit amont, soit
une courbe de fonctionnement.
2.3 Enregistrement
L’ensemble des résultats est stocké dans un fichier .mat qui sert de rappel au programme
en cas de consultation des résultats ainsi que des caractéristiques du projet. Un fichier .txt est
également disponible dans le cas ou l’on veut consulter et traiter les résultats grâce à un autre
logiciel tel que Excel. Le fichier .txt est formaté pour être facilement consultable sur Excel. Il
suffit de choisir le séparateur « point-virgule » (;) lors de l’ouverture dans Excel.
212
Annexes
Annexe 9
Tableaux de valeurs concernant la validation du logiciel
déversoir dans le cas d’une crête déversante pour une
hauteur de crête w=7 cm.
• La configuration considérée dans laquelle la pente du collecteur amont est de 3 ‰ et des
pentes de 1 ‰ pour le déversoir et le canal aval permet de reproduire une ligne d’eau
caractéristique du régime fluvial pour un débit injecté en continu à l’amont du système de
24,3 m3/h. Le tableau ci-dessous regroupe les valeurs de hauteur d’eau simulées et
mesurées ainsi que les valeurs des débits injectés et conservés par le déversoir.
213
Annexes
214
Annexes
Annexe 10
Tableaux de valeurs concernant la validation du modèle
pour les déversoirs autres que w=7cm.
• Déversoir à crête simple w=5cm
215
Annexes
216
Annexes
217
Annexes
218
Annexes
219
Annexes
220
Annexes
221
Annexes
222
Annexes
Annexe 11
Présentation du logiciel confluence
Le logiciel confluence vise à reproduire le comportement du regroupement de deux
collecteurs en un seul. Ce modèle est très simplifié puisque la confluence est un ouvrage de
nature bidimensionnelle et dont le comportement, pour être précis, doit être décrit par un
logiciel de simulation 2D voire 3D.
¾ Longueur.
¾ Largeur.
¾ Pente du collecteur.
¾ Coefficients de Strickler.
¾ Pente de la conduite.
¾ Nombre de pas d’espace (nombre de cellules de calcul) :
10 → valeur moyenne
30 → grande précision sur la ligne d’eau
5 → Rapidité de calcul mais faible précision
¾ Angle par rapport au collecteur aval.
¾ Débit injecté dans les collecteurs amont 1 et amont 2.
223
Annexes
Les résultats que l’on obtient grâce à ce logiciel, pour la configuration présentée sur la
fenêtre d’accueil, sont représentés ci-dessous. Les tronçons amont sont tous deux à pente
fortes et le collecteur aval à pente faible. On observe après convergence la formation de
ressauts hydrauliques dans les canaux amonts.
Cette fenêtre visualise les lignes d’eau (en rouge) obtenues dans les deux collecteurs
amont et le collecteur aval. On représente en plus les hauteurs normales (en jaune) et critiques
(en vert) pour chaque tronçon.
224
Annexes
Cette figure permet de visualiser l’évolution du débit dans chaque canal. On observe en
s’approchant de la confluence une légère augmentation du débit. Ceci est du d’une part à la
présence des ressauts hydrauliques dans les canaux amonts car dès qu’il y a ressaut
hydraulique on observe un léger pic de débit engendré par le schéma numérique. Le schéma
numérique est tel que lorsqu’il y a présence d’une discontinuité, celui-ci joue sur le débit pour
assurer la convergence. D’autre part, cette erreur sur les débits est sans doute due à une
mauvaise gestion des conditions aux limites au niveau de la confluence.
225