# sandbox/M1EMN/Exemples/floodwave.c

# Resolution of flood wave numericaly and with characteristics

## Problem

how a shock forms with the transport equation for a flood wave

# Model equations

The flood wave or kinetic wave, see Whitham p 82, Fowler p 76, or Chanson, is a long wave approximation of shallow water equations were inertial and pressure gradient are neglected. THe waves go only downstream. We have only mass conservation, aband momentum equationb reduces to a balance between friction and mild slope only: \displaystyle \frac{\partial h}{\partial t} + \frac{\partial Q}{\partial x} =0\;\;\;\ \;\;\;0 = - g h Z'_b - c_f \frac{Q^2}{h^2} the slope -Z'_b>0 is constant, hence the flux is Q\sim h^{3/2}.

This is the turbulent case, in the laminar one Q= -gZ_b' h^3/(3 \nu) and mass equation is \displaystyle \frac{\partial h} {\partial t} + \frac{\partial }{\partial x}((\frac{-gZ_b'}{3\nu}) h_{\;}^3) =0, See example from viscous collapse on a slope (Huppert’s problem).

Then, for the turbulent case, without dimension we have to solve: \displaystyle \frac{\partial}{\partial \bar t} \bar h + \frac{\partial}{\partial \bar x}(\bar h^{3/2})= 0

the method is similar to advection which explains the notions of advection, testing the flux, coded with Basilisk. The flux is here \bar Q=\bar h^{3/2}, and we solve: \displaystyle \frac{\partial}{\partial \bar t} \bar h + \frac{\partial}{\partial \bar x} \bar Q = 0 \text{ solved as } \frac{\partial}{\partial \bar t} \bar h + \bar c \frac{\partial}{\partial \bar x} \bar h = 0 \text{ with } \bar c = \partial \bar Q/\partial \bar h

## Code

mandatory declarations:

```
#include "grid/cartesian1D.h"
#include "run.h"
```

definition of the field h, the flux Q, Boundary conditions

```
scalar h[];
scalar Q[];
h[left] = neumann(0);
h[right] = neumann(0);
```

the flux is \bar Q =\bar h^{3/2}. Let wfrite it \bar Q =\bar h^{m} so that we can test m=1 (pure advection) and m=2 (Burgers).

```
double m=3./2; // change for tests, m=1 advection, m=2 Burgers
double flux(double z)
{
return pow(fabs(z),m);
}
```

the velocity \bar c = \partial \bar Q/\partial \bar h:

```
double celerity(double z)
{
return fmin(m*pow(fabs(z),m-1.),10000);
}
```

Main with definition of parameters, note that time step is small

initial elevation: a constant level plus a gaussian perturbation

begin the time loop, print data, in practice a max time of 5 is enough.

```
event printdata (t += 1; t <=5)
{
foreach()
fprintf (stdout, "%g %g %g \n", x, h[], t);
fprintf (stdout, "\n\n");
}
```

integration step, at each time step

```
event integration (i++) {
double dt = DT;
double cDelta;
```

finding the good next time step

` dt = dtnext(dt);`

the algorithm is based on the flux. Approximation of the numerical flux taking into account \displaystyle Q_i = \frac{h_i^{3/2}+h_{i-1}^{3/2}}2 - c \frac{(h_i-h_{i-1})}2

```
foreach()
{
cDelta = (celerity(h[])+celerity(h[-1]))/2.;
Q[] = (flux(h[])+flux(h[-1]))/2. - cDelta *(h[]-h[-1])/2;
}
boundary ({Q});
```

explicit step update \displaystyle h_i^{n+1}=h_i^{n} -{\Delta t} \dfrac{Q_{i+1}-Q_{i}}{\Delta x}

` Boundary condition`

```
boundary ({h});
}
```

## Run

Then compile and run, either by `qcc`

either with `make`

:

```
qcc -g -O2 -DTRASH=1 -Wall floodwave.c -o floodwave ;./floodwave > out
make floodwave.tst
make floodwave/plots
make floodwave.c.html
source ../c2html.sh floodwave
```

## Results

`Time evolution with splot:```
reset
set xlabel "x"
set ylabel "t"
set zlabel "h"
sp [-5:10][0:5][0:2]'out' u 1:3:2 w l
```

The equation is of advection type, with \bar c = \partial \bar Q/\partial \bar h we have \displaystyle \frac{\partial}{\partial \bar t} \bar h + \frac{\partial}{\partial \bar x} \bar Q = 0 \text{ is } \frac{\partial}{\partial \bar t} \bar h + \bar c \frac{\partial}{\partial \bar x} \bar h = 0 In the x,t plane, along the lines \bar c = (d \bar x/d \bar t) the value of the solution is constant. This constant is the value in t=0, at the position x=\xi, given by function h_0. We write x=ct+\xi, and the value of the solution is constant along this line. This constant value along the line x-ct, it is h(x,t=0)=h(\xi,0). Given at time t=0 an initial value h(x,0)=h_0(x), then we construct the solution. Every x as an initial \xi, so that the solution is \displaystyle x=\xi +t c(\xi,0)\;\;\;\text{and}\;\; h(x,t)=h_0(x-t c(h_0(\xi))).

On the graph we plot the numerical result and the analitycal one.

```
reset
set xlabel "x"
set ylabel "h"
h(x)=.5+exp(-x*x)
c(x)=3./2*sqrt(h(x))
set parametric
set dummy x
p[:][:]'out' w l ,x+1*c(x),h(x) not w l lc -1,\
x+2*c(x),h(x) not w l lc -1,\
x+3*c(x),h(x) not w l lc -1,\
x+4*c(x),h(x) not w l lc -1,\
x+5*c(x),h(x) not w l lc -1, x ,0 not w l lc -1
```

We see the formation of the shock (red computed values) when the solution is no more a function (analytical solution in black).

## Links

## Bibliography

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

G. B. Whitham “Linear and Nonlinear Waves” Wiley-Interscience, p 82