A Shallow Water Analogue for the Standing Accretion Shock Instability

    We use the Saint-Venant solver to reproduce the experimental setup of Foglizzo et al. 2012. We use radial coordinates and Basilisk View to visualise the results.

    #include "grid/multigrid.h"
    #include "radial.h"
    #include "saint-venant.h"
    #if dimension > 1
    #  include "view.h"

    We use seven levels of refinement (i.e. 128^2 grid points) by default. The geometrical parameters and other inputs are those of the experimental setup. Note that we chose centimetres and seconds as length and time units.

    int LEVEL = 7;
    #define Ri    4.     // inner radius (cm)
    #define Ro    32.    // outer radius (cm)
    #define Ho    0.074  // injection layer thickness (cm)
    #define Q     1e3    // flow rate (cm^3/s)
    #define R_45  5.6    // hyperbolic surface scaling (cm)

    We double the equivalent viscosity compared to that given in Foglizzo et al. 2012 (0.03 cm^2/s), to obtain a more stable spiral mode. The “tube height” is set (by trial and error using 1D runs) so that the shock radius obtained numerically is close to 20 cm.

    #define NU    0.06   // equivalent viscosity (cm^2/s)
    #define STEP  6.98   // height of exit tube (cm)

    The injection velocity is the flow rate divided by the surface area of the outer injection ring.

    #define Uo (Q/(2.*pi*Ro*Ho))

    The maximum runtime (seconds).

    #define TMAX 150.

    Boundary conditions

    The outer injection ring corresponds to the right boundary of the radial coordinate system.

    Note that this code can also be run in one (radial) dimension by changing the first line to “grid/multigrid1D.h”.

    u.n[right] = - Uo;
    #if dimension > 1
    u.t[right] = dirichlet(0.);
    h[right] = Ho;

    The inner outflow is the left boundary.

    u.n[left] = dirichlet(- Ro*Uo*Ho/(STEP*Ri));

    Main program

    The level of refinement can be given as a command-line argument.

    int main (int argc, char * argv[])
      if (argc > 1)
        LEVEL = atoi (argv[1]);
      dtheta = 2.*pi;
      G = 981.; // acceleration of gravity (cm.s^-2)
      size (Ro - Ri);
      origin (Ri, 0.);
    #if dimension > 1
      periodic (top);
      init_grid (1 << LEVEL);

    Initial conditions

    We setup the hyperbolic bottom profile and set an initial uniform fluid layer and radial velocity identical to the injection conditions.

    event init (i = 0)
      foreach() {
        zb[] = - sq(R_45)/r;
        h[] = Ho;
        u.x[] = - Uo;
      zb[left]  = dirichlet(- sq(R_45)/r);
      zb[right] = dirichlet(- sq(R_45)/r);


    At t=40 seconds the shock is close to being stationary, we perturb the velocity field to trigger the instability.

    We add a systematic mode 1 perturbation of the radial component of the velocity, tapering off at R_i and R_o. This gives a clean growth of the mode and allows a simple estimate of the frequency and growth rate (see below).

    Optionally we can also add a small, constant azimuthal perturbation which will trigger the “spiral” instability (see second movie below).

    event perturb (t = 40) {
      foreach() {
        u.x[] += 1e-2*Uo*((r - Ri)/(Ro - Ri))*(1. - (r  - Ri)/(Ro - Ri))*sin(theta);
    #if 0
        u.y[] += 1e-2*Uo;


    We use a time-implicit discretisation of the linear friction.

    event friction (i++) {
      if (NU > 0.)
        foreach() {
          double a = h[] < dry ? HUGE : 1. + NU*dt/sq(h[]);
    	u.x[] /= a;


    We log the evolution of the fluid height and norms of the azimuthal (in 2D) or radial (in 1D) component of the velocity as well as the coordinates of the center-of-mass of the fluid.

    event logfile (i += 10; t <= TMAX) {
      stats s = statsf (h);
    #if dimension > 1  
      norm n = normf (u.y);
      norm n = normf (u.x);
      double sumx = 0., sumy = 0.;
      foreach (reduction(+:sumx) reduction(+:sumy))
        sumx += r*cos(theta)*h[]*dv(),
        sumy += r*sin(theta)*h[]*dv();
      fprintf (stderr, "%g %d %g %g %g %g %g %g %g\n", t, i, s.min, s.max, s.sum,
    	   n.rms, n.max, sumx/s.sum, sumy/s.sum);

    The evolution with time of the y-coordinate of the center of mass is illustrated below. It clearly shows the exponential growth of the oscillation of this position followed by non-linear saturation (see also the animation below).

    set xlabel 'Time (sec)'
    set ylabel 'Position (cm)'
    plot [0:150]'log' u 1:9 w l t ''
    Evolution of the y-coordinate of the center-of-mass of the liquid (script)

    Evolution of the y-coordinate of the center-of-mass of the liquid (script)

    This can be used to estimate the period of oscillation (2.74 s) and the initial growth rate (0.07 s^{-1}) as illustrated below

    set ylabel 'Absolute value of position (cm)'
    set logscale y
    set key top left
    set grid
    plot [40:150][1e-3:10]'log' u 1:(abs($9)) w l t 'numerical', \
         '' u 1:(exp(0.07*$1)/1050*abs(sin(2.*pi*$1/2.74+0.55))) w l t 'exp(0.07*x)*|sin(2.*pi*x/2.74)|'
    Same as above but using a logscale (script)

    Same as above but using a logscale (script)

    While we run, we display the evolution of a vertical cross-section of the flow.

    event profiles (i += 100) {
      static FILE * fp = popen ("gnuplot", "w");
      if (i == 0)
        fprintf (fp,
    	     "set term x11 noraise\n"
    	     "set grid\n"
    	     "set xrange [%g:%g]\n", X0, X0 + L0);
      fprintf (fp,
    	   "set title 't = %.2f'\n"
    	   "p '-' u 1:3:2 w filledcu lc 3 t ''\n", t);
      foreach (serial)
    #if dimension > 1
        if (y < Delta)
          fprintf (fp, "%g %g %g\n", x, eta[], zb[]);
      fprintf (fp, "e\n\n");
      fflush (fp);

    We save the stationary profile at t=50.

    event stationary (t = 50) {
      FILE * fp = fopen ("stationary", "w");
      foreach (serial)
    #if dimension > 1
        if (y < Delta)
          fprintf (fp, "%g %g %g %g\n", x, eta[], zb[], u.x[]);
      fclose (fp);
    set xlabel 'radius (cm)'
    set grid
    p 'stationary' u 1:3:2 w filledcu lc 3 t '', \
      '' u 1:(-$4/20.) w l lt 1 lw 2 t '-u_r/20'
    Stationary profile at t=50 (script)

    Stationary profile at t=50 (script)


    We create an animation of the free surface height \eta. We first define a coordinate mapping function which will be used by Basilisk View.

    #if dimension > 1
    void radial (coord * p) {
      double r = p->x, theta = p->y*dtheta/L0;
      p->x = r*cos(theta), p->y = r*sin(theta);
    event movie (t = 60; t += 0.1) {
      view (map = radial, fov = 45, width = 600, height = 600, samples = 1);
      squares ("eta", min = -1.4, max = -0.2, linear = true);
      save ("eta.mp4");

    SWASI instability. Growth and non-linear saturation of mode 1.

    SWASI instability. Spiral mode produced by adding a small azimuthal perturbation.



    Thierry Foglizzo, Frédéric Masset, Jérôme Guilet, and Gilles Durand. Shallow water analogue of the standing accretion shock instability: Experimental demonstration and a two-dimensional model. Physical Review Letters, 108(5):051103, 2012. [ http ]