seawall.c

Solitary wave overtopping a seawall

Solitary wave overtopping a seawall

This test case seeks to reproduce the experimental data of Hsiao and Lin, 2010 for a solitary wave overtopping a seawall. It was also studied by Lannes and Marche, 2014. Both nonlinear and dispersive effects are important. To illustrate this, we model the wave both with non-dispersive Saint-Venant equations and with the dispersive Green-Naghdi equations or layered model in one dimension.

#include "grid/multigrid1D.h"
#if SAINT_VENANT
# include "saint-venant.h"
#elif ML
# include "layered/hydro.h"
# include "layered/nh.h"
scalar h;
vector u;
#else
# include "green-naghdi.h"
#endif

The domain is 15 metres long. The acceleration of gravity is set to 9.81 m/s². The problem is discretised with 2048 grid points.

int main()
{
  L0 = 15.;
  G = 9.81;
  N = 1 << 11;
#if ML
  breaking = 0.15;
  TOLERANCE = 1e-4 [*];
  CFL_H = 0.5;
#endif
  run();
}

The initial conditions are for a Green-Naghdi soliton in a water depth of 0.2 metres and a relative soliton amplitude of 0.35 (type 1 wave of Hsiao and Lin, 2010).

double h0 = 0.2, A = 0.35;

double sech2 (double x) {
  double a = 2./(exp(x) + exp(-x));
  return a*a;
}

double soliton (double x, double t)
{
  double c = sqrt(G*(1. + A)*h0), psi = x - c*t;
  double k = sqrt(3.*A*h0)/(2.*h0*sqrt(h0*(1. + A)));
  return A*h0*sech2 (k*psi);
}

event init (i = 0)
{
  double c = sqrt(G*(1. + A)*h0);
  foreach() {
    double eta = soliton (x - 5.9, t);

Here we define the bathymetry (Fig. 1 of Hsiao and Lin). We need to shift the origin by 3 metres to match the shift in origin introduced by Hsiao and Lin when reporting computation results.

    x += 3.;
    zb[] = (x < 10. ? 0. :
	    x > 13.6 && x <= 13.9 ? 3.6/20. + (x - 13.6)*0.076/(13.9 - 13.6) :
	    x > 13.9 && x <= 13.948 ? 3.6/20. + 0.076 :
	    x > 13.948 && x <= 14.045 ? 
	    3.6/20. + 0.076 - (x - 13.948)*(0.076 - 0.022)/(14.045 - 13.948)
	    : (x - 10.)/20.)
      - h0;
    h[] = max (0., eta - zb[]);
    u.x[] = c*eta/(h0 + eta);
  }
}

We output timeseries of free surface elevation for the first 12 seconds for comparison with experimental records.

event timeseries (i++; t <= 12) {

We also add a quadratic bottom friction term of the form \displaystyle \partial_tu_x = - 10^{-2}|\mathbf{u}|u_x/h

  foreach() {
    double a = h[] < dry ? HUGE : 1. + 1e-2*dt*norm(u)/h[];
    foreach_dimension()
      u.x[] /= a;
  }

The coordinates of the various gauges are given in the legend of Figure 8 of Hsiao and Lin.

  static FILE * fp0 = fopen ("g0", "w");
  fprintf (fp0, "%g %g\n", t, interpolate (eta, 5.9, 0));
  fprintf (stderr, "%g %g\n", t, interpolate (eta, 7.6, 0));
  static FILE * fp10 = fopen ("g10", "w");
  fprintf (fp10, "%g %g\n", t, interpolate (eta, 9.644, 0));
  static FILE * fp22 = fopen ("g22", "w");
  fprintf (fp22, "%g %g\n", t, interpolate (eta, 10.462, 0));
  static FILE * fp28 = fopen ("g28", "w");

These gauges are on “dry” land and seem to use as reference the local height of the topography (rather than sealevel), presumably so that they measure “inundation depth” rather than wave elevation.

  fprintf (fp28, "%g %g\n", t, (interpolate (eta, 10.732, 0)
				- interpolate (zb, 10.732, 0)));
  static FILE * fp37 = fopen ("g37", "w");
  fprintf (fp37, "%g %g\n", t, (interpolate (eta, 11.005, 0)
				- interpolate (zb, 11.005, 0)));
  static FILE * fp40 = fopen ("g40", "w");
  fprintf (fp40, "%g %g\n", t, (interpolate (eta, 11.12, 0)
				- interpolate (zb, 11.12, 0)));
}

set term svg enhanced font ",10" size 640,800
set multiplot layout 7,1
set tmargin 1.
set bmargin 1.
set yrange [-0.02:0.1]
set ytics 0,0.05,0.1
set grid
set pointsize 0.5

