# sandbox/M1EMN/Exemples/ressaut_mascaret.c

## 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
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+ 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)

 set xlabel 't'
set ylabel 'hmax'
plot [0:][0:1.5]'maxht.OUT' u 1:2  w l 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)

 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)

# Bibliography

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

OK v2