sandbox/M1EMN/Exemples/belangerturb.c

Saint- Venant (Shallow Water) Bélanger relation with friction and slope

Scope

We test the Bélanger relation for velocity and water height before and after a steady standing jump. We study the influence of friction (turbulent u^2) and driving angle of topography (-Z'_b>0) (gravity is tilted)

 set arrow from 2,4 to 3,1.5 front
 set label  "g" at 2.5,2 front
 set label  "gravity tilted" at 2.5,4 front
 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.6,1.5 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+x/25 :3) w filledcurves x1 linec 3 t'free surface'
 reset

hydrolic jump configuration (script)

Equations

Saint-Venant

\displaystyle \frac{\partial }{\partial t} h \; +\; \frac{\partial }{\partial x} uh=0 \displaystyle \frac{\partial }{\partial t} hu + \frac{\partial }{\partial x} \dfrac{(hu)^2}{h} + \frac{\partial }{\partial x}g\dfrac{h^2}{2} = - gh \frac{d }{d x} Z_b-C_f {|u|}u + \frac{\partial }{\partial x}(\nu_e h \frac{\partial }{\partial x}u) over an inclined plane Z_b=-Z'_b x with the -Z'_b positive angle of the tilted plane.

Démonstration

Bélanger

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}.

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}

Chézy

the equilibrium velocity \displaystyle 0= - g h Z'_b - C_fu^2 is the Chézy friction velocity u_{chezy}=\sqrt{-g h Z'_b/C_f}, which arises after the jump, where we have a long calm river.

As the flow raet is constant Q=u_{chezy}h, so h=\left(\frac{C_f Q}{g(-Z'_b)} \right)^{2/3}

Code

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

double tmax,Fr,deltah,x_rs,W,Zpb,Cf,uchezy,hmax,xF1=0,hF1=1;

int main(){
  X0 = 0;
  L0 = 5;
  N = 128*2;
  G = 1;
  tmax = 100;
  Fr = 2.5;    // Nombre de Froude initial
  Zpb = -.1;

  W = 0;  //vitesse du ressaut, ici 0
  x_rs=1.8;    // position nitiale du saut dans [X0;X0+L0]
  DT = 0.002;
    
 double CF[6]={.15, .17 , .18 , 0.2  , .25 , .3};  // loop on Froude
    
 for (int iF=0;iF< sizeof(CF) / sizeof(CF[0]);iF++){
        Cf = CF[iF];
        run();
        fprintf(stderr,"%lf %lf \n",Cf, xF1);
  }

}

sortie libre

h[left] = dirichlet(1);
h[right] = h[];
u.n[left] =dirichlet(Fr);
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+deltah*(1+tanh((x-x_rs)/.01))/2;
    // vitesse associé par conservation du flux Fr/h + translation 
     u.x[]=Fr/h[]+W;
      
      deltah = pow(sq( Fr*sqrt(G*1)*1)/(G*(-Zpb)/Cf),1./3);
      h[]= ((x-x_rs)<0?  1: deltah) ;
      
       //Q^2 = (hmax^3*G*(-Zpb)/Cf);
    }
  boundary({zb,h,u});
}

At each timestep

We use a simple explicit scheme for imposed slope (-Z'_b ) implicit scheme to implement quadratic bottom friction with C_f. i.e. \displaystyle \frac{\partial \mathbf{u}}{\partial t} = -Z'_b g h - C_f|\mathbf{u}|\frac{\mathbf{u}}{h}

event friction (i++) {
    foreach() {
        double a = h[] < dry ? HUGE : 1. + Cf*dt*norm(u)/(h[]);
        foreach_dimension()
        u.x[] = (u.x[] - G*Zpb*dt)/a;
    }
    boundary ((scalar *){u});
}

we put a NON pysical longitudinal dissipation (as in Razis et al. paper) \displaystyle \frac{\partial h \mathbf{u}}{\partial t} = \nu_e \frac{\partial }{\partial x } (h \frac{\partial u}{\partial x } )

 event friction_long (i++) {
   foreach_dimension() {
     face vector g[];
     scalar a = u.x;
     double nue=.0;
       foreach_face()
           g.x[] = nue*((h[] + h[-1,0])/2)*(a[] - a[-1,0])/Delta ;
       foreach ()
           if(h[]>dry){ u.x[] += dt/Delta*(g.x[1,0] - g.x[] + g.y[0,1] - g.y[])/((h[] + h[-1,0])/2);}
      boundary ((scalar *){u});
   }}

compute Chezy velocity associated to the slope and position of the jump for plots and analysis

event plot (t<tmax;t+=0.05) {
    hmax=0;
    foreach(){
        hmax=(h[]>hmax?h[]:hmax);
        xF1=(sq(u.x[])/(G*h[])<1?xF1:x);
        hF1=(sq(u.x[])/(G*h[])<1?hF1:h[]);
    }
    uchezy=sqrt(hmax*G*(-Zpb)/Cf);
   
#ifdef gnuX
    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 1:( $2*(-1 + sqrt(1+8*$3*$3/%g/$2))/2) t'hj' w l lt 6,\\\n"
     "'+' u (%lf):(%lf) t'jump theo' , %lf , %lf \n",
           t,X0,X0+L0,G,G,xF1,hF1,uchezy,Fr*1/hmax);
    foreach()
    printf ("%g %g %g %g %g\n", x, h[], u.x[], zb[], t);
    printf ("e\n\n");
#endif
}

final plot

#ifdef gnuX
#else
event end(t=tmax ) {
    foreach()
    fprintf (stdout,"%g %g %g %g %g\n", x, h[], u.x[], zb[], t);
}
#endif

Run and Results

Run

To compile and run:

 qcc  -DgnuX=1  belangerturb.c  -lm; ./a.out 2>log | gnuplot
 
 make belangerturb.tst
 make belangerturb/plots;
 
 
 source  ../c2html.sh belangerturb
 

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}

 set xlabel "x"
 p [:1][1:1.2] 'out' t'free surface' w l lc 3,1+ 0.2*x lc 2,1+.25*x lc 1,1+.3*x
 
 

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

 set xlabel "x"
 p [:][-1:5] 'out' t'free surface' w l lc 3,'' u 1:3 t'speed' w l, \
   '' u 1:4 t'zb' w l lc -1
 

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

Position of shock

 set xlabel "Cf"
 p [:][0:] 'log' t' x Jump' w lp
 

(script)

Liens

• http://basilisk.fr/sandbox/M1EMN/Exemples/belanger.c Cas normal

• http://basilisk.fr/sandbox/M1EMN/Exemples/belangerdisp.c Cas dispersif

• http://basilisk.fr/sandbox/M1EMN/Exemples/belangerturb.c Cas “turbulent”

Bibliographie

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

• Razis, Kanellopoulos, van der Weele, “Continuous hydraulic jumps in laminar channel flow” J. Fluid Mech. (2021), vol. 915, A8,