src/test/rising.c

    Rising bubble

    A two-dimensional bubble is released in a rectangular box and raises under the influence of buoyancy. This test case was proposed by Hysing et al, 2009 (see also the FEATFLOW page).

    We solve the incompressible, variable-density, Navier–Stokes equations with interfaces and surface tension. We can solve either the axisymmetric or planar version. We can used standard or “reduced” gravity. We also test levelset interface tracking and a momentum formulation.

    #if AXIS
# include "axi.h" // fixme: does not run with -catch
#endif
#if MOMENTUM
# include "momentum.h"
#else
#include "navier-stokes/centered.h"
#if CLSVOF
# include "two-phase-clsvof.h"
#elif LEVELSET
# include "two-phase-levelset.h"
#else
# include "two-phase.h"
#endif
#endif
#if LEVELSET
# include "integral.h"
#else
# include "tension.h"
#endif
#if REDUCED
# include "reduced.h"
#endif

#ifndef LEVEL
# define LEVEL 8
#endif

    The boundary conditions are slip lateral walls (the default) and no-slip on the right and left walls.

    #if MOMENTUM
q.t[right] = dirichlet(0);
q.t[left]  = dirichlet(0);
#else
u.t[right] = dirichlet(0);
u.t[left]  = dirichlet(0);
#endif

int main() {

    The domain will span [0:2]\times[0:0.5] and will be resolved with 256\times 64 grid points.

      size (2 [1]);
  DT = 1. [0,1];
  init_grid (1 << LEVEL);

    Hysing et al. consider two cases (1 and 2), with the densities, dynamic viscosities and surface tension of fluid 1 and 2 given below.

      rho1 = 1000.[0], mu1 = 10.;  // works also with rho1 = [-3,0,1]
#if CASE2
  rho2 = 1., mu2 = 0.1;
#else
  rho2 = 100., mu2 = 1.;
#endif

#if LEVELSET
  #if CASE2
  const scalar sigma[] = 1.96;
  #else
  const scalar sigma[] = 24.5;
  #endif
  d.sigmaf = sigma;
#else // !LEVELSET
  #if CASE2
  f.sigma = 1.96;
  #else
  f.sigma = 24.5;
  #endif
#endif // !LEVELSET

    We reduce the tolerance on the Poisson and viscous solvers to improve the accuracy.

      TOLERANCE = 1e-4 [*];
#if REDUCED
  G.x = -0.98;
  Z.x = 1.;
#endif
  run();
}

event init (t = 0) {

    The domain is a rectangle. We only simulate half the bubble.

      mask (y > 0.5 ? top : none);

    The bubble is centered on (0.5,0) and has a radius of 0.25.

    #if LEVELSET
  foreach()
    d[] = sqrt (sq(x - 0.5) + sq(y)) - 0.25;
#else
  fraction (f, sq(x - 0.5) + sq(y) - sq(0.25));
#endif
}

    We add the acceleration of gravity.

    #if !REDUCED
event acceleration (i++) {
  face vector av = a;
  foreach_face(x)
    av.x[] -= 0.98;
}
#endif

    A utility function to check the convergence of the multigrid solvers.

    void mg_print (mgstats mg)
{
  if (mg.i > 0 && mg.resa > 0.)
    printf ("%d %g %g %g %d ", mg.i, mg.resb, mg.resa,
	    mg.resb > 0 ? exp (log (mg.resb/mg.resa)/mg.i) : 0.,
	    mg.nrelax);
}

    We log the position of the center of mass of the bubble, its velocity and volume as well as convergence statistics for the multigrid solvers.

    event logfile (i++) {
  double xb = 0., vb = 0., sb = 0.;
  foreach(reduction(+:xb) reduction(+:vb) reduction(+:sb)) {
    double dv = (1. - f[])*dv();
#if MOMENTUM
    vb += q.x[]*dv/rho(f[]);
#else
    vb += u.x[]*dv;
#endif
    xb += x*dv;
    sb += dv;
  }
  static double sb0 = 0.;
  if (i == 0) {
    printf ("t sb -1 xb vb dt perf.t perf.speed\n");
    sb0 = sb;
  }
  printf ("%g %g %g %g %g %g %g %g ", 
	  t, (sb - sb0)/sb0, -1., xb/sb, vb/sb, dt, perf.t, perf.speed);
#if !MOMENTUM
  mg_print (mgp);
  mg_print (mgpf);
  mg_print (mgu);
#endif
  putchar ('\n');
  fflush (stdout);
}

    At t=3 we output the shape of the bubble.

    event interface (t = 3.) {
  output_facets (f, stderr);
}

