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.

     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)

    hydrolic jump configuration (script)

    Démonstration

    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}

    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
      W = 0;  //vitesse du ressaut, ici 0
      x_rs=L0/2;    // position nitiale du saut dans [X0;X0+L0]
     
      run();
    }

    sortie libre

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

    initialisation: Belanger with h_1=1 and h_2 computed \displaystyle \big(\frac{h_2}{h_1}\big) = \frac{-1+\sqrt{1+8 F_1^2}}{2}

    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é par conservation du flux Fr/h + translation 
         u.x[]=Fr/h[]+W;
        }
      boundary({zb,h,u});
    #ifdef gnuX
      printf("\nset grid\n");
    #endif
    }
    
    #ifdef gnuX
    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
    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  -DgnuX=1  belanger.c -o belanger -lm; ./belanger | gnuplot
     
     make belanger.tst
     make belanger/plots;
     source  c2html.sh belanger
     

    Plot of jump

    Plot of jump comparision with analytical result \displaystyle \big(\frac{h_2}{h_1}\big) = \frac{-1+\sqrt{1+8 F_1^2}}{2} for Fr=2.5 and h1=1 we have h2=(((-1 + sqrt(1+8*Fr*Fr))/2))
    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)

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

    Liens

    Cas dispersif

    Bibliographie

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