sandbox/yonghui/smalltest/bubble3d.c
Bubble rising 3D
(to get a precise results, it would caused a long simulation time, since we need MAX_LEVEL = 10+) it’s a simple version edit from the exemples on basilisk, intend to compared the results with axi.h: Bubble.c & rising.c
#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "navier-stokes/perfs.h"
#include "two-phase.h"
#include "tension.h"
#include "reduced.h"
#include "view.h"
#define Z0 0.5
#define tmax 3.
#define bubble(x,y,z) (sq(x) + sq(y - Z0) + sq(z)) - sq(0.25)
The main function .
int LEVEL = 7;
int MAX_LEVEL = 10;
int main() {
//We set the domain geometry and initial refinement.
size (2);
origin (-L0/2., 0, -L0/2.);
init_grid (64);
We set the physical parameters: densities, viscosities and surface tension.
rho1 = 1000., mu1 = 10.;
rho2 = 1., mu2 = 0.1, f.sigma = 1.96;
TOLERANCE = 1e-4;
run();
}
//cylinder bc, slip laterial wall
bid tube;
uf.n[tube] =0.;
pf[top]=dirichlet(0.);
p[top]=dirichlet(0.);
//no slip flow at two side wall
u.t[top] = dirichlet(0.);
u.t[bottom] = dirichlet(0.);
u.r[top] = dirichlet(0.);
u.r[bottom] = dirichlet(0.);
event init (t = 0) {
//we refine the mesh around the bubble,in order to get a better fraction field
refine (sq(x) + sq(y - Z0) + sq(z) - sq(0.3) < 0 && level < LEVEL);
fraction (f, bubble(x, y, z));
//we also set a cylinder field
mask (sq(x)+sq(z) > sq(0.5) ? tube :none);
view(fov = 30.65,ty= - 0.1, width = 640, height = 640);
clear();
box();
draw_vof("f", edges = true, lw = 0.5);
save("f.png");
}
We add the acceleration of gravity (unity) in the downward (-y) direction.
event acceleration (t < tmax;i++) {
face vector av = a;
foreach_face(y){
av.y[] -= 1.;
}
}
Outputs
Every 100 iteration, we output a fraction field.
int j=0;
event bviewer (i += 5; t <= tmax){
j += 1;
char legend[1000];
char photo[40];
sprintf(legend, "t = %0.2g", t);
sprintf (photo, "pic-%d", j);
clear();
view (fov = 20, ty = -0.5, width = 720, height = 720);
draw_string(legend, 1, size = 30., lw = 2.);
squares ("f", min = 0, max = 1);
if ((j % 100) == 0)
save (file = photo);
save ("field.mp4");
clear();
view (fov = 20, ty = -0.5, width = 720, height = 720);
squares ("u.y", min = -0.5, max = 0.5);
save ("axivelo.mp4");
clear();
view (fov = 20, ty = -0.5, width = 720, height = 720);
squares ("u.x", min = -0.5, max = 0.5);
save ("planvelo.mp4");
}
int k=0;
event snapshot (t += 0.5)
{
k += 1;
char name[40];
char velofield[40];
sprintf (velofield, "ve-%d", j);
sprintf (name, "dump-%d", i);
dump (file = name);
output_field((scalar *){u}, fopen (velofield, "w"));
}
Every ten timesteps, we output the time, volume, position, and velocity of the bubble.
event loggfile (i++) {
double xb = 0., yb = 0., zb = 0., sb = 0.;
double vbx = 0., vby = 0., vbz = 0.;
foreach(reduction(+:xb) reduction(+:yb) reduction(+:zb)
reduction(+:vbx) reduction(+:vby) reduction(+:vbz)
reduction(+:sb)) {
double dv = (1. - f[])*dv();
xb += x*dv;
yb += y*dv;
zb += z*dv;
vbx += u.x[]*dv;
vby += u.y[]*dv;
vbz += u.z[]*dv;
sb += dv;
}
fprintf (ferr,
"%.8f %.8f %.8f %.8f %.8f %.8f %.8f %.8f\n",
t, sb, xb/sb, yb/sb, zb/sb,
vbx/sb, vby/sb, vbz/sb);
fflush (ferr);
}
We adapt the mesh by controlling the error on the volume fraction and velocity field.
event adapt (i++) {
double uemax = 1e-2;
// adapt_wavelet ({f,u}, (double[]){1e-3,uemax/10.,uemax/10.,uemax/10.}, MAX_LEVEL, 1 );
adapt_wavelet ({f,u}, (double[]){5e-3,uemax,uemax,uemax}, LEVEL, 3);
}
output
set grid
set xlabel 'Time'
set key bottom right
plot 'log' u 1:7 w l t'real 3D velocity'