sandbox/M1EMN/Exemples/ressaut_mascaret.c

    From hydrolic jump to undular bores, “du ressaut au mascaret”

    The problem

    how an hydrolic jump moves and is deformed by dispersion. This is done with the Green-Naghdi equations and with SV with Boussinesq correction.

    The left level of water is
    h_1, velocity is u_1, the froude number is F_1=u_1^2/(gh_1), this jump goes on a level h_2 of still water (u_2=0).

    We want to see how a hydrolic jump changes into an ondular bore.

    Dispersive terms

    The full system with Boussinesq correction reads: \displaystyle \frac{h^{n+1}-h^n}{\Delta t}+ \frac{\partial Q}{\partial x}=0 \displaystyle \frac{Q^{n+1}-Q^n}{\Delta t}+\frac{\partial }{\partial x}(\frac{Q^2}{h}+\rho g \frac{h^2}{2})=\alpha h^3 \frac{\partial }{\partial t}\frac{\partial^2 u}{\partial x^2}

    we do here an extra splitting, we estimate first Saint Venant estimate h^* and Q^* without dispersion, next we add the extra dispersive term. We compare with the Serre Green-Nagdhi, we will compare with multilyer latter.

    Code

    if SV defined, then we use the saint-venant model. else we use serre-green-naghdi model.

    We define an extra DISP, when defined, we solve the saint-venant model plus an extra term (Boussinesq)

    #undef SV
    #define SV
    #define DISP
    
    
    #include "grid/multigrid1D.h"
    #ifdef SV
    #include "saint-venant.h"
    #else
    #include "green-naghdi.h"
    #endif

    extrapolation of Int Airy

    double IntgAiry(double x)
    { double ai;
        if((x>-4)&&(x<3)){
        ai=-0.47800749642926166 + (4 + x)*(0.11541826749089747 + (-3 + x)*
         (-0.0013612022054726482 + x*(-0.03341069433229331 + (2 + x)*
        (0.006480518707412483 + (-2 + x)*(0.003011122308321213 - 0.0008866964376086442*(3 + x))))));}
        else{ai=0/0.1;}
      return ai;
    }

    Neumann conditions

    u.n[left] = neumann(0);
    h[left] = neumann(0);
    u.n[right] = neumann(0);
    h[right] = neumann(0);
    double h0,h1,h2,W,u1,u2,U1,U2,tmax;
    scalar F[],f[];
    F[left] = neumann(0);
    f[left] = neumann(0);

    Main and parameters The domain is XXX long, centered on the origin. The problem is without dimensions. The problem is solved in one dimension with 2048 grid points.

    int main()
    {
      X0 = -20;
      L0 = 75;
      G = 1.;
      h0 = 1;
      h1 = h0;
      h2 = h0-.1;
    #ifdef SV
       N = 1024;
    #else
       gradient = NULL;
       N = 1024;
    #endif
        tmax = 45;  // L0=75 512 45
      run();
    }

    The initial conditions are h(x<0,t=0)=h_1>h_2 and u(x<0,t=0)=u_1 and for positive x: h(x>0,t=0)=h_2<h_1 and u(x>0,t=0)=u_2. we define U_1= -\sqrt{ g (\frac{h_1+h_2}{2})\frac{h_2}{h_1}} and U_2= -\sqrt{ g (\frac{h_1+h_2}{2})\frac{h_1}{h_2}}. they correspond to the standing jump: if u_1=U_1<0 (subcritical |u_1/\sqrt{g h_1}|<1) and u_2=U_1<0 (supecritical |u_2/\sqrt{g h_2}|>2), then W=0.

    If we do a moving jump, we impose any value to u_1, then the velocity after the shock is u_2=U_2 + (u_1 -U_1) and the shock velocity is W=(u1 -U1). Here we take as value of u_1=U_1 - U_2 in order to have a flow at rest after the shock u_2=0, the shock velocity is W=-U_2>0.

    event init (i = 0){
       foreach(){
           h[] = (h1+h2)/2 + (h2-h1)/2*tanh(x) ; // x < 0. ? h1 : h2;
         U1 = -sqrt(G/2 *(1+h2/h1)*h2);
         U2 = -sqrt(G/2 *(1+h1/h2)*h1);
         u1 = U1 + (- U2) ;
    //      u1 = U1;   // uncomment for standing jump
         W = 0 + (u1 -U1);
         u2 = U2 + (u1 -U1);
         u.x[]= (u1+u2)/2 + (u2-u1)/2*tanh(x) ; // x < 0. ? u1 : u2;
         zb[] = 0.;
       }
    }

    old value usefull for dispersion

    scalar un[];
    event init_u (i = 0) {
        foreach() {
            un[] = u.x[];
        }
    }
    
    #ifdef DISP

    Poor s man dispersive model

    Let say that pressure is the sum of the hydrostatic pressure plus a perturbation of it: p= \rho g h + \Pi. This perturbation is linked to the up to now neglected transverse momentum: -\frac{\partial \Pi}{\partial y} \simeq \rho \frac{\partial v}{\partial t}, as we linearise around a mean level of water say h_0. Hence, as v \simeq - y \frac{\partial u}{\partial x}, we subsitute in transverse momentum and integrate -\frac{\partial \Pi}{\partial y} to obtain \Pi \simeq \rho \frac{h_0^2}{2} \frac{\partial^2 u}{\partial t \partial x } we estimate teh extra pressure gradient that will modiy the longitudinal momentum: \displaystyle -\frac{\partial \Pi}{\partial x} \simeq \rho \frac{h_0^2}{2} \frac{\partial }{\partial t}\frac{\partial^2 u}{\partial x^2} which is integrated again across the depth \displaystyle - \int_0^h\frac{\partial \Pi}{\partial x} dy \simeq \alpha \rho h^2 \frac{\partial^3 u}{\partial t \partial x^2} the coefficient \alpha takes into acount the approximations. It may be tuned to match the dispersion relation, or to match with Boussinesq model, it gives then \alpha=1/3. The non hydrostatic part of the pressure gives an extra contribution to Saint Venant solver \displaystyle \alpha h^3 \frac{\partial }{\partial t}\frac{\partial^2 u}{\partial x^2} we do here an extra splitting we estimate first Saint Venant estimate h^* and Q^*: \displaystyle \frac{h^{*}-h^n}{\Delta t}+ \frac{\partial Q}{\partial x}=0 \displaystyle \frac{Q^{*}-Q^n}{\Delta t}+ \frac{\partial }{\partial x}(\frac{Q^2}{h}+ \rho g \frac{h^2}{2}) = 0 then as u^*=Q^{*}/h^* and u^n=Q^{n}/h^n, the longitudical momentum equation may be written as u^{*}-u^n = F and in the new source term (\alpha h^3 \frac{\partial }{\partial t}\frac{\partial^2 u}{\partial x^2})\Delta t the variation \frac{\partial u}{\partial t}\Delta t is replaced by the total increase (u^{n+1}-u^n) so that the source is: \alpha h^2 \frac{\partial^2 }{\partial x^2}(u^{n+1}-u^n) then, the splitting reads: \displaystyle u^{n+1}-u^n = F + \alpha h^2 \frac{\partial^2}{\partial x^2}(u^{n+1}-u^n) it allows to compute implicitely the increase in time of u \displaystyle (1 - \alpha h^2 \frac{\partial^2 }{\partial x^2}) (u^{n+1}-u^n) = F

    event disperssse(i++,i>1){

    Store old value and estimate of explicit variation of velocity after Saint Venant splitting u^{*}-u^n = F, the corrected velocity increase is then u^{n+1}-u^n =f where f =\left((1-\frac{h^2}{3}\frac{\partial^2 }{\partial x^2})\right)^{-1}F

        scalar uo[];
        foreach(){
          uo[] = un[];
          F[] =  (u.x[] -uo[]);
          f[] = F[];
        }
        boundary({F,f});
        
        double dh = change (u.x, un);
        if ( (i > 1 && dh < 1e-7) || i==10000000) {
            foreach()
            fprintf (stderr, "%g %g\n", x, h[]);
            return 1; /* stop */
        }

    compute total increase of u by inversion of (1 - \frac{h^2}{3} \frac{\partial^2 }{\partial x^2 } ) f = F

        double eps=.0,fn=0.,omega=.25,s;
        do
        { eps=0;
            foreach()
            {   s =   sq(h[])/3/sq(Delta);
                fn = (F[] + s*(f[-1] + f[1]))/(1+ 2*s);
                eps = eps + sq(f[]-fn);
                f[] = omega * fn + (1-omega) * f[] ;
            }
        }while(sqrt(eps)>1.e-15);
        boundary ({f});

    we update the velocity u(x, t+\Delta t) = u (x, t) + f

        foreach(){
         u.x[] = uo[] + (f[]);
         un[] = u.x[];
        }
      boundary({u.x,un});
    
    }
    #endif

    We use a simple implicit scheme to implement quadratic bottom friction by splitting i.e. \displaystyle \frac{d\mathbf{u}}{dt} = - C_f|\mathbf{u}|\frac{\mathbf{u}}{h} with C_f=.05.

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

    output

    event output1(t+=.5){
        static FILE * fp = fopen("maxht.OUT","w");
        double hmax;
        hmax=0;
        foreach()
         if(h[]>=hmax)hmax=h[];
        fprintf (fp, " %g %g \n",t,hmax);
    
    }
    
    
    #ifdef gnuX
    event output  (t += .1; t <= tmax) {
        fprintf (stdout, " reset \n set title \"t=%g\"\n",t);
        fprintf (stdout, " set label \"s\" at %g,%g \n",W*t,h0);
        fprintf (stdout, "p[:][:]  '-' u 1:2  not  w l linec 3, '' u 1:3  w l\n");
    #else
      event output (t += 1; t <= tmax) {
    #endif 	
    	
      foreach()
        fprintf (stdout, "%g %g %g %g %g \n", x , h[],u.x[], (x- W*t)/pow(t+dry,1./3),t);
    #ifdef gnuX
        fprintf (stdout, "e\n\n");
    #else    
        fprintf (stdout, "\n");     
    #endif 
    }
    
    event outputf (t = end) {
        foreach()
          fprintf (stdout, "%g %g %g %g %g \n", x , h[],u.x[], (x- W*t)/pow(t+dry,1./3),t);
        foreach()
          fprintf (stderr, "%g %g %g %g %g \n", x , h[],u.x[], (x- W*t)/pow(t+dry,1./3),t);
        }

    compilation

    qcc -g -O2 -DTRASH=1 -Wall  -DgnuX=1  -o ressaut_mascaret ressaut_mascaret.c -lm
    ./ressaut_mascaret 2> log | gnuplot
     
     
    make ressaut_mascaret.tst
    make ressaut_mascaret/plots
    make ressaut_mascaret.c.html
    
    
    source c2html.sh ressaut_mascaret
     

    Results

    The standing jump is dispersed, we see the formation of the ondular bore.

    set xlabel 'h'
    set ylabel 'z'
     plot [0:][0:1.5]'out' u 1:($5==5?$2:NaN)  w l t 't=05',\
     ''  u 1:($5==10?$2:NaN) w l t 't=10',\
     ''  u 1:($5==20?$2:NaN) w l t 't=20',\
     ''  u 1:($5==30?$2:NaN) w l t 't=30',\
     ''  u 1:($5==40?$2:NaN) w l t 't=40'
    Fluid depth profile. (script)

    Fluid depth profile. (script)

     set xlabel 't'
     set ylabel 'hmax'
     plot [0:][0:1.5]'maxht.OUT' u 1:2  w l
    max wave elevation (script)

    max wave elevation (script)

    iai(x)=((x<3)&&(x>-4)?-0.478+(4+x)*(0.115418+(-3+x)*(-0.00136120+x*(-0.033411+(2+x)*(0.006480+(-2+x)*(0.0030111-0.0008867*(3+x)))))) : NaN)
    iai(x)=((x<3)&&(x>-4)?-0.47800749642926166 + (4 + x)*(0.11541826749089747+(-3+x)*(-0.0013612022054726482+x*(-0.03341069433229331+(2+x)*(0.006480518707412483 + (-2+x)*(0.003011122308321213-0.0008866964376086442*(3+x)))))): NaN)
    p[-3:2]'out' u ($4):($5>10?$2:NaN) w lp,''u ($4):(-iai($4*4)+.33)*3/2.*.1+.9
    (script)

    (script)

     reset
     set pm3d map
     set palette gray negative
     unset colorbox
     set tmargin at screen 0.95
     set bmargin at screen 0.15
     set rmargin at screen 0.95
     set lmargin at screen 0.15
     set xlabel "x"
     set ylabel "t"
     unset key
     splot [][][.85:1.1]'out' u 1:5:2
    (script)

    (script)

    Bibliography

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

    OK v2