sandbox/M1EMN/Exemples/belangerdisp.c

    Saint Venant Shallow Water with surface tension effects

    Scope

    We test the Bélanger relation for velocity and wtaer height before and after a steady standing jump. Here we add Laplace pressue due to surface tension.

     set arrow from 4,0 to 4,1.2 front
 set arrow from 4,1 to 4,0 front
 set arrow from 8,0 to 8,2. 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 label  "L0" at 6,.15 front
 set label  "depth h(x,t)" at 2.9,1.4 front
 set label  "depth h2" at 8.2,2.2 front
 set label  "water" at 1.,.8 front
 set label  "air" at 1.4,2.05 front
 set xlabel "x"
 set ylabel "h"
 p [0:10][0:3]0 not,(x<5? 1+1*exp(1*(x-5))*cos(6*(x-5)) :2) w filledcurves x1 linec 3 t'free surface'
 reset
    configuration (script)

    Cases

    ideal fluid, no surface tension, no viscosity: the jump

    Conservation de la quantité de mouvement U_1^2h_1+g\frac{h_1^2}{2}= U_2^2h_2+g\frac{h_2^2}{2} Conservation de la masse: U_1h_1=U_2h_2

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

    viscous fluid, no surface tension: Burgers

    bla bla

    slightly viscous fluid, with surface tension: dispersive jump

    That is the case we look at

    Code

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

double tmax,Fr,deltah,x_rs,W,sigma,nu;
scalar K[],hxc[],hxcm[];
int main(){
  X0 = 0;
  L0 = 256;
  X0 = -200;
  N = 1024*2;
  Fr = 1.5;    // Nombre de Froude initial
  G = 1;
  tmax = 150;
  
  W = 0;       //vitesse du ressaut, ici 0
  sigma = 0.5;  // surface tension
  nu = 0.0625;      // additional viscosity
  x_rs=0;    // position nitiale du saut dans [X0;X0+L0]
  DT = 1.e-3;   // 0.25 5.e-3;
  run();
}

    sortie libre

    u.n[left]  = dirichlet(1.5);
h[left]    = dirichlet(1);

h[right]   = neumann(0);
u.n[right] = neumann(0);//radiation(Fr/((-1 + sqrt(1+8*Fr*Fr))/2));

K[right]   = dirichlet(0);//neumann(0);
hxc[right]  = dirichlet(0);//neumann(0);
hxcm[right] = dirichlet(0);//neumann(0);

K[left]    = neumann(0);
hxc[left]   = neumann(0);//dirichlet(0);
hxcm[right] = neumann(0);//dirichlet(0);

    initialisation: Belanger with h_1=1 and h_2 computed

    event init (i = 0)
{
  foreach(){
    zb[]=0.;
    //  Bélanger formula with 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)/2))/2;
    // vitesse associée par conservation du flux Fr/h + translation
     u.x[]=Fr/h[]+W;
    }
  boundary({zb,h,u});
#ifdef gnuX
  printf("\nset grid\n");
#endif
}

    Compute the approximate curvature:

    K = \frac{\partial^2 h}{\partial x^2} wigh gives Laplace pressure: p = -\sigma \frac{\partial^2 h}{\partial x^2}

    event courv(i++)
{

    Compute the a centered derivative:

    foreach()
    hxc[] = (eta[+1]-eta[-1])/2/Delta;
boundary ({hxc});

    mean value gives an estimation of the derivative with an error O(\Delta^2) (2 \frac{h(1) - h(-1)}{2 \Delta} + \frac{h(0) - h(-2)}{2 \Delta}+ \frac{h(2) - h(0)}{2 \Delta})/4 = h_x + \frac{5 \Delta^2}{12}h_{xxx}

    foreach()
    hxcm[] =  (hxc[1] + 2* hxc[]+ hxcm[-1])/4 ;
boundary ({hxcm});

    derivative h_{xx} + 1/2 \Delta h_{xxx} + (7/12) \Delta^2 h_{xxxx}+ (1/4) \Delta^3 h_{xxxxx}+... the curvature, second order derivative is:

    foreach()
    K[] = (hxcm[0] -  hxcm[-1])/(Delta);
boundary ({K});

    Curvature contribution to pressure \frac{\partial u}{\partial t}= + \sigma \frac{\partial K}{\partial x} \text{ or } \frac{\partial u}{\partial t}= + \sigma \frac{\partial^3h}{\partial x^3} 

    foreach()
    u.x[] += dt*(sigma*(K[1] - K[0])/Delta) ;
 boundary ({u.x});

    all the differences all together give: h_{xxx}+ (1/2) \Delta^2 h_{xxxxx}+... precision is O(\Delta^2).

    }

    Some viscous dissipation; at the entrance and at the output…. \frac{\partial u}{\partial t}= + \nu \frac{\partial^2 u}{\partial x^2} 

    event visq(i++)
{
foreach()
    K[] = ((x<X0+8)||(x>X0+L0-32))*nu*(u.x[0]-u.x[-1])/Delta;
boundary ({K});
    foreach()
    u.x[] += dt*((K[1] - K[0])/Delta) ;
 boundary ({u.x});
}

    Plots

    #ifdef gnuX
event plot (t<tmax;t+=0.1) {
    printf("set title 'Ressaut dispersif 1D  --- t= %.4lf '\n"
	 "p[%g:%g][-.25:2]  '-' 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 (%lf):(1+%lf) t'jump theo' \n",
           t,X0,X0+L0,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");
}
#endif
event end(t=tmax ) {
    foreach()
    fprintf (stderr,"%g %g %g %g %g\n", x, h[], u.x[], zb[], t);
}

    Run and Results

    Run

    To compile and run:

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

    Plot of jump

    set xlabel "x"
 set size ratio .33333
 p [:56][0:3] 'log' t'free surface' w l lc 3
    result, free surface (blue) (script)

    Plot of jump

    Plot of jump comparision with analytical result for Fr=1.5 and h1=1 we have h2=(((-1 + sqrt(1+8*Fr*Fr))/2)) 
     set xlabel "x"
 Fr=1.5
 h2=(((-1 + sqrt(1+8*Fr*Fr))/2))
  set size ratio .4
 p [-150:56][-1:3] 'log' t'free surface'w l lc 3,'' u 1:3 w l t'speed',\
'' u 1:4 t'zb' w l lc -1,(x>0?h2:NaN) t'h2'  w l lc 1,(x<0?1:NaN)  t'h1' w l lc 1
    result, free surface (blue) and bottom (black) (script)

    Liens

    Bibliographie

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