floodwave.c

Resolution of flood wave numericaly and with characteristics
- Problem
Model equations
- Code
- Run
- Results
- Links
- Bibliography

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

int main() {
    L0 = 12.;
    X0 = -2;
    N = 256*2;
    DT = (L0/N)/16/4;
    run();
}

initial elevation: a constant level plus a gaussian perturbation

event init (t = 0) {
    foreach(){
        h[] = .5+exp(-x*x) ;
    }
    boundary ({h});
}

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 = 1;

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}

foreach()
    h[] +=  - dt* ( Q[1] - Q[] )/Delta;
boundary ({h});

 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

(script)

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