src/examples/swasi.c

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"
#endif

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.);
#endif
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);
#endif
  init_grid (1 << LEVEL);
  run();
}

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);
}

Perturbation

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;
#endif
  }
}

Friction

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[]);
      foreach_dimension()
	u.x[] /= a;
    }
}

Outputs

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);
#else
  norm n = normf (u.x);
#endif
  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)

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)

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)
#endif
      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)
#endif
      fprintf (fp, "%g %g %g %g\n", x, eta[], zb[], u.x[]);
  fclose (fp);
}

reset
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)

Animation

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);
  clear();
  squares ("eta", min = -1.4, max = -0.2, linear = true);
  save ("eta.mp4");
}
#endif

References

• Videos
• PRL supplementary material