sandbox/M1EMN/Exemples/damb_dressler.c

    Turbulent fluid “dam break problem”,

    The classical Shallow Water Equations in 1D with turbulent friction. \displaystyle \left\{\begin{array}{l} \partial_t h+\partial_x Q=0\\ \partial_t Q+ \partial_x \left[\dfrac{Q^2}{h}+g\dfrac{h^2}{2}\right] = -C_f u^2 \end{array}\right. with at t=0, a lake on the left h(|x|<0,t=0)=1 and zero water h(|x|>0,t=0)=0 for x>0. As the flow is a bit viscous, the convective term is important.

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

double tmax;

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

    start by a lake on the left and dry on the right, the soil is flat z_b=0

    event init (i = 0)
{
  foreach(){
    zb[] =  0;
    h[] =  (x<0);
    u.x[] =0;
   }
   boundary({zb,h,u});

#ifdef gnuX
  printf("\nset grid\n");
#endif    
}

    position of domain X0, length L0, no dimension G=1,
    run with 512 points (less is not enough)

    int main() {
  X0 = -5.;
  L0 = 10.;
  G = 1;
  N = 512;
  tmax=2;
  
  run();
}

    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. + .05*dt*norm(u)/h[];
    foreach_dimension()
      u.x[] /= a;
  }
}

    Output in gnuplot if the flag gnuX is defined

    #ifdef gnuX
event plot (t<tmax;t+=0.01) {
    printf("set title 'rupture de barrage en 1D --- t= %.2lf '\n"
      "h(x)=(((x-0)<-t)+((x-0)>-t)*(2./3*(1-(x-0)/(2*t)))**2)*(((x-0)<2*t));t= %.2lf  ; "
      "p[%g:%g][-.5:2]  '-' u 1:($2+$4) t'free surface' w l lt 3,"
      "'' u 1:($2*$3) t'Q' w l lt 4,\\\n"
      "'' u 1:4 t'topo' w l lt -1,\\\n"
      "'+' u (0):(0.296296296296296) t'theo Qmax ' , h(x) t 'theo h(x)'\n",
           t,t,X0,X0+L0);
    foreach()
    printf ("%g %g %g %g %g\n", x, h[], u.x[], zb[], t);
    printf ("e\n\n");
}

    else put the outputs in out

    #else
event plot(t<tmax;t+=.5 ) {
    foreach()
    printf ("%g %g %g %g %g\n", x, h[], u.x[], zb[], t);
}
#endif

    Run

    To compile and run with gnuplot:

      qcc -DgnuX=1   -O2   -o damb damb_dressler.c -lm;  ./damb_dressler | gnuplot

    To generate the plot with out

     qcc  -O2   -o damb_dressler damb_dressler.c -lm;  ./damb_dressler > out

    Plots

    Plot of the surface and flux at time t=1,2:

    comparison of heights with and without friction

    set xlabel "x"
h(x,t)=(((x-0)<-t)+((x-0)>-t)*(2./3*(1-(x-0)/(2*t)))**2)*(((x-0)<2*t))
u(x,t) =2./3*(1+x/t)+0.05*F1t(x,t)
F1t(x,t)=t*(-108./7/(2-x/t)**2. +12/(2-x/t)-8./3+8.*sqrt(3)/189*(2-x/t)**1.5)


F2t(x,t)= ((6./5)/(2-x/t)-2./3)*t+4*sqrt(3)/189.*((2-x/t)**(3./2))*t
 
#h(x,t)=((2./3*(1-(x-0)/(2*t))) +0.05*F2t(x,t) )**2.
p [-2:3][0:1] 'out' u 1:(($5==.5)||($5==1)||($5==1.5)? $2 :NaN)w p, h(x,.5),h(x,1.),h(x,1.5)
    result for h, comparison between computation in red and ideal fluid theory blue and green (script)

    result for h, comparison between computation in red and ideal fluid theory blue and green (script)

    using the Dressler 1952 solution for velocities u ans h: \displaystyle u(x,t)= 2./3*(1+x/t)+C_fF_1(x,t)t \displaystyle c= 1./3*(2-x/t) + C_fF_2(x,t)t \displaystyle h = c^2 with \displaystyle F_1(x,t)=(-\frac{(108./7)}{(2-x/t)^2} +\frac{12}{(2-x/t)}-8./3+8.\sqrt(3)/189*(2-x/t)^{3/2}) \displaystyle F_2(x,t)= ((6./5)/(2-x/t)-2./3)t+4(\sqrt(3)/189)((2-x/t)^{(3./2)}) See Dressler, see as well Chanson and Delestre (who put a 1./135 in F_2).

    set xlabel "x"
h(x,t)=(((x-0)<-t)+((x-0)>-t)*(2./3*(1-(x-0)/(2*t)))**2)*(((x-0)<2*t))
u(x,t) =2./3*(1+x/t)+0.05*F1t(x,t)
F1t(x,t)=t*(-108./7/(2-x/t)**2. +12/(2-x/t)-8./3+8.*sqrt(3)/189*(2-x/t)**1.5)
p[-2:4][0:2] 'out' u 1:(($5==.5)||($5==1)||($5==1.5)? $3 :NaN)w p, u(x,.5),u(x,1.),u(x,1.5)
    result for velocity, comparison between computation in red and Dressler theory blue and green (script)

    result for velocity, comparison between computation in red and Dressler theory blue and green (script)

    set xlabel "x"
h(x,t)=(((x-0)<-t)+((x-0)>-t)*(2./3*(1-(x-0)/(2*t)))**2)*(((x-0)<2*t))
c(x,t) =  sqrt(h(x,t)) +0.05*F2t(x,t)
F1t(x,t)=t*(-108./7/(2-x/t)**2. +12/(2-x/t)-8./3+8.*sqrt(3)/189*(2-x/t)**1.5)
F2t(x,t)=t*(6./5/(2-x/t) -2./3+4.*sqrt(3)/89*(2-x/t)**1.5)
hd(x,t)=c(x,t)*c(x,t)*((x-0)<2*t)
p[-2:4][0:2] 'out' u 1:(($5==.5)||($5==1)||($5==1.5)? $2 :NaN)w p, hd(x,.5),hd(x,1.),hd(x,1.5)
    result for h, comparison between computation in red and Dressler theory blue and green (script)

    result for h, comparison between computation in red and Dressler theory blue and green (script)

    Links

    see the same with friction and see non viscous dam break with standard C and with Basilisk

    Bibliographie

    Version 1: Paris 2018