Équations de Saint Venant et Écoulements
Équations de Saint Venant et Écoulements
Résumé
Cours salle M19 CESTIA, Abidjan L 24/09/18 au V28/09/18 + L01/10/18 08 :00-12 :00
planning du cours MU4MEF04 2019-2020 : 55-65 101 13 :45
[Link]
cours, notes... [Link]
voir page Moodle pour l’année 2020-2021 [Link] voir page Moodle pour l’année 2021-2022
[Link]
2021-2022 :
MU4MEF04 - Ondes et Ecoulements en milieu naturel 55-66 101 Wednesday, January 19 8 :30-10 :30am Weekly on Wednesday, until Apr 6, 2022
1 Introduction
1.1 Généralités
Source de toute vie sur Terre, l’eau est un enjeu évident pour l’humanité.
La compréhension des écoulements en milieu naturel est un enjeu fondamental pour la société humaine.
Enjeu scientifique (modélisation, simplification pour compréhension des phénomènes et du monde qui nous entoure), industriel (par l’énergie que l’on en tire depuis
les anciens moulin à marée, aux actuelles usines marémotrices et aux fermes à courant en cours, ainsi que toute l’énergie des barrages, appelée houille blanche dans
la région de Grenoble...) et enfin et surtout humain car la majeure partie des humains vivent le long des fleuves ou des côtes et que l’eau est indispensable à la
vie. (70 % des côtes sont en érosion, 80 % de la population mondiale habite à basse altitude et plus de 20 % à proximité d’un océan ou d’un estuaire. (source
[Link]
1
Saint-Venant
De nombreux lieux sont en péril par ensablement, sédimentation : Mont St Michel, Lagune Venise, Lagune Cocody (Afrique de l’Ouest), ensablement des ports,
comblement par sédimentation des barrages, grands deltas
érosion côtière, reculée des terres...
Au contraire on veut faire des Polders pour augmenter la surface habitable : Hollande Danemark, ou les alluvions du Nil sont à l’origine de la civilisation des Pharaons
Ce sont des problèmes d’aménagement du territoire BRGM ([Link] et de nombreuses sociétés de service, Veolia, Eau de Paris, Artelia, etc ;
Comme déjà signalé, on tire de l’énergie de l’eau, la parties liées à l’énergie intéressent EDF (arrivées d’eaux aux centrales nucléaires, hydrauliques) et les autres
acteurs de l’énergie plus ou moins verte...
Les écoulements ne sont pas toujours maı̂trisés : inondation par débordement de fleuves, identification de zones inondables,
Tsunamis, ressaut, mascaret....
la pollution peut s’y déverser, prolifération d’algues
prévention, simulation des accidents
Codes EDF utilisés pour la protection de l’environnement, le calcul des tempêtes, de marées des courants de marées, pour estimer les contraintes exercées sur les
ouvrages (ponts barrages, digues), aide à la navigation (prédiction des marées).
• inondations
- Un problème antédiluvien : l’épopée de Gilgamesh dont les racines viennent d’il y a 3200 ans [Link] L’Épopée
raconte plusieurs exploits du roi Gilgamesh, c’est le plus vieux récit écrit de l’humanité, il a été retrouvé à Ninive sur des tablettes transcrites vers -700. L’histoire du
déluge et du navire de Uta-napishti qui a sauvé quelques hommes et les animaux y tient une grande place à la fin de l’épopée.
- En 709, soixante ans avant l’avènement de Charlemagne, un coup de mer a détaché Jersey de la France. D’autres sommets des terres antérieurement submergées
sont, comme Jersey, visibles. Ces pointes qui sortent de l’eau, sont des ı̂les. C’est ce qu’on nomme l’archipel normand. Victor Hugo ([Link]
wiki/Mythe_du_raz-de-mar%C3%A9e_de_mars_709)
- Inondations de Grenoble 1219 et 1859
- Mac Mahon : dans la nuit du 23 au 24 juin 1875, une importante crue de la Garonne se produit. Visitant des villes et des villages dévastés, ne sachant que dire,
il déclara le célèbre que d’eau. . . que d’eau !. . . Le préfet du département lui répondit alors : Et encore, Monsieur le Président, vous n’en voyez que le des-
sus. . . ! [Link]
- Grande inondation de Paris 1910
- ...
- Les inondations en 2013 ont fait environ 10.000 morts
- Les inondations en 2013 ont couté 50.000.000 Euros en dégâts.
- inondation Grande Bretagne décembre 2015
- Cannes inondation octobre 2015 ([Link]
- Brague inondation octobre 2015, 250 mm de pluie en 3 heures, 20 victimes 600ME dégats assurés
- 22/01/2018 alerte inondation, 1m d’eau à Ornans (Doubs), niveau de la Seine et du Rhin inquiétant, le Zouave a les pieds dans l’eau
- mercredi 24 janvier 2018, la présence d’eau a été constatée dans des galeries techniques situées au 2ème sous-sol du campus Pierre et Marie Curie. Des pompages y
- ESV . 2-
Saint-Venant
sont réalisés.
- Inondations meurtrières dans la nuit du 18 au 19 juin 2018 à Cocody Abidjan, environ 20 morts.
- Inondation meurtrières dans l’Aude nuit du 14 octobre 2018
- 8 avril 2019 inondations en Iran.
- 23 octobre 2019 inondation dans l’Hérault, 3 mois en 48h, inondation en Catalogne et en iIalie la veille
- 11 mars 2020 Apilly : commune inondée depuis 1 mois ; le barrage de Manicamp déverse l’eau de l’Oise sur cette commune.
- Pluies fortes en Afrique de l’Est, indondation et coulée de Boues, 120 morts au Kenya le 06/12/19 (Le Monde)
- 2020 inondations au japon pendant le confinement
- 2020 inondations en Chine
- 2020 ”La Chine est confrontée à des crues exceptionnelles, Plus de 140 personnes sont mortes (région de Wuhan), L’état de certains barrages inquiète. 20 millions d
personnes affectées, pertes 7.76 Milliards d’Euros. Barrage des trois Gorges Le monde 12-13/07/20
- Les crues du 2 octobre 2020 ont littéralement dévasté certaines parties des vallées de la Vésubie et de la Roya dans l’arrière-pays niçois. Des travaux colossaux de
reconstruction des routes, de ponts, des réseaux d’eau et d’électricité, sont nécessaires désormais pour un coût d’au moins un milliard d’euros (France Bleue). - Crues
2021 a Garonne a atteint 10,22 m à Marmande.
- février 2021 Rupture d’un glacier dans l’Himalaya, au moins dix-huit morts et deux cents disparus : La rupture d’une partie d’un glacier dans le nord de l’Inde a pro-
voqué dimanche des flots torrentiels qui ont tout emporté sur leur passage, y compris un barrage hydroélectrique. [Link]
- Journée de l’eau 22 mars 2021
- 15 juillet 2021 At least 20 people dead in destructive flooding that hit Germany and Belgium this week Heavy rainfall across parts of western Europe has resulted
in destructive and unfortunately [Link]
- Crues 14 et 15 juillet 2021, en Allemagne et Europe de l’Ouest, la crue est due à des précipitations record pour la saison. Une des pires catastrophes naturelles du
début du siècle en nombre de victimes.
- novembre d’importantes inondations ont lieu den Colombie-Britannique (Canada), provoquant des glissements de terrain. Plusieurs personnes sont portées disparues.
- Inondation en Australie mars 22 [Link]
php
• Tsunami
- 1 novembre 1755 tremblement de terre et raz de marée à Lisbonne : ”À peine ont-ils mis le pied dans la ville, en pleurant la mort de leur bienfaiteur, qu’ils sentent la
terre trembler sous leurs pas ; la mer s’élève en bouillonnant dans le port, et brise les vaisseaux qui sont à l’ancre. Des tourbillons de flammes et de cendres couvrent
les rues et les places publiques ; les maisons s’écroulent, les toits sont renversés sur les fondements, et les fondements se dispersent ; trente mille habitants de tout âge
et de tout sexe sont écrasés sous des ruines.” Candide 1759, Voltaire
- 1958 Mégatsunami de 1958 de la baie Lituya [Link] (mettre un accent ”é” à Méga
pour aller sur wiki)
- [Link] 22 mai 1960 Chili
- tsunami 26/12/2004
- tsunami Fukushima 11/03/2011
- 2015 The Tyndall Glacier landslide occurred at about 8 :19 pm local time on 17th October 2015 at the toe of the Tyndall Glacier. The landslide is enormous
– Stark and Ekstrom estimate from the seismic data that it had a mass of about 180 million tonnes, which would give a volume in the order of 72 million cubic
metres. The landslide flowed into Taan Fjord, triggering a localised tsunami that was detected 155 km away [Link]
tyndall-glacier-landslide-1/
- Tsunami aux Célèbes sept 2018.
- Tsunami aux Iles Tonga janvier 2022 [Link]
- Suivi des tremblements de terre [Link] générateurs
possibles
• ruptures de barrage
- ESV . 3-
Saint-Venant
- 02/12/1959 rupture du barrage EDF de Malpasset (Var) 423 victimes. La rupture du barrage de Malpasset en 1959 a été très documentée à l’époque (les gendarmes
ont mesuré l’extension maximale, des transformateurs EDF en disjonctant ont donné l’heure de passage de la vague), cette documentation en fait un cas de simulation
à re simuler.
- rupture du barrage de Vajont, 09/10/1963, 3 ans après sa mise en eau, 1900 morts dans la vallée du Piave.
- La rupture du barrage de Banqiao (Chine du nord), 1975, aurait causé plus de victimes que n’importe quelle autre rupture de barrage.
- le barrage d’Oroville aux US qui risquait de se déverser en février 2017..[Link] à no-
ter un an après en Colombie : un barrage hydroélectrique en construction menace de céder : 25 000 personnes évacuées 200 000 habitants en danger https:
//[Link]/meteo/inondations/alerte-rouge-en-colombie-un-barrage-menace-de-ceder_2763047.html.
- Le plus haut barrage des Etats-Unis, situé dans la région d’Oroville et plein à ras bord après des semaines de fortes pluies, menace de céder. fv 2017
- 25/01/19 rupture du barrage de Brumadinho 150 morts 200 disparus.
• glissements de terrain
- Glissements de terrain et de boue suite à des pluies (Rio Brésil 1988)
- Glissement de terrain En Éthiopie, au moins 113 morts dans l’éboulement d’une décharge.
- Un glissement de terrain a fait au moins 113 victimes à Addis-Abeba samedi 11 mars 2017.
- 15 avr. 2017 - Au moins quinze personnes, dont quatre enfants, ont été tuées, vendredi 14 avril, au Sri Lanka, dans l’effondrement d’une montagne d’ordures
([Link] Monde > Asie-Pacifique), ...
- Pluies fortes en Afrique de l’Est, inondation et coulée de Boues, 120 morts au Kenya le 06/12/19 (Le Monde)
- Suite aux incendies de Los Angeles, les pluies ont provoqué des coulées de boue Janvier 2018, 13 morts 47 disparus
• érosion
- 10 mai 2017 - Depuis 1984, la plage de Dooagh (Irlande) n’était que rochers. Le sable est revenu en quelques semaines, comme en atteste la photo ... https:
//[Link]/environnement/[Link]
- ”Pacifica” érosion de la côte aux US, recul en 20 ans, 28 janvier 2016 Le Monde [Link]
4855511_3244.html
- ”The frontline of climate change” The Norfolk village falling into the sea January 30, 2021 [Link]
- 41◦ 51’25” N 1◦ 13’17” O Construit à 200m de la mer en 1965 à Soulac-sur-Mer l’immeuble ”le Signal” est devenu inhabitable à cause de l’érosion du littoral, 21 Jnv
2021 les occupants reçoivent un dédommagement de 70 ◦ /o(Le Monde 31/01/21)
- ESV . 4-
Saint-Venant
- ESV . 5-
Saint-Venant
3 Plan approximatif
Nous allons poser les équations de Saint-Venant et montrer d’abord comment elle s’écrivent en supposant un écoulement constant par tranche, puis en faisant un
bilan sur une tranche à la section §5. C’est la manière la plus simple, elle a de plus l’intérêt d’être la forme utile pour une résolution numérique avec une approche
volumes finis (et c’est la démarche historique de Saint-Venant lui même). Puis nous prenons le temps de réétablir ces équations en partant des Equations de Navier
Stokes complètes écrites en couche mince. Cela fait ressortir l’équilibre hydrostatique. En intégrant sur l’épaisseur cela fait apparaı̂tre le frottement au fond et les flux
de masse et de quantité de mouvement. On discute quelques cas d’écoulement stationnaires invariants §6.7.1 et §6.7.2. Nous discutons ensuite suivant le profil choisi
quel est le lien entre le frottement et le débit et la hauteur dans §7, c’est là qu’apparaissent Chézy et Manning dont les noms sont associés aux hypothèses sur le frottement.
Ayant fixé le système et compris d’où il venait, nous pouvons attaquer sa résolution dans des cas particuliers dans §II. D’abord sur des échelles très très longues en
x : les deux équations se réduisent à une équation : notamment le cas de l’onde de crue §8.5 pour lequel on introduit la notion de transport et de caractéristique. Puis,
avant que ces notions ne soient reprises pour les Equations de Saint-Venant proprement dites : c’est le cas de la rupture de barrage, puis de la formation de ressaut,
nous regardons les cas linéarisés à faible vitesse et faible variation de hauteur d’eau (cas des marées par exemple). Enfin, l’écoulement sur un obstacle est aussi abordé
dans les cas linéaires et non linéaires.
- ESV . 6-
Saint-Venant
Un chapitre spécial est consacré aux résolutions en volumes finis en C (adaptable facilement en python). Une fois introduites les notions de base, des exemples
réalistes de calculs de marée, de tsunamis sont ensuite présentés en utilisant Basilisk.
Ce cours s’inscrit dans le triptyque Théorie, Expérience & Simulation. Ce fichier est la partie théorique où l’on essaye de comprendre la modélisation et les équations,
ma résolution numérique est présentée et des exemples pratiques seront étudiés avec le code Basilisk. Enfin, une cuve à vague à laquelle nous aurons accès permet
d’observer in vitro les écoulements abordés.
Bibliographie
Ces notions sont dans différents ouvrages de mes estimés collègues Olivier Thual (Hydrodynamique de l’environnement), Hubert Chanson (The Hydraulics of Open
Channel Flow) Christophe Ancey (Ondes de crue et de rupture de barrage [Link]
Barrages/[Link]) et bien d’autres que l’on trouve en bibliothèque ou sur le ouaibe.
- ESV . 7-
Saint-Venant
Figure 1 – Quart haut droit : un ressaut dans la rivière le Frémur causé par la montée de marée dans la baie de la Fesnaye, modélisé par Saint-Venant écrit sous forme
compacte (haut gauche), puis sous forme discrétisée, elle même codée en bas à gauche. Le résultat numérique (simulation numérique) et en bas à droite la simulation expérimentale
(modélisation expérimentale) dans une cuve à vagues de 5 m de long à Sorbonne-U (salle Galillée).
- ESV . 8-
Saint-Venant
Première partie
Le Fluide : dérivation des équations de Saint-Venant
4 Introduction
Nous allons établir les équations de Saint-Venant et présenter quelques exemples d’applications dans le cas d’un écoulement d’eau sur une surface donnée. La clef
principale est que la distribution de pression est reliée principalement aux variations de la surface de l’eau par la relation d’équilibre hydrostatique et que l’écoulement
est en couche mince..
Figure 2 – A gauche état de repos, pas de vitesse, le fond est plat, la pression varie avec la hauteur z, la pression est P0 dans l’air au dessus, la pression varie de manière
hydrostatique dans l’eau : p(z) = P0 + ρg(h0 − z) pour 0 < z < h0 . A droite, on suppose une petite perturbation de vitesse, de forme de fond (f ) et/ou de surface libre (η), la
hauteur d’eau étant h = η − f , la pression est supposée rester hydrostatique dans l’eau : p(z) = P0 + ρg(η − z) pour f < z < f + h = η.
traduisant la forme du fond de la rivière. Le fond varie lentement en x, y dans la description de Saint-Venant. Nous noterons aussi z = Zb (x, y) la forme du fond,
”b” pour bottom. Les deux notations sont équivalentes (et seront employées indifféremment). Le fond est fixe dans la théorie de base. Dans les cas d’érosion, il varie
très lentement par rapport au temps caractéristique du mouvement de l’eau. On définit z = η(x, y, t) comme étant la cote de la surface libre. La hauteur d’eau est la
distance du fond à la surface libre : h(x, y, t) = η(x, y, t) − f (x, y). Au dessus de la surface libre, l’air assure une pression atmosphérique P0 que l’on suppose constante
et uniforme (on néglige la densité de l’air (1.2kgm−3 par rapport à celle de l’eau 103 kgm−3 ). On néglige donc complètement les mouvements de l’air supérieur. On a
- ESV . 9-
Saint-Venant
Cette formulation hydrostatique est une sévère limitation, en effet, les vagues, par exemple, ne pourront pas être décrites das le cadre de Saint-Venant (voir
[Link] Nous corrigerons en §12 cette approximation hydrostatique.
∂h ∂(hu) ∂(hv)
+ + = 0
∂t ∂x ∂y
∂(hu) ∂(hu2 ) ∂(huv) ∂η τx
+ + = −hg − + Fx , (1)
∂t ∂x ∂y ∂x ρ
∂(hv) ∂(huv) ∂(hv 2 ) ∂η τy
+ + = −hg − + Fy .
∂t ∂x ∂y ∂y ρ
τx et τy sont les contraintes de frottement au fond, Fx et Fy sont des forces (volumiques) comme l’accélération de Coriolis ou une force d’entraı̂nement due au vent.
Toutes ces forces seront explicitées, discutées, détaillées plus loin paragraphes 7.2.3, 7.3.
On remarque que dans cette partie le fond est fixe f (x, y), c’est une donnée du problème, on cherche à établir comment l’eau s’écoule sur une topographie donnée.
Les conditions aux limites sont spéciales, on va supposer que le bord du domaine (la berge) est un talus très abrupt, et on ne va pas se soucier des vaguelettes
qui mouillent la berge. Il y a ”glissement au bord” (ce bord étant les berges). Ensuite, la vitesse est imposée à l’entrée (c’est le flux d’entrée), puis on se donne une
condition de sortie libre.
∂h ∂(hu)
+ =0
∂t ∂x
- ESV . 10-
Saint-Venant
∂(hu) ∂(hu2 ) ∂η τx
+ = −hg − + Fx ,
∂t ∂x ∂x ρ
Les équations seront démontrées par la suite à partir des équations de Navier Stokes. La démarche sera la plus rigoureuse possible en tenant compte d’hypothèses
η(x, t) ∆x
h −
→ u(x) u(x + ∆x)
g
f (x) f (x)
Figure 3 – volume de contrôle si on veut établir des équations 1D.
simplificatrices successives. On peut cependant commencer la démarche en partant d’une vision simple des écoulements sans faire appel à Navier Stokes. La manière
simple de présenter les équations est de se donner le dessin de la figure 3. On a une tranche de fluide d’épaisseur ∆x, les bornes x et +∆x de cette tranche sont fixes,
dans cette tranche, il circule un écoulement de vitesse constante par tranche, u(x, t). Avec cette hypothèse on met en défaut le non glissement à la paroi (effet de la
viscosité). L’effet de la viscosité est cependant réintroduit dans l’équation de quantité de mouvement par le frottement à la paroi. On dit que u(x, t) est une vitesse
moyenne par tranche.
• Il entre un flux h(x, t)u(x, t) en x et il sort un flux en x + ∆x qui est h(x + ∆x, t)u(x + ∆x, t), mais il ne faut pas oublier qu’il y a un flux dû au déplacement de la
surface libre vers le haut : ∆x∂η/∂t, donc, puisque f ne dépend pas du temps, le bilan de masse s’écrit :
∂
(h∆x) = h(x, t)u(x, t) − h(x + ∆x, t)u(x + ∆x, t).
∂t
On peut dire aussi que la variation de volume est égale au flux qui rentre moins celui qui sort. Insistons sur le fait que cette écriture, que l’on peut réécrire avec le flux
Q = hu est très générale, elle traduit une loi de conservation, elle est le fondement de la méthode des volumes finis qui va nous servir pour la résolution numérique :
∂
(h∆x) = Q(x, t) − Q(x + ∆x, t)
∂t
en faisant apparaı̂tre la dérivée u(x + ∆x, t)h(x + ∆x, t) = u(x, t)h(x, t) + ∆x ∂uh
∂x + ... on retrouve la relation de flux que l’on a annoncée :
∂h ∂(hu)
=− .
∂t ∂x
• Passons à la quantité de mouvement, on a une tranche de fluide d’épaisseur ∆x, il y circule un écoulement de vitesse constante en z qui est u(x, t) et la surpression
p(x, z, t) = ρg(η − z) qui varie du haut au bas de la tranche. Le bilan de quantité de mouvement est fait dans la tranche. La masse dans la tranche ρh∆x, l’accélération
d ∂ ∂
de la tranche si elle bouge dans son mouvement est du dt , (la dérivée totale est bien notée dt = ∂t + u ∂x ). Mais ici on a une tranche d’épaisseur ∆x, les bornes x et
x + ∆x de cette tranche sont fixées, la tranche est soumise à des forces de pression que l’on écrit dans le cadre hydrostatique précédent. Les forces de pression poussent
vers la droite en x, soit sur toute la hauteur h = η − f on a :
Z η Z η
(η − z)2 η (η − f )2 h2
p(x, z, t)dz = ρg(η − z)dz = −ρg f
= 0 − (−ρg ) = ρg
f f 2 2 2
- ESV . 11-
Saint-Venant
et poussent vers la gauche en x + ∆x. Le bilan des forces de pression est donc :
η(x,t) η(x+∆x,t)
h(x, t)2 h(x + ∆x, t)2
Z Z
p(x, z, t)dz − p(x + ∆x, z, t)dz = ρg( − )
f (x) f (x+∆x) 2 2
∂ h2
= −∆xρg + ...
∂x 2
Attention à ce stade à ne pas omettre un terme un peu subtil de poids dû à la pente faible du fond f 0 , qui projette la pesanteur en une très faible force le long de
l’écoulement :
df
−ρgh(∆x) .
dx
La topographie f varie très lentement, ce terme est donc très très petit. Ce terme introduit ainsi, il y a une petite contradiction car les coordonnées sont bien telles
que z est suivant la gravité →−
g (en fait il faudrait légèrement corriger g), mais ce terme traduit la non connaissance exacte de l’horizontale car la pente est très très
faible. Une fois réintroduit nous allons voir que ce terme préserve l’équilibre du lac. Enfin, il ne faut pas oublier la contrainte de frottement visqueux τx qui agit sur
la surface en z = f (x), et qui dépend de u, sa résultante sera −τx ∆x, en mettant tous ces termes bout à bout et en factorisant ∆x on trouve la résultante globale
(pression + pente + frottement) intégrale qui s’exerce sur la tranche :
∂ h2 df
F = (−ρg ( ) − ρgh − τx )∆x.
∂x 2 dx
Pour calculer l’accélération de la tranche il faut faire un bilan de la variation de quantité de mouvement dans la tranche (∆x) ∂(ρuh) ∂t . Il faut tenir compte du flux de
quantité de mouvement qui sort à droite −ρu(x + ∆x, t)2 h(x + ∆x, t) et de celui qui rentre à gauche +ρu(x, t)2 h(x, t) (en supposant encore que le profil est plat, sinon,
on ne pourrait pas simplement multiplier hu et u). Ce bilan de flux suit la même démarche que pour la masse. On trouve donc un bilan en écrivant le développement
2
de Taylor : ρu(x, t)2 h(x, t) − ρu(x + ∆x, t)2 h(x + ∆x, t) = −(∆x)[ ∂(ρu ∂x
h)
+ ...]. Dans le cas de la masse, il n’y a pas de terme source, ici il y en a un c’est F que l’on
vient de calculer (pression, pente et forttement), le bilan :
∂(ρuh) ∂(ρu2 h)
(∆x) = −(∆x)[ + ...] + F.
∂t ∂x
L’équation totale s’écrit bien alors
∂(ρuh) ∂(ρu2 h) ∂ h2 df
+ = −ρg ( ) − ρgh − τx .
∂t ∂x ∂x 2 dx
∂(ρuh) ∂(ρu2 h)
On peut aussi l’écrire avec d
dt = ∂
∂t
∂
+ u ∂x et compte tenu de la conservation de la masse : ∂t + ∂x = ( ∂(ρu) ∂(ρu) ∂h
∂t + u ∂x ) + ρu( ∂t +
∂(hu)
∂x ) = ( ∂(ρu) ∂(ρu)
∂t + u ∂x ) + 0
on a ainsi car h + f = η :
du ∂ h2 df ∂u ∂u ∂η τx
ρh = −ρg ( ) − ρgh − τx , donc ρ( + u ) = −ρg −
dt ∂x 2 dx ∂t ∂x ∂x h
∂p ∂η
on aurait pu directement écrire cette forme à partir de l’expression de la pression − ∂x = −ρg ∂x , sans passer par la moyenne sur la tranche.
∂η
Remarquons que ces équations vérifient ”l’équilibre du lac”, à savoir, si la surface de l’eau est à hauteur constante : η = f + h,i.e. ∂x = 0 alors le fluide est bien au
repos u = 0 et η = 0. Il faudra toujours bien faire attention à ce problème d’horizontalité qui est récurrent. Les différents utilisateurs négligent souvent la variation de
g avec l’angle.
Remarquons enfin que c’est quasiment la présentation originale d’Adhémar Jean-Claude Barré de Saint-Venant lui même en 1871 au paragraphe ”théorie et
équations générales du mouvement non permanent des eaux courantes”. Il termine en écrivant ”Je donnerai dans un autre article l’intégration de ces équations pour
un cas étendu”. Le fichier est ici [Link]
- ESV . 12-
Saint-Venant
∂h ∂(hu)
+ = 0
∂t ∂x
∂(hu) ∂ g df τx (2)
+ (hu2 + h2 ) = −gh − + Fx ,
∂t ∂x 2 dx ρ
h = η − f.
En première lecture on saute la suite (qui correspond à l’établissement de ces équations à partir de Navier Stokes) §7 et on va directement à la section ”Hypothèses
pour obtenir Saint-Venant” à la sous section ”Hypothèses de fermeture, choix d’un profil” ; on y discute le fait que la vitesse n’est peut être pas constante
dans l’épaisseur mais que l’on peut considérer effectivement la vitesse moyenne, quitte à rajouter un coefficient correctif tenant compte de la forme du profil de vitesse.
De même on y estime le frottement (qui lui aussi est fortement dépendant de la forme exacte du profil de vitesse).
Sans trop divulgâcher la suite (§7.2.6), on pose Fx = Cg2 u2 avec un coefficient de frottement de ”Chézy” coefficient empirique Cc , nous verrons que les équations
c
de Saint-Venant admettent plusieurs dégénérescences suivant les cas : bosse f courte ou longue et petite perturbation, frottement ρ Cg2 u2 faible (§9) ou fort... toute
c
cette zoologie sera examinée après dans les cas ”longue échelle”, que nous l’appellerons aussi ”Onde d’inondation”ou ”flood wave ” (§8.5) dans les cas d’échelles
intermédiaires nous l’appellerons aussi ”Onde de crue diffusante” (§8.6), etc etc. Tout ceci sera vu dans la suite.
- ESV . 13-
Saint-Venant
- ESV . 14-
Saint-Venant
Quand on regarde un fleuve, il est clair que la vitesse transverse U2 (vers le haut ou le bas) est faible par rapport à la vitesse longitudinale U1 . Il est aussi clair
que les variations longitudinales sont sur des échelles plus longues que la hauteur d’eau dont la valeur caractéristique est disons h0 . Si U0 est l’ordre de grandeur de la
vitesse et ε le rapport entre la profondeur et la longueur de variation longitudinale, alors un bon choix pour U1 = U0 Ū1 , U2 = εU0 Ū2 , z = z̄h0 et x = x̄h0 /ε. Ce choix
est excellent pour que l’équation d’incompressibilité reste entière :
∂U1 ∂U2 U0 ∂ Ū1 εU0 ∂ Ū2 ∂ Ū1 ∂ Ū2
+ = 0. devient + = 0 puis + = 0.
∂x ∂z h0 /ε ∂ x̄ h0 ∂ z̄ ∂ x̄ ∂ z̄
On dit aussi que le choix de ces échelles est guidé par ”moindre dégénérescence” : on garde le maximum de termes dans l’équation. Dans la dérivée totale, on mesure
∂ ∂ ∂
le temps t = τ t̄, et on note que ( + U1 + U2 ) est donc
∂t ∂x ∂z
U0 h0 ∂ ∂ ∂
( )(( ) + Ū1 + Ū2 )
h0 /ε U0 τ ε ∂ t̄ ∂ x̄ ∂ z̄
tous les termes restent aussi si on mesure le temps avec la longueur longitudinale h0 /ε/U0 . Remarquons au passage que l’on peut définir un Strouhal St = ( Uh0 0τ ε ). Si
on mesure le temps avec h0 /U0 alors St = 1/ε, mais si on mesure le temps avec la longueur longitudinale, St = 1 et on plus de termes (moindre dégénérescence), la
dérivée totale est entière :
εU0 ∂ ∂ ∂
( )( + Ū1 + Ū2 )
h0 ∂ t̄ ∂ x̄ ∂ z̄
Dans l’équation de conservation de quantité de mouvement longitudinale, dans le cas Newtonien laminaire
∂τ11 ∂τ12 ∂ 2 U1 ∂ 2 U1
+ =µ + µ ,
∂x ∂z ∂x2 ∂z 2
2 2
il est clair que l’on peut négliger le terme µ ∂∂xU21 par rapport à µ ∂∂zU21 (le premier est d’ordre ε2 par rapport au second). Il ne reste que :
De même substituons les mêmes échelles dans l’équation suivant U2 , le terme de dérivée longitudinale de contrainte est négligé :
∂ Ū2 ∂ Ū2 ∂ Ū2 (δP ) ∂ p̄ gh0 ((µU0 /h0 )) ∂ 2 Ū2
( + Ū1 + Ū2 )=− 2 2 − 2 + + ...
∂ t̄ ∂ x̄ ∂ z̄ ρU0 ε ∂ z̄ ρU0 ε2 ρU02 ε3 ∂ z̄ 2
U0
il paraı̂t naturel de choisir (δP ) = ρU02 (en fait on peut aussi prendre (δP ) = ρgh0 = ρU02 /F 2 le nombre de Froude F = √ est nous le verrons d’ordre 1). On récrit
gh0
l’équation en multipliant par ε2
∂ Ū2 ∂ Ū2 ∂ Ū2 ∂ p̄ 1 ν ∂ 2 Ū2
ε2 ( + Ū1 + Ū2 )=− − 2 + (ε−1 ) + ...
∂ t̄ ∂ x̄ ∂ z̄ ∂ z̄ F U0 h ∂ z̄ 2
- ESV . 15-
Saint-Venant
Le terme convectif est alors négligeable, tant que ε2 << 1. Dans le cas laminaire, le terme visqueux est aussi négligeable si Uν0 h est assez petit et (ε−1 ) pas trop
grand. Il n’est pas très petit si la viscosité est forte (on fera attention dans le cas turbulent il peut y avoir une contribution à ce terme dans la pression, de même dans
les cas non newtoniens, ce terme peut revenir dans les équations). En pratique, si le Reynolds Re = U0 h/ν est assez grand, il ne reste que
∂ p̄ 1
0=− − 2
∂ z̄ F
il n’y a pas de terme de pression hors le terme hydrostatique.
µ
Les équations s’écrivent maintenant puisque l’ordre de grandeur de ρU0 hε est aussi 1/(εRe) qui peut être quelconque.. :
• Transformation de Prandtl
On vient de voir que l’on garde x̄ mais que l’on change z̃ = z̄ − f¯, par la dérivation en chaı̂ne :
∂ ∂ ∂ ∂ ∂ f¯ ∂
= et = −
∂ z̄ ∂ z̃ ∂ x̄ ∂ x̄ ∂ x̄ ∂ z̃
ainsi
∂ Ū1 ∂ Ū2 ∂ Ū1 ∂ f¯ ∂ Ū1 ∂ Ū2
+ = 0. devient − + = 0.
∂ x̄ ∂ z̄ ∂ x̄ ∂ x̄ ∂ z̃ ∂ z̃
- ESV . 16-
Saint-Venant
ou encore
∂ Ū1 ∂ ∂ f¯ ∂ Ū1 ∂ Ũ2
+ (Ū2 − Ū1 ) = 0 devient + = 0.
∂ x̄ ∂ z̃ ∂ x̄ ∂ x̄ ∂ z̃
¯
en définissant une vitesse transformée Ũ2 = Ū2 − ∂∂ fx̄ Ū1 associée à la transformation de Prandlt, z̃ = z̄ − f¯. On constate alors que
∂f
z → z − f, et U2 → U2 − U1
∂x
• conservation de la masse
∂U1 ∂U2
+ =0
∂x ∂z
au fond, U1 = U2 = 0, l’interface, définie de manière implicite par la fonction F (x, z) = (h(x, t) − z). Cette interface (transformée par Prandtl) est conservée dans le
mouvement du fluide donc la valeur de la vitesse transverse à l’interface en fonction du déplacement de l’interface :
∂h ∂h
U2 (x, h, t) = + U1 .
∂t ∂x
- ESV . 17-
Saint-Venant
• conservation de la quantité de mouvement écrite sous forme conservative (compte tenu de l’incompressibilité et toujours de la transformation de Prandtl)
avec τ12 (h) = 0 pas de contrainte en haut. Dans le cas laminaire que nous avons vu τ12 = µ ∂U
∂z .
1
Remarque :
Ces équations pourraient se résoudre numériquement avec un schéma numérique adapté. Bien entendu, cette résolution dépasse le cadre de ce cours. Elle passe
par la méthode dite ”multilayer” de Audusse et Sainte Marie codée dans Basilisk (cf F. De Vita, P.-Y. Lagrée, S. Chibbaro, S. Popinet ”Beyond Shallow Water :
appraisal of a numerical approach to hydraulic jumps based upon the Boundary Layer Theory.” European Journal of Mechanics / B Fluids 79 (2020) 233-246)
[Link]
Suite :
Nous allons par la suite intégrer sur l’épaisseur h ces équations, et nous allons voir que l’on obtient les équations de Saint-Venant. Auparavant, nous avons besoin
de rappeler la dérivation des intégrales puisque que l’on va intégrer ces équations sur une hauteur qui varie.
- ESV . 18-
Saint-Venant
6.4 Rappel sur les dérivées d’intégrales, lien avec la formule de la divergence
Pour l’équation de la masse, nous avons besoin de dérivées en x, mais ensuite, pour la quantité de mouvement nous dérivons en temps t.
Z x2 Z x2
∂t gdx = ∂t g(x, t)dx + g(x2 , t)∂t x2 − g(x1 , t)∂t x1 .
x1 x1
La variation de la quantité intégrée, est égale à l’intégrale de la variation à laquelle on ajoute le flux qui sort g(x2 , t)∂t x2 et le flux qui entre −g(x1 , t)∂t x1 .
Utilisons cette relation, intégrons l’équation de conservation de la masse du fond f à la surface η (sans faire la transformation de Prandtl)
Z η
∂U1 ∂U2
( + )dz = 0
f ∂x ∂z
cela donne la relation suivante où sont apparues les vitesses transversales en bas et en haut de la couche liquide :
Z η
∂U1
( )dz + U2 (x, η, t) − U2 (x, f, t) = 0
f ∂x
- ESV . 19-
Saint-Venant
Remarque, la formule est miraculeusement vraie en fluide parfait. Dans le cas d’un fluide parfait au fond, U2 (x, f, t), la vitesse transverse n’est pas nulle, mais
U2 (x, f, t) = U1 (x, f, t) ∂f
∂x , la vitesse à la surface est toujours la même et le résultat global est toujours valide.
intégrons l’équation de conservation de la masse du fond 0 à la surface h (en ayant fait la transformation de Prandtl)
Z h
∂U1 ∂U2
( + )dz = 0
0 ∂x ∂z
cela donne la relation suivante où sont apparues les vitesses transversales en bas (nulle) et en haut de la couche liquide :
Z h
∂U1
( )dz + U2 (x, h, t) − 0 = 0
0 ∂x
Remarque, la formule est plus simple à démontrer avec la transformation de Prandtl puisque l’on fait disparaı̂tre la position du fond.
et les termes suivants, dont certains sont peut être plus grand que d’autres , l’équation de conservation de la masse (incrompressibilité)
∂U1 ∂U2
+ =0
∂x ∂z
- ESV . 20-
Saint-Venant
et l’équation de conservation de la quantité de mouvement dont les principaux termes sont le long de x, avec seulement le frottement transverse, et remarquant qu’il
n’y a pas d’équation de quantité de mouvement suivant z (comme en couche limite).
2
∂U ∂ h
( 0 U12 dz) − U12 (x, h, t) ∂h
R
• ∂x1 devient ∂x ∂x + 0
∂U1 U2
• ∂z devient U1 (x, h, t)U2 (x, h, t) − U1 (x, 0, t)U2 (x, 0, t)
En facteur de U1 (x, h, t), on retrouve la condition de vitesse à la surface : ∂h ∂h
∂t + U1 (x, h, t) ∂x = U2 (x, h, t). Puis l’expression des conditions d’adhérence au fond
permettent d’annihiler U1 (x, 0, t) et U2 (x, 0, t). Il ne reste donc pour le membre de gauche
Z h Z h Z h !
∂U1 ∂U1 U1 ∂U1 U2 ∂ ∂
( ρ( + + ) dz) = ρ ( U1 dz) + ( U 2 dz)
0 ∂t ∂x ∂z ∂t 0 ∂x 0 1
- ESV . 21-
Saint-Venant
R 2
ΓQ2
Rh Γ h
On va donc faire l’approximation ( 0 U12 dz) = h 0
U1 dz = h , et poser une vitesse moyenne u = Q/h pour que :
Z h Z h
∂ ∂ ∂ ∂ ΓQ2 ∂ ∂
( ( U1 dz) + ( U12 dz)) s’écrive ( Q + ( )) ou encore ( (hu) + (Γhu2 ))
∂t 0 ∂x 0 ∂t ∂x h ∂t ∂x
Rh Rh
Ce terme Γ = (h 0 U12 dz)/Q2 est appelé le coefficient de Boussinesq. Le rapport (h2 0 U13 dz)/Q3 sera le coefficient de Coriolis... Ce sont des facteurs de forme (c.f.
Chanson HOCF p28).
On va ensuite montrer que supposer Γ ' 1 permet bien de retrouver le terme d’accélération du modèle 1D. Mais auparavant nous devons aussi nous concentrer sur
le terme de frottement au fond τ12f qui doit lui aussi être écrit en fonction de h et Q.
- ESV . 22-
Saint-Venant
- ESV . 23-
Saint-Venant
• conservation de la masse
∂U1 ∂U2
+ =0
∂x ∂z
au fond, U1 = U2 = 0, l’interface, définie de manière implicite par la fonction F (x, z) = (h(x, t) − z). Cette interface est conservée dans le mouvement du fluide donc
la valeur de la vitesse transverse à l’interface en fonction du déplacement de l’interface (on néglige d’éventuels termes en < u01 h0 >) :
∂h ∂h
U2 (x, h, t) =
+ U1 .
∂t ∂x
• conservation de la quantité de mouvement écrite sous forme conservative (compte tenu de l’incompressibilité)
∂U1 ∂U12 ∂U1 U2 ∂h ∂τ12
ρ( + + ) = −ρg + ρgα +
∂t ∂x ∂z ∂x ∂z
avec τ12 (h) = 0 pas de contrainte en haut.
On a gardé que le terme dominant (somme de la viscosité normale plus la viscosité turbulente, on verra que l’on néglige souvent la viscosité µ) :
∂U1 ∂U1 ∂U1
τ12 = µ( ) − ρ < u01 u02 >= (µ + ρνt )( ) = (µ + µt )( )
∂x2 ∂x2 ∂x2
Il nous faut maintenant simplifier ces équations pour tenir compte de diverses hypothèses physiques sur la turbulence.
- ESV . 24-
Saint-Venant
∂
0 = (ρu2∗ /h)(z − h) + (µ + ρνt ) U (z),
∂z
ou encore
z ∂ z ∂
u2∗ (1 − ) = ν U (z)− < u0 w0 > ou u2∗ (1 − ) = (ν + νt ) U (z).
h ∂z h ∂z
En mesurant z avec h et U avec la vitesse moyenne U0
νU0 dū νU0 U0 /u∗
− < ū0 w̄0 >= (1 − z̄), = on pose R∗ = u∗ h/ν
hu2∗ dz̄ hu2∗ R∗
comme le Reynolds est grand, il reste
− < ū0 w̄0 >= (1 − z̄).
les fluctuations sont nécessairement linéaires en z. Or à la paroi les fluctuations son nulles, ce qui n’est pas le cas pour cette équation, il faut donc introduire une
couche limite pour assurer − < ū0 w̄0 >= 0 en z̄ = 0... On va trouver l’échelle qui convient dans la suite et surtout on va modéliser µt pour l’instant inconnu.
Viscosité turbulente longueur de mélange pour écoulement 2D sur une paroi lisse
Ludwig Prandtl a introduit la longueur de mélange pour construire la viscosité turbulente, par définition, la viscosité turbulente de Prandtl dépend du cisaillement
(la rotation, ou la vorticité) et d’une longueur qui est la taille des tourbillons les plus grands (arguments heuristiques de Prandtl), au plus simple cette longueur est
proportionnelle à la distance à la paroi avec une constante κ (qui a été appelée constante de Kármán, on prend κ = 0.41 = 1/2.439 ) donc on écrit grâce à cette
hypothèse :
∂
νt = `2 U (z), avec ` = κz.
∂z
- ESV . 25-
Saint-Venant
∂
On a vu plus haut la relation d’équilibre pente frottement intégrée 0 = (u2∗ )(z/h − 1) + (ν + νt ) ∂z U (z). On y substitue νt , elle devient :
z ∂ ∂
0 = (u2∗ )( − 1) + (ν + κ2 z 2 U (z)) U (z),
h ∂z ∂z
on veut garder le maximum de termes dans cette expression. Par moindre dégénérescence on voit que u2∗ joue un rôle important, on voit qu’il faut prrendre u2∗ ∼
∂ ∂
κ2 z 2 ( ∂z U (z))2 donc l’échelle de vitesse est u∗ , ensuite u2∗ ∼ ν ∂z U (z)), donc l’échelle en z est ν/u∗ . On voit que si utilise ν/u∗ comme unité de longueur, et u∗ comme
∂
échelle de vitesse, on a le PMD et on choisit historiquement d’écrire y + = zu∗ /ν sans dimension et aussi u+ = u/u∗ sans dimension, donc ∂z U (z) = (u2∗ /ν) ∂y∂+ u+ et
2 +2 ∂ + 2
ainsi comme νt = ν(κ y ∂y+ u ) en simplifiant par u∗ , on voit que l’on a obtenu les bonnes échelles pour notre équation
ν + ∂ ∂
0=( y − 1) + (1 + κ2 y +2 + u+ ) + u+ .
hu∗ ∂y ∂y
On peut la résoudre, en supposant toujours que hu∗ /ν 1, on voit que si y + 1 alors la vitesse est linéaire
∂ + 1
pour en pratique y + > 10, +
u = , l’intégration donne un logarithme
∂y κy +
C’est la ”loi log” de la turbulence. En fait, on remarque que l’équation 1 = (1 + κ2 y +2 ∂y∂+ u+ ) ∂y∂+ u+ s’intègre exactement en
√ 2 κ2 +1
q q
4y+
1
y+ − y+ + 2κ sinh−1 (2y+ κ) 1− 2 κ2 + 1
4y+ log 2 κ2 + 1 + 2y κ
4y+ +
u+ = , ou aussi +
2κ2 2y+ κ2 κ
on retrouve le comportement linéaire et logarithmique, voir le tracé figue 5. On a bien u+ = y + pour y + petit on a bien la loi linéaire et pour y + grand on retrouve
une loi log u+ = κ1 ln (y + ) + Log(κ)−1+Log(4)
κ = 2.439Log(y+ ) − 1.23 (car 1/κ = 1/0.41 = 2.439). Retenons que la vitesse est
1
u+ = ln y + + B, avec une constante B.
κ
Expérimentalement, il a été observé depuis 100 ans que B ' 5.5, et 2.4 varie des fois en 2.5, le 5.6 est des fois 5.75 donc
A partir d’ici, on néglige dans cette expression µ [Attention en aérodynamique, on ne néglige pas cette contribution qui joue un rôle car elle redonnera la ”sous
couche” laminaire sous la ”couche logarithmique”]. On néglige µ car la sous couche visqueuse en ν/u∗ va être négligeable par rapport à la taille des grains de sable ou
des cailloux au fond de la rivière.
Pour les tuyaux lisses, on revient aux équations de Navier Stokes en axisymétrique (il y a des r en plus par rapport au cas de la fente considéré auparavant) ;
l’équilibre entre la pression et le frottement dans les équations de Navier Stokes turbulentes donne
∂p 1 ∂τ ∂p ∂p D
0=− + , intégré sur la section cela donne 0 = − πR2 + 2πR(0 − τ (0)), soit τ (0) = − car D = 2R.
∂x r ∂r ∂x ∂x 4
- ESV . 26-
Saint-Venant
yν
Figure 5 – Profil de vitesse adimensionné u∗ en fonction de νy/u∗ en semi Log dans le cas de la fermeture de longueur de mélange de Van Driest λ = κy(1 − e− u∗ /26 )
utilisée en aérodynamique et correspondant aux expériences en rouge et dans le cas longueur de Prandtl en bleu. La courbe linéaire est en vert. Les pointillés sont les
lois Log approchées. A droite extrait de Schlichting, expériences montrant la loi log (u+ = 2.5 ln (y + ) + 5.5 ou en décimal, u+ ' 5.75 log10 (y + ) + 5.5), ici principalement
de Nikuradse.
∂p 2
Historiquement on a posé λ = (−D ∂x )/(ρ u2 ) (avec le diamètre D et la vitesse moyenne u, ce choix est historique). Ce coefficient est 4 fois le coefficient de frottement.
∂p
(−D ∂x ) τ (0) ρu2 2 ρu2
λ= u2
=8 , donc par définition de ce lambda τ (0) = λ , donc ρu∗ = λ ,
(ρ 2 ) ρu2 8 8
Dans le cas des tuyaux de rayon R (ou des couches limites turbulentes), une dernière observation expérimentale indique que près du centre (r ∼ 0) la vitesse a la forme
déficitaire U1 (r) = U1 (r = 0) − u∗ F ((R − r)/R) plus loin du centre en venant près de la paroi z = (R − r), on doit retrouver le profil log, donc −F ((R − r)/R) ∼ κ1 ln R
z
,
donc près de la paroi, la vitesse déficitaire se comporte
u∗ R
U1 (r) = U1 (r = 0) + ln( )
κ z
RR
on en prend le flux sur la section, r = R − z donc comme 0 2π(R − z) ln(R/z) = 3/2, et avec u la valeur moyenne de la vitesse, (Ryhming p264). On a donc
u = U1 (r = 0) − 3u approximation la forme logarithmique U1 (r = 0) = u∗ ( κ1 ln Ru
2κ en utilisant en première + B), donc on a une relation en tre la vitesse moyenne
∗ ∗
ν
et la vitesse de frottement u ∼ u∗ ln Ru
∗
ν
1 Ru∗ 3
u = u∗ ( ln +B− ).
κ ν 2κ
√
On a dit que l’on pose historiquement λ = 8(u∗ /u)2 , donc Ru ν
∗
= ( 4√λ2 ) uD
ν , on a donc la ”loi de friction universelle des tuyaux de Prandtl”, c’est une meilleure
approximation sous forme implicite pour les tuyaux lisses (elle est toujours utilisée de nos jours) :
√ !
uD √
r
8 1 λ uD 3 1
= ( ln ( √ ) +B− ) réécrit en √ = 2.0 log10 λ − 0.8,
λ κ 4 2 ν 2κ λ ν
- ESV . 27-
Saint-Venant
voir Schlichting p 611, les constantes ln 10, 1/κ et B sont ajustées expérimentalement pour trouver cette expression...
Nikuradse a montré expérimentalement qu’elle est correcte jusqu’à Re = 106 . En fait, la loi de Blasius (on en aprlera plus loin avec Chézy) λ = 0.316(uD/ν)1/4
et la ”loi de friction universelle des tuyaux de Prandtl” (relation implicite précédente) sont confondues jusqu’à Re = 105 , avec toujours une meilleure validation
expérimentale pour Prandtl.
Nous venons de voir que la vitesse est dans une large zone de la forme (assez loin du fond)
1
u+ = ln y + + B, avec une constante B qui dépend du modèle exact de longueur de mélange
κ
toute théorie de couche limite turbulente doit retrouver cette expression qui est bien connue depuis que Prandtl et Von Kármán se sont livrés dans les années 1920 à
une compétition amicale mais rude sur ce problème. Ces résultats sont dans le livre de Schlichting. Rappelons que pour un écoulement turbulent dans un tuyau (et de
nombreuses expériences depuis 100 ans l’on montré) on estime que l’on a en bonne approximation (avec 2.4 ou 2.5 pour 1/κ)
U zu∗
= 2.5 ln( )+B
u∗ ν
avec par définition u∗ la vitesse de frottement, c’est à dire que par définition τ0 = ρu2∗ . On a on estime donc que pour une rivière à fond lisse ou pour une couche
limite de vent, on est comme dans un tuyau et donc que la loi logarithmique est valable (c’est une approximation usuelle). Pour un fond rugueux, on s’attend à la
même forme, mais avec une autre échelle que ν/u∗ pour dimensionner la profondeur, soit z0 et alors U = uκ∗ ln( zz0 ) (attention, on peut aussi écrire U = uκ∗ ln( z+z z0 )
0
suivant la position de z = 0 et suivants les auteurs..., mais z0 est petit). Or, dans un tuyau rugueux de rugosité ks (taille des rugosités ou mini bosses plus ou moins
régulièrement réparties sur la paroi) Nikuradse a décomposé la vitesse du tuyau lisse en faisant apparaı̂tre ks artificiellement (dans le cas lisse il n’y a pas de ks !) :
U u∗ z z ks u∗
= 2.5 ln( ) + 5.5 = 2.5 ln( ) + 5.5 + 2.5 ln( )
u∗ ν ks ν
et défini ainsi une fonction de rugosité B, cette fonction est exactement 5.5 + 2.5ln(ks u∗ /ν) dans le cas lisse (ks → 0), et sinon elle a été mesurée expérimentalement
et tabulée (voir figure 6). Cette fonction B, est donc telle que pour un écoulement rugueux quelconque
U z
= 2.5 ln( ) + B
u∗ ks
avec en construisant un Reynolds sur les rugosités Rep
z
U = 2.5 ln( ) + B avec cas lisse Rep < 5, B = 5.5 + 2.5ln(ks u∗ /ν) avec cas rugueuxRep > 70, B = 8.5
ks
u∗ z
pour les cas très rugueux, B est constante=8.5, et comme e8.5/2.5 ' 30 d’où en pratique U = κ ln( z0 ) avec z0 ∼ d/30.
Dans le cas des écoulements turbulents en rivière, en canaux, ou en fleuves, on s’attend à retrouver un profil établi de forme logarithmique qui est caractéristique
d’un écoulement turbulent [attention, les profils turbulents ont été observés et mesurés dans les tuyaux, dans les fleuves c’est moins clair, cette hypothèse de forme
logarithmique est une hypothèse usuelle, il est vraisemblable qu’elle est fausse, le coefficient 1/κ est ajusté.... des fois on prend une loi de puissance plutôt qu’un log] :
u∗ z
U (z) = ln( )
κ z0
z = 0 est le fond où sont présentes des granulosités de taille z0 , u∗ est appelée la vitesse de frottement, elle est telle que le frottement au fond τ0 = ρu2∗ .
- ESV . 28-
Saint-Venant
Figure 6 – Gauche :Diagramme de la fonction Roughness B de Nikuradse (Schlichting), on voit bien le comportement lisse B = 5.5 + 2.5ln(ks u∗ /ν) et rugueux
B = 8.5 Centre : profil de vitesse en log et définitition du z0 , droite deux exemples de profils théroriques en log fittant des données de terrain (Graeme Smart ”turbulent
velocity profiles and boundary shear in gravel bed rivers” JOURNAL OF HYDRAULIC ENGINEERING / FEBRUARY 1999
Schlichting a introduit la longueur de mélange pour construire la viscosité turbulente, par définition, la viscosité dépend du cisaillement (la rotation, ou la vorticité)
et d’une longueur qui est la taille des tourbillons les plus grands (arguments heuristiques de Prandtl) :
∂
νt = `2 U (z)
∂z
∂
or nous avons vu que l’on néglige µ pour obtenir 0 = (ρu2∗ /h)(z − h) + ρνt ∂z U (z) et on veut obtenir comme solution un profil logarithmique dont la dérivée est
∂U (z)/∂z = u∗ /(κ(z)), ces deux expressions nous donnent la longueur de mélange ad hoc (puisque la couche sous couche a été enlevée)
p
` = κ(z 1 − z/h).
Le frottement au fond est bien :
∂ ∂
τ0 = (ρ`2 U (z)) U (z)|z=0 soit τ0 = ρu2∗ .
∂z ∂z
La contrainte en z = 0 est ρu2∗ bien que l’on ait négligé le µ∂U/∂z|0 .
En résumé dans le cas d’un écoulement turbulent établi, on a donc l’équilibre pesanteur/ frottement dans le régime turbulent :
∂ ∂
0 = ρgα + [ρνt ( U (z))],
∂z ∂z
∂ ∂
la loi de viscosité turbulente µt /ρ = νt = `2 ∂z U (z) est une fonction donnée de z et ( ∂z U (z)) pour que l’on puisse résoudre le profil numériquement (par exemple
p u∗ z
` = (κ(z 1 − z/h)), qui redonne la forme logarithmique U (z) = κ ln( z0 )).
- ESV . 29-
Saint-Venant
ρ(Q/h0 )2
Q = u∗ h0 (ln(h0 /z0 )), et le frottement τ0 =
(ln(h0 /z0 ))2
Le frottement au fond compense exactement le terme moteur ρgα. La vitesse de frottement, qui est aussi la vitesse caractéristique est
p
u∗ = ghα
√
Dans le cas turbulent, la vitesse caractéristique n’est plus proportionnelle à α comme en laminaire, mais à α.
Nous retrouverons cette dépendance dans la loi empirique de Chézy (§7.2.3) et dans les discussions de §7.2.5.
- ESV . 30-
Saint-Venant
Soit on a sauté cette partie et on a avec u(x, t) la vitesse moyenne dans la couche d’épaisseur h posée sur un fond f ,
∂h ∂
+ (hu) = 0
∂t ∂x
∂ ∂ df ∂h
ρ( (hu) + (hu2 )) = −ρgh − ρgh − τ12f (6)
∂t ∂x dx ∂x
h = η − f.
Q = hu
Rh
U1 dz
Pour l’équation de la masse, pas de surprise, Il s’agit bien de la même équation, la vitesse moyenne u = 0
h et l’équation de la masse intégrée sur l’épaisseur
Z h
∂η ∂
+ ( U1 dz) = 0 (le fond est fixe ∂t f = 0 et h = η − f.) est bien l’équation :
∂t ∂x 0
∂h ∂
+ (hu) = 0
∂t ∂x
Rh
pour la quantité de mouvement, on a un problème car le carré de la vitesse intervient : ( 0 U12 dz), et bien entendu
Z h Z h
1
( U12 dz) 6= ( ( U1 dz)2 ).
0 h 0
Il faut écrire cette intégrale en fonction de u la vitesse moyenne dans la tranche et h la hauteur d’eau. Ce que l’on ne peux pas faire si on n’a pas la solution, à moins
d’inventer une forme à la solution U1 (x, z). Cette solution aura une forme qui ressemble à la réalité. Cette idée est liée aussi aux solutions ”semblables” qui sont
des solutions simples des équations aux dérivées partielles et dans lesquelles les deux variables x et z sont assemblées en une seule. On dira donc que la vitesse suit
l’Ansatz :
U1 (x, z) = u(x, t)ū(z/h)
où on a décomposé la dépendance en x et t dans u(x, t) la dépendance transverse est comprise dans ū, en fait ū est une fonction de z̄ = z/h, qui varie de 0 à 1 quand
Rh
z̄ varie de 0 à 1. Le u est tel qu’il est la moyenne de la vitesse : ( 0 U1 dz) = hu.
- ESV . 31-
Saint-Venant
3
ū = z̄(2 − z̄)
2
R1 R1 6
où z̄ est z/h. La forme de vitesse choisie elle est telle que le flux 0
ūdz̄ = 1 et on peut calculer le carré inconnu 0
ū2 dz̄ = 5 ainsi que le frottement au fond ∂ ū/∂ z̄|0 = 3
donc la fermeture laminaire de nos quantités à modéliser est :
Z h
6 2 u
( U12 dz) = u h, τ = µ∂u/∂z|0 = 3µ
0 5 h
∂(hu) ∂ 6 2 ∂η u ∂η ∂
+ ( hu ) = gαh − gh − 3ν avec toujours + (hu) = 0 et h = η − f.
∂t ∂x 5 ∂x h ∂t ∂x
• lorsque l’écoulement est bien établi, on constate que l’équilibre est
∂(hu) ∂ ∂η 1 ∂η ∂
+ (hu2 ) = gαh − gh − cf u2 , avec toujours + (hu) = 0 et h = η − f.
∂t ∂x ∂x 2 ∂t ∂x
En s’inspirant des écoulements dans les tuyaux (voir Schlichting chapitre XX, et voir plus loin), on peut l’estimer en se fondant sur les mesures dans les tuyaux de
Nikuradse, les variations de ce coefficient avec la taille des rugosité et le Reynolds sont indiqués au §6.7.3.
En réalité, la science de l’hydrologie est plus ancienne que celle de l’aérodynamique (du début du XXème), d’autres simplifications historiques ont été proposées.
Le plus ancien est celui de Chézy 1778 que nous allons voir (puis Manning et les autres) dans les paragraphes suivants (§7.2.3).
- ESV . 32-
Saint-Venant
Lorsque l’écoulement est bien établi, pour une rivière quelconque avec une
pente α faible, on a un équilibre entre le terme moteur lié à la pente de la rivière
qui va agir sur la surface mouillée, zet le terme de friction au sol qui agit sur tout
le périmètre mouillé :
0 = ρgαA − τ12 P
On retrouve en 2D pour une largeur infinie l’équation 7 (pour une rivière large
A/P ∼ h). En utilisant le coefficient de frottement
√ turbulent cf défini au para-
graphe précédent on voit que la vitesse est en α :
r
g A
r
u= α,
cf /2 P
Figure 7 – Périmètre mouillé et surface d’une rivière : pour une rivière large PA ∼ h.
c’est à dire que u est reliée à la pente α par la relation de proportionna- Attention, la définition de rayon hydraulique n’est pas fixée, certains auteurs définissent
RH = A/P , donc pour une gouttière A = πR2 /2 P = πR, RH = R/2 ), et pour une rivière
r
√ A
lité u ∝ α, on pose le coefficient CC tel que u = CC α . Par définition, large c’est A ∼ h (et non le double comme le donnerait la première définition).
P P
q
A
le coefficient de proportionnalité entre la vitesse moyenne et P α est noté
q
CC = cfg/2 , il est appelé le coefficient de Chézy, il n’est pas sans dimension, il
est en m1/2 /s. Cela provient d’un choix historique et expérimental de Antoine Chézy (1718-1798), il mesurait la hauteur d’eau, la vitesse et la pente. Il a remarqué
que la vitesse est en racine de A/P et de la pente α.
Historiquement en 1776, pour le canal du Courpalet, Chézy trouva Cc = 31m1/2 /s et 44 pour la Seine. Les valeurs vont de 30 pour un canal petit et rugueux à 90
m /s pour un canal large et lisse. Dans le cas de la cuve FC80 utilisée en TP, le manuel de la cuve donne une valeur du Chézy de 51m1/2 s−1 , des mesures effectuées
1/2
Malheureusement, ce coefficient historique dimensionné continue à être utilisé de nos jours. On écrit donc pour le frottement :
g 2 g
τ =ρ u = ρ 2 2 Q2
Cc2 Cc h
Il est clair qu’il faut le redéfinir pour extraire un coefficient de frottement sans dimension cf .
- ESV . 33-
Saint-Venant
Robert Manning (1816 1897), irlandais, autodidacte en mécanique des fluides est associé à une formule plus précise que celle de Chézy, mais il semble que c’est le
français Philippe Gauckler (1826-1905) qui l’a établie grâce aux expériences de Henry Bazin. Le Suisse Albert Strikler (1887-1963) a aussi popularisée cette formule
√
avec de nombreuses expérimentations : u ∝ (A/P )2/3 α . Strickler, Gauklert et Manning sont donc associés à cette formule pour l’avoir développée et popularisée.
On définit aussi nGM le coefficient de Gauckler Manning (sm−1/3 ) tel que
1 −2 nGM 2 2
CC = (A/P )1/6 pour un fleuve large on a CC = nGM 2 (h)−1/3 et le frottement est τ = ρg Q
nGM h7/3
On définit aussi le coefficient de Strickler qui est l’inverse du Manning. Dans la suite, on va souvent négliger cette variation assez lente en 1/6 et garder une
1
valeur constante pour CC . On remarque que pour un grand W (une rivière large est telle que A/P ∼ h) on a CC = nGM (h)1/6 , et donc le coefficient de frottement
g 2 −1/3 2 3
cf /2 = C 2 = gnGM h . On peut donc interpréter (gnGM ) comme une rugosité équivalente du fond
C
Dans le cas Gauckler-Manning-Strickler, frottement τ s’écrit (à grande largeur) τ = ρgn2 Q2 /h7/3 à l’équilibre τ = ρghS0 avec S0 la pente du lit, donc ρgn2 u2 /h1/3 =
ρghS0 donc on trouve la vitesse de frottement de Manning √
S0 2/3
u= h
n
valeurs en sm−1/3 : béton lisse 0.01 bois 0.014 sable assez fin 0.02 sable grossier, gravier 0.03, Galets 0.05, herbe végétation 0.2 (voir tableau [Link]
edu/geowater/FX3/help/8_Hydraulic_Reference/Mannings_n_Tables.htm
Débit de la Seine dans le XIIIème environ 300m3 /s, profondeur environ 1.5m. En cas de crue, janvier 2011, il passe à 985m3 /s (niveau 3.38m)
7.2.5 Autres coefficient de frottement (inspirés du tuyau) : Blasius, Prandtl, Fanning et Darcy-Weisbach
• frottement en ρu2
Chézy, élève des Ponts et Chaussées, fut le directeur ce cette école, et se rajouta une particule ”de Chézy”, Prony le remplaça, et après Darcy, Bazin... Pendant que
- ESV . 34-
Saint-Venant
la science de l’hydrologie se transmettait d’ingénieur des Ponts en ingénieurs des Ponts, la mécanique des fluides d’inspiration plus mathématique se développait dans
les autres branches. L’impulsion vient de l’Université Allemande fondée sur une interaction forte entre Mathématiques Appliquées et Expériences. Nous faisons ici un
rappel sur les tuyaux lisses, voir la référence H. Schlichting ”Boundary Layer Theory” Chapitre XX ”turbulent flow through pipes” qui synthétise les recherches des
années 1920 de Blasius, Prandtl, von Kármán sur la compréhension de la turbulence, leurs idées sont restées justes 100 ans après. Globalement, pour un frottement
turbulent, ρu2 est l’ordre de grandeur du frottement, on introduit comme nous l’avons déjà dit plus haut un coefficient de friction (introduit aussi par Gustave Eiffel
1832-1923 qui utilisait sa tour comme une énorme soufflerie verticale) cf et parfois appelé f coefficient de Fanning (l’américain John Fanning (1837-1911) cf ≡ f ). Le
frottement est écrit (en conservant l’anachronisme du 1/2 héritage des Suisses Bernoulli dont Daniel Bernoulli 1700-1782) avec la vitesse moyenne u au carré :
u2
τ = cf ρ .
2
• frottement dans ls tuyaux
Pour les tuyaux, si on revient aux équations de Navier Stokes (nous avons déjà vu ceci au §6.7.3) ; l’équilibre entre la pression et le frottement dans les équations de
Navier Stokes turbulentes donne
∂p 1 ∂τ ∂p ∂p D
0=− + , intégré sur la section donne 0 = − πR2 + 2πR(0 − τ (0)), soit τ (0) = − car D = 2R.
∂x r ∂r ∂x ∂x 4
∂p 2
Historiquement on a posé λ = (−D ∂x )/(ρ u2 ) (avec la vitesse moyenne u). Il faut dire que lors de la révolution industrielle, à partir de 1840, de nombreux scientifiques
se sont posé des questions sur les arrivées d’eau et le écoulements dans les tuyaux (bien sûr, le tout premier d’entre eux attesté date de l’époque romaine, c’est Sextinus
Frontin vers 70, il s’intéressait au diamètre des conduites des aqueducs et aux pertes d’eau par des voleurs). C’est vraisemblablement Henry Bazin (1829-1917) qui
est initiateur d’expériences fines, meilleures que celles originelles de Gaspard de Prony (1755-1839), le coefficient s’appelant coefficient de Henry Darcy (1803-1858),
dont Bazin était l’élève). La formule a ensuite été modifiée par Julius Weisbach (Allemand 1806-1871) en 1845. Toujours est il que de nos jours la relation entre ce λ
historique et cf le coefficient de friction utile est telle que cf = λ/4, la friction est définie par
u2 1
τ = cf ρ = λρu2 .
2 8
En laminaire on vérifie facilement que λ = 64/Re avec Re = uD/ν (Reynolds sur la vitesse moyenne u et le diamètre D).
Pour les tuyaux lisses, en turbulent, Blasius a montré expérimentalement (Re = uD/ν et D = 2R) une dépendance en Re−1/4 :
−1/4
uD
Formule de ”Blasius turbulente” (à ne pas confondre avec la laminaire) : λ = 0.3164
ν
voir Schlichting ”Boundary Layer Theory” p 596. Le frottement turbulent dans un tuyau suit la loi suivante (ré obtenue par Ivan Nikuradse 1894 -1979 et d’autres,
écrite avec D ou R = D/2 et comme 0.3164 /8=0.03955 et 21/4 = 1.19) on a les expressions suivantes avec u vitesse moyenne :
2
−1/4 −1/4 2 −1/4
ρu uR 2 2 uR ρu uD
τ0 = 0.03955
ν 1/4 = 0.03325 ν
ρu = ρ(u∗ ) = 0.0665 et aussi λ = 100
( ) ν 2 ν
uD
p
en ayant défini u∗ = τ0 /ρ la vitesse de friction, on a alors à partir de cette expression la dépendance en rayon en 1/7 pour a vitesse moyenne u :
1/7
u u∗ R
= 6.99 .
u∗ ν
- ESV . 35-
Saint-Venant
Cette dépendance a été utilisé longtemps comme une première bonne approximation de la vitesse turbulente. Comme déjà dit en §6.7.3, pour faire mieux que Blasius
et que Kármán, car avec Prandtl ils étaient en compétition, Prandtl a proposé une autre formule ”La loi de friction universelle des tuyaux de Prandtl”, c’est une
meilleure approximation pour les tuyaux lisses (Re = uD/ν vitesse moyenne u et le diamètre D) :
uD √ Re √
1 1
√ = 2.0 log10 λ − 0.8 ou √ = 2.0 log10 λ
λ ν λ 2.51
voir Schlichting p 611. [Pour l’établir il faut savoir que vers le centre d’un tuyau, le profil de vitesse est de la forme déficitaire U1 (r) = Ucentre − u∗ F (r/R), donc et
raccorder les deux]. Nikuradse a montré expérimentalement qu’elle est correcte jusqu’à Re = 106 . En fait, la loi de Blasius et la loi de Prandtl sont confondues jusqu’à
Re = 105 , avec toujours une meilleure validation pour Prandtl.
La formule de Nikuradse a été généralisée par Colebrook (1917-1997) et White 1937 puis Colebrook 1939 dans le cas d’un tuyau avec des rugosités, k taille de la
rugosité, Re = uD/ν (Reynolds sur la vitesse moyenne u et le diamètre D) :
1 2, 51 k
√ = −2 log( √ + )
λ Re λ 3, 7D
et la ”harpe de Nikuradse” a été incluse dans le diagramme de Moody.
Par la conversion subtile le ”diamètre hydraulique” D est le diamètre du tuyau, mais le rayon hydraulique Rh n’est pas le rayon mais Rh = D/4, pour les
écoulements en canaux on utilise donc
1 2, 51 k
√ = −2 log( √ + )
λ Re λ 14.8Rh
k 2.51
si on pose b = 14.8(D/4) , a = Re on peut résoudre
1
λ= 2
ln 10 b
2W 2a 10 2a b
−
ln 10 a
Dans le cadre des rivières, le coefficient cf est parfois appelé fF coefficient de Fanning, on définit aussi fDW = 4fF coefficient de Darcy-Weisbach (c’est le λ).
Lorsque le fond est rugueux, la partie en Re est négligeable, il ne reste que le morceau avec kh taille de la rugosité importante. Comme g/Cc2 = λ/8 par définition, de
ces coefficients, on peut employer la formule suivante pour le Chézy, formule modifiée de celle de Colebrook
p
CC = 8g(2 log(14.8Rh /kh ))
avec kh taille de la rugosité, Rh . on trouve dans la littérature des fois avec 12.7 au lieu du 14.8 (comme plus haut dans Colebrook).
- ESV . 36-
Saint-Venant
u|u| cf cf
Sf = cf soit τ = ρ u2 = 2 Q2
2gh 2 2h
Fx = 2ρhωsinλv, Fy = −2ρhωsinλu.
Pour les calculs pratiques, les domaines ne sont pas très étendus et le coefficient fc est constant avec fc = 2ωsinλ, à la latitude de la France 48 ˚, fc = 1.08 10−4 , (un
jour sidéral est de 23h56mn04s).
- ESV . 37-
Saint-Venant
∂ ∂ ∂ ∂ 1 ∂ ∂2 ∂2 ∂2
u + (u + v )u + w u = − ps + (νx 2 + νy 2 + νz 2 )u + f v
∂t ∂x ∂y ∂z ρ ∂x ∂x ∂y ∂z
∂ ∂ ∂ ∂ 1 ∂ ∂2 ∂2 ∂2
v + (u + v )v + w v = − ps + (νx 2 + νy 2 + νz 2 )v − f u
∂t ∂x ∂y ∂z ρ ∂y ∂x ∂y ∂z
vraisemblablement νx = νy 6= νz . On s’autorise des viscosités différentes dans les différentes directions. Ces équations sont légèrement différentes par l’existence de
cette surface ”rigide”, mais de la résolution du système sous la contrainte d’incompressibilité, on obtient ps donc la perturbation de vague.
- ESV . 38-
Saint-Venant
Si on ne l’a pas fait, on peut maintenant revenir à la section ” Saint Venant à partir de Navier Stokes”.
Dans la suite nous examinons des solutions de ces équations dans des cas pratiques, inondation, crue, tsunamis, marées, rupture de barrage.
Figure 9 – Photo PYL. La Yamunâ fait partie des Sept rivières sacrées de
l’Inde, sa longueur est 1 370 kilomètres, elle se jette dans le Gange. Environ
57 millions de personnes dépendent des eaux de la Yamunâ. Avec un volume
annuel d’environ 10 000 millions de m3 et l’utilisation de 4 400 millions de
m3 - dont 96 % pour l’irrigation - la rivière participe pour plus de 70 % à
l’approvisionnement en eau de Delhi. La Yamunâ est, après le Gange et la
Sarasvati aujourd’hui disparue, le cours d’eau le plus sacré en Inde. Elle est
considérée comme la fille de Sûrya, le dieu du soleil, et la sœur de Yama,
le dieu de la mort et, selon la tradition, ceux qui prennent un bain dans les
eaux saintes du fleuve ne craignent pas la mort. (source wiki)
- ESV . 39-
Saint-Venant
Deuxième partie
Le Fluide : Equations de Saint-Venant, grandes classes de
simplifications
On a maintenant établi les équations, rappelons leur différentes formes, puis montrons quelques exemples de dégénérescences caractéristiques de ces équations.
C’est à dire quelques exemples de simplifications de certains des termes pour obtenir des équations encore plus simples. Bien entendu, ces simplifications seront liées à
des hypothèses comme par exemple la forme du terrain, s’il est par exemple en pente forte ou faible, si le frottement est fort ou faible... On comprend bien aussi que
l’on a des formes un peu différentes suivant que l’on considère un canal de largeur constante, ou infinie (c’est pareil !), ou de largeur variable.
On fait donc attention à la définition du flux, des fois, il s’agit du flux par unité de largeur (cas invariant transverse : infini, ou largeur constant), ou il s’agit du
vrai flux (Q ou uh ou uA...), on fait alors attention à la surface mouillée et à l’aire A. Des fois la position de la surface libre η intervient, des fois non.
cf
Chacune de ses formes a une utilité pour résoudre un problème donné. Nous verrons ensuite les applications. On met ici le frottement g C12 u2 plutôt que 2 u2
C
∂h ∂(uh)
+ =0
∂t ∂x
∂(hu) ∂ ∂η u2
+ hu2 = −gh −g 2
∂t ∂x ∂x CC
η =f +h
- ESV . 40-
Saint-Venant
• Forme non conservative issue de la précédente, Q = hu, on ne se préoccupe toujours pas de la largeur du canal, f forme du fond (pas de η)
∂h ∂(hu)
+ =0
∂t ∂x
∂u ∂ ∂h u2
+ u u = −gf 0 − g −g 2
∂t ∂x ∂x CC h
• Forme conservative avec prise en compte d’une section variable Si A est la section mouillée, Q est le flux à travers cette surface
∂A ∂Q
+ =0
∂t ∂x
∂Q ∂ Q2 df ∂h Q2
+ = −gA − gA − Ag 2 2
∂t ∂x A dx ∂x CC A (A/P )
( ∂h ∂(hu)
+ =0
∂t ∂x (8)
∂u ∂ ∂h u2
+ u u = +gα − gfb0 − g −g 2 .
∂t ∂x ∂x CC h
On cherche les échelles pertinentes de manière à avoir des coefficients sans dimension et une dégénérescence significative. Il s’agit bien entendu de la même analyse
phénoménologique (voir 6.7.1) qu’auparavant, mais effectuée sur les équations intégrales de Saint-Venant en turbulent.
Une première idée consiste à utiliser la profondeur caractéristique h0 de l’écoulement pour mesurer h = h0 h̄ et écrire x̄ = x/(h0 /ε) = εx/h0 . Cet ε est là pour nous
rappeler qu’intrinsèquement les équations de Saint-Venant sont valide pour des couches minces (longueur h0 /ε h0 ). En fait, il faut utiliser les équations de Navier
Stokes si x̄ = x/h0 .
Le deuxième problème est dans l’estimation de la vitesse, celle ci peut être imposée par un moyen extérieur (déversoir, source...) ou être la réponse à la pente du
fleuve α qui génère l’écoulement. Notons U0 la vitesse, la conservation de la masse est alors adimensionnée par moindre dégénérescence (tous les termes sont de même
ordre de grandeur), car on prend pour le temps t̄ = tU0 /(h0 /ε) : ∂t̄ h̄ + ∂x̄ (ūh̄) = 0. Dans la quantité de mouvement
∂u ∂ ∂ ū ∂ U 2 ∂ ū ∂
+ u u = εU02 /(h0 )( + u u) = g(ε 0 ( + u u))
∂t ∂x ∂ t̄ ∂x gh0 ∂ t̄ ∂x
U02
on fait apparaı̂tre gh0 devant l’accélération, c’est F 2 le nombre de Froude au carré, les équations sans dimension sont donc :
2
ū2
∂ h̄ ∂(ūh̄) ∂ ū ∂ U0 ∂
+ = 0, & εF 2 ( + ū ū) = α(1 − √ ( )) − ε (h̄ + f¯b ).
∂ t̄ ∂ x̄ ∂ t̄ ∂ x̄ CC αh0 h̄ ∂ x̄
- ESV . 41-
Saint-Venant
on a fait apparaı̂tre le rapport C U
√0 , en effet à la pente α qui globalement emporte l’écoulement sur une hauteur h0 on associe la vitesse de Chézy en observant
C αh0
u2 √
gα ∼ g 2 , d’où CC αh0 (est la vitesse de Chézy), si on mesure avec cette vitesse,
CC h
∂ h̄ ∂(ūh̄) ∂ ū ∂ ū2 ∂
+ = 0, & εF 2 ( + ū ū) = α(1 − ) − ε (h̄ + f¯b ).
∂ t̄ ∂ x̄ ∂ t̄ ∂ x̄ h̄ ∂ x̄
U
Comme α est très petit, et ε est aussi très petit, et que F peut parfois être très petit, de même CC
√0
αh0
est grand ou petit.... il faut voir au cas par cas.
échelle d’espace, imposée par exemple par la topographie f¯ = f¯b , les équations précédentes deviennent :
( ∂ h̄ ∂(ūh̄)
+ =0
∂ t̄ ∂ x̄ (9)
∂ ū ∂ ∂
F 2( + ū ū) = − (h̄ + f¯).
∂ t̄ ∂ x̄ ∂ x̄
avec F d’ordre un. Nous parlerons d’écoulement super critique (ou torrentiel) si F > 1 ou sub critique (ou fluvial) si F < 1. Il s’agit des équations de Saint-Venant
sans frottement (cf plus loin). Il s’agit des équations inviscid Shallow water.
√
Remarquons qu’un choix judicieux de l’échelle de vitesse est alors U0 = gh0 qui fait disparaı̂tre le Froude :
( ∂ h̄ ∂(ūh̄)
+ =0
∂ t̄ ∂ x̄ (10)
∂ ū ∂ ∂
( + ū ū) = − (h̄ + f¯).
∂ t̄ ∂ x̄ ∂ x̄
∂ h̄ ∂(ūh̄)
+ =0
∂t ∂ x̄
∂ ū ∂ ∂ ū2
F 2( + ū ū) = 1 − (h̄ + f¯b ) −
∂ t̄ ∂ x̄ ∂ x̄ h̄
Si même F est d’ordre 1, soit plus grand soit plus petit que 1 (nous parlerons d’écoulement super critique si F > 1 ou sub critique si F < 1), tout est d’ordre 1. On a
alors le sytème le plus riche et le plus compliqué à résoudre.
- ESV . 42-
Saint-Venant
8.2.4 Saint-Venant, équations adimensionnées, échelle intermédiaire, mais vitesse lente (fluvial)
Si maintenant, on considère un écoulement lent, F 1, la vitesse s’ajuste toujours au frottement (on mesure la vitesse avec Chézy) et aux variationx de surface
libre et de fond mais l’inertie est faible : On obtient alors le problème de l’onde de crue diffusante :
∂ h̄ ∂(ūh̄)
+ =0
∂t ∂ x̄
∂ ū2
0=1− (h̄ + f¯b ) −
∂ x̄ h̄
¯
On va voir que l’on peut exprimer ū avec h̄ et fb dans l’équation de quantité de mouvement et tout mettre dans la première équation. Ce cas sera traité avec plus de
détails au paragraphe §8.6.
∂ h̄ ∂(ūh̄) ū2
+ = 0, et seulement 0 = (1 − )
∂ t̄ ∂ x̄ h̄
Pour des échelles longues par rapport, on voit que l’écoulement est toujours un équilibre entre le frottement au fond et la pente qui le provoque. Ce cas ”longue échelle”
sera appelé ”Onde d’inondation” ou ”onde de crue” ”flood wave” ou onde cinématique (kinematic wave). Ce cas sera traité avec plus de détails au paragraphe §8.5.
ū2
= 1.
h̄
On a débranché tous les termes les uns après les autres pour la simplification ultime..
- ESV . 43-
Saint-Venant
∂ h̄ ∂(ūh̄) ∂ ū ∂ ∂ ū2
+ = 0, et F 2( + ū ū) = 1 − (h̄ + f¯b ) −
∂t ∂ x̄ ∂ t̄ ∂ x̄ ∂ x̄ h̄
dans le cas d’une petite bosse de taille fb = O(), on pose fb = f1 . ū = 1 + ū1 , h̄ = 1 + h̄1 , (ūh̄) = 1 + (ū1 + h̄1 ) + O(2 ) donc
- ESV . 44-
Saint-Venant
Troisième partie
Le Fluide : Equations de Saint-Venant, quelques propriétés
La suite de ce cours consiste à reprendre toutes les simplifications proposées et les étudier plus précisément une à une pour voir comment elles se répondent et quel
modèle simplifié sera utilisé pour quel cas simplifié.
8.5 Onde d’inondation ou ”onde de crue” (simplification à longue échelle), onde cinématique
8.5.1 Equation en onde longue
Dans cette section, nous allons regarder le cas où le domaine est très grand, et où on est en ”longue échelle”. Il s’agit du cas où l’on peut négliger les termes non
linéaires de convection, on est dans la description ”onde longue”, on a toujours l’équilibre entre le frottement et une pente très douce −α. On suppose que le fond
est en fait f = −αx sans perturbation supplémentaire de cette surface plane. Dans la littérature anglaise, on parle de Flood wave ou de kinetic wave Whitham [27] p
82, ou Fowler [8] p76 ou Chanson [2], voir aussi Daubert [5] On néglige donc à la fois les termes non linéaires de convection, les termes de pression, (on est dans la
description onde longue), on a principalement l’équilibre entre le frottement et une pente très douce, et bien sûr, la conservation de la masse :
∂A ∂Q Q2
+ =0 & 0=α− 2 A2 (A/P ) ,
∂t ∂x CC
où, si on préfère, il n’y a qu’une équation liant les dérivées partielles de A (ou de Q puisque l’un est directement fonction de l’autre) :
r
∂A ∂ A3
+ (CC α ) = 0.
∂t ∂x P
Pour des rivières lors d’inondation P le périmètre varie peu par rapport à A, on a Q ∼ A3/2 .
- ESV . 45-
Saint-Venant
4
3.5
3
2.5
2
1.5
1
0.5
-3 -2 -1 0 1 2 3 4
• Examinons maintenant le cas où la célérité c n’est plus constante. Dans ce cas, le flux Q est constant le long des trajectoires caractéristiques x(t) définies par
dx
= c(x(t), t),
dt
cette vitesse c est une fonction de Q uniquement (car h est fonction de Q). Le long de ces caractéristiques, on a bien Q(x(t), t) constant
d ∂Q ∂x ∂Q ∂Q ∂Q
Q(x(t), t) = + = c+
dt ∂x ∂t ∂t ∂x ∂t
∂Q d
or comme par définition ∂t + c ∂Q
∂x = 0, on a bien Q(x(t), t) = 0.
dt
• Réciproquement, le long des droites x = ct + ξ du plan (x, t) la valeur de la solution est constante, c’est la valeur en t = 0, c’est à dire en ξ qui vaut ensuite x − ct
qui est conservée.
x = ξ + t c(ξ, 0) et Q(x, t) = Q(ξ, 0)
Donc, si on a une distribution au temps t = 0 de débit Q(x, 0) = Q0 (x) initial, chacun de ces points x est en fait un des points initiaux ξ, alors la solution est
Q = Q0 (x − t c(Q0 (ξ)))
(voir par exemple ”Riemann Solvers and Numerical Methods for Fluid Dynamics A Practical Introduction” 3rd edition E. Toro page 67 [26]). ou voir Chapitre 2 de
Whitham p 19 [27]).
On peut vérifier que c’est vrai en posant Q00 = dQ0 /dξ et ξ = x − ct, c est une fonction de Q donc de ξ, par dérivation on a puisque Q(x, t) est une fonction de ξ :
∂Q ∂ξ ∂Q ∂ξ
= Q00 , = Q00
∂t ∂t ∂x ∂x
∂ξ ∂(ct) ∂ξ ∂c ∂ξ ∂c ∂ξ ∂ξ ∂c ∂ξ
avec : =0− , =1− t ce qui donne aussi : =− t − c, =1− t, on regroupe pour écrire
∂t ∂t ∂x ∂x ∂t ∂ξ ∂t ∂x ∂ξ ∂x
∂ξ −c ∂ξ 1
= ∂c
, = ∂c
∂t 1 + ∂ξ t ∂x 1 + ∂ξ t
donc en substituant dans l’équation ∂t Q + c∂x Q on voit qu’elle est bien vérifiée :
∂Q ∂Q ∂ξ ∂ξ (−c + c) 0
+ c(Q) =( + c(Q) )Q00 = ∂c
Q0 = 0,
∂t ∂x ∂t ∂x 1 + ∂ξ t
- ESV . 46-
Saint-Venant
La référence [18] Miller 1983 ”Basic Concepts of Kinematic-Wave Models” [Link] donne des applications pratiques.
La simulation numérique de cette ondede crue (flood wave) est présentée dans les codes
[Link] en C
[Link] en Basilisk,
[Link]
!
"+
3 !"
" =0
en Python.
!T 2 !#
pour simplifier les notations on pose t=T et x=(2/3)#. L'équation d'advection
3/2 3
On remarque que si Q ∼ A , alors c = !"
!t
+"
!"
!x
=0
2 Q/A, donc c est supérieur à la vitesse moyenne de l’écoulement Q/A. Le système forme un choc : une crue va avoir
tendance à former un mur d’eau.
peut se résoudre connaissant "(x,t=0)="0(x) la distribution initiale. Le long des droites:
x=" (#)t+# (passant en t=0 en x=# la où la hauteur était " (x=#)) la hauteur est conservée:
0 0
"(x,t)="0(#) = "0(x-"0(#)t)
On vérifie !u/!t="0'(#)(!#/!t) et !u/!x="0'(#)(!#/!x), or 2 2 2
!#/!x=1/("0'(#)t+1) et (!#/!t)= -"0(#)/("0'(#)t+1). OK. La queue de la vague s'étalle de plus
en plus tandis que le haut dépasse l'avant. 1.5 1.5 1.5
1 1 1
t3
x
0.5 0.5 0.5
0 0 0
-0.5 -0.5 -0.5
t2 -1 -1 -1
x
-1.5 -1.5 -1.5
t1
x -3 -2 -1 0 1 2 3 4 -3 -2 -1 0 1 2 3 4 -3 -2 -1 0 1 2 3 4
Lorsque le haut dépasse le bas, le choc se produit lorsque l'on résout nos équation (car " est
une fonction de x, la méthode des caractéristique introduit 3 valeurs).
Les équations SW ne sont en réalité plus valides: il se passe des phénomènes à des échelles
courtes. un bricolage; la viscosité qui va lisser le choc. On ajoute un terme de viscosité
Figure 11 – gauche De bas en haut les hauteurs pour différents temps, construction de la solution par conservation le long des caractéristqiuers : on a tracé deux
artificielle en $!2u/!x2:
!" !" !2"
hauteurs (initialement la max et la demi hauteur à droite) qui vont se conserver le long de la courbe caractérisqtique de pente 1/c, le sommet va plus vite (dt/dx = 1/c)
!t
+"
!x
=$
!x2
Cela fait apparaître une échelle courte $.
que la demi hauteur. Raidissement du profil : le haut qui va plus vite dépasse le bas. La solution prédit un haut qui continue et surplombe la partie droite. centre
initial 4 . 3et
. 2 . final du
note sur la film de droite Raidissement et déferlement d’une onde d’inondation, [click to launch the movie, Adobe Reader required]
marée
* Amplification de l'onde de marée dans un bassin.
dans un bassin rectangulaire, de longueur L, en x=0 on donne "0cos(%t)
"="0cos(&x)/cox&L cos(%t); &=%/c.
Il y a donc résonance pour &L='/2... Marées importantes dans la baie de Fundy (Canada) et en
Manche
Cette sous section est l’occasion d’une réflexion sur les dérivées Eq. 3, ou aussi Eq. 13, rappelons que nous avons établi l’équation de conservation de la masse en
partant du bilan sur une tranche fixe x1 et x2 (de longueur ∆x),
Z x2
∂
hdx = Q(x1 , t) − Q(x2 , t)
∂t x1
Rx ∂
que l’on a écrit avec Q(x2 , t) − Q(x1 , t) = [Q]xx21 = x12 ∂x Qdx donc la masse sous forme intégrale (on peut permuter l’intégration et la dérivation car le domaine est
fixe : Z x2
∂ ∂ ∂ ∂
h+ Q dx = 0, vrai pour tous x1 et x2 donc : h + Q = 0.
x1 ∂t ∂x ∂t ∂x
- ESV . 47-
Saint-Venant
Figure 12 – volume de contrôle, cas normal et cas traversé par une discontinuité, q est le Q du texte.
Si la tranche n’est pas fixe, la relation de Leibniz (3) pour h dans une tranche x1 (t) et x2 (t) s’écrit :
Z x2 Z x2
∂t hdx = ∂t h(x, t)dx + h(x2 , t)∂t x2 − h(x1 , t)∂t x1 .
x1 x1
La variation de la quantité intégrée, est égale à l’intégrale de la variation à laquelle on ajoute le flux qui sort h(x2 , t)∂t x2 et le flux qui entre −h(x1 , t)∂t x1 .
On re suppose maintenant que x1 et x2 sont fixes, mais qu’une discontinuité traverse le domaine, à l’instant t est en xc (t) avec x1 < xc < x2 . On va supposer que
la discontinuité se propage à la vitesse W = ∂t xc (ne pas confondre avec la largeur W infinie du canal !). Donc, puisque x1 et x2 sont fixes, et que xc bouge, il est
associé à ce mouvement un flux sortant pour le domaine (x1 , x− +
c ) et un flux rentrant pour le domaine (xc , x2 ) :
Z x2 Z x−
c
Z x2
∂t hdx = ∂t hdx + ∂t hdx
x1 x1 x+
c
En additionnant
Z x2 Z x2 Z x2 Z x2
− −
∂t hdx = ∂t h(x, t)dx − W (h(x+
c , t) − h(xc , t)) que l’on note ∂t hdx = ∂t h(x, t)dx − W [|h|], avec [|h|] = h(x+
c , t) − h(xc , t)
x1 x1 x1 x1
- ESV . 48-
Saint-Venant
R x− x− R x2
chacun des deux termes x1
∂t hdx + [Q]xc1 et x+
∂t hdx + [Q]xx2+ est nul car il traduit la conservation de la masse dans une tranche fixe normale. On en déduit
c
[|Q|]
[| − W h + Q|] = 0 donc la vitesse de la discontinuité est W = .
[|h|]
remarquons que [| − W h + Q|] = 0 s’écrit, puisque Q = uh, aussi [|(−W + u)h|] = 0.
Nous reverrons plus loin ce type de discontinuité dans le cas du ressaut (avec deux équations cette fois, une pour la masse, que l’on vient de voir à l’instant, et
l’autre pour la quantité de mouvement).
Bien entendu, la discontinuité traduit le fait que le détail de la physique n’est pas clair sur une échelle bien plus petite que cette de développement de l’onde de
crue. Dans le paragraphe suivant, nous réintroduisons justement ∂h ∂x qui n’est plus négligeable.
∂F ∂F ∂Q ∂F ∂h
= +
∂t ∂Q ∂t ∂h ∂t
on trouve l’équation de l’onde de crue diffusante
∂Q 3Q ∂Q h3 Cc2 ∂ 2 Q
+ =
∂t 2h ∂x 2Q ∂x2
Le choc en onde longue, se fait sur une distance courte qui est donc réintrduite ici par l’intermédiaire du gradient de pression. La diffusion introduite permet de
lisser le choc et lui donner une épaisseur. Bien entendu le détail du choc n’est encore pas bien décrit....
- ESV . 49-
Saint-Venant
- ESV . 50-
Saint-Venant
2
∂h
Figure 13 – Formation d’une discontinuité dans l’équation de Burgers haut, visqueuse ∂t
+ h ∂h
∂x
= ν ∂∂xh2 , la solution se déplaçant à la vitesse c = (h1 + h2 )/2, passant de h1 à
∂h [|Q|]
h2 sur une épaisseur ν/(h1 − h2 ) , bas inviscide ∂t
+ h ∂h
∂x
= 0 le choc se déplace à la vitesse c = [|h|] = (h1 + h2 )/2, passant de h1 à h2 brusquement
- ESV . 51-
Saint-Venant
δQ (x − U0 t)2
Q = Q0 + √ exp(− )
A 4πDt 4Dt
en supposant les coefficients constants. Voir Chanson [2] p 341 pour d’autres exemples de solution du diffusive wave problem et notamment la fameuse équation de
Cunge-Muskingum.
4
3.5
3
2.5
2
1.5
1
0.5
-3 -2 -1 0 1 2 3 4
Figure 14 – Déplacement et diffusion d’une onde de diffusion [click to launch the movie, Adobe Reader required]
∂V
+ Qout − Qin = 0
∂t
puis, on dit de manière empirique que le volume instantané est fonction des débits : V = kQin − k(1 − )Qout avec k coefficient de stockage et coefficient de
pondération obtenus par essai erreur et calibration. C’est la méthode de Muskingum. Cunge a réinterprété l’équation discrétisée de Muskingum en écrivant l’équation
cinématique discrétisée.
- ESV . 52-
Saint-Venant
Nous allons linéariser les équations de Saint-Venant pour obtenir la propagation des ondes.
Dans ce cas le frottement est négligeable, l’inertie redevient importante, pour simplifier on prend un fond plat. La perturbation se déplace sur de très grandes
distances, l’équation de Saint-Venant
∂h ∂Q ∂Q ∂ Q2 ∂η Q2
+ = 0, + = ghα − gh −g 2 2
∂t ∂x ∂t ∂x h ∂x h CC
α pente faible du fond, η = f + h. Avec f = 0 et sans frottement, en supposant η petit η = εη1 alors on pose u = εu1 (et Q = εu1 h0 ), les équations deviennent
8.7.1 Tsunamis
Une source (tremblement de terre, faille, explosion nucléaire) crée au temps t = 0 une perturbation de l’eau qui va se propager ensuite. Le premier problème est de
déterminer la source : c’est un problème de géophysique, c’est aussi un problème de mécanique car ce sont les équations de la mécanique qui décrivent le mouvement et
la déformation des plaques (sur lesquelles sont les continents) et du magma terrestre. Ce domaine est passionnant mais n’est pas l’objet ici. En fait le fond bouge très
vite : un faille glisse, ou le fond remonte brusquement... Cela déplace la masse d’eau initialement, cette première étape cruciale est en fait assez compliquée puisqu’elle
mêle la description du fond et des premières épaisseurs d’eau qui commencent à se déformer alors que le fond n’a pas fini de se déformer...
Une fois la condition initiale générée, le reste n’est que propagation sur le fond donné. Bien entendu, le bon modèle est Saint-Venant puisqu’il contient (presque)
tous les termes. La théorie la plus simple ensuite est de considérer l’onde linéarisée simple, la perturbation se déplace sur de très grandes distances en eau peu profonde
(relativement à l’étendue du Tsunami), l’équation de Saint-Venant précédente devient la fameuse équation d’onde (ou équation de d’Alembert 1747) :
∂2η ∂2η ∂2 1 ∂2
2
= c20 2 avec c20 = gh0 . Que l’on note parfois à l’aide de l’opérateur appelé d’Alembertien : η = 0, avec = 2
− 2 2
∂t ∂x ∂x c0 ∂t
∂2η 2
2∂ η
− c0 =0
∂t2 ∂x2
√
avec c = gh0 .
La solution de cette équation peut s’obtenir de diverses manières, repartons des équations linéarisées autour d’une grande profondeur h0 , d’un fond plat et du repos
(vitesse nulle) :
∂η ∂u ∂u ∂η ∂(gη) ∂(c0 u) ∂(c0 u) ∂(gη)
+ h0 = 0, = −g ou + c0 = 0, = −c0
∂t ∂x ∂t ∂x ∂t ∂x ∂t ∂x
- ESV . 53-
Saint-Venant
∂ ∂ ∂ ∂
( + c0 )(c0 u + gη) = 0, ( − c0 )(c0 u − gη) = 0
∂t ∂x ∂t ∂x
ce qui permet de dire que c0 u + gη = F (x − c0 t) et c0 u − gη = −G(x + c0 t) ou, la solution s’écrit sous la forme de deux fonctions F (x − c0 t) et G(x + c0 t) :
1 1
u= (F (x − c0 t) − G(x + c0 t)) η = (F (x − c0 t) + G(x + c0 t)).
c0 g
∂2η 2
2∂ η 4
= c0 avec c20 = gh0 .
∂t2 ∂x2 3.5
3
En posant ξ = x − c0 t, et ζ = x + c0 t , on voit que comme 2.5
2
∂ ∂ξ ∂ ∂ζ ∂ ∂ ∂ξ ∂ ∂ζ ∂
= + et = + 1.5
∂x ∂x ∂ξ ∂x ∂ζ ∂t ∂t ∂ξ ∂t ∂ζ
1
∂ ∂ ∂ ∂ ∂ ∂ 0.5
= + et = c0 (− + )
∂x ∂ξ ∂ζ ∂t ∂ξ ∂ζ -3 -2 -1 0 1 2 3 4
∂2η
l’équation devient alors par ce changement de variables 4c20 ∂ζ∂ξ = 0. Figure 15 – Advection à célerité c0 d’une perturbation à 4 temps différents . La
Remarque, on peut ainsi écrire vague se déplace sans changer de forme, ici la solution f (x − c0 t) est tracée ainsi
∂2 2 que quelques droites ”caractéristiques” t = x/c0 + cst).
2 ∂ ∂ ∂ ∂ ∂
− c 0 = ( − c0 )( + c0 )
∂t2 ∂x2 ∂t ∂x ∂t ∂x
2
∂ η
Les solutions de ∂ζ∂ξ = 0 sont par intégration en ξ : ∂η
∂ζ = cst, cette constante
est constante par rapport à ξ mais elle dépend de ζ, appelons la g 0 (ζ), ) puis η = g(ζ) + cst, mais cette constante est une fonction de ξ, appelons la f (ξ). Donc les
solutions sont de la forme : f (ξ) + g(ζ) :
η = f (x − c0 t) + g(x + c0 t)...
La perturbation de forme f se déplace vers la droite sans changer de forme à la vitesse c0 (g se déplace vers la gauche).
Application numérique pour une profondeur de l’océan d’environ h0 = 5km, la vitesse est de 800km/h.
- ESV . 54-
Saint-Venant
Figure 17 – exemple de calcul du Tsunami de décembre 2004. Elévation calculée à 7200s de l’onde rouge : +2 m bleu -2. A droite observation du Tsunami par le
satellite Jason-1 comparé et calcul de la simulation. [Link]
Remarque, un Mile est une minute d’angle (60 minutes = 1 degrés, le tour
de la terre est 360 degrés, et la circonférence de la terre est 2πRterre ) : la
longeur du Mille marin est 6378000 ∗ 2 ∗ π/(360 ∗ 60) = 1852m la distance
mesurée au bout de 5 heures est de 32 degrés, donc 32*60=1920 Miles, soit
en 5 heures : 1920/5= 384 noeuds, 1 noeud = 0.5m/s, on retrouve bien les
200m/s.
- ESV . 55-
Saint-Venant
Figure 18 – exemple de calcul du Tsunami du 11 mars 2011. Elévation calculée à 7200s de l’onde rouge : 0.5 m bleu -0.5m. A droite observations du Tsunami par
les bouées DART dans le pacifique, calcul G erris. [Link]
Le calcul complet des Tsunamis nécessite non seulement des calculs de Saint-Venant non linéaires avec frottement, mais surtout la connaissance fine de la forme
des fonds. Ce sont de grands effets de variation de fond qui déforment le signal initial. Les réflexions sur les côtes produisent des signaux qui paraissent dispersifs (la
perturbation initiale est suivie de nombreuses vagues qui correspondent aux réflexions, et non pas le signal de l’accident initial ; voir la figure avec la comparaison avec
le satellite Jason-1 en 2004).
Plus récemment, on a enrichi (particulièrement dans Basilisk) les calculs de tsunamis en rajoutant les termes dispersifs (équations de Boussineq ou de Serre Green
Naghdi). Ce sont les termes que nous verrons pour le Mascaret.
- ESV . 56-
Saint-Venant
Jusqu’à présent les termes de dérivée en temps et en espace étaient de même ordre dans la dérivée totale (Strouhal unitaire). Dans le cas de la marée, il s’agit d’un
forçage extérieur avec une échelle de temps imposée. On adimensionne le temps avec la période imposée par l’alternance des marées, on obtient donc en continuant de
négliger les termes non linéaires, et le frottement au fond ainsi que la pente moyenne (qui n’intervient plus puisque l’on est à l’échelle d’une mer comme la Manche)
∂η ∂Q ∂Q ∂η
+ = 0; et = −gh
∂t ∂x ∂t ∂x
donc, comme la variation de hauteur de vague est faible, l’équation de
∂2η ∂ ∂η
= (g(h0 (x)) )
∂t2 ∂x ∂x
c’est l’équation de type tidal wave, onde de marée (aussi long gravity wave.
Il faut la résoudre connaissant le fond, en fait, h0 (x) désigne la hauteur d’eau compte tenu de la forme du fond f (la bathygraphie), En conditions aux limites, on se
donne la marée océanique au loin provoquée par la lune et le soleil, donc une superposition de ces deux modes. Rappel sur les marées : il y en a deux par jour, liées
au jour lunaire 24h50min28s qui donne deux marées plus la contribution du soleil toutes les 12h00.
• La baye de Fundy, ou ”baie des Français” (Landau p 40, ENS StCLoud 80, Ellipses, Walker 4.55
Cet exemple est peu général, mais se fait à la main, il traite de l’endroit où les marées sont les plus grandes au monde (avant la Manche). Il s’agit en première
approximation d’une onde linéaire simple en fond peu profond, la forme de la baie est assez rectangulaire, on l’approxime à un rectangle et on fait du 1D. On se donne
(et c’est une approximation, la hauteur à l’entrée de la baie η0 e−iωt , donc un seul mode en première approximation). La marée importante s’explique par le fait que
- ESV . 57-
Saint-Venant
cos(kx − L) −iωt p
η = η0 e avec k = ω/ hg,
coskL
√
la marée est extrémale au fond de la baie, si k = (2n+1)π/2 on a résonance. A.N. η0 =2.2m, période 12h25 min, profondeur de la baie h = 73m, donc c0 = hg = 27m/s,
on a k = 5.2610−6 m, on trouve η(L, t) = 14.4m..
Ce cas correspond à une résonance due à la forme rectangulaire de la baie, dans le cas de la marée à Saint Malo, à une éventuelle résonance due au fait que l’on se
trouve dans une cavité coincée au creux de la Bretagne et de la Normandie, l’effet s’additionne à celui de la force de Coriolis.
baie
L
276km
mer
Figure 20 – La baie de Fundy au Canada est connue pour avoir les marées les plus hautes du monde.
• A titre d’exemple de résolution des équations sous un forçage périodique on peut traiter ici le cas de faible profondeur linéarisé en 2D et en tenant compte de la
force de Coriolis, avec le paramètre de Coriolis fc = 2ω sin λ. On va résoudre dans le cas de l’écoulement dans un canal infini en x et borné en largeur y (largeur W ,
ce qui peut être le modèle de l’écoulement dans la Manche qui est longue et assez rectangulaire).
Nous allons voir ”l’Amphidromie de Kelvin” : on fait apparaı̂tre par la résolution les points de hauteur de marée nulle, ce sont les points amphidromiques cf Hervouet
[12]. Ce cas est en fait générique. On se donne un canal de profondeur h0 et de largeur W , les équations linéarisées autour de h0 et du repos, en négligeant les frottements,
mais en mettant la force de Coriolis :
∂η ∂u ∂v
+ + =0
∂t ∂x ∂y
∂u ∂η ∂v ∂η
− fc v = −g , + fc u = −g ,
∂t ∂x ∂t ∂y
cherchons des solutions de la forme (σ pulsation) :
F = F (x, y)eiσt ,
on obtient la solution pour le champ de vitesses
−g ∂η ∂η g ∂η ∂η
u= (iσ + fc ), & v = 2 (−fc + iσ )
fc2 − σ 2 ∂x ∂y fc − σ 2 ∂x ∂y
- ESV . 58-
Saint-Venant
Figure 21 – Onde de Kelvin se propageant dans un canal vers la droite [click Figure 22 – 1 point amphidromique de marée nulle [click to launch the movie,
to launch the movie, Adobe Reader required]. Adobe Reader required].
On appelle onde de Kelvin la solution particulière, se propageant ici vers la droite (dans le cas d’une largeur infinie a priori, mais en fait comme v = 0, il y a glissement
automatique sur al paroi supérieure) : r
g − √fc y iσ(t− √ x )
u= η, v = 0, η = η0 e gh0 e gh0
.
h0
L’onde se propage le long d’une des frontières, (vers le pôle pour une frontière est, vers l’équateur pour une frontière ouest). Le maximum est sur la frontière.
L’amplitude varie de manière√ exponentielle d’un bord à l’autre du canal. La marée est plus haute d’un côté que de l’autre. L’onde décroı̂t vers le large avec une
longueur caractéristique gh0 /fc qui s’appelle rayon de Rossby. Les lignes d’égale amplitude sont parallèles aux bords du canal et les lignes cotidales (x =constante)
sont normales à l’axe du canal.
Les ondes de marées se réfléchissant et se réfractant sur les côtes, superposons deux ondes ainsi calculées se déplaçant en sens inverse dans le canal. La surélévation
résultante s’écrit alors :
( −iσx−f
√ cy ) ( iσx+f
√ cy )
η = eiσt [η1 e gh0
+ η2 e gh0
].
soit a = η2 /η1 qui caractérise le rapport des amplitudes des deux ondes. C’est vrai pour toute largeur de canal car v = 0. Sur chaque bord, l’onde se propage en sens
inverse. L’amplitude de la marée : r
2fc y 4fc y
√ p √ p
|η| = (1 + a2 e gh0 cos(2σx/ gh0 ))2 + a2 e gh0 sin(2σx/ gh0 ))2
√
Cette
√ expression nous montre aussi qu’il existe des points à marnage nul (d’amplitude nulle), ils sont situés sur la droite y = gh0 ln(1/a)/(2fc ), séparés de
x = gh0 (2k + 1)π/(2σ),
En ces points où la marée est toujours nulle, la mer est pleine (ou basse) à chaque instant ; toutes les lignes cotidales passent par ces points autour desquels elles
paraissent tourner : ce sont les points amphidromiques.
- ESV . 59-
Saint-Venant
Figure 23 – Calcul ”TELEMAC” d’EDF des courants de marées [click to launch the movie, Adobe Reader required].
les courants ont été mesurés par des bouées. Et les hauteurs par des marégraphes. La simulation numérique permet de tracer des champs continus utiles aux navigateurs.
Un exemple de mise en oeuvre est sur [Link] on y utilise les équations de Saint-Venant complètes
avec Coriolis et frottement au fond dans la bathymétrie de la Manche.
- ESV . 60-
Saint-Venant
Figure 24 – La marée est plus élevée en France de par la force de Coriolis, document SHOM. A droite points Amphidromiques, image issue de
[Link] document SHOM ?
- ESV . 61-
Saint-Venant
Quatrième partie
Le Fluide : Equations de Saint-Venant à frottement faible, échelle
assez courte, quelques exemples de solutions
Nous avons vu que les équations de Saint-Venant contenaient un terme moteur lié à la pente α. Dans le cas des fleuves, ce terme est équilibré par le terme de
frottement de Chézy. Cette description équilibrant pente et frottement est pertinente pour les problèmes caractérisés par une variation lente longitudinale. Dans cette
partie, on étudie plutôt les écoulements à variation rapide (mais pas trop pour rester dans l’hypothèse couche mince de Saint-Venant) qui sont caractérisés par un
frottement faible dans lequel on peut négliger la pente (voir §8.2 pour la discussion générale). On s’intéresse en plus à des écoulements établis, stationnaires en temps
dans une première partie pour introduire la notion de ”charge”. On examine ensuite les cas instationnaires pour introduire ”la rupture de barrage” et le ”ressaut”.
h f u2 u2 u2 u2
=1− + 0 (1 − 2 ) ou encore h+f + = h0 + 0 .
h0 h0 2gh0 u0 2g 2g
2
Notons qie h + u2g est classiquement appelée ”la charge spécifique” en hydrodynamique. Comme h/h0 = u0 /u, on écrit une relation implicite entre f et u. Connaissant
u20
la vitesse u et F 2 = gh0 on trouve la topographie
f F2 u2 u0
=1+ (1 − 2 )) −
h0 2 u0 u
Cette relation est utilisée en pratique pour mesurer la topographie en fonction de la vitesse de surface mesurée, mais elle est délicate car il ne faut pas que la vitesse
soit trop faible.
- ESV . 62-
Saint-Venant
u2 Q2 Q2
Hs = h + =h+ , on a dHs /dh = 1 − ,
2g 2gh2 gh3
On voit que pour h → ∞ et pour h → 0 on a Hs → ∞. Plus généralement, pour deux hauteurs différentes, nous avons le même Hs , mais il existe une valeur minimale,
2 2
appelée valeur critique hc ; telle que Hs est minimale. On a dHs /dh = 0 telle que hc = ( Qg )1/3 et Hsc = 23 hc = 32 ( Qg )1/3
√u Q
Le nombre de Froude F = gh
c’est aussi √
gh3/2
. On a alors
Q dHs
F =√ et = 1 − F 2.
gh3/2 dh
2
Le nombre F est tel que en h = hc on a Fc = √g(hQc )3/2 = ( g(h
Q
c)
3)
1/2
= 1. Le froude est unité à la valeur critique. Pour une valeur de la charge spécifique, il existe
donc deux valeurs possibles pour la hauteur, une valeur petite à nombre de Froude supérieur à l’unité et une valeur plus haute à Froude inférieur à l’unité.
3/2
u2 F2
hc Hs h
Hs = h + s’écrit aussi puisque F = sous la forme = (1 + ).
2g h hc hc 2
En pratique, on appelle ”écoulement graduellement variés” des écoulements sur une échelle suffisamment longue pour que l’on soit toujours à l’équilibre. On peut
alors analyser l’écoulement avec la notion de charge.
- ESV . 63-
Saint-Venant
— Si F 2 (x) > 1 les variations de f et de u sont de sens opposé, une augmentation du fond se traduit par une décélération.
— Si F 2 (x) < 1 les variations de f et de u sont de même sens, une augmentation du fond se traduit par une accélération. Ce résultat généralise le résultat linéaire
précédent.
— Si F 2 (x) 6= 1, la vitesse est extrémale là où la bosse atteint sont sommet.
— Si F 2 (x) = 1 au sommet, la vitesse continue d’augmenter, ∂u 2
∂x ne change pas de signe après le sommet, mais F (x) − 1 a changé de signe, donc la bosse peut
∂f
décroı̂tre ∂x < 0.
Le passage sub critique supercritique peut se produire exactement au sommet d’une bosse. On a alors la figure à droite ou la figure . La vitesse est accélérée, la
surface se creuse. Au col l’écoulement est critique F = 1, après le col il est supercritique. La surface ne remonte pas.
l’interprétation physique de ces calculs est liée au fait que des ondes (des vagues) se déplacent sur la surface et transmettent ainsi des informations à l’écoulement
en amont et en aval. La vitesse des ondes dépend de la hauteur d’eau.
Dans le cas subcritique, la vitesse du courant est inférieure à la vitesse des ondes (u < c0 , F < 1), les ondes peuvent remonter le courant et ”prévenir” l’eau en
amont qu’il va y avoir bientôt une bosse et qu’il faut adapter la vitesse : accélérer pour conserver le débit. Par Bernoulli, cette accélération provoque une dépression,
le niveau de l’eau baisse au dessus de la bosse.
Dans le cas supercritique, la vitesse du courant est supérieure à la vitesse des ondes (U > c0 , F > 1), les ondes ne peuvent plus remonter le courant et ne
”préviennent” plus l’eau en amont qu’il va y avoir bientôt une bosse et qu’il faut adapter la vitesse : ”surprise” par la bosse, l’eau s’accumule devant elle et est
brusquement freinée. Par la loi de Bernoulli, cette décélération provoque une surpression, le niveau de l’eau augmente au dessus de la bosse.
Le cas subtil est celui subcritique tel que l’accélération au niveau du sommet est telle qu’elle permet à l’écoulement de devenir supercritique.
Bien entendu, cette relation est comparable à la relation dans les tuyères, avec le nombre de Mach.
f¯
f¯ − F 2 ū1 + ū1 = 0, donc ū1 = .
1 − F2
- ESV . 64-
Saint-Venant
Ou, si on repasse en valeur dimensionnée, mais sans oublier que la perturbation du fond est faible :
(f /h0 )
u = u0 (1 + ).
1 − F2
2
F 2 f¯
Bernoulli sans dimension (F 2 (u/u2 0 ) + η/h0 ) = (F 2 12 + 1), linéarisé devient F 2 ū1 + η1 = 0, on peut ainsi exprimer la déviation relative η1 = 1−F 2 position de la surface
libre
F 2f
η = h0 + 2 .
F −1
ou de la hauteur d’eau :
F 2f f
h = h0 + 2 − f = h0 + 2 .
F −1 F −1
Suivant le régime, la réponse de la hauteur d’eau (supercritique) sera une élévation ou un creusement (subcritique).
∂f
Si on part de u1 ∂u 2
∂x (F (x) − 1) + ∂x = 0, on retrouve bien sûr les mêmes résultats.
9.6 Appilcations
Exprimons l’expression de la déviation de la surface libre en fonction de la forme de la topographie et la varaition de vitesse, respectivement :
F 2f f
resp u0 .
F2 − 1 h0 (1 − F 2 )
Pour un régime subcritique, F < 1, on constate que la surface libre de l’eau est creusée et que le courant est accélérée. Inversement, pour un écoulement supercritique
la surface libre est déviée vers le haut et le fluide est ralenti. La figure 26 montre à gauche un écoulement subcritique et la figure 27 montre un écoulement supercritique.
Figure 26 – A gauche F < 1 la surface libre de l’eau est creusée. à Droite F < 1 à F > 1 au col, puis ressaut, Photo issue du livre de Hulin Guyon Petit [10]
De exemples de résolution numérique Basilisk à différents Froude pour l’écoulement sur une bosse sont dans
[Link]
- ESV . 65-
Saint-Venant
Figure 27 – Ecoulement supercritique F > 1, la surface libre est déviée vers le haut et le fluide est ralenti. Cuve de l’ENSTA Palaiseau, photo PYL
dh(x) h df (x)
= 2 .
dx Q (x)/(gh(x)3 ) − 1 dx
On peut aussi rajouter le frottement dans cette équation qui permet de calculer la hauteur par intégration. Remarquons que dans le cas où le fond est faible, le vitesse
varie peu Q2 (x)/(gh(x)3 ) est à peu près constant c’est F 2 . et on peut intégrer h = h0 + F 2f−1 , on retrouve bien la solution linéarisée.
– la courbe débit / h
les courbes classiques etc etcx
- ESV . 66-
Saint-Venant
10 Caractéristiques et chocs
10.1 Encore et encore ∂’Alembert ?
Reprenons le système de base que nous allons linéariser pour en obtenir un système simple (menant à l’équation des ondes)
introductif au cas non linéaire. Cherchons les invariants de l’équation linéaire pour u(x, t) et h(x, t).
∂h ∂(uh) ∂u ∂ u2 ∂h
+ = 0 et + = −g
∂t ∂x ∂t ∂x 2 ∂x
linéarisons la vitesse et
√ la hauteur autour de l’état de base (u, h) = (0, h0 ), donc u = 0 + εu1 + ... et h = h0 + εh1 + ... etc bien
entendu on pose c0 = gh0 , on a donc au premier ordre
∂h1 ∂(u1 ) ∂u1 ∂h1 ∂ h1 ∂u1 ∂u1 ∂ h1
+ h0 = 0 et = −g on ré arrange en [c0 ] + c0 = 0 et + c0 [c0 ] = 0
∂t ∂x ∂t ∂x ∂t h0 ∂x ∂t ∂x h0
on peut les additionner et les soustraire
∂ h1 ∂ h1 ∂ h1 ∂ h1
[u1 + c0 ] + c0 [u1 + c0 ] = 0 et [u1 − c0 ] − c0 [u1 − c0 ] = 0
∂t h0 ∂x h0 ∂t h0 ∂x h0 Figure 28 – Caractéristiques, ici
des droites. Le long dx/dt = c0
On rappelle qu’en posant ξ = x − c0 t, et ζ = x + c0 t , on voit que comme
on a u1 + c0 hh10 constant, Le long
∂ ∂ξ ∂ ∂ζ ∂ ∂ ∂ξ ∂ ∂ζ ∂ ∂ ∂ ∂ ∂ ∂ ∂ dx/dt = −c0 on a u1 − c0 hh10
= + et = + et donc = + et = c0 (− + )
∂x ∂x ∂ξ ∂x ∂ζ ∂t ∂t ∂ξ ∂t ∂ζ ∂x ∂ξ ∂ζ ∂t ∂ξ ∂ζ constant.
ce qui donne pour ∂t ± c0 ∂x les expressions suivantes :
∂ ∂ ∂ ∂ ∂ ∂
+ c0 = 2c0 puis − c0 = −2c0
∂t ∂x ∂ζ ∂t ∂x ∂ξ
∂ h1 ∂ h1 ∂ h1 ∂ h1
[u1 ± c0 ] ± c0 [u1 ± c0 ] = 0 deviennent 2c0 [u1 + c0 ] = 0 et 2c0 [u1 − c0 ] = 0
∂t h0 ∂x h0 ∂ζ h0 ∂ξ h0
h1 2c0
- la première équation porte sur [u1 + c0 ] qui ne varie pas en ζ = x + c0 t, donc une fonction ξ = x − c0 t donc manifestement (le coefficient est introduit pour
h0 h0
simplifier) :
h1 2c0 dx
[u1 + c0 ]= f (x − c0 t) constant le long de = c0
h0 h0 dt
h1
- la seconde équation porte sur [u1 − c0 ] qui est indépendante de ξ, c’est une fonction de ζ on écrit de même :
h0
h1 −2c0 dx
[u1 − c0 ]= g(x + c0 t) constant le long de = −c0
h0 h0 dt
En recombinant, on a
c0
h1 (x, t) = f (x − c0 t) + g(x + c0 t) et u1 (x, t) =(f (x − c0 t) − g(x + c0 t))
h0
le cas linéarisé est très simple et reconduit à d’Alembert, on retrouve ainsi les deux solutions en x ± c0 t.
- ESV . 67-
Saint-Venant
h1 2c0
• Par exemple, si une onde se propage sur la droite exclusivement, c’est que g est identiquement nul, u1 − c0 hh01 = 0. Il ne reste que u1 + c0 = f (x − c0 t) donc que
h0 h0
h1
u1 = c0 , ou u1 = c0 f (x − c0 t), et h1 = h0 f (x − c0 t).
h0
• Si au temps t = 0, on se donne u1 = U (x) et h1 = H(x), on peut résoudre à tout temps car on en déduit
1 h0 U 1 −h0 U
f= + H et g = +H .
2 c0 2 c0
Nous allons reprendre ce calcul linéarisé en non linéaire sur les équations non linéaires. Nous allons le faire de deux manières et montrer que ce résultat est plus
général. Mais la compréhension de ce cas linéarisé est fondamentale.
∂h ∂(uh) ∂u ∂ u2 ∂h
+ = 0, et + = −g
∂t ∂x ∂t ∂x 2 ∂x
est de la forme
∂U ∂F (U ) ∂U ∂F ∂U
+ = 0, ou + · = 0.
∂t ∂x ∂t ∂U ∂x
√
h uh ∂F u h
avec U = et F (U ) = u2 et A = ∂U = . La matrice A est diagonalisable, ses valeurs propres sont λ1,2 = u ± gh, soit
u 2 + gh
g u
√
u − gh 0√
Λ= ,
0 u + gh
p p
−1 − h/g h/g
ses vecteurs propres Vi=1,2 tels que : A · Vi = λi Vi et V = R · U , où la matrice de changement de base R est telle que l’on a R = , et telle que
1 1
p
− g/h/2 1/2
R−1 = p et R · R−1 = I. On a R · Λ · R−1 = A et A · R = R · Λ ainsi que R−1 · A = Λ · R−1 . Dans le cas d’une équation linéaire, en multipliant à
g/h/2 1/2
gauche par R−1 l’équation matricielle, et en remplaçant R−1 · A par : Λ · R−1
∂U ∂U
R−1 · + Λ · R−1 =0
∂t ∂x
On obtient donc la forme suivante :
∂V ∂V
+Λ =0
∂t ∂x
la solution est donc de la forme Vi (x − λi t)
Vérifions
p dans ce cas non linéaire si par hasard
√ ça marche aussi. Calculons R−1 · ∂t U : sur la première ligne :
on a − g/h/2∂t h + (1/2)∂t u = (1/2)∂t (u − 2 gh), etc pour la seconde ligne, et la dérivée ∂x donnera une expression similaire, donc après calculs, on trouve la forme
suivante :
∂ p ∂ p ∂ p ∂ p
[ + (u + gh) ](u + 2 gh) = 0, [ + (u − gh) ](u − 2 gh) = 0
∂t ∂x ∂t ∂x
- ESV . 68-
Saint-Venant
dx
dx
=u−c C+ C− dt
=u−c
dt
t C− t dx
=u−c C+
dx
dt
=u+c dx
=u+c
dt
t C−
dt dx
=u+c
dt
C+
x x
Figure 29 – à gauche régime fluvial, la vitesse est plus petite que la vitesse des x
ondes ; l’écoulement en un point dépend de l’amont et de l’aval. A droite ; régime
torrentiel, les ondes vont moins vite que le courant ; l’écoulement en un point ne Figure 30 – Des caractéristiques se croisent : il y a un choc, voir §10.6.
dépend que de l’amont, il n’influence que l’aval. u + 2c constant le long des C +
et u − 2c constant le long des C −
La fonction
√ qui est constante le long de la caractéristique est appelée invariant
√ de Riemann. Pour les ondes non linéaires on a conservation des invariants de Riemann
(u ± 2 gh) le long des caractéristiques aux vitesses caractéristiques (u ± gh).
On a : √
d dx
dt (u + 2√gh) = 0 le long des courbes C+ d’équation dt =u+c
d dx
dt (u − 2 gh) = 0 le long des courbes C− d’équation dt =u−c
- ESV . 69-
Saint-Venant
Remarque
Nous faisons le lien ici avec l’équation des ondes. Si on linéarise autour d’une hauteur d’eau fixe et d’une vitesse nulle : u = 0 + εu1 + ... et h = h0 + εh1 + ... donc
h1
c = c0 (1 + εh1 /h0 + ...)1/2 = c0 + ε + ...
2h0
dx
• u + 2c = 2c0 + ε(u1 + c0 h1 /h0 ) + ... est constant le long de la courbe = c0 ,
dt
h1 dx
on retrouve ce que l’on avait en linéaire : [u1 + c0 ] fonction de (x − c0 t) constant le long de = c0
h0 dt
dx
• u − 2c = −2c0 + ε(u1 − c0 h1 /h0 ) + ... est constant le long de la courbe = −c0 ,
dt
h1 dx
on retrouve ce que l’on avait en linéaire : [u1 − c0 ] fonction de (x + c0 t) constant le long de = −c0
h0 dt
Une mise en oeuvre de ces solutions est le problème de rupture de barrage. Ce genre d’accidents est rare, un cas tristement célèbre est celui de la rupture du
barrage de Malpasset (Var) le 2/12/59. Ce jour là, la rupture du barrage fait déferler une vague de 5 m de haut à la vitesse de 70 km/h (voir le film du lien suivant à
- ESV . 70-
Saint-Venant
dx
=u−c
dt
C− t C+
η(x, t) dx
=u+c
dt
Figure 31 – Rupture x
de √barrage, entre √x =
−t gh0 et x = 2t√ gh0 , Figure 32 – Les
u = (2/3)(x/t +√ gh0 ) caractéristiques pour
et gh = ((2/3)( gh0 − l’onde simple centrée Figure 34 – Calcul EDF [14]
x/(2t))2 . lors de la rupture de de la rupture du barrage. Ce cas
barrage. u + 2c constant est un cas test numérique car la
Figure√33 – Rupture√de barrage, entre
le long des C+ et u − 2c rupture de chaque ligne électrique
x = −t gh0 et x = 2t gh0 .
constant le long des C− croisée par l’inondation a été re-
levée temporellement.
√
la minute 04 :41)) (effectivement 3.6 5 ∗ 9.81 = 50). Balayant tout sur son passage sur une distance d’une douzaine de kilomètres, elle débouche sur Fréjus 20 minutes
plus tard, avant de se jeter dans la mer faisant 423 victimes.
- ESV . 71-
On vérifie !u/!t="0'(#)(!#/!t) et !u/!x="0'(#)(!#/!x), or !" !" !2" Saint-Venant
+ " - # =0
!t !x !x2
!#/!x=1/("0'(#)t+1) et (!#/!t)= -"0(#)/("0'(#)t+1). OK. La queue de la vague s'étalle de plus Si maintenant on ne met pas de dispersion mais uniquement des termes non linéaires.
et en s'atténuant, le choc est étalé sur quelques points
exemple propagation d'un pulse donné au temps t=0:
en plus tandis que le haut dépasse l'avant. [Link]
t3
x il se déplace en se raidissant
un soliton:
[Link]
!" !" !2"
+ " - # =0
8. transport de sédiments.... !t !x !x 2
On lira le livre
Si maintenant onde
neYalin...
met pas On consulteramais
de dispersion les documents
uniquement EDF. L'idée
des termes à retenir,
non est que les
linéaires.
sédiments sont emportés
exemple propagation comme
d'un des scalaires
pulse donné passifs
au temps t=0: par la couche limite turbulente...
...
et en s'atténuant, le choc est étalé sur quelques points
9. films visqueux
t2 [Link]
9.1. ressaut très visqueux
il se déplace en se raidissant
Lorsque le haut dépasse le bas, le choc se produit lorsque l'on résout nos équation (car " est duit sur une échelle courte et qui met en jeu des phénomènes physiques
8. transport de sédiments....
Figure 35 – De bas en haut les hauteurs de vague pour différents temps. Rai- absents
On lira le deslivre
équations
de Yalin... de Saint-Venant.
On consultera L’enroulement
les documents
et en s'atténuant, le choc est étalé sur quelques points
deà retenir,
EDF. L'idée la surface libre
est que les
une fonction de x, la méthode des caractéristique introduit 3 valeurs). sédiments sont emportés comme des scalaires passifs par la couche limite turbulente...
dissement du profil : le haut qui va plus vite dépasse le bas. la solution prédit un sur elle même est remplacée par une discontinuité.
[Link]
... En pratique, cette dis-
Les équations SW ne sont en réalité plus valides: il se passe des phénomènes à des échelles continuité estvisqueux
lissée sur quelques points ∆x par les schémas numériques,
enroulement. Cet enroulement est détruit par rupture de la surface en gouttes 9. films
etcourtes.
bulles un
: il bricolage;
s’agit en lafait
viscosité qui va lisser le choc. On ajoute un terme de viscosité
du déferlement. 9 . 1 . ressaut très visqueux
[Link]
artificielle en $!2u/!x2: pyl@[Link] ondes fev 98 page 64 le 18 Octobre 1999 à 14:57
4 . 3 . 2 . note
Alternative : sur la marée 9. films visqueux
* Amplification
Bien de l'ondeon
entendu, comme de marée dans undes
a considéré bassin.
ondes allant vers la droite, avec u − 2c constant9 .u1 −
. ressaut
2c = −2c0très visqueux
partout, on aurait aussi pu partir de l’autre équation :
pyl@[Link] ondes fev 98 page 64 le 18 Octobre 1999 à 14:57
dans un bassin rectangulaire, de longueur L, en x=0 on donne "0cos(%t)
∂ ∂
"="0cos(&x)/cox&L cos(%t); &=%/c. ( + (u + c) )(u + 2c) = 0
∂t ∂x
Il y a donc résonance pour &L='/2... Marées importantes dans la baie de Fundy (Canada) et en
qui donne pusique u = 2c − 2c0
Manche
∂ ∂
( + (3c − 2c0 ) )(4c + 2c0 ) = 0
pyl@[Link] ondes fev 98 page 27 le 18 Octobre 1999 à 14:57 ∂t ∂x
et comme c(h), donc on retrouve bien l’équation précédente.
Note historique
Ce calcul a été fait par J. Proudman en 1957
q ”On the series that represent tides and surges in an estuary” avec un point de vue proche, avec en plus un contre courant
−u0 . Il trouve que la vitesse est alors ((3 1 + hη0 − 2 − uc00 )c0 ) (exercice : vérifiez le).
Mais en fait il a été proposé par Adhémar Jean-Claude Barré de Saint-Venant lui même en 1871 dans la première partie de son article de 1871 Le fichier est ici
[Link]
- ESV . 72-
Saint-Venant
∂ η̄ 3 ∂ η̄
+ (1 + ε η̄) =0
∂ t̄ 2 ∂ x̄
Plaçons nous dans le repère qui se déplace ξ = x − t et utilisons pour le temps la variable τ = t/ε, c’est un temps long : ∂∂x̄ = ∂ξ
∂
et ∂
∂ t̄
∂
= − ∂ξ ∂
+ ε ∂τ donc l’équation
aux temps longs et dans le repère de l’onde qui se déplace devient l’équation dite de Bürgers (on enlève la barre au dessus du η) :
∂η 3 ∂η
+ η = 0.
∂τ 2 ∂ξ
Cette équation non linéaire a tendance à créer un choc car les vagues de plus grande hauteur rattrapent les plus basses.
Équation de Bürgers
Dans les équations de Navier Stokes, on a toujours négligé le terme de Rdérivée seconde (∂x2 ) en prétendant qu’il était négligeable. Une autre difficulté est aussi que
l’on ne peut pas écrire simplement l’intégrale moyenne sur l’épaisseur de ∂x2 u1 dz en fonction ∂x2 (hu). Cependant, ce terme est parfois rajouté de manière ad hoc. Il
permet justement de lisser les chocs, sur une distance peu physique, mais utile en pratique :
∂η 3 ∂η ∂2η
+ η =ν 2
∂τ 2 ∂ξ ∂ξ
Pour résoudre cette équation (équation de Bürgers) on fait la transformation dite de Hopf en posant η = −2ν ∂log(F
∂ξ
)
d’où
∞
∂2F (x − ζ)2
Z
∂F 1
− ν 2 = g(t)F, On obtient : F (x, t) = exp(− )F (ζ, 0)dζ.
∂τ ∂ξ 4πνt −∞ 4νt
Une autre solution de cette équation, en posant ζ = (ξ − cτ )/ν pour faire disparaı̂tre la viscosité artificielle et en cherchant une solution F (ζ) :−cF 0 + F F 0 + F 00 =
∂2h
0, soit − cF + F 2 /2 + F 0 = cst on a vu déjà que ∂h ∂h
∂t + h ∂x = ν ∂x2 s’il y a mouvement la solution est de la forme :
h2 + h1 h2 + h1 h2 − h1 (x − ct)
c= , η= + tanh((h1 − h2 ) )
2 2 2 4ν
ce qui donne la solution se déplaçant à la vitesse c = (h1 + h2 )/2, passant de h1 à h2 sur une épaisseur ν/(h1 − h2 ) :
- ESV . 73-
Saint-Venant
11 Ressaut
11.1 Système
Le système qui nous occupe est toujours le système de Saint-Venant sur fond plat sans frottement. Rappelons que nous l’avons établi avec la démarche de
conservation sur des petites tranches. On a donc sous établi sous la forme ”conservative” en 1D
∂h + ∂(hu)
=0
∂t ∂x 2
∂(hu) + ∂ (hu2 + g h ) = 0.
∂t ∂x 2
Les équations de Saint-Venant s’écrivent en 2D
∂h ∂(hu) ∂(hv)
+ + =0
∂t ∂x ∂y
h2
∂(hu) ∂ 2 ∂(huv) ∂U ∂Fx ∂Fy
+ (hu + g ) + =0 de la forme + + =0
∂t ∂x 2 ∂y ∂t ∂x ∂y
h2
∂(hv) ∂(huv)
∂
+ + (hv 2 + g ) =0
∂t ∂x ∂y 2
En fait, pour établir les équations 1D, nous avons fait un bilan sur une tranche de longueur ∆x que l’on a fait tendre vers 0. Mais nous aurions pu nous arrêter à la
tranche sans la faire tednre vers 0. Voire, considérer un domaine quelconque, pour l’écrire sur un domaine (D de frontière ∂D) qui n’est pas forcément infinitésimal :
h2
Z Z Z
d d
ρhdS = 0 et ρuhdS + ρg d` = 0.
dt D dt D ∂D 2
Nous allons montrer que ce système permet d’écrire des relations de part et d’autre d’une discontinuité. Cette ”discontinuité” est le ressaut. Le vrai ressaut physique
a une certaine longueur, mais on suppose l’extension suffisamment faible que vu de loin, c’est une variation brusque des valeurs. On a vu plus haut comment se forme
un choc et que sa structure dépend de la viscosité. En fait, la description est plus compliquée car l’écoulement devient turbulent et donc le ressaut est plus épais, mais
encore assez mince.
- ESV . 74-
Saint-Venant
Figure 37 – Discontinuité traversant un volume, notations pour la dérivation d’intégrales, ne pas confondre V le volume et v la vitesse relative.
S1 et Σ entourent V1 , le flux est pris sur la surface S1 et sur Σ, on note C1 la valeurs prises sur σ en venant de l’intérieur de V1 . On obtient la même expression pour
−−→ −
V2 , mais on va faire apparaı̂tre dSΣ = → n Σ dSΣ car la normale à Σ en venant de 1 est opposée à celle venant de 2 :
−
→ −−→
Z
C2 W · (−dSΣ )
Σ
→ −−→
−
car on garde la normale précédente et que l’on suppose que le choix est tel que W · dSΣ > 0, la noramle est dirigée de l’extérieur de 1 vers 2. On regroupe
−
→ → −−→
−
Z Z Z Z Z Z
d d d ∂ →
−
Cdv = Cdv + Cdv = ( C)dv + C u · dS + (C1 − C2 )W · dSΣ
dt V dt V1 dt V2 V ∂t S Σ
−
→ −
on utilise souvent la notation de saut |C| = C2 − C1 et W = W · →
n Σ , el la vitesse relative, la relation fondamentale de saut est donc :
−
→
Z Z Z Z
d ∂ →
−
Cdv = ( C)dv + C u · dS − |C| W dSΣ
dt V V ∂t S Σ
il y a trois contributions ; les variations locales à l’intérieur du volume figé, le flux lié à la convection à la surface S et les variations dues à la discontinuité.
- ESV . 75-
Saint-Venant
−
→ −
→
On peut réécrire cette dernière relation, en faisant intervenir la vitesse relative →
−
v =→ −u − W , et en réécrivant la divergence : en évaluant Σ C →
−
R
u · dS avant la
discontinuité, cela permet d’écrire un flux pour 1 (puis pour 2), avec les vitesses matérielles et
→
− −
→ −−→
Z Z Z
→
−
∇ · (C u )dv = →
−
C u · dS + (C → −u )1 · dSΣ
V1 S1 Σ
−
→
idem avec 2, une autre forme utile de dérivation des intégrales est donc (→
−
v =→
−
u − W ) au final
→
−
Z Z Z
d ∂ →
−
Cdv = ( C) + ∇ · (C u ) dv + |Cv| dSΣ
dt V V ∂t Σ
Cette forme fait intervenir la vitesse relative, et les équations sous forme locale dans le volume (ce qui va nous servir par la suite).
On a besoin aussi d’un terme sur la périphérie, Z Z
d
Cdv + KdS = 0
dt V S
R −→ R −→ R −
→
Le second terme S K dS va se décomposer en S1 K dS + S2 K dS, on va ajouter et soustraire de part et d’autre de la discontinuité :
−
→ −−→ −→ −−→ −−→
Z Z Z Z Z
K dS + K1 dSΣ + K dS − K2 dSΣ + (K2 − K1 )dSΣ
S1 Σ S2 Σ Σ
d x2 d x2 h2
Z Z
hdx = 0, et hudx + [g dx]xx21 = 0
dt x1 dt x1 2
donc comme on connaı̂t le résultat d’une dérivation d’intégrale lorsqu’il y a une discontinuité, que l’on écrit ici en 1D
Z Z
d ∂C ∂(Cu)
Cdx = ( + )dx − |C(w − u)|
dt ∂t ∂x
d
R x2
ici C = h donc ici dt x1
hdx = 0, s’écrit
Z
∂h ∂(hu)
( + )dx − |h(w − u)| = 0
∂t ∂x
∂(Cu)
L’équation locale ( ∂C
∂t + ∂x ) = 0 fait disparaı̂tre l’intégrale de ”volume”, il ne reste que que les discontinuités au bord.
|h(w − u)| = 0
- ESV . 76-
Saint-Venant
d
R x2 2
De même si C = hu donc dt x1
hudx + [g h2 dx]xx21 = 0 devient
∂ h2
Z
∂h ∂(hu)
( + + ( ))dx − |hu(w − u)| = 0
∂t ∂x ∂x 2
où C = h pour la masse et C = uh pour la quantité de mouvement, les discontinuités sont :
1
|ρh(w − u)| = 0 et |ρuh(w − u) − ρgh2 | = 0.
2
Remarque
On aurait pu utiliser la démarche aussi de §8.5.2, qui donne bien sûr la même chose. En effet, rappelons que l’on avait montré
∂ ∂
h+ Q = 0 donnait en discontinuité [| − wh + Q|] = 0
∂t ∂x
de même on identifie
∂ ∂ Q2 h2 Q2 h2
Q+ ( + g ) = 0 qui donne en discontinuité [| − wQ + ( + g )|] = 0
∂t ∂x h 2 h 2
les discontinuités sont donc bien en remplaçant Q = uh :
1
|ρh(w − u)| = 0 et |ρuh(w − u) − ρgh2 | = 0.
2
Pour un ressaut fixe la Conservation de la masse et la Conservation de la quantité de mouvement donnent respectivement les relations de saut suivantes ;
h21 h2
U1 h1 = U2 h2 , resp. U12 h1 + g = U22 h2 + g 2 . (11)
2 2
Pour un ressaut mobile, la vitesse de la discontinuité W n’est pas nulle :
h22 h2
conservation de la masse (W − U2 )h2 = (W − U1 )h1 , Conservation de la quantité de mouvement (W − U2 )h2 U2 − (W − U1 )h1 U1 = g − 1 (12)
2 2
Nous allons traiter ces deux cas, qui sont en fait les mêmes dans les sections suivantes.
Nous remarquons au passage que la modélisation nous donne l’impression que le ressaut est infiniment fin, on voit bien sur l’image 38 qu’il n’en est rien. Le ressaut
2
est constitué de bulles et d’écume sur une épaisseur d’environ la profondeur. Même si on avait tenu compte des termes en ∂∂xu21 dans la partie visqueuse de Navier
Stokes, ils n’auraient pas expliqué la structure intime du ressaut.
- ESV . 77-
Saint-Venant
Figure 38 – Un ressaut créé par la montée de la marée à Port à La Duc baie de la Fresnaye, 22, Photo PYL 22/08/09, en dessous sa simplification en deux niveaux.
Figure 39 – ”The Lord of the Ring”, un ressaut au cinéma (délirant : on voit des chevaux d’écume dans le ressaut)
- ESV . 78-
Saint-Venant
Figure 40 – Un ressaut créé par l’eau sortant d’un robinet au fond d’un évier, Photo PYL, on a 1 litre/20s le rayon est d’environ 7cm, la hauteur d’eau environ 0.2
mm ce qui fait un Froude d’environ 12, et une hauteur après le ressaut 16 fois plus haute environ.
h2 h21 h2
U12 h1 h2 + g = U12 h21 + g 2 h2
2 2
d’où en mettant les vitesses à gauche et la gravité à droite et comme on reconnaı̂t une identité remarquable avec (h22 − h21 ) :
- ESV . 79-
Saint-Venant
1 U2 h2
U1 h1 (h2 − h1 )[− 2 − g( − 1)]
2 h1 h1
- ESV . 80-
Saint-Venant
montée de marée
ut
ressa
rivière
Figure 41 – Un ressaut créé par la montée de la marée à Port à La Duc baie de la Fresnaye, Côtes du Nord, Photo PYL 22/08/09. Dans le cas de la rivière avec un
ressaut de marrée qui remonte U2 > 0 et U1 > 0, la vitesse de remontée est ici W > 0.
h1 − h2
(W − U2 )h2 [U2 − (W − (h2 /h1 )(W − U2 )] = (W − U2 )h2 [U2 − W ]( )
h1
d’où en simplifiant par (h1 − h2 )
h1 + h2 h1
(W − U2 )2 = g( )
2 h2
Ce qui donne r
h1 + h2 h1
W = U2 + g( )
2 h2
Reprenant l’équation donnant (W − U2 )2 , on multiplie par 2 et on la divise par gh2
(W − U2 )2 h1 2 h1
2 = +
gh2 h2 h2
cette équation du second degré a une racine positive qui est la relation de Bélanger avec un ressaut mobile :
q
(W −U2 )2
h1 −1 + 1 + 8( gh2 )
=
h2 2
On a donc maintenant toutes les quantités du ressaut mobile.
Remarque : q q
Si on se donne par exemple U2 et h1 et h2 , on voit que W = U2 + g( h1 +h 2 ) h2 , et U1 = gh1 (h
2 h1 1 +h2 )
2h2 (1 − h2 /h1 ) + U2 . On retrouve bien que pour tout U2 donné,
alors U1 et W sont de la forme U2 + F (h1 , h2 ), ce qui montre que l’on peut changer de repère Galiléen le ressaut sans changer sa forme.
- ESV . 81-
Saint-Venant
11.7 Application en instationnaire des invariants de Riemann et des ressauts : Rupture de Barrage sur fond non sec
Considérons maintenant une rupture de barrage sur un lit déjà rempli d’eau. La solution précédente fait apparaı̂tre un ressaut...
....
....
....
Application numérique :
à Port à La Duc h1 = 0.30m, h2 = 0.1m, on a W = 3.4km/h la vitesse de la rivière est de -1.4 m/s et l’eau monte à environ 0.15m/s.
à Saint Pardon h1 = 3.0m, h2 = 2.5m, on a W = 17km/h la vitesse de la rivière est de -.9 m/s et l’eau monte à environ 0.03m/s.
Cuve FC 80
Débit (3,30E-04 m3 /s) h1= 9,50E-03 m, h2=1,40E-02 Fr∼ 1.4 h2/h1= 1,47±0,5 la relation de Bélanger donne 1,57 (à 6%)
Débit (20E-03 m3 /s) h1= 0.019 m, h2=0.064, Fr∼ 3.04 h2/h1= 3,37±0,5 la relation de Bélanger donne 3.8 (14% d’errur)
- ESV . 82-
Saint-Venant
Le débit est écrit sur l’écran de la cuve. Afin de vérifier que la valeur sur l’écran du débitmètre est correcte on a aussi mesuré par ailleurs le temps nécessaire pour
remplir un Bécher de 5L.
L’existence des deux types de ressauts, fixe et dynamique, nécessite deux méthodes. Pour ceux du premier type il faut deux vannes positionnées dans la cuve de la
manière suivante : en premier lieu une vanne amont qui sera fixe pour toutes les mesures, puis à débit fixé, on joue avec la vanne aval afin de rendre le ressaut statique.
Pour chaque débit il faut faire des petits ajustements avec la vanne aval.
Une fois le ressaut fixé, on réalise les mesures de hauteur :
- première méthode, analogique, on utilise le pointeau pour mesurer les hauteurs d’eau.
- seconde méthode, on prend une photo du ressaut puis on mesure les hauteurs d’eau avec le logiciel imageJ. Connaissant l’espace entre deux vis (distance de 15 cm)
qui sert pour créer une échelle, on mesure la distance du fond de la cuve à la surface de l’écoulement. L’illustration a été réalisée sur Paint, d’où l’imprécision des
traits pour les hauteurs h1 et h2.
Pour visualiser le second type de ressaut, le ressaut dynamique, il suffit d’utiliser une seule vanne, de définir un débit puis de fermer brusquement cette vanne. On
filme alors le ressaut qui parcourt la cuve puis on utilise le logiciel imageJ pour mesurer la position du ressaut pour deux images successives. Connaissant le nombre
d’images par seconde de la caméra on trouve ainsi la vitesse du ressaut.
Exemples de tracés de la relation de Bélanger liant les hauteurs et le Froude relatif pour un ressaut comparé à des expériences.
Figure 42 – Données et procédures de L. Gengembre et D. Ruiz 2021 pour vérifier la relation de Bélanger.
- ESV . 83-
Saint-Venant
Lorsque la profondeur est du même ordre de grandeur que la longueur d’onde, nous savons modéliser et décrire le phénomène, dans le cas des petites perturbations :
c’est le problème de la houle de Airy : ([Link] Nous avons vu l’équation de dispersion pour la houle
simple, on a vu que pour des ondes de surface, il faut chercher la solution du problème de Euler linéarisé
∂u ∂v ∂u ∂p ∂v ∂p
+ = 0, ρ =− , ρ =− ,
∂x ∂y ∂t ∂x ∂t ∂y
sous forme d’onde exp(i(ωt − kx)) avec à la surface pour la perturbation de pression p = ρgη.
Nous allons partir de cette description et regarder ce qu’il se passe lorsque λ = 2π/k devient plus grand que h0 (ou si on veut kh0 de plus en plus petit). On espère
ainsi être dans ce régime intermédiaire où la profondeur est faible mais pas trop. Cela va nous donner des idées pour corriger les équations de Saint-Venant lorsque la
profondeur devient pas trop faible.
On a trouvé pour le problème de la houle d’Airy des solutions en ondes telles que ω 2 = gk tanh(kh0 ). Si on développe dans le cas peu dispersif de l’eau peu profonde,
à grande longueur d’onde kh0 → 0, on fait un développement limité en kh0 de la tangente puis de la racine, (et c20 = gh0 ) :
(kh0 )3 (kh0 )2 (kh)2
ω 2 = gk(kh0 − + ....) = (gh0 )k 2 (1 − + ....), dont les racines sont ω = ±c0 k(1 − + ....).
3 3 6
2
On en déduit en prenant la valeur + des ondes qui se déplacent vers la droite. On la forme d’onde η = η0 exp(i(ωt − kx)), avec iω = (ik)c0 (1 + (h60 ) (ik)2 + ....) ce qui
veut dire que puisque ∂t η = iωη et que ∂x η = −ikη l’équation de dispersion linéarisée est celle correspondant au problème suivant (elle s’appelle équation de KdV
linéarisée) :
∂η ∂η c0 h20 ∂ 3 η
+ c0 + = 0.
∂t ∂x 6 ∂x3
(h0 )2 2
Ce la veut bien dire que cette équation a pour relation de dispersion : iω = c0 (ik)(1 + 6 (ik) ).
R∞
Nous allons résoudre cette équation en supposant que la surface déplacée −∞ ηdx est une donnée. La première idée est de se déplacer avec la vitesse c0 et de poser
ξ = x − c0 t et α = c0 h20 , l’équation devient :
∂η α ∂3η
=− .
∂t 6 ∂ξ 3
Cette équation, se résout par la technique des solutions semblables. Par invariances par dilatations on cherche des solutions semblables... Consulter
t = T t̂
le changement d’échelle ξ = X ξˆ (17)
η = H η̂
R∞ R∞ R∞ R∞
La conservation de la masse totale −∞ ηdξ devient HX −∞ η̂dξˆ mais comme on veut l’invariance −∞ ηdξ = −∞ η̂dξ, ˆ donc HX = 1 préserve la conservation de la
R∞ 3
surface déplacée −∞ ηdx = 1. De même pour l’équation, si T = X 3 cela préserve l’invariance de l’équation qui s’écrit identiquement ∂∂η̂t̂ = − α6 ∂∂ ξ̂η̂3 . La variable de
- ESV . 84-
Saint-Venant
ξ ξ
similitude est ζ = t1/3 et la surface est de la forme : η = t−1/3 f ( t1/3 ). Par substitution et dérivation la fonction f (ζ) vérifie −α/6f 000 = −ζf 0 /3 − f /3 . En intégrant,
00
et comme f est nulle à l’infini, on a : αf = 2ζf . R∞
La solution de y 00 (x) = xy(x) avec y(∞) = 0 est y = Ai(x) la fonction d’Airy, de plus −∞ Ai(x)dx = 1. On a donc la solution pour f = (2/α)1/3 Ai((2/α)1/3 ζ),
puisque ξ = x − c0 t la solution est au final :
2 1/3 2 1/3 (x − c0 t)
η(x, t) = ( ) Ai[( ) ].
c0 h20 t c0 h20 t1/3
Nous allons retrouver ce résultat autrement dans la suite.
∂η ∂η 3c0 ∂η c0 h20 ∂ 3 η
+ c0 + η + = 0.
∂t ∂x 2h0 ∂x 6 ∂x3
2
Cette équation possède la bonne non linéarité et a aussi été construite de manière a obtenir pour une onde eiωt−kx la relation de dispersion : iω = (ik)c0 (1 + (h60 ) (ik)2 )
qui est la bonne approximation à l’ordre (h0 k)2 .
- ESV . 85-
Saint-Venant
Figure 44 – Soliton observé dans une cuve à Vague au Palais de la Découverte [Link] lagree/SIEF/SIEF97/[Link]
(f 0 )2 = 6cf 2 − 3f 3
Il est alors facile de voir qu’une solution de la forme α/ cosh(βs)2 convient ! ! avec α = 2c et β = (3c/2)1/2 .
La solution de (f 0 )2 = 6cf 2 − 3f 3 est donc
2c
f (x − ct) =
cosh[(3c/2)1/2 (x − ct)]2
C’est l’onde solitaire :
a
η= a
cosh[(3a/4h3 )1/2 (x − (1 + ( 2h ))t)]2
ou, on peut aussi trouver des ondes cnoı̈dales... succession de bosses...
The extra term comes from the fact that the pressure is no more hydrostatic due to vertical acceleration term. We linearise around a mean level of water say h0 .
Let say that pressure is the sum of the hydrostatic pressure plus a small perturbation of it : p = ρg(h − y) + Π. This perturbation is linked to the up to now neglected
- ESV . 86-
Saint-Venant
y 2 −h20 ∂ 2 u
transverse momentum : − ∂Π ∂v ∂u ∂Π
∂y ' ρ ∂t , Hence, as from continuity : v ' −y ∂x , we substitute in transverse momentum and integrate − ∂y to obtain Π(y) ' ρ 2 ∂t∂x
(we have Π(h) = 0). This is integrated again across the depth (St Venant equation are integrated across the depth) :
y h0
y3 yh2 ∂2u h30 ∂ 2 u
Z Z
Πdy ' ρ[( − 0 )]h0 0 , so Πdy ' −ρ( )
0 6 3 ∂t∂x 0 3 ∂t∂x
We guess that we can estimate the extra integrated again across the depth pressure gradient that will modify the longitudinal momentum :
h0
∂3u
Z
∂Π 1
− dy ' ρh20
0 ∂x 3 ∂t∂x2
12.5 Mascaret
Lorsque la mer monte dans une embouchure de rivière, la marée étant de plus en haute, elle peut créee par accumulation des non linéarités un ressaut qui remonte
l’écoulement, voir figure 41. Si la profondeur du fleuve est adéquate, le ressaut se casse par la dispersion et un train d’onde apparaı̂t. C’est ce que l’on observe sur la
figure 48. Bien plus en amont, ce mascaret finit par être détruit, un seul soliton peut être éventuellement observé.
[Link]
- ESV . 87-
Saint-Venant
- ESV . 88-
Saint-Venant
Cinquième partie
Résolution numérique avec Basilisk
13 Résolution numérique
13.1 Résolution numérique : schémas pour les équations 1D
Les équations ressemblent beaucoup aux équations de la dynamique des gaz. De nombreuses méthodes ont été proposées pour ces équations. On consultera princi-
palement les ouvrages de Eleuterio Toro ou de Randall LeVeque.
Le plus important est d’avoir un oeil critique sur le résultat : bien comprendre le résultat numérique obtenu, pouvoir l’interpréter, voir ses limites.
- ESV . 89-
Saint-Venant
Figure 49 – Calcul EDF [14] effectués à l’aide du code TELEMAC d’EDF, dispersion d’eau chaude dasn un port et crue de la Severn. Cette dernière rivière est
célèbre pour son mascaret, voir plus loin.
- ESV . 90-
Saint-Venant
cd basilisk/src
export BASILISK=$PWD
export PATH=$PATH:$PWD
cp [Link] config
make
La page d’installation [Link] la suivre scrupuleusement. Noter que make -k force la compilation sans s’arrêter sur des erreurs
intermédiaires.
15 Exemples
15.1 Equation d’advection
Examinez le code [Link] en Basilisk ou le code [Link] en
C standard.
- ESV . 91-
Saint-Venant
- ESV . 92-
Saint-Venant
15.6 Tas Visqueux s’effondrant sur une surface plane (tas de Huppert), onde de crue diffusante
Cet exemple est extrait de Huppert ”The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface” J . Fluid Mech.
(1982), vol. 121, p p . 43-58
/Users/pyl/macintoshHD/DOKUMENTS/ENSEIGN/ENSTA/PC1213/[Link]/huppert/[Link]
Il s’agit d’étudier avec les équations de Saint-Venant la chute d’un tas initialement au repos et rectangulaire sous l’action de la gravité et avec un frottement
visqueux de Poiseuille
∂h ∂(Q)
+ =0
∂t ∂x
∂Q ∂ Q2 g Q
+ ((5/4) + (h2 )) = −3ν 2
∂t ∂x h 2 h
pour les temps longs, on a équilibre entre la pression hydrostatique et la friction 3νQ/h2 = −gh∂x h. On substitue Q dans l’équation de la masse
∂h g ∂ 3 ∂h ∂h ∂ ∂h ∂h ∂ 2 h4
− h = 0 qui est de la forme − k h3 = 0 ou aussi − (k/4) =0
∂t 3ν ∂x ∂x ∂t ∂x ∂x ∂t ∂x2
Par dilatations des échelles l’invariance par changement d’échelle équivaut à H/T = H 4 /X 2 pour cette équation, et comme la masse conservée HX = 1, cela donne
H = T −1/5 et X = T 1/5 donc une solution de la forme autosemblable
h = t−1/5 H(xt−1/5 )
on substitue et on a
−5kH(η)3 H00 (η) − 15kH(η)2 H0 (η)2 − ηH0 (η) − H(η) = 0
intégrale première
(−5kH(η)3 H0 (η))0 − (ηH(η))0 = 0, par symétrie H0 (0) = 0
d’où
−(5k)(H(η)2 H0 (η)) = η or (5k/3)H(η)3 = b2 /2 − η 2 /2
e tas arrête quand H = 0 au point η = b
(5k/3)H(η)3 = b2 /2(1 − (η/b)2 )
soit
H(η) = (3b2 /(10k)(1 − (η/b)2 ))1/3
La surface totale Z b
A= H(η)dη
−b
On a √
b
2 πΓ 13
Z
2 1/3
(1 − (η/b) )) dη = b
5Γ 56
−b
- ESV . 93-
Saint-Venant
Figure 50 – Gauche Plot at t=100,200,300...1000 of h(x, t) with Basilisk. Droite Plot at t=100,200,300...1000 of t1/5 h(x/t1/5 , t) with Basilisk, with a FV code order 1
HLL and the analytical solution (3b2 /(10k)(1 − (η/b)2 ))1/3 .
- ESV . 94-
Saint-Venant
∂h gZ 0 ∂ 3
− h =0
∂t 3ν ∂x
de la forme
∂h ∂h
− kh2 =0
∂t ∂x
par changement d’échelle, l’invariance par dilatation donne H/T = H 3 /X, la masse conservée donne HX = 1, ce qui permet de trouver H = T −1/3 et X = T 1/3
donc une solution de la forme t−1/3 H(x/t1/3 ) soit
−ηH0 (η) − H(η) + 3kH(η)2 H0 (η) = 0
si H(0) = H0
(−ηH(η) + kH(η)3 ) = kH03
on a la forme implicite
H03
η = k(H(η)2 − )
H(η)
On part d’une masse localisée entre x0 et x1 , on va donc prendre H0 = 0 donc η = kH(η)2 , simplement, on va avoir
p
H(η) = (η)/k,
la solution de
∂h ∂ ∂h
− k h2 =0
∂t ∂x ∂x
est simplement pour x < xf r
x
q
−1/3 −1/3
h=t (x)t )/k =
kt
R x1
On se donne une masse initiale A = 0
h(x, 0)dx qui se déplace entre xf t, position du front et 0 le front arrière qui ne bouge pas, la masse est donc ultérieurement
xf
Z r
x q
dx = (2/3) x3f /(kt)
0 kt
2
qui est A par conservation de la masse soit xf = ( 9A4 kt )1/3
Remarque : on peut aussi le résoudre avec les caractéristiques.
- ESV . 95-
Saint-Venant
La solution est assez subtile, en temps que caractéristique, la solution est constante le long des caractéristiques
dx
= c(x, t) avec c = kh2
dt
pour une distribution donnée, h(x, t = 0) = h0 (x) alors le long des droites x = ct + ξ
Le long des droites x = ct + ξ la valeur de la solution est constante, c’est la valeur en ξ qui vaut ensuite x − ct qui est conservée.
Donc, si on a une distribution au temps t = 0 de hauteur h0 (x) initiale, chacun de ces points x est en fait un des points initiaux ξ, alors la solution est h = h0 (x−tc(h0 (ξ)))
avec c = kh2 donc
x = ξ + tkh2 (ξ, 0)
p
donc h = (x − ξ)/(kt) et pour de grands x, on a s
(x)
h=
(kt)
comme présenté.
- ESV . 96-
Saint-Venant
Figure 52 – Gauche Plot at t=100,200,300...1500 of h(x, t) with Basilisk. Droite Plot at t=100,200,300...1500 of t1/3 h(x/t1/3 , t) with Basilisk.
- ESV . 97-
Saint-Venant
event init (i = 0)
{
terrain (zb, "/Users/pyl/basilisk/basilisk/src/test_pyl/DS_OK/WORD/etopo2", NULL);
définir des jauges
Gauge gauges[] = {
// file lon lat description
{"coco", 96.88, -12.13, "Cocos Islands, Australia"},
{NULL}
};
Le premier problème est de déterminer la source : c’est en fait le fond qui bouge très vite et qui déplace la masse d’eau,
- ESV . 98-
Saint-Venant
Figure 53 – exemple de calcul du Tsunami de décembre 2004. Elévation calculée à 7200s de l’onde rouge : +2 m bleu -2. A droite observation du Tsunami par le
satellite Jason-1 comparé et calcul de la simulation. [Link]
- ESV . 99-
Saint-Venant
Figure 54 – exemple de calcul du Tsunami du 11 mars 2011. Elévation calculée à 7200s de l’onde rouge : 0.5 m bleu -0.5m. A droite observations du Tsunami par
les bouées DART dans le pacifique, calcul G erris. [Link]
- ESV . 100-
Saint-Venant
- ESV . 101-
Saint-Venant
Sixième partie
Granulaires, Neige, Glace, boue...
17 Introduction
Nous venons de voir des écoulements d’eau en couche mince, ils sont décrits par les équations de Saint-Venant en bonne première approximation (tant que la
pression reste hydrostatique et tant que la profondeur est assez faible par rapport au développement longitudinal et que la pente n’est pas trop forte). Nous avons
considéré le cas du fluide Newtonien laminaire puis turbulent. Souvent dans la nature, la pente peut être plus forte, et surtout ce n’est pas uniquement de l’eau qui
coule mais un mélange complexe de débris variés. la rhéologie de ces fluides n’est pas encore bien connue ni modélisée. Pour ces fluides plus complexes, on va cependant
utiliser les mêmes hypothèses et essayer d’écrire le frottement τ à la paroi en fonction de Q et h et de coefficients pertinents.
Ce fluides étant ”collants”, il peut y a voir un seuil τseuil tel que si les contraintes sont sous ce seuil, alors le fluide ne bouge pas (cas de la boue, des rochers, du
sable....). Il faut que τ > τseuil pour que le mouvement se produise. Dans le cas où il y a beaucoup de roches, le caractère ”solide” va l’emporter dans le frottement.
On utilisera la loi d’Amontons-Coulomb :
où µ est le coefficient de friction et où la contrainte normale est prise égale à la pression hydrostatique ρgh en première approximation (il existe certaine controverse
à ce sujet liée au fait que l’on parle d’un solide et non d’un fluide, mais nous allons au plus simple). Pour les granulaires secs, ce sera en effet une bonne description,
pour la neige, il faudra y rajouter une partie dynamique.
18 Equations de Saint-Venant pour les fluides plus complexes : granulaires & neige, ...
Le système que nous considérons est celui des Equations de Saint-Venant (CRAS 1871 Adhémar Jean-Claude Barré de Saint- Venant, Shallow Water) :
∂t h + ∂x Q = 0
Q2 h2
τ , (18)
∂t Q + ∂x Γ +g = −gh∂x Z −
h 2 ρ
avec h(x, t) est la hauteur d’eau, u(x, t) la vitesse moyenne, Q(x, t) = h(x, t)u(x, t) le flux, Γ facteur de forme, τ frottement au fond et Z(x) est la ”topographie” ou
forme du fond, g la gravité...
On va garder ces équations où en fait les hypothèses de rhéologie n’ont pas été appliquées (seulement couche mince, faible angle, moyenne en épaisseur). On a
besoin du facteur de forme Γ, mais on prendra 1, on a surtout besoin de modéliser τ le frottement à la paroi. On remarque ici, que le détail de τ dans l’épaisseur est
sans importance. On a besoin de τ à la paroi uniquement.
- ESV . 102-
Saint-Venant
En pratique pour un glacier n = 1/3 c’est la loi de Glen. L’écoulement est tellement lent que l’on ne se soucie pas de Γ, on est dans la description de l’onde
d’inondation : n
∂h ∂Q Q ∂h
+ = 0; avec cn µn 2
= −gh∂x Z − gh
∂t ∂x h ∂x
- ESV . 103-
Saint-Venant
Un modèle correspondant à la lente avancée des glaciers a été proposé par M. W. Mahaffy 1976, l’image suivante montre un exemple de répartition de la glace (A
three-dimensional numerical model of ice sheets : Tests on the Barnes Ice Cap, Northwest Territories, JGR). La solution autosemblable a été écrite par P. Halfar 1981
(”On the Dynamics of the Ice Sheets” JGR).
voir [Link]
18.4 La boue
u n
La boue peut être considérée aussi comme un fluide en loi de puissance τ = cn µn avec 0.1 < n < 0.4. C’est la description de Ng et Mei. (1994), ”Roll waves
h
on a shallow layer of mud modelled as a power-law fluid” Fluid Mech. C’est aussi envisagé dans Lervolino Vacca & Di Cristo. Simplified Wave Models Applicability
to Shallow Mud Flows Modeled as Power- Law Fluids Mountain Science · (2015) DOI : 10.1007/s11629-014-3065-6.
Cependant, une modélisation plus fine de la boue doit tenir compte du fait qu’il existe un seuil en contrainte pour l’écoulement τseuil . Les nombreux articles de
Balmforth évoqués dans [Link] montrent alors qu’il faut analyser plus précisément le profil de vitesse
pour écrire l’onde cinétique diffusante.
19 Conclusion
En première approximation, les équations de Saint-Venant restent valides pour la plupart des fluides géophysiques qui s’écoulent en couches mince. Ces équations
(ou des simplifications ”onde cinétique”) sont utilisées effectivement de plus en plus dans des codes dédiés pour reproduire les évènements passés, et pour enuite r
prévenir les risques sur le terrain.
En toute première approximation, une loi newtonienne suffit (lave...), mais parfois il faut rajouter des détails de transfert de chaleur et de solidification (lave...),
ou changer la rhéologie avec des lois en puissance (écoulement de glaciers, boue....) L’aspect friction est important pour la Neige, les avalanches de roches, les dunes
de sable.... Cependant, dans le cas des fluides de type ”boue” si on veut une description précise de l’arrêt on aura besoin du détail sur l’épaisseur de l’écoulement pour
résoudre effectivement. Cela sort donc un peu du cadre Saint-Venant.
Plus de détails voir [Link] qui précise les écoulements de type boue et milieux granulaires qui sont
des écoulements à seuil.
- ESV . 104-
Saint-Venant
Septième partie
Les Sédiments : Problèmes de sédimentation érosion
20 Introduction
Il s’agit du problème général de transport de sédiments (de particules, de sable, d’alluvions, de boues), on parle aussi de courant de turbidité, mais c’est aussi
le transport de déchets, de polluants dans un écoulement de fluide (d’eau ou d’air). Les écoulements envisagés sont tous ceux liés aux fleuves et à la mer, ou au
déplacement des dunes de sable. Les ordres de grandeurs vont des rides de sable sur la plage qui sont de quelques centimètres aux bancs de sable (méga ripples) de
plusieurs centaines de mètres sous la Manche. Les dunes du Sahara vont de 5m à plusieurs centaines de mètres.
La complexité vient ici en plus du fait que la présence des particules modifie la viscosité de l’écoulement, mais surtout, que le déplacement des particules modifie
le fond de l’écoulement et modifie donc aussi l’écoulement lui même.
On simplifiera en supposant que les sédiments sont emportés comme des scalaires passifs par la couche limite (turbulente) avec en plus une vitesse de chute de
sédimentation. On n’étudiera pas ici le problème des avalanches qui est un problème de transport dans un milieu granulaire.
F = 6πµ(d/2)v, donc mdv/dt = −mg + meau g + 6πµ(d/2)v donc (4/3π(d/2)3 )ρp dv/dt = (4/3π(d/2)3 )(ρp − ρ) g − 6πµv(d/2)
2
la vitesse terminale est la vitesse de Stokes Vf = 18dρ ν (ρp − ρ)g on a aussi un temps caractéristique outre cette vitesse caractéristique (vitesse de Stokes, ou vitesse de
chute (settling velocity) ou vitesse finale). On réécrit donc
dv Vf − v ρp d2 (ρs − ρ)gd2
= avec TStokes = , et Vf = (19)
dt TStokes ρ ν 18µ
On a ainsi déterminé la vitesse de sédimentation, dans le cas très visqueux pour un grain sphérique de diamètre d. Le temps caractéristique de déposition, dans le
cas visqueux est le temps de chute d’une hauteur d à la vitesse de Stokes : Vds (cf. Charru et Mouilleron-Arnould [?]).
Si la vitesse est trop rapide, on quitte le régime laminaire. Une estimation simpliste de Vf est souvent présentée sous la forme d’un équilibre entre la traı̂née
aérodynamique (attention ici prise en Vf2 , Stokes n’étant qu’une approximation à faible vitesse), la poussée d’Archimède et le poids de la particule (d diamètre, cd
coefficient de traı̂née..).
L’équation de mise à l’équilibre :
dv 4 cd
(4/3πd3 )ρp = ( πd3 )(ρp − ρ) g − πd2 v 2 (20)
dt 3 2
Cette fois comme noté par Andreotti [?], c’est une longueur caractéristique qui apparaı̂t naturellement :
s
dv Vf2 − v 2 ρp 8(ρs /ρ − 1)gd
' avec lsat = d. et Vf = . (21)
dt lsat ρ 3cd
- ESV . 105-
Saint-Venant
On peut construire une vitesse avec simplement la taille des grains et l’accélération de la pesanteur, c’est un peu une version appauvrie du cas précédent :
p
V = gd.
Enfin, on peut construire une vitesse à partir de la vitesse du fluide elle même, la vitesse pertinente sera liée au cisaillement de vitesse. En laminaire, entre le haut
et le bas d’une particule dans un écoulement cisaillé γ = ∂u ∂u
∂y , la différence de vitesse est d ∂y . C’est donc avec cette vitesse que l’on va construire le nombre de Reynolds
2
associé à la particule ( dν ) ∂u 2
∂y . Dans le cas turbulent, on utilise la vitesse de frottement construite sur la contrainte à la paroi : τ = ρu∗ . Le nombre de Reynolds associé
sera donc ρu∗ d/µ.
Par la suite nous prendrons une de ces vitesses comme un des paramètres du problème. Ces vitesses, longueurs et temps sont ensuite systématiquement utilisés
pour dimensionner les phénomènes. On peut cependant imaginer que cette vitesse n’est pas constante, qu’elle varie avec la densité... Il s’agit d’une résolution hyper
simplifiée des équations de la mécanique pour les grains qui traduit notre ignorance des équations exactes.
Expérimentalement on mesure une valeur du paramètre de Shields critique (noté Sc ) à partir duquel il peut y avoir mise en suspension et déplacement des particules.
(ρs −ρ)gd2 µ ∂u d ∂u
Comme τ = µ ∂u
∂y et que la vitesse de Stokes est Vf = 18µ , on a S = ∂y
(ρs −ρ)gd = ∂y
18Vf , rapport de deux vitesses d ∂u
∂y et Vf . .
Ce phénomène de mise en mouvement lorsqu’un seuil a été dépassé a déjà été remarqué par du Boys en 1879 “ un caillou posé au fond d’un courant liquide, peut
être déplacé par l’impulsion des filets qui le rencontrent : le mouvement aura lieu si la vitesse est supérieure à une certaine limite (...) [qui] dépend de la densité, du
volume et de la forme du caillou ; elle dépend aussi de la densité du liquide.”
La figure 55 présente la valeur du nombre de Shields critique pour différents grains en fonction du nombre de Reynolds construit ici sur le diamètre et la vitesse de
frottement (définie pour les profils turbulents au paragraphe 6.7.2), rappelons que la valeur du frottement au fond τ0 permet d’introduire une vitesse fictive u∗ mais
qui est quand même un bon ordre de grandeur de la fluctuation de vitesse telle que le frottement au fond τ0 = ρu2∗ . Le Shields dépend du nombre de Reynolds du
grain :
∂u d2 ∂u
∂y u∗ d
dans le cas visqueux, on utilise d comme vitesse et Re = ρ , dans le cas turbulent Re = ρ .
∂y µ µ
On remarque sur la figure 55 que S ∼ 0.04 sur une grande plage. Le nombre de Shields critique est pris souvent constant en première approximation car il varie
faiblement en Re.
On peut aussi définir un nombre de Rouse, rapport entre la vitesse de sédimentation et la vitesse u∗ : par définition R = Vf /(κu∗ ). La définition historique y
rajoute la constante de Kármán 0.41. (cf ??).
- ESV . 106-
Saint-Venant
Mass: Threshold,
18 The Shield criteria
Mass: Threshold, The Shield criteria
Figure 55 – Diagramme de Shields : valeur critique Sc du frottement adimen- mation le surcoût est proportionnel à la pente Λ ∂f
vitesse est supérieure à une certaine limite qu’il (S. Gras) nomme vitesse d’entraı̂nement. Cette vitesse limite dépend de la densité, du
τ ∂x .
sionné par le poids déjaugé (S = (ρs −ρ)gd ) à partir duquel les grains bougent volume et de la forme du caillou; elle dépend aussi de la densité du liquide et de la profondeur du courant.”
∂f
S − Sc − Λ .
∂x
- ESV . 107-
Saint-Venant
La notation x+ veut dire x si x > 0 et 0 si x est négative. Mais il n’existe pas de vrai consensus, voir figure 57 avec la loi de Nielsen. et des lois de la forme suivante
sont utilisées (V On peut écrire le flux en fonction du débit qs = AQ (Q − Qs )B B
+ ou en fonction de la vitesse moyenne qs = AV (V − Vs )+ ou en fonction du frottement
B B
qs = AF (τ − τs )+ voire en fonction d’une puissance qs = AP (τ V − τs Vs )+ (voir la figure 58 extraite du livre de Yang [29]) Au final on retiendra que qs est proportionnel
à une vitesse et une longueur qui sont environ la vitesse caractéristique Vf et au diamètre du grain disons d et à l’écart entre le frottement réduit et le seuil critique
élevé à une puissance β déterminée empiriquement, on écrira donc
µ ∂u
∂y
q/(Vf d) = S 3 avec Vf = ∆ρgd2 /(18µ) avec le Shields S =
∆ρgd
Charru Mouilleron-Arnould J. Fluid Mech. (2002), ”Instability of a bed of particles sheared by a viscous flow” 2002 en pratique (Charru Mouilleron-Arnould Eiff
”Erosion and deposition of particles on a bed sheared by a viscous flow” JFM 04) proposent :
2
µ ∂u
∂y
q/(Vf d) = 0.45S(S − 0.12) avec Vf = ∆ρgd /(18µ) avec le Shields S =
∆ρgd
- ESV . 108-
Saint-Venant
- ESV . 109-
Saint-Venant
avec par définition u∗ la vitesse de frottement, c’est à dire que par définition τ = ρu2∗ . on estime donc que pour une rivière à fond lisse ou pour une couche limite de
vent, on est comme dans un tuyau et donc que la loi logarithmique est valable (c’est une approximation usuelle)
U zu∗
= 2.5 ln( )
u∗ ν
Pour un fond rugueux, on s’attend à la même forme, mais avec une autre échelle que ν/u∗ pour dimensionner la profondeur, soit z0 et alors U = uκ∗ ln( zz0 ). Or, dans
un tuyau rugueux de rugosité ks (taille des rugosités ou mini bosses plus ou moins régulièrement réparties sur la paroi) Nikuradse a décomposé la vitesse du tuyau
lisse en faisant apparaı̂tre ks artificiellement (dans le cas lisse il n’y a pas de ks !) :
U u∗ z z ks u∗
= 2.5 ln( ) + 5.5 = 2.5 ln( ) + 5.5 + 2.5 ln( )
u∗ ν ks ν
et défini ainsi une fonction de rugosité B, cette fonction est exactement 5.5 + 2.5ln(ks u∗ /ν) dans le cas lisse (ks → 0), et sinon elle a été mesurée expérimentalement
et tabulée. Cette fonction B, est donc telle que pour un écoulement rugueux quelconque
U z
= 2.5 ln( ) + B
u∗ ks
avec
z
U = 2.5 ln( ) + B avec cas lisse Rep < 5, B = 5.5 + 2.5ln(ks u∗ /ν) avec Rep > 70, B = 8.5
ks
u∗ z
pour les cas très rugueux, B est constante=8.5, et comme e8.5/2.5 ' 30 d’où en pratique U = κ ln( z0 ) avec z0 ∼ d/30.
22.3.2 Flux
Cette analyse est aussi valide pour l’air, mais en fait, pour l’air, Bagnold a observé que la hauteur z0 dépendait de la force du vent. Cette longueur est l’ordre de
grandeur de la hauteur et longueur de saut des grains qui est environ ` = u2∗ /g en faisant un bilan de quantité de mouvement, avec u− et u+ vitesses de chute et de
montée :
`
τ ` = q(u− − u+ ) ou si on intuite la possibilité d’un seuil τt threshold tel que q = − (τ − τt )
(u − u+ )
ordre de grandeur u∗ et ` ∼ u2∗ /g et τ ∼ ρu2∗ le flux est donc en ordre de grandeur
ρf 3
q= u
g ∗
Les formules suivantes ont été alors proposées (avec la vitesse de frottement τ = ρu2∗ )
p ρf 3 ρf
Bagnold 1936 : q = Cb d/D u avec Cb = 18, D = 0.25mm, Lettau et Lettau 1978 : q = CL u2∗ (u∗ − u∗t )+ avec CL = 4.1
g ∗ g
Cependant, des campagnes de mesures plus récentes de Ho et al. ont constaté que le flux dépendait du sol, s’il est érodable on non.
S’il est érodable ils n’ont pas observé la dépendance de ` ni de (u− − u+ ) en u∗ , ils pensent donc que q ∼ 27(ρf /g)(u2∗ − u2∗t ) avec u∗t = 0.15m/s. Si le fond n’est
pas érodable, ils ont observé que ` augmente en (u∗ − u∗t )2 et (u− − u+ ) en (u∗ − u∗t ) soit au final q ∼ (u∗ − u∗t )(u2∗ − u2∗t )+ avec u∗t = 0.2m/s :
- ESV . 110-
a major role in determining the dynamics and morphology of desert dunes. It is therefore necessary to
understand the principles governing the entrainment and transport of sand by the wind before Saint-Venant
examining
dune processes and dynamics.
Transport of sediment by the wind involves interactions between the wind and the ground surface.
Understanding the operation of these processes requires a knowledge of relevant surface characteristics
(e.g. sediment texture, vegetation cover, degree of cohesion and crusting) as well as the dynamics of airflow
over the surface. There are three distinct modes of aeolian transport (Figure 2.1). These depend primarily on
the grain size of the available sediment (Bagnold 1941). Very small particles (<60–70 µm) are transported
in suspension and kept aloft for relatively long distances by turbulent eddies in the wind. Particles of this
size range play a minor or
q
& g ! z " ! (u g ! z " . !10"
l
In the following we restrict the discussion to the grain
FIG. 1. Sketch of the horizontal sand flux q caused by saltating borne shear stress & g0 at the ground that is given by the
grains passing the vertical rectangle. The dashed rectangle shows momentum transfer from the grains to the bed during their
the surface area l times a that is used to calculate the flux ' of impacts as sketched in Fig. 1. We denote quantities that are
grains impacting onto the surface, where l is the length of the sal- taken at the ground with an index 0, e.g. & g0 ! & g (z 0 ), where
tation trajectory.
z 0 is the height at the ground. We need the impact and ejec-
tion velocities of the grains or moreover the change in hori-
61 –leads
of Eq. !5"
IntegrationFigure Fluxfinally
issu deto[24]
the et inspiré
well known de Bagnold
law of zontal velocity (u g0 at the ground to calculate the momen-
the wall for turbulent flow and therefore to the logarithmic tum transfer. In order to keep the discussion simple we will
profile of the atmospheric boundary layer directly formulate the model in terms of mean values for the
trajectory length l and its time T instead of writing every-
u z thing in terms of not well known distribution functions. - ESV . 111-
v! z " ! * ln , !6"
Saint-Venant
∂f ∂q
φ =− .
∂t ∂x
C’est la relation pratique de conservation de la masse des sédiments. Le paramètre φ représente la compacité. C’est le rapport entre le volume des grains et le volume
qu’il occupent. Ce coefficient est environ 0.6.
23.2 hypothèses
Il s’agit d’une interaction entre un écoulement et le fond érodable sur lequel il s’écoule : la forme du fond gouverne l’écoulement, ce dernier modifie la forme du
fond par érosion et sédimentation.
Pour simplifier, on suppose que l’écoulement est quasistationnaire (la variation de la forme du fond en fonction du temps est supposée lente).
Se donnant une première forme de fond f (x, t = 0), on doit suivre la démarche suivante :
• On travaille à l’échelle de temps lente, à chaque instant (dans cette échelle lente), la forme du fond est fixée, l’écoulement est calculé. Comme les phénomènes
d’interaction fluide/sol se produisent par définition près de la paroi, le frottement pariétal est évalué par résolution des équations de Navier Stokes incompressibles
stationnaires (en fait en pratique Saint-Venant), à viscosité constante. On obtient de τ le frottement au fond, donc la valeur du Shields S :
- ESV . 112-
Saint-Venant
S. Cordier, M.H. Le, T. Morales de Lun ”Bedload transport in shallow water models : why splitting (may) fail, how hyperbolicity (can) help”
M. J. Castro Dı́az. D. Fernández-Nieto. M. Ferreiro ” Sediment transport models in Shallow Water equations and numerical approach by high order finite volume
methods” Computers & FluidsMarch 2008...
Kouakou Lagrée, avec une perturbation d’un écoulement cisaillé (pas Saint-Venant )
- ESV . 113-
Saint-Venant
17 18
Masse : loi de conservation Masse : loi de conservation Masse : loi de conservation
f f f
∂R ∂R ∂R ∂q
= ... + Γ = ... + Γ = − + Γ c’est en fait une
Figure 62 – Dans un écoulement avec ∂t érosion et sédimentation, considérons un petit volume
∂t de longueur ∆x de grains. L’épaisseur est ici
∂t arbitraire,
∂x
hauteur de grains roulants (notée d1 ensuite). Le nombre de grains en mouvement augmente par érosion ṅe ∆x (à gauche), il tombe par sédimentation ṅd ∆x ) (au
∂f
centre). Pour les grains mobiles, le bilan ∂f
de =l’érosion
−Γ sédimentation est ṅe ∆x − ṅd ∆x, ces grains ∂f flux qui rentre q(x) et
= −Γen mouvement ont un flux noté q (à droite), le
= −Γ
∂t ∂t ∂t
sort q(x + ∆x)apporte ∂x q∆x grains au bilan. Pour les grains fixes, le bilan est l’opposé de l’érosion sédimentation −ṅe ∆x + ṅd ∆x. Les grains dans le sol sont fixes.
∂q
ls + q = q0 ((S − Ss )+ )β
∂x
L’interprétation de cette équation est qu’il faut une certaine distance ls = γdV2fd1 avant que le flux saturé ne soit atteint. Plus la vitesse de chute est grande, plus ls
est petit, plus le cisaillement est grand (γ) plus ls est grand. Ces résultats ont été établis pour un liquide.
ρgrain
D’autres auteurs (voir les travaux de Andreotti Claudin & Douady) estiment qu’en fait fait ls = ρf luide d, cette estimation tient plutôt pour l’air.
Cette formule peut se résumer par les phrases de du Boys, 1879. Commentant la vitesse limite d’entraı̂nement évoquée par ces prédécesseurs, il évoque la notion
de saturation de flux en ces termes “ une fois une certaine quantité de matières en mouvement sur le fond du lit, la vitesse des filets liquides devient trop faible
pour entraı̂ner davantage : le cours d’eau est alors saturé. L’existence d’une longueur de saturation est traduite par : un cours d’eau non saturé tend à le devenir en
entraı̂nant une partie des matériaux qui composent son lit, et en choisissant de préférence les plus petits.”
On remarque incidemment que l’équation de l’évolution du fond s’écrit aussi (équation d’Exner 1925) en éliminant les taux d’érosion et de sédimentation :
∂f ∂q
φ =− .
∂t ∂x
C’est la relation pratique de conservation de la masse des sédiments. Le paramètre φ représente la compacité. C’est le rapport entre le volume des grains et le volume
qu’il occupent. Ce coefficient est environ 0.6.
- ESV . 114-
Saint-Venant
Devauchelle Malverty Lajeunesse Lagrée Josserand KD Nguyen Thu Lam ” Stability of bedforms in laminar flows with free-surface : from bars to ripples”JFM
2010 mais en linéarisant les équations de laminaires NS (pas Saint Venant)
Kouakou Lagrée, dans un écoulement laminair cisaillé
Lagrée avec un écoulement sans surface libre
Andreotti Claudin avec un écoulement sans surface libre
Revue Courrech Dupond
idée initiale de Kroy
- ESV . 115-
Saint-Venant
26 Transport en suspension
Nous venons de voir le cas où les grains sont assez lourds, ils restent près du fond, roulent et glissent les un sur les autres en étant toujours en contact. Si les grains
deviennent plus petit, tout petits, on parle de sédiments. Le courant sera assez fort pour les faire s’envoler. Ils perdent le contact et peuvent se trouver dans toute
l’épaisseur du fluide. Ils sont en suspension. Bien entendu en pratique, on a un mélange des deux cas ”charriage” et ”suspension”, et il faudra réadditioner les deux
contributions. ici nous décomposons pour mieux modéliser les phénomènes.
Les sédiments sont maintenant supposés très légers et ils sont en suspension dans l’eau. Ils ne roulent plus au fond. Le mélange est constitué d’eau à la densité ρ
ρ −ρ
et de solide à la densité ρs . La densité totale est ρ = ρs c + (1 − c)ρf , ou ρ = ρf (1 + sρf f c) . On fait une approximation de type ”Boussinesq” comme en thermique
qui suppose toujours l’incompressibilité pour l’eau
∂u ∂v ∂h ∂Q
+ = 0, qui donne + =0
∂x ∂y ∂t ∂x
Dans le cadre simplifié choisi, la conservation de la quantité de mouvement pour les particules est résolue de manière simplifiée en supposant que cette vitesse est la
composition de deux phénomènes, une chute constante à la vitesse de chute Vf vers le bas et une diffusion aléatoire des sédiments autour de la vitesse d’entrainement
u, v + Vf . La conservation de la masse des particules devient
∂c ∂(cu) ∂(c(v − Vf )) ∂ ∂c ∂ ∂c
+ + = (( D ) + ( D ))
∂t ∂x ∂y ∂x ∂x ∂y ∂y
Dire que les particules sont transportées à la vitesse d’entrainement u, v + Vf et diffusent autour de cette trajectoire revient à avoir résolu l’équation de quantité de
mouvement des particules.
La conservation de la masse des particules va donner en hypothèse couche mince, et en intégrant sur l’épaisseur
Z Z
∂ ∂ ∂c
cdy + cudy = −Vf c(x, f ) − D |f
∂t ∂x ∂y
on pourra ainsi définir C = h1 cdy et réécrire cette équation en supposant que cudy = CQ (avec un facteur de forme pris égal à 1) et en posant E = −D ∂y ∂c
R R
|f et
comme c(x, f ) ∝ C on peut écrire D = Vf C (avec un facteur de forme encore pris égal à 1) ,
∂Ch ∂CQ ∂C ∂C 1
+ = (E − D), ou avec Q = hu, on a aussi +u = (E − D).
∂t ∂x ∂t ∂x h
l’excès de masse (qui intervient par la concentration c, c = 0 pas d’excès, pas de solide, et attention c << 1) qui joue le rôle de l’excès de température (on sait que
ρ −ρ
l’approximation de Boussinesq thermique est à α(T − T0 ) 1) : l’analogue de la dilatabilité ρ = ρ0 (1 − α(T − T0 )) est donc ici ρf (1 + sρf f c). On conserve donc
l’incompressibilité pour l’eau. Et la contrepartie de l’équation de la chaleur est l’équation de conservation de la masse de sédiment...
Les équations dynamiques sont calquées sur Boussinesq. quantité de mouvement pour le fluide :
∂u ∂ ∂ ∂p ∂2u ∂2u
ρ( + u u + v u) = − + µ( 2 + 2 )
∂t ∂x ∂y ∂x ∂x ∂y
∂v ∂ ∂ ∂p ρs − ρf ∂2v ∂2v
( + u v + v v) = (− − g) − cg + ν( 2 + 2 )
∂t ∂x ∂y ρf ∂y ρf ∂x ∂y
En première approximation, on n’en garde que l’équilibre hydrostatique, mais avec le surpoids des sédiments
h
ρs − ρf ρs − ρf
Z
∂p
0 ' (− − g) − cg donc p(x, y, t) = ρf g(h(x, t) − y) − ρf cdy
ρf ∂y ρf ρf y
- ESV . 116-
Saint-Venant
• besoin d’un modèle pour le flux au fond, les sédiments sont resuspendus par une érosion E = q0 ((S − Sc )+ )β et se déposent (D) avec une vitesse de chute −Vf
• calcul du nouveau fond après érosion et sédimentation par actualisation du fond,
∂f
φ = −(E − D)
∂t
- ESV . 117-
Saint-Venant
On peut mélanger avec du charriage, Provenzale Balmforth ”Pattern of Dirt” chapitre 15 p 377 :
∂f ∂qc
φ + =D−E
∂t ∂x
R ∂C
en posant q = qc + cudy, ou en supposant le facteurde forme égal à 1 : q = qc + CQ/h et en négligeant −
∂t
∂f ∂q
φ + =0
∂t ∂x
- ESV . 118-
Saint-Venant
Huitième partie
Exemples de résolution
27.5 Principe de la Résolution couplée simplifiée
Ce type de formulation permet le calcul du déplacement de dunes éoliennes ou sous marines. Pour fixer les idées nous proposons ici un mécanisme très simplifié de
mouvement de dune sous marine.
Il s’agit d’une interaction entre un écoulement et le fond érodable sur lequel il s’écoule : la forme du fond gouverne l’écoulement, ce dernier modifie la forme du
fond par érosion et sédimentation. Ce problème interactif couplé, très compliqué, est simplifié en supposant que l’écoulement est quasistationnaire (la variation de la
forme du fond en fonction du temps est supposée lente). Donc, on travaille à l’échelle de temps lente, à chaque instant (dans cette échelle lente), la forme du fond est
fixée, l’écoulement est calculé. Il faut ici distiguer deux types d’écoulements, le cas subcritique à nombre de Froude inférieur à 1, et le cas supercritique. Les deux cas
produisent des structures différentes !
Comme les phénomènes d’interaction fluide/sol se produisent par définition près de la paroi, le frottement pariétal est évalué à partir de l’expression de la vitesse.
C’est le frottement pariétal qui provoque l’érosion, les sédiments sont emportés par l’écoulement dans la couche limite. Ils se redéposent ensuite plus loin, modifiant
ainsi la forme du fond... et à l’instant suivant, le fond est donc différent... on continue...
En pratique, c’est lorsque l’écoulement est accéléré qu’il érode et lorsqu’il est freiné qu’il sédimente. Pour la dune, il y a accélération puis décélération, donc érosion
puis déposition : la dune se déplace dans le sens de l’écoulement. Pour l’antidune, c’est l’inverse, il y a freinage, donc déposition, puis accélération donc érosion.
L’anti-dune se déplace donc en remontant l’écoulement.
- ESV . 119-
Saint-Venant
- ESV . 120-
Saint-Venant
Figure 65 – Dune qui descend un courant Figure 66 – AntiDune qui remonte le courant
Nous avons déjà vu que l’on peut calculer la vitesse par la conservation du flux d’eau, et par Bernoulli (on suppose le fluide faiblement visqueux) la pression
hydrostatique (hypothèse de couche mince, soit
u(h0 + f − η) = u0 h0 , u2 + gη = u20 .
On en déduit ainsi en linéarisant au premier ordre (en supposant que la hauteur du fond est faible par rapport à la hauteur d’eau ) :
u = u0 (1 + (−1/(F 2 − 1))f ) + ... et η = (F 2 /(F 2 − 1))f ) + ... où F 2 = u20 /(gh0 )
On se souvient qu’un écoulement subcritique (F < 1) provoque un creux de la surface libre et une accélération de la vitesse. En revanche, un écoulement supercritique
provoque un bourrelet, et la vitesse diminue.
Puis, on calcule le flux de matière dont les variations sont proportionnelles à celles de τ en première approximation, et comme les variations de τ sont proportionnelles
à celles de la vitesse, on a au final que les variations de q sont proportionnelles aux variations de u. Soit K la constante de proportionnalité qui englobe tous les
paramètres. La conservation de la masse de sédiments (Relation d’Exner) : par perte par sédimentation et par gain par érosion,
∂f ∂q
φ=− .
∂t ∂x
compte tenu de l’expression de q en fonction de u et en fonction de η cela nous donne l’équation d’évolution du fond sous la forme d’une équation d’advection :
∂f ∂f
= −v .
∂t ∂x
avec pour la vitesse de déplacement du fond : v = K(−1/(F 2 − 1))/φ Dans ce cadre simplifié, la forme de f ne varie pas !
Interprétation :
si F < 1, la vitesse du fluide augmente quand la bosse croı̂t puis diminue quand elle décroı̂t, donc q fait de même. Donc, avant le sommet, il y érosion ( ∂q/∂x
est positif) ; après, il y a sédimentation ( ∂q/∂x est négatif). Donc, la bosse diminue avant le sommet, augmente après, le résultat global est un déplacement vers la
droite. Pour un écoulement subcritique F < 1, les dunes ont une vitesse positive,
C’est la conclusion inverse si F > 1 : la vitesse diminue avant le sommet, il y a donc déposition du sable provenant de l’amont, après le sommet, la vitesse
réaugmente et provoque donc une érosion. Le sable se déposera plus loin, on a alors un train d’antidunes qui remontent le courant d’eau (F > 1, elles ont une vitesse
négative).
Films d’antidunes :
[Link]
[Link]
- ESV . 121-
Saint-Venant
∂f ∂q 1
τ /(∆ρgd) = Ku(x + `), ls ∂x q + q = Ku(x + `), φ = − , et enfin u = u0 (1 + f)
∂t ∂x 1 − F2
comme on cherche des solutions en fk eσt+ikx , le déphasage est donc inférieur à la longueur d’onde on trouve alors σ(k). Le fond est instable lorsque la partie réelle de
σ est positive
−ikKeik`
σfk = −ikKeik` fk /(ikls + 1) donc σ =
ikls + 1
Pour 0 < k < kc le fond est donc instable car Re(σ) > 0. L’instabilité vient de l’existence du déphasage entre le frottement et la vitesse. On a ensuite atténuation.
L’atténuation est causée par ls , l’effet de pente Λ ∂f
∂x donnerait le même résultat d’atténuation. Ce modèle est un peu simple mais donne une idée de ce qui se passe
dans la réalité.
28 Conclusion
Nous avons vu de manière sommaire les différents concepts utiles pour comprendre le déplacement d’un fond érodable : le nombre de Shields permet de mettre en
mouvement les particules, il y a alors un flux de matériaux, la loi de conservation d’Exner permet de modifier le fond. La mécanique des fluides par les approximations
de type Saint-Venant permet de calculer l’écoulement. Les modèles sont ici un peu simplifiés mais il semble que dans la Nature, les outils proposés permettent de
commencer à comprendre les phénomènes.
- ESV . 122-
Saint-Venant
écoulement cisaillé pur
fluid
! 1
f 0
Σ
-1
-2
Figure 67 – Sur un fond de rides, le fluide va de gauche à droite, le frottement Figure 68 – taux d’amplification (partie réelle de σ = −ikKe
ikls +1 ) des rides en
est en avance par rapport au sommet de la bosse. C’est ce qui produit le flux fonction de la longueur d’onde. Dans le cas du sable dans l’eau kmax correspond
représenté par les flèches sur les rides. On voit qu’il y a un apport de matière au à environ 10cm.
92
Carry Le Rouet / 15 juin de
sommet 2005la ride
... ...venant de sa face gauche, et que la face doitr est en revanche
u=y
Figure 69 – Un exemple de ”dune” numérique. Solution 3D des équations simplifiées dans le cas d’un écoulement cisaillé. On observe la formation des cornes
caractéristique des barkhanes. animation
- ESV . 123-
Saint-Venant
Références
[1] Bagnold R. A. ”The physics of blown sand and desert dunes” Dover 1954 red. 2005
[2] H. Chanson ”The hydraulics of open channel flow : an introduction ; basic principles” Elsevier Butterworth-Heinemann, Amsterdam, The Netherlands. 2004. 2nd
ed. ISBN 0 7506 5978 5. 585 p google book
[3] F. Charru (2007) Instabilités hydrodynamiques EDP Sciences google book
[4] Darrozès JS, François C. Monavon A. ”Exercices de Mécanique des fluides avec Solutions” Ed ENSTA
[5] Daubert 1964,Quelques aspects de la propagation des crues - La Houille Blanche [Link]
[6] Pierre-Gilles de Gennes ”Avalanches, torrents, rivières”, Cours Collège de France, 01/99 - 02/99
[7] Feynman (1979) ”Mécanique 2” Interédition 391p (Voir chapitre 51)
[8] A. C. Fowler (2004) : ”Mathematics and the environment”. Mathematical Institute lecture notes. (Revised edition, September 2006.)
[Link]
[9] Paul Germain, Mécanique - Tome 1 ISBN : 9782729886530
[10] Guyon E., Hulin J.-P. & Petit L. (1991) ”Hydrodynamique Physique”, InterEditions, Ed. CNRS 506p, (voir page 208 et suivantes)
EDP Science 2001 pour la deuxième édition
[11] Guyon E., Hulin J.-P., Petit L. & Mitescu C.D. (2001) : ”Physical hydrodynamics” Oxford University Press (version anglaise)
[12] J.-M. Hervouet (1991) ”Une présentation des équations de Saint-Venant” document technique EDF HE 43/91-20 40p
[13] Jean-Michel Hervouet, ”Hydrodynamique des écoulements à surface libre : Modélisation numérique avec la méthode des éléments finis” (Broché) presses ponts et
chaussées
[14] J.-M. Hervouet (1991) ”Présentation du Système TELEMAC” document technique EDF HE 43/96/039/1p
[15] Ho Valance Dupont Ould El Moctar Sacling laws in eolian transport PRL 106, 094501 (2011)
[16] Lancaster Nicholas (2005) Geomorphology of Desert Dunes
[17] M. J. Lighthill and G. B. Whitham On Kinematic Waves. I. Flood Movement in Long Rivers Proc. R. Soc. Lond. A 1955 229, doi : 10.1098/rspa.1955.0088,
published 10 May 1955.
[18] . Miller 1983 USGS ”Basic Concepts of Kinematic-Wave Models” [Link]
[19] A. Monavon 2009 ”Cours de Mécanique des Fluides M1” UPMC.
[20] Nielsen ”Coastal bottom boundary layers and sediment transport” sur googb
[21] R.R. Long ”Stratified flows” [Link] et le film en streaming ou en [Link] youtube
[22] Paterson A.R. (1983) : ”A first course in fluid dynamics”, Cambrige, 528p (voir chapitre XV)
[23] A. Petitjean (1998) Recent Progresses in Dam-Break Modelling in France, [Link] Case Studies : Dam Break,
Dam Break Modelling
[24] Gerd Sauermann Klaus Kroy, and Hans J. Herrmann Continuum saltation model for sand dunes PHYSICAL REVIEW E, Vol 64, 031305
[25] O. Thual [Link]
[26] E. Toro ”Riemann Solvers and Numerical Methods for Fluid Dynamics A Practical Introduction” 3rd edition, Springer
[27] Whitham ”Linear and Nonlinear Waves” google book
[28] K. Whipple [Link]
lecture-notes/4_sediment_transport_edited.pdf
[29] C.T. Yang, Sediment Transport : Theory and Practice McGraw-Hill, NewYork, 1995, p. 480.
- ESV . 124-
Saint-Venant
- ESV . 125-