fricfluv_rain_darcy.c

Rain in saint-venant : Fluvial Test case with rain and Darcy friction

Rain in saint-venant : Fluvial Test case with rain and Darcy friction

In this example, we reproduce the MacDonald benchmark^{[1]} and we add rain. We test the friction ans the rain in a fluvial flow. The leading equations can be written as : \displaystyle \begin{array}{cc} \partial_t h + \partial_x q = R(x,t) \\ \partial_t q + \partial_x \left(\frac{q^2}{h}+ \frac{1}{2} g h^2 \right) = -gh (S_0 + S_f) \end{array} S_0 is the slope, S_f is the friction term in its kinematic form and R(x,t) the rain intensity.

Declarations

We call the saint venant solver in 1D and we add the Darcy friction term and the rain.

#include "grid/cartesian1D.h"
#include "saint-venant.h"
#include "saintvenant/darcy.h"
#include "saintvenant/rain.h"

int LEVEL;

scalar e[];
double e1 = 0., e2 = 0., emax = 0.;
int ne = 0;
double pause, tmax, q0 = 1, z0, zf, h0, tb = 2000, intens = 1e-3;

// Q is not cst
double q(double x){
  return q0+x*intens;
}
// Analytical solution for h and dh/dx
double hex(double x) {
  return pow(4/G,1/3.)*(1 + 0.5*exp(-16*pow(-0.5 + x/1000.,2)));
}

double dhex(double x) {
  return -0.016*pow(4/G,1/3.)*exp(-16*pow(-0.5 + x/1000.,2))*(-0.5 + x/1000.);
}

// Darcy friction term in kinematic formulation
double sf(double x){
  return  -f/(8*G)*q(x)*q(x)/pow(hex(x),3);
}

// Z and dz/dx
double dzex(double x) {
  return (q(x)*q(x)/(G*pow(hex(x),3))-1)*dhex(x)-2*q(x)*intens/(G*sq(hex(x)))+sf(x);
}
double zex(double x, double z){
  double dx = L0/N;
  return z + dx/4.*(dzex(x-dx)+2*dzex(x-0.5*dx)+dzex(x));
}

Parameters

Definition of parameters and calling of the saint venant subroutine run().

int main()
{
  f = 0.093;
  pause=0.2;
  L0 = 1000.;
  X0 = 0;
  G = 9.81;
  tmax = 2500;
  for( LEVEL = 4; LEVEL <= 8; LEVEL++){  
    e1 = e2 = emax = 0.;
    ne = 0;
    N = 1 << LEVEL;
    run();
    fprintf (stderr, "%d %g %g %g\n", N, e1/ne, sqrt(e2)/ne, emax/ne);
  }
}

Boundary condition

We fix h and q (u) at both boundaries (fluvial case).

h[left] = dirichlet(max(hex(0),0));
eta[left] =  dirichlet(max(hex(0)+zb[],zb[]));
u.n[left] = dirichlet(max(q0/hex(0),0));

h[right] = dirichlet(max(hex(1000),0));
eta[right] =  dirichlet(max(hex(1000)+zb[],zb[]));
u.n[right] = dirichlet(max(q(L0)/hex(1000),0));

Initial conditions

event init (i = 0)
{
  // Because the slope is initially dry, we fix an artificial time-step. 
  DT = 1e-2;
  z0=0;
  foreach(){
    rain[] = intens;
    zb[] = zex(x,z0);
    z0=zb[];
    u.x[] = 0;
    h[]=0;
    zf=z0;
  }
  boundary(all);
}

Computing error

Noticing tb the time when the flow is already stationary, we define the following relative error norms [2]:

\displaystyle |h|_1 = \frac{\int_{tb}^T \int_0^L |h() - hex()| dx dt}{(T-tb)*L} \displaystyle |h|_2 = \frac{\sqrt{\int_{tb}^T \int_0^L (h() - hex()^2 dx dt}}{(T-tb)*L} \displaystyle |h|_{max} = \frac{\int_{tb}^T max(h() - hex())}{T-tb}

event error (i++; t<=tmax) {
  foreach()
      e[] = (h[] - hex(x));
  norm no = normf (e);
  if(t > tb){
    e1 += no.avg;
    e2 += no.rms*no.rms;
    ne++;
    if (no.max > emax)
      emax = no.max;
  }
   if( N == 128){
    static FILE * fp1 = fopen("ErrorN128.dat","w");
    fprintf(fp1,"%g \t %g \t %g \t %g \n",t,no.avg,no.rms,no.max);
  }
}

Gnuplot output

We can use gnuplot to produce an animation of the water surface.

/*
event plot ( t <= tmax; t += 10 ) {
    if( N == 32 || N == 128 ){
      printf("set title 'Manning friction fluvial ----- t= %.3g , N = %i '\n"
	     "p[%g:%g][-5:1]  '-' u 1:($2+$4) t'free surface' w p pt 1,"
	     "'' u 1:4 t'topo' w l lt 4,"
	     "'' u 1:5 t'Analytical' w l lt -1 \n",t,N,X0,X0+L0);
      foreach()
	printf (" %g %g %g %g %g %g\n", x, h[], u.x[], zb[],hex(x)+zb[], t);
      printf ("e\n"
	      "pause %.4g \n\n",pause);
    }
}
*/

Print the water profile along the channel at final time.

event printprofile ( t = tmax-1 ){
    char name[100];
    FILE * fp;
    sprintf(name,"profil-%i.dat",N);
    fp=fopen(name,"w");
    foreach() {
      fprintf(fp,"%g  \t %g  \t %g \t %g \t %g   \n"
	      ,x,h[],zb[],hex(x),u.x[]);
    }
    fclose(fp);
}

Results

Error convergence Error time Water depth profiles

References

[1] I. MacDonald, M. Baines, N. Nichols, and P. G. Samuels, “Analytic Benchmark Solutions for Open-Channel Flows,” no. November, pp. 1041–1045, 1997.

[2] S. Popinet, “Quadtree-adaptive tsunami modelling,” Ocean Dyn., vol. 61, no. January, pp. 1261–1285, 2011.