sandbox/M1EMN/Exemples/belanger.c

Saint Venant Shallow Water Bélanger relation

Scope

We test the Bélanger relation for velocity and wtaer height before and after a steady standing jump. The file is in french.

Nous vérifions les relation de Bélanger de part et d’autre d’un ressaut. Les relations sont établies pour un ressaut stationnaire de vitesse W=0.

 set arrow from 4,0 to 4,1 front
 set arrow from 4,1 to 4,0 front
 set arrow from 8,0 to 8,3. front
 set arrow from 0.1,.1 to 9.9,0.1 front
 set arrow from 9.1,.1 to 0.1,0.1 front
 set arrow from 1,0.75 to 3,0.75 front
 set arrow from 1,0.5 to 3,0.5 front
 set arrow from 1,0.25 to 3,0.25 front
 set arrow from 6,0.5 to 7,0.5 front
 set arrow from 6,1.50 to 7,1.5 front
 set arrow from 6,0.25 to 7,0.25 front
 set label  "L0" at 6,.15 front
 set label  "depth h1" at 3.2,1.1 front
 set label  "depth h2" at 8.2,3.2 front
 set label  "water" at 1.,.9 front
 set label  "air" at 1.4,2.05 front
 set xlabel "x"
 set ylabel "h"
 p [0:10][0:6]0 not,(x<5? 1 :3) w filledcurves x1 linec 3 t'free surface'
 reset

hydrolic jump configuration (script)

Démonstration

Dans le repère qui se déplace avec le ressaut à vitesse constante, le ressaut est fixe. On écrit les équations de conservation de part et d’autre du ressaut. Conservation de la quantité de mouvement \displaystyle U_1^2h_1+g\frac{h_1^2}{2}= U_2^2h_2+g\frac{h_2^2}{2}. Conservation de la masse: \displaystyle U_1h_1=U_2h_2,

à partir de ces relations de conservation nous allons écrire l’expression de U_2/U_1 et h_2/h_1 en fonction du nombre de Froude F_1^2=\frac{U_1^2}{gh_1}. Partant de la première multipliée par h_2 dans lequel U_2h_2 est remplacé par U_1 h_1 compte tenu de la deuxième \displaystyle U_1^2h_1h_2+g\frac{h_2h_1^2}{2}= U_1^2h_1^2+g\frac{h_2^2}{2}h_2 d’où en mettant les vitesses à gauche et la gravité à droite et comme on reconnaît une identité remarquable avec (h_2^2-h_1^2): \displaystyle U_1^2h_1(h_2- h_1) = g\frac{(h_2-h_1)(h_2+h_1)}{2} \;\;\text{ ou }\;\; U_1^2 = g\frac{ (h_2+h_1)h_2}{2h_1}. Le nombre de Froude est donc \displaystyle F_1^2=\frac{U_1^2}{gh_1} = \frac{ (h_2+h_1)h_2}{2h_1^2} si h_2>h_1 le nombre F_1 est supérieur à 1. On pourrait faire de même pour l’indice 1, on trouve alors \displaystyle F_2^2=\frac{U_2^2}{gh_2} = \frac{ (h_2+h_1)h_1}{2h_2^2}, cette fois, si h_2>h_1 le nombre F_2 est inférieur à 1. On constate que \displaystyle \frac{F_2}{F_1}=\big(\frac{h_1}{h_2}\big)^{3/2}. Reprenons l’équation de F_1^2 et exprimons la comme un polynôme en h_2/h_1 soit: \displaystyle \big(\frac{h_2}{h_1}\big)^{2} + \big(\frac{h_2}{h_1}\big) - 2 F_1^2=0 cette équation du second degré a une racine positive qui est: \displaystyle \big(\frac{h_2}{h_1}\big) = \frac{-1+\sqrt{1+8 F_1^2}}{2}. On a donc maintenant toutes les quantités avant et après le ressaut fixe: \displaystyle \big(\frac{h_2}{h_1}\big)= \big(\frac{U_1}{U_2}\big)= \big(\frac{F_1}{F_2}\big)^{2/3} = \frac{-1+\sqrt{1+8 F_1^2}}{2}.

Remarque: ces équations classiques sont établies pour un ressaut fixe. On peut faire bouger le ressaut en ajoutant une vitesse d’ensemble W qui sera la vitesse d’avancée du ressaut. En effet les équations de Saint Venant sans frottement sur fond plat sont invarianets par changement de référentiel Galiléen.

Code

#include "grid/cartesian1D.h"
#include "saint-venant.h"

double tmax,Fr,deltah,x_rs,W;

