# 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.

``````#if AXIS
# include "axi.h" // fixme: does not run with -catch
#endif
#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 $\left[0:2\right]×\left[0:0.5\right]$ and will be resolved with $256×64$ grid points.

``````  size (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.σ = 1.96;
#else
rho2 = 100., mu2 = 1., f.σ = 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 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.

``````  fraction (f, sq(x - 0.5) + sq(y) - 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``````

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();
xb += x*dv;
vb += u.x[]*dv;
sb += dv;
}
printf ("%g %g %g %g %g %g %g %g ",
t, sb, -1., xb/sb, vb/sb, dt, perf.t, perf.speed);
mg_print (mgp);
mg_print (mgpf);
mg_print (mgu);
putchar ('\n');
fflush (stdout);
}``````

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

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

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