belangerdisp.c

Saint Venant Shallow Water with surface tension effects
- Scope
- Cases
Code
- Compute the approximate curvature:
Run and Results
Liens
Bibliographie

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

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

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:

\displaystyle K = \frac{\partial^2 h}{\partial x^2} wigh gives Laplace pressure: \displaystyle 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) \displaystyle (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 \displaystyle 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 \displaystyle \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: \displaystyle h_{xxx}+ (1/2) \Delta^2 h_{xxxxx}+... precision is O(\Delta^2).

Some viscous dissipation; at the entrance and at the output…. \displaystyle \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