sandbox/piersonj/hysing3D.c

    Rising bubble

    This case is almost exactly the same as the one proposed by S. Popinet in 2D. All credit to him !

    A three-dimensional bubble is released in a box and raises under the influence of buoyancy. This test case is similar to the case 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 used standard or “reduced” gravity.

    #define REDUCED 1
#define CASE2 1
#define ADAPT 0

#include "grid/multigrid3D.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
#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.

    u.t[right] = dirichlet(0);
u.t[left]  = dirichlet(0);

    We make sure there is no flow through the top and bottom boundary, otherwise the compatibility condition for the Poisson equation can be violated.

    uf.n[bottom] = 0.;
uf.n[top] = 0.;

int main() {

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

      size (2);
  origin (0., -0.5, -0.5);
  dimensions(nx=4,ny=2,nz=2);
  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., mu1 = 10.;
#if CASE2
  rho2 = 1., mu2 = 0.1, f.sigma = 1.96;
#else
  rho2 = 100., mu2 = 1., f.sigma = 24.5;
#endif

    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 bubble is centered on (0.5,0,0) and has a radius of 0.25.

      fraction (f, sq(x - 0.5) + sq(y) + sq(z) - sq(0.25));
}

    We add the acceleration of gravity.

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

    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 extract (i++) {
  double xb = 0., vb = 0., sb = 0.;
  foreach(reduction(+:xb) reduction(+:vb) reduction(+:sb)) {
    double dv = (1. - f[])*dv();
    xb += x*dv;
    vb += u.x[]*dv;
    sb += dv;
  }
  double s = interface_area(f);

  char name[80];
  sprintf (name, "velocity_pos.dat");
  static FILE * fp = fopen (name, "w");
  fprintf (fp,"%g %g %g %g %g %g %g\n", 
	  t, sb, -1., xb/sb, vb/sb, s, dt);
  //putchar ('\n');
  //fflush (stdout);
  fflush (fp);
}

event logfile(i++) {
  printf ("i = %d t = %g\n", i,t);
  //printf ("utc = %g", utc);
  fflush(stdout);
}



event interface (t += 1; t <= 3) {

    We only output the interface in this function.

      char name[80];
  sprintf (name, "interface-%f-%d.txt", t,pid());
  FILE* fp = fopen (name, "w");
  output_facets (f, fp);
  fclose (fp);
   
  char named[80];
  sprintf (named, "dump-%g", t);
  dump (named);

}

    If gfsview is installed on the system, we can also visualise the simulation as it proceeds.

    #if 0
event gfsview (i += 10) {
  static FILE * fp = popen("gfsview2D rising.gfv", "w");
  scalar vort[];
  vorticity (u, vort);
  output_gfs (fp);
}
#endif

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