plot [-2:2]'../wg1' pt 7 t 'WG1', \
           '../seawallsv/g0' w l lt 4 lw 1 t '', \
           '../seawall-ml/g0' w l lt 6 lw 1 t '', \
           'g0' w l lt 3 lw 2 t ''
plot [0:10]'../wg3' pt 7 t 'WG3', \
           '../seawallsv/log' w l lt 4 lw 1 t '', \
           '../seawall-ml/log' w l lt 6 lw 1 t '', \
           'log' w l lt 3 lw 2 t ''
plot [2:12]'../wg10' pt 7 t 'WG10', \
           '../seawallsv/g10' w l lt 4 lw 1 t '', \
           '../seawall-ml/g10' w l lt 6 lw 1 t '', \
           'g10' w l lt 3 lw 2 t ''
plot [2:12]'../wg22' pt 7 t 'WG22', \
           '../seawallsv/g22' w l lt 4 lw 1 t '', \
           '../seawall-ml/g22' w l lt 6 lw 1 t '', \
           'g22' w l lt 3 lw 2 t ''
plot [2:12]'../wg28' pt 7 t 'WG28', \
           '../seawallsv/g28' w l lt 4 lw 1 t '', \
           '../seawall-ml/g28' w l lt 6 lw 1 t '', \
           'g28' w l lt 3 lw 2 t ''
plot [2:12]'../wg37' pt 7 t 'WG37', \
           '../seawallsv/g37' w l lt 4 lw 1 t '', \
           '../seawall-ml/g37' w l lt 6 lw 1 t '', \
           'g37' w l lt 3 lw 2 t ''
plot [2:12]'../wg40' pt 7 t 'WG40', \
           '../seawallsv/g40' w l lt 4 lw 1 t '', \
           '../seawall-ml/g40' w l lt 6 lw 1 t '', \
           'g40' w l lt 3 lw 2 t ''
unset multiplot

Comparison between experimental (dots) and numerical timeseries (script)

The agreement is very satisfactory for this difficult problem and can be compared with the 2D Navier-Stokes VOF simulations of Hsiao and Lin (Figure 8 left column), as well as with the results of Lannes and Marche, 2014 (Figure 10). Dispersive effects are particularly clear for gauge 3 after 6 seconds and are well reproduced by the model.

The results for the Saint-Venant solver are the thin magenta lines and the results for the layered model are the thin dark blue lines.

We also output wave profiles at times corresponding to those of Figure 2 of Hsiao and Lin.

event profile (t = {2.63,2.89,3.01,3.19,3.35,3.71}) {
  char name[20]; sprintf (name, "p-%g", t);
  FILE * fp = fopen (name, "w");
  foreach()
    fprintf (fp, "%g %g %g %g\n", x, h[], u.x[], zb[]);
  fclose (fp);
}

reset

set term svg enhanced font ",10" size 640,800
set multiplot layout 6,2
set rmargin 2.
set lmargin 4.
set tmargin 0.
set bmargin 0.
set yrange [-0.1:0.2]
set ytics -0.1,0.1,0.2
set size ratio -1
unset key

plot [9.6:10.6]'p-2.63' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [9.6:10.6]'../seawallsv/p-2.63' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [9.6:10.6]'p-2.89' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [9.6:10.6]'../seawallsv/p-2.89' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'p-3.01' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'../seawallsv/p-3.01' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'p-3.19' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'../seawallsv/p-3.19' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'p-3.35' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'../seawallsv/p-3.35' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'p-3.71' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'../seawallsv/p-3.71' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1

unset multiplot

Green-Naghdi (left column) versus Saint-Venant (right column) (script)

The results are displayed both for the Green-Naghdi solver (left column) and for the Saint-Venant solver (right column). The Green-Naghdi results are again remarkably similar to the experimental snapshots of Hsiao and Lin (Figure 2) despite the complexity of the wave breaking process.

reset

set term svg enhanced font ",10" size 640,800
set multiplot layout 6,2
set rmargin 2.
set lmargin 4.
set tmargin 0.
set bmargin 0.
set yrange [-0.1:0.2]
set ytics -0.1,0.1,0.2
set size ratio -1
unset key

plot [9.6:10.6]'p-2.63' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [9.6:10.6]'../seawall-ml/p-2.63' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [9.6:10.6]'p-2.89' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [9.6:10.6]'../seawall-ml/p-2.89' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'p-3.01' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'../seawall-ml/p-3.01' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'p-3.19' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'../seawall-ml/p-3.19' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'p-3.35' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'../seawall-ml/p-3.35' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'p-3.71' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1
plot [10.2:11.2]'../seawall-ml/p-3.71' u 1:($2+$4) w l lc -1 lw 2, \
     '' u 1:4 w filledcu x1 lc -1

unset multiplot

Green-Naghdi (left column) versus layered (right column) (script)