int main(){
  X0 = 0;
  L0 = 10;              
  N = 512;
  G = 1;
  tmax = 5;
  Fr = 2.5;      // Nombre de Froude initial à gauche avec indice 1
  W = 0;         //vitesse du ressaut, ici 0, on peut faire varier W
  x_rs= L0/2;    // position nitiale du saut dans [X0;X0+L0]
  DT = HUGE;
  FILE * fp = fopen ("WT.txt", "w");
  fclose(fp);    // fichier pour la position du ressaut 
  run();
}

sortie et entrée libres

h[left] = h[];
h[right] = h[];
u.n[left] = neumann(0);
u.n[right] = neumann(0) ;

initialisation: Bélanger avec h_1=1 et h_2 est calculé (attention compiler avec -disable-dimensions !) \displaystyle \big(\frac{h_2}{h_1}\big) = \frac{-1+\sqrt{1+8 F_1^2}}{2} Fr est F_1

event init (i = 0)
{
  foreach(){
    zb[]=0.;
    // formule de Bélanger avec h1=1;
    deltah=(((-1 + sqrt(1+8*Fr*Fr))/2)-1);
     h[]=1+(((-1 + sqrt(1+8*Fr*Fr))/2)-1)*(1+tanh((x-x_rs)/.01))/2;
    // vitesse associée par conservation du flux Fr/h + translation 
     u.x[]=Fr/h[]+W;
    }
 
#ifdef gnuX
  printf("\nset grid\n");
#endif
}

position du choc x_c

event findS(t<tmax;t+=0.2) {
   FILE * fp = fopen ("WT.txt", "a");
   double xc=0.;
// si la hauteur dépasse un peu h1=1, c'est le ressaut   
   foreach(){
     if(h[]>1.0001){  xc=xc;} else { xc=x;};
   }
    fprintf (fp, "%g %g   \n", t, xc );
    fflush(fp);
}
#ifdef gnuX
// plot à la volée avec gnuplot
event plot (t<tmax;t+=0.01) {
    printf("set title 'Ressaut en 1D --- t= %.2lf '\n"
   "p[%g:%g][-.5:5]  '-' u 1:($2+$4) t'free surface' w l lt 3,"
   "'' u 1:3 t'velocity' w l lt 4,\\\n"
   "'' u 1:4 t'topo' w l lt -1,\\\n"
     "'' u 1:(sqrt($3*$3/($2*%g))) t'Froude' w l lt 2,\\\n"
     "'' u 1:($2*%g+$3*$3/2) t'Charge' w l lt 5,\\\n"
     "'+' u (%lf):(1+%lf) t'jump theo' \n",
           t,X0,X0+L0,G,G,x_rs+W*t,deltah/2);
    foreach()
    printf ("%g %g %g %g %g\n", x, h[], u.x[], zb[], t);
    printf ("e\n\n");
}
#else
// affichage uniquemeny du dernier temps
event end(t=tmax ) {
    foreach()
    printf ("%g %g %g %g %g\n", x, h[], u.x[], zb[], t);
}
#endif

Run and Results

Run

To compile and run:

 qcc -disable-dimensions -DgnuX=1  belanger.c -o belanger -lm; ./belanger | gnuplot
 
 make belanger.tst
 make belanger/plots;
 source  c2html.sh belanger
 

Plot of jump

set xlabel "x"
 Fr=2.5
 h2=(((-1 + sqrt(1+8*Fr*Fr))/2))
 p [:][-1:5] 'out' t'free surface' lc 3,'' u 1:3 w p t'speed', '' u 1:4 t'zb' w l lc -1,(x>5?h2:NaN) t'h2'  w l lc 1,(x<5?1:NaN)  t'h1' w l lc 1
 

result, free surface (blue) and bottom (black) (script)

Plot de la position du choc en fonction du temps. Dans le cas de Bélanger W=0.
En jouant avec W dans le code, on voit que les ressauts mobiles (W\ne 0) sont par translation en W des ressauts fixes…

set xlabel "t"
set ylabel "jump position"
p[0:5][0:10] 'WT.txt' u 1:2 w lp

(script)

Liens pour aller plus loin

Cas dispersif, tension de surface explicite

http://basilisk.fr/sandbox/M1EMN/Exemples/belangerdisp.c

cas avec frottement turbulent

http://basilisk.fr/sandbox/M1EMN/Exemples/belangerturb.c

Cas dispersif, mascaret à partir des équations de Boussinesq (implicite) avec frottement turbulent

http://basilisk.fr/sandbox/M1EMN/Exemples/ressaut_mascaret.c

en multicouche laminaire….

http://basilisk.fr/sandbox/M1EMN/Exemples/higuera.c

Bibliographie

• Lagrée P-Y “Equations de Saint Venant et application, Ecoulements en milieux naturels” Cours MSF12, M1 Sorbonne-Université