src/test/seawall.c

    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 in one dimension.

    #include "grid/multigrid1D.h"
    #if SAINT_VENANT
    # include "saint-venant.h"
    #else
    # include "green-naghdi.h"
    #endif

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

    int main()
    {
      L0 = 15.;
      G = 9.81;
      N = 1 << 11;
      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 '', \
               '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 '', \
               '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 '', \
               '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 '', \
               '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 '', \
               '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 '', \
               '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 '', \
               'g40' w l lt 3 lw 2 t ''
    unset multiplot
    Comparison between experimental (dots) and numerical timeseries (script)

    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.

    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
    (script)

    (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.