sandbox/bugs/rising_bug.c
Rising bug
This is a small variation on the rising bubble example to illustrate the problem. An inviscid drop freely falls under the action of gravity in a closed box. If the drop is initially tangenting grid cells, an instability develops as soon as the timestep is too small. In the following example, only half of the problem is simulated. If the flag -DNO_MIRROR is used, the full drop is simulated, but the problem persists. A workaround consists in introducing a small shift in the initial drop position (via -DSLIGHT_SHIFT or -DNO_MIRROR_SHIFTED).
Here are illustrated typical evolutions of the pressure field for the first 30 iterations:
#include "navier-stokes/centered.h"
#include "vof.h"
#include "tension.h"
#define rholiq 1000.
#define rhogas 10.
#define SIGMA 24.5
#define LEVEL 9
scalar f[];
scalar * interfaces = {f};
face vector alphav[];
scalar rhov[];
u.t[right] = dirichlet(0);
uf.t[right] = dirichlet(0);
u.t[left] = dirichlet(0);
uf.t[left] = dirichlet(0);
uf.n[top] = 0;
uf.n[bottom] = 0;
p[top] = neumann(0);
p[bottom] = neumann(0);
int main() {
size (2);
init_grid (1 << LEVEL);
alpha = alphav;
rho = rhov;
f.sigma = SIGMA;
TOLERANCE = 1e-4;
DT = 1e-6;
run();
}
event init (t = 0) {
vertex scalar phi[];
float small_amount = sqrt(3.0) / pi;
#ifdef SLIGHT_SHIFT
foreach_vertex()
phi[] = sq(x - 0.5 + small_amount*Delta) + sq(y) - sq(0.25);
#elif NO_MIRROR
foreach_vertex()
phi[] = sq(x - 0.5) + sq(y - 0.5) - sq(0.25);
#elif NO_MIRROR_SHIFTED
foreach_vertex()
phi[] = sq(x - 0.5 + small_amount*Delta) + sq(y - 0.5 + small_amount*Delta) - sq(0.25);
#else
foreach_vertex()
phi[] = sq(x - 0.5) + sq(y) - sq(0.25);
#endif
fractions (phi, f);
}
event acceleration (i++) {
face vector av = a;
foreach_face(x)
av.x[] -= 0.98;
}
#define rho(f) ((f)*rhogas + (1. - (f))*rholiq)
event properties (i++) {
foreach_face() {
double ff = (f[] + f[-1])/2.;
alphav.x[] = fm.x[]/rho(ff);
}
foreach()
rhov[] = cm[]*rho(f[]);
}
event movies (i++; i < 30) {
static FILE * fp = popen ("ppm2gif > closeup.gif", "w");
output_ppm (p, fp, box = {{0.15,0},{0.85,0.4}}, n = 256);
}