sandbox/wmostert/turbulent_bubble.c
This is a standard file to reproduce the simulations presented in:
Deformation and breakup of a bubble in a turbulent environment.
A fully developed turbulent flow is generated as a precursor simulation, based on the example case isotropic.c, itself based on the case investigated by Rosales & Meneveau, 2005; once fully developed, a gas bubble is inserted.
This simulation is in two parts. First, an equivalent version of isotropic.c is run. Second, the simulation is dumped and restarted, and the bubble is initialized in the flow. The turbulent flow may be forced or unforced in the liquid phase after bubble insertion (usually forced), but in any case there is no forcing inside the bubble.
Both parts are incorporated into this one problem file using preprocessor commands.
Include the Output_facets to have a good description of the interface
Include a way to restart the files. Should have a file dump file called dump.
Level in the bulk fixed at 8.
Write several dump files to avoid problems when you want to restart your simulations.
Initially flow inside the bubble put to 0.
#include "grid/octree.h"
#include "navier-stokes/centered.h"We use the \lambda_2 criterion and Basilisk View for visualisation of vortices.
Include surface tension and two-phase schemes.
#include "two-phase.h"
#include "navier-stokes/conserving.h"
#include "tension.h"We monitor performance statistics and control the maximum runtime.
#include "navier-stokes/perfs.h"
#include "maxruntime.h"Include tag functionalities for droplet/bubbles calculation.
#include "tag.h"Include density and viscosity ratios.
#define RHOR 850.0
#define MUR 100.0Define domain size.
#define WIDTH 120.0We need to store the variable forcing term.
face vector av[];The code takes the level of refinement as optional command-line argument (as well as an optional maximum runtime).
int maxlevel = 8; //only for the compilation
double MAXTIME = 20; //idem
int FORCED = 0; //equals 0 if it is not forced, 1 otherwise
double amp_force = 0.1; //amplitude of the forcing
double SIGMA = 18.;
double R0 = 8.0;
int main (int argc, char * argv[]) {
if (argc > 1)
{
maxlevel = atoi(argv[1]); //maximum refinement level
MAXTIME = atoi(argv[2]); //maximum runtime
SIGMA = atof(argv[3]); //surface tension strength, which sets the Weber number
R0 = atof(argv[4]); //bubble diameter, if inserted
FORCED = atoi(argv[5]); //equals 1 if forced, 0 otherwise
amp_force = atof(argv[6]); //amplitude of the forcing
}
size (WIDTH);
foreach_dimension()
periodic (right);The acceleration is defined below. The level of refinement is maxlevel.
mu1 = 0.01 * sq(WIDTH/(2.0*pi))/(2.0);
mu2 = mu1/MUR;Use reference density for fluid, and define gas by density ratio.
rho1 = 1.0;
rho2 = rho1/RHOR;
f.sigma = SIGMA;
a = av;
#if TREE
N = 1 << (maxlevel-3);
#else
N = 1 << maxlevel;
#endifSet tolerance on flow divergence.
TOLERANCE = 1e-4;
run();
}Initial conditions
The initial condition is “ABC” flow. This is a laminar base flow that is easy to implement in both Basilisk and a spectral code. The init event is structured so that if there are no dump files from which to restart, a precursor simulation (i.e. without a bubble) is begun. If a dump file is present, its name will determine what kind of simulation is to be done. If it is called “restart”, then a bubble is injected into the already turbulent flow. If it is called “dump”, then it is merely a checkpoint from which to restart, and can be a precursor or not.
If there are no existing dump files from which to restart, initialize the flow (without a bubble).
if (!restore (file = "restart") && !restore(file="dump"))
{
double waveno = WIDTH/(2.0*pi);
double amp = WIDTH/(2.0*pi);
foreach() {
f[] = 1.0;
double xw = x/waveno;
double yw = y/waveno;
double zw = z/waveno;
u.x[] = amp*( cos(yw) + sin(zw) );
u.y[] = amp*( sin(xw) + cos(zw) );
u.z[] = amp*( cos(xw) + sin(yw) );
}
}If there is a dump file named “dump”, do not do anything.
If there is a dump file named “restart”, then insert the bubble.
else
{
fraction (f, sq(x-WIDTH/2.0) + sq(y-WIDTH/2.0) + sq(z-WIDTH/2.0) - sq(R0));
foreach(){
foreach_dimension(){
u.x[] = f[]*u.x[];
}
}
}
}Include accelerations:
event acceleration (i++) {Linear forcing term
We compute the average velocity and add the corresponding linear forcing term, if desired at runtime by FORCED variable.
if (FORCED)
{
coord ubar;
foreach_dimension() {
stats s = statsf(u.x);
ubar.x = s.sum/s.volume;
}
foreach_face()
av.x[] += f[]*amp_force*((u.x[] + u.x[-1])/2. - ubar.x);
}
}Outputs
We log the evolution of viscous dissipation, kinetic energy, and microscale Reynolds number. This is as in the isotropic.c example case.
event logfile (i++) {
coord ubar;
foreach_dimension() {
stats s = statsf(u.x);
ubar.x = s.sum/s.volume;
}
double ke = 0., vd = 0., vol = 0.;
foreach(reduction(+:ke) reduction(+:vd) reduction(+:vol)) {
vol += dv();
foreach_dimension() {
// mean fluctuating kinetic energy
ke += dv()*sq(u.x[] - ubar.x);
// viscous dissipation
vd += dv()*(sq(u.x[1] - u.x[-1]) +
sq(u.x[0,1] - u.x[0,-1]) +
sq(u.x[0,0,1] - u.x[0,0,-1]))/sq(2.*Delta);
}
}
ke /= 2.*vol;
vd *= mu1/vol;
static FILE * fd = fopen("stats.dat","a");//"w" before
if (i == 0) {
fprintf (ferr, "t dissipation energy Reynolds\n"); //screen
fprintf (fd, "t dissipation energy Reynolds\n");
}
fprintf (ferr, "%g %g %g %g\n",
t, vd, ke, 2./3.*ke/mu1*sqrt(15.*mu1/vd));
fprintf (fd, "%g %g %g %g\n",
t, vd, ke, 2./3.*ke/mu1*sqrt(15.*mu1/vd));
}We also want to count the drops and bubbles in the flow.
event countDropsBubble(i++)
{
scalar m1[]; //droplets
scalar m2[]; //bubbles
foreach(){
m1[] = f[] > 0.5; //i.e. set m true if f[] is close to unity (droplets)
m2[] = f[] < 0.5; //m true if f[] close to zero (bubbles)
}
int n1 = tag(m1);
int n2 = tag(m2);Having counted the bubbles, now we find their size. This example is similar to the jet atomization problem. We are interested in the volumes and positions of each droplet/bubble.
double v1[n1]; //droplet
coord b1[n1]; //droplet
double v2[n2]; //bubble
coord b2[n2]; //bubbleWe initialize:
for (int j=0; j<n1; j++)
{
v1[j] = b1[j].x = b1[j].y = b1[j].z = 0.0;
}
for (int j=0; j<n2; j++)
{
v2[j] = b2[j].x = b2[j].y = b2[j].z = 0.0;
}We proceed with calculation.
foreach_leaf() //droplets
{
if (m1[] > 0) {
int j = m1[] - 1;
v1[j] += dv()*f[]; //increment the volume of the droplet
coord p = {x,y,z};
foreach_dimension()
b1[j].x += dv()*f[]*p.x;
}
}
foreach_leaf() //bubbles
{
if (m2[] > 0) {
int j = m2[] - 1;
v2[j] += dv()*(1.0-f[]);
coord p = {x,y,z};
foreach_dimension()
b2[j].x += dv()*(1.0-f[])*p.x;
}
}Reduce for MPI.
#if _MPI
MPI_Allreduce (MPI_IN_PLACE, v1, n1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce (MPI_IN_PLACE, b1, 3*n1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce (MPI_IN_PLACE, v2, n2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce (MPI_IN_PLACE, b2, 3*n2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#endifOutput the volume and position of each droplet to file.
static FILE * fdrop = fopen("droplets.dat","w");
static FILE * fbubb = fopen("bubbles.dat","w");
for (int j=0; j<n1; j++)
{
fprintf (fdrop, "%d %g %d %g %g %g %g\n", i, t,
j, v1[j], b1[j].x/v1[j], b1[j].y/v1[j], b1[j].z/v1[j]);
}
for (int j=0; j<n2; j++)
{
fprintf (fbubb, "%d %g %d %g %g %g %g\n", i, t,
j, v2[j], b2[j].x/v2[j], b2[j].y/v2[j], b2[j].z/v2[j]);
}
}We generate a movie of the vortices. Commented lines are optional alternative outputs.
#define POPEN(name, mode) fopen (name ".ppm", mode)
event movie (t += 0.25)
{
view (fov = 44, camera = "iso", ty = .2,
width = 600, height = 600, bg = {1,1,1}, samples = 4);
clear();
squares ("u.y", linear = true, n = {0,1,0});
squares ("u.x", linear = true, n = {1,0,0});
squares ("u.z", linear = true, n = {0,0,1});
//scalar omega[];
//vorticity (u, omega);
//squares ("omega", linear = true, n = {0,1,0});
//scalar l2[];
//lambda2 (u, l2);
//isosurface ("l2", -1);
draw_vof ("f");
{
static FILE * fp = POPEN ("movie", "w");
save (fp = fp);
}
}Gas-liquid interface data is also outputted.
event getInterface(i++)
{
char fixname[100];
char fiyname[100];
char fizname[100];
char fiintname[100];
sprintf(fixname, "fix_%d.dat", pid() );
sprintf(fiyname, "fiy_%d.dat", pid() );
sprintf(fizname, "fiz_%d.dat", pid() );
sprintf(fiintname, "fint_%d.dat", pid() );
FILE * fix = fopen(fixname, "a");
FILE * fiy = fopen(fiyname, "a");
FILE * fiz = fopen(fizname, "a");
FILE * fint = fopen(fiintname, "a");
fprintf (fix, "%g ", t);
fprintf (fiy, "%g ", t);
fprintf (fiz, "%g ", t);
fprintf (fint, "%g\n", t);
output_facets(f,fint);
foreach()
{
if (f[] > 0.4 && f[] < 0.6)
{
fprintf (fix, "%g ", x);
fprintf (fiy, "%g ", y);
fprintf (fiz, "%g ", z);
}
}
fprintf (fix, "\n");
fprintf (fiy, "\n");
fprintf (fiz, "\n");
fprintf (fint, "\n");
fclose(fix);
fclose(fiy);
fclose(fiz);
fclose(fint);
return 0;
}We refine on the velocity field in the precursor, and on the VOF field in the full simulation.
double femax = 0.01;
adapt_wavelet ((scalar *){u}, (double[]){uemax,uemax,uemax}, maxlevel);
//adapt_wavelet ((scalar *){f, u}, (double[]){femax, uemax, uemax, uemax},maxlevel);
}
#endifWe output a full snapshot every time unit.
int j = 0;
event snapshot (t=1; t <=MAXTIME; i+=10)
{
//dump(NAME);
//
char namedump[80];
sprintf(namedump,"dump_%d",j);
dump(file = namedump);
j++;
if (j==2)
j=0;
}
event allsnapshots (t=60; t<=MAXTIME; t +=1)
{
char step[100];
p.nodump = false;
sprintf(step, "time=%g",t);
dump(step);
}End the simulation after 300 time units.
