multiriverinflow.c

Flow rates for multiple rivers

Flow rates for multiple rivers

In this example, we use the discharge() function in order to impose different flow rates on different rivers situated on the same boundary of a Saint-Venant simulation. We use a cartesian grid (non-adaptative).

#include "grid/cartesian.h"
#include "saint-venant.h"
#include "discharge.h"

The domain is 100 metres squared, centered on the origin. Time is in seconds.

#define LEVEL 7

int main(){
  L0 = 10.;
  X0 = - L0/2.;
  Y0 = - L0/2.;
  G = 9.81;
  N = 1 << LEVEL; // = pow(2,LEVEL);
  run();
}

Initial conditions

We chose a reasonably complicated river bed with two “valleys” (see Figure below). We allocate a new field river which is set to different values (1 and 2) for each of the rivers.

event init (i = 0){

We start with a dry riverbed, therefore the problem does not have a natural timescale the Saint-Venant solver can use. We set a maximum timestep to set this timescale.

  DT = 1e-2;

  foreach() {
    zb[] = 0.05*pow(x,4) - x*x + 2. + 0.2*(y + Y0);

The scalar river[] is now defined by default in the discharge.h library .We use it to store the location of each river. Its value is 1 on the river 1 (west side) and 2 on the river 2 (east side).

    river[] =  x < 0 ? 1 : 2;
  }
  boundary ({zb, river});
}

Boundary conditions

The tangential velocity on the top boundary u.t is set to zero and a neumann condition is fixed on the normal velocity.

u.n[top] = neumann(0);
u.t[top] = dirichlet(0);

u.n[bottom] = neumann(0);

To impose a given flow rate for the two rivers on the top boundary, we compute the elevation \eta of the water surface necessary to match this flow rate. This gives two elevation values eta1 and eta2, one for each river. The flow rates are set to 4 and 2 m³/sec for river 1 and 2 respectively.

double eta1, eta2;

event inflow (i++) {
  eta1 = discharge (4, top, 1);
  eta2 = discharge (2, top, 2);

Once we have the required values for the water surface elevations at the top boundary of both riverbeds, we impose them on both h and \eta=h+z_b.

  h[top] = max ((river[] == 1 ? eta1 : eta2) - zb[], 0);
  eta[top] = max ((river[] == 1 ? eta1 : eta2), zb[]);
}

Outputs

We compute the evolution of the water volumes in both riverbeds.

event volume (i += 10) {
  double volume1 = 0, volume2 = 0;
  foreach() {
      double dv = h[]*sq(Delta);
      if (x < 0) volume1 += dv;
      else volume2 += dv;
  }
  fprintf (stderr, "%g %g %g %g %g\n",
	   t, volume1, volume2, eta1, eta2);
}

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

event animation (t <= 1; t += 0.05) {
  double dx = L0/(1 << LEVEL), dy = dx;
  printf ("set title 't = %.3f'\n"
	  "sp [%g:%g][%g:%g][-5:5] '-'"
	  " u 1:2:($3+$4-.05) t 'free surface' w l lt 3,"
	  " '' u 1:2:4 t 'topography' w l lt 2\n",
	  t, X0, -X0, Y0, -Y0);
  for (double x = X0;  x <= X0 + L0; x += dx) {
    for (double y = Y0; y <= Y0 + L0; y += dy)
      printf ("%g %g %g %g\n",
	      x, y, interpolate (h, x, y),  interpolate (zb, x, y));
    putchar ('\n');
  }
  printf ("e\n"
	  "pause %.5lf \n\n", 0.);
}

Results

set key top left
set xlabel 'Time'
set ylabel 'Volume'
f(x)=a*x + b
fit [0.2:1] f(x) './log' u 1:2 via a,b
g(x)=c*x + d
fit [0.2:1] g(x) './log' u 1:3 via c,d
title1=sprintf("f(x) = %1.3f*t + %1.3f", a,b)
title2=sprintf("g(x) = %1.3f*t + %1.3f", c,d)
plot [0:1] './log' u 1:2 t 'river 1', \
     f(x) t title1, \
     './log' u 1:3 t 'river 2', \
     g(x) t title2

Evolution of water volumes in both rivers (script)

reset
set view 80,03
set xlabel 'x'
set ylabel 'y'
set zlabel 'z'
set hidden3d; unset ytics ; unset xtics
set term gif animate
set output 'movie.gif'
load './out'

Animation of the free surface. (script)