sandbox/explosion.c

Two-dimensional explosion

We solve the Euler equations for a compressible gas.

#include compressible.h

We make boundary conditions free outflow.

w.y[top]    = neumann(0);
w.y[bottom] = neumann(0);
w.x[left]   = neumann(0);
w.x[right]  = neumann(0);

The domain spans [1:1]×[1:1].

#define LEVEL 7

int main() {
  origin (-1, -1);
  size (2.);
  init_grid (1 << LEVEL);
  run();
}

Initial conditions come from Toro’s book (Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd Edition, Springer Ed.) Chapter 17 section 17.1.1 are given in terms of density (ρ), pression (p), velocity (u) both at the left and right side of the discontinuity placed at R=0.4.

event init (t = 0)
{
  double R = 0.4 ;
  double rhoL = 1., rhoR = 0.125 ;
  double VmL = 0.0, VmR = 0.0 ;
  double pL = 1.0,  pR = 0.1 ;

Left and right initial states for ρ, w and energy E=ρu2/2+p/(γ1).

  foreach() {
    double r = sqrt(sq(x) + sq(y));
    double p;
    if (r <= R) {
      ρ[] = rhoL;
      w.x[] = w.y[] = VmL;
      p = pL;
    }
    else {
      ρ[] = rhoR;
      w.x[] = w.y[] = VmR;
      p = pR;
    }
    E[] = ρ[]*sq(w.x[])/2. + p/(gammao - 1.);
    w.x[] *= x*ρ[]/r;
    w.y[] *= y*ρ[]/r;
  }
}

event print (t = 0.25)
{

At t=0.25 we output the values of ρ and the normal velocity un as functions of the radial coordinate.

  foreach() {
    double r = sqrt(sq(x) + sq(y));
    double wn = (w.x[]*x + w.y[]*y)/r;
    printf ("%g %g %g\n", r, ρ[], wn/ρ[]);
  }

For reference we also output a cross-section at y=0.

  for (double x = 0; x <= 1; x += 1e-2)
    fprintf (stderr, "%g %.4f %.4f\n", x,
             interpolate (ρ, x, 0.),
             interpolate (w.x, x, 0.));
}
 web trial

On quadtrees, we adapt the mesh by controlling the error on the density field.

#if QUADTREE
event adapt (i++) {
  adapt_wavelet ({ρ}, (double[]){1e-5}, LEVEL + 1);
}
#endif

Results

Results are presented in terms of ρ and normal velocity un for Cartesian (7 levels) and adaptive (8 levels) computations. The numerical results compare very well with Toro’s numerical experiments.

Radial density profile

Radial density profile

Normal velocity

Normal velocity