#if ADAPT
event adapt (i++) {
  adapt_wavelet ({f,u}, (double[]){5e-4,1e-3,1e-3}, LEVEL);
}
#endif

    Results

    The final shape of the bubble is compared to that obtained with the MooNMD Lagrangian solver (see the FEATFLOW page) at the highest resolution. We also display the shape of the axisymmetric version of the test. The axisymmetric bubble moves much faster.

    set term push
set term @SVG size 640,320
set size ratio -1
set grid
plot [][0:0.4]'../c1g3l4s.txt' u 2:($1-0.5) w l t 'MooNMD', \
              'log' u 1:2 w l t 'Basilisk', \
              '../rising-levelset/log' u 1:2 w l t 'Basilisk (levelset)', \
              '../rising-clsvof/log' u 1:2 w l t 'Basilisk (CLSVOF)', \
              '../rising-axi/log' u 1:2 w l t 'Basilisk (axisymmetric)', \
              '../rising-axi-clsvof/log' u 1:2 w l t 'Basilisk (axi + CLSVOF)', \
              '../rising-axi-momentum/log' u 1:2 w l t 'Basilisk (axi + momentum)'
    Bubble shapes at the final time (t=3) for test case 1. (script)

    Bubble shapes at the final time (t=3) for test case 1. (script)

    For test case 2, the mesh in Basilisk is too coarse to accurately resolve the skirt.

    set key bottom left
plot [][0:0.4]'../c2g3l4s.txt' u 2:($1-0.5) w l t 'MooNMD', \
              '../rising2/log' u 1:2 w l t 'Basilisk', \
              '../rising2-levelset/log' u 1:2 w l t 'Basilisk (levelset)', \
              '../rising2-clsvof/log' u 1:2 w l t 'Basilisk (CLSVOF)'
    Bubble shapes at the final time (t=3) for test case 2. (script)

    Bubble shapes at the final time (t=3) for test case 2. (script)

    The agreement for the bubble rise velocity with time is also good.

    set term pop
reset
set grid
set xlabel 'Time'
set key bottom right
plot [0:3][0:]'../c1g3l4.txt' u 1:5 w l t 'MooNMD', \
              'out' u 1:5 w l t 'Basilisk', \
              '../rising-levelset/out' u 1:5 w l t 'Basilisk (levelset)', \
              '../rising-clsvof/out' u 1:5 w l t 'Basilisk (CLSVOF)',     \
              '../rising-axi/out' u 1:5 w l t 'Basilisk (axisymmetric)',  \
              '../rising-axi-clsvof/out' u 1:5 w l t 'Basilisk (axi + CLSVOF)',  \
              '../rising-axi-momentum/out' u 1:5 w l t 'Basilisk (axi + momentum)'
    Rise velocity as a function of time for test case 1. (script)

    Rise velocity as a function of time for test case 1. (script)

    reset
set grid
set xlabel 'Time'
set ylabel '(vb - vb_0)/vb_0'
set key bottom left
plot [0:3]'out' u 1:2 w l t 'Basilisk', \
          '../rising-levelset/out' u 1:2 w l t 'Basilisk (levelset)',   \
          '../rising-clsvof/out' u 1:2 w l t 'Basilisk (CLSVOF)',	\
          '../rising-axi/out' u 1:2 w l t 'Basilisk (axisymmetric)',	\
          '../rising-axi-clsvof/out' u 1:2 w l t 'Basilisk (axi + CLSVOF)',	\
          '../rising-axi-momentum/out' u 1:2 w l t 'Basilisk (axi + momentum)'
    Relative volume difference as a function of time for test case 1. (script)

    Relative volume difference as a function of time for test case 1. (script)

    reset
set grid
set xlabel 'Time'
set key bottom right
plot [0:3][0:]'../c2g3l4.txt' u 1:5 w l t 'MooNMD', \
              '../rising2/out' u 1:5 w l t 'Basilisk', \
              '../rising2-levelset/out' u 1:5 w l t 'Basilisk (levelset)', \
              '../rising2-clsvof/out' u 1:5 w l t 'Basilisk (CLSVOF)'
    Rise velocity as a function of time for test case 2. (script)

    Rise velocity as a function of time for test case 2. (script)

    reset
set grid
set xlabel 'Time'
set ylabel '(vb - vb_0)/vb_0'
set key top left
plot [0:3]'../rising2/out' u 1:2 w l t 'Basilisk',		       \
          '../rising2-levelset/out' u 1:2 w l t 'Basilisk (levelset)', \
          '../rising2-clsvof/out' u 1:2 w l t 'Basilisk (CLSVOF)'
    Relative volume difference as a function of time for test case 2. (script)

    Relative volume difference as a function of time for test case 2. (script)