sandbox/farsoiya/viscous_droplet_breakup.c
Droplet breakage in turbulent flow
Role of viscous droplet breakup for the Farsoiya et. al. 2023.
//~ #define MIN_LEVEL 5
//~ #define LEVEL 10
int MAX_LEVEL = 8; //only for the compilation
//~ #define dR_refine (2.*L0/(1 << LEVEL))
#define F_ERR 1e-10
#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "navier-stokes/conserving.h"
#include "tension.h"
#include "tag.h"
#include "view.h"
We monitor performance statistics and control the maximum runtime.
#include "navier-stokes/perfs.h"
#include "maxruntime.h"
#include <sys/stat.h>
Include density and viscosity ratios.
#define RHOR 1.
//#define MUR 0.004
Defined domain size is
#define WIDTH 120.0
The code takes the level of refinement as optional command-line argument (as well as an optional maximum runtime).
double MAXTIME = 500; //end time default
int FORCED = 1; //equals 0 if it is not forced, 1 otherwise
double amp_force = 0.1; //amplitude of the forcing
double SIGMA = 500.0; // surface tension coeff
double R0 = 8.0; //Radius of the bubble
double SC1 = 1.0;
double mur = 1.;
double We =1.;
double visp = 1.;
double gravity = 0;
double tst = 1;
double conc_liq1 = 0;
//double conc_gas1 = 1.0;
int main (int argc, char * argv[]) {
if (argc > 1)
{
MAX_LEVEL = atoi(argv[1]);
MAXTIME = atoi(argv[2]);
SC1 = atof(argv[3]);
FORCED = atoi(argv[4]); //equals 1 if forced, 0 otherwise
amp_force = atof(argv[5]); //amplitude of the forcing
//visp for Re = 38, visp=2, Re = 55, visp = 1 and Re = 77, visp = 0.5
visp = atof(argv[6]);
We = atof(argv[7]);
mur = atof(argv[8]);
tst = atof(argv[9]);
}
size (WIDTH);
foreach_dimension()
periodic (right);
Use reference density for fluid, and define gas by density ratio.
double epsilon = 10.;
rho1 = 1.0;
rho2 = rho1/RHOR;
SIGMA = 2.*rho1*pow(epsilon,2./3)*pow(2.*R0,5./3)/We;
mu1 = visp* 0.01 * sq(WIDTH/(2.0*pi))/(2.0)/2.0;
mu2 = mu1*mur;
And the surface tension.
f.sigma = SIGMA;
#if TREE
N = 1 << (MAX_LEVEL-2);
#else
N = 1 << MAX_LEVEL;
#endif
Set tolerance on flow divergence.
TOLERANCE = 1e-4;
NITERMAX = 20;
struct stat st = {0};
char name[80];
sprintf (name, "dat");
char newname[80];
sprintf (name, "dat");
if (stat(name, &st) == -1)
{
mkdir(name, 0700);
}
sprintf (name, "Images");
if (stat(name, &st) == -1)
{
mkdir(name, 0700);
}
Preserve old stats restoring the simulation with a timestamp suffixed.
sprintf (name, "stats.dat");
if (stat(name, &st) == 0)
{
sprintf (newname, "stats-%ld.dat",time(0));
rename(name, newname);
}
sprintf (name, "perfs");
if (stat(name, &st) == 0)
{
sprintf (newname, "perfs-%ld",time(0));
rename(name, newname);
}
sprintf (name, "gas_transfer_stats.dat");
if (stat(name, &st) == 0)
{
sprintf (newname, "gas_transfer_stats-%ld.dat",time(0));
rename(name, newname);
}
sprintf (name, "bubbles.dat");
if (stat(name, &st) == 0)
{
sprintf (newname, "bubbles-%ld.dat",time(0));
rename(name, newname);
}
sprintf (name, "droplets.dat");
if (stat(name, &st) == 0)
{
sprintf (newname, "droplets-%ld.dat",time(0));
rename(name, newname);
}
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.
#include <dirent.h>
event init (i = 0) {
if (!restore (file = "restart") && !restore(file="dump"))
{
double waveno = WIDTH/(2.0*pi);
double amp = WIDTH/(2.0*pi);
foreach() {
f[] = 1.;
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) );
}
}
else if(restore (file = "dump")){
// do nothing
}
else
{
refine(sq(2.0*R0)-sq(x-WIDTH/2.0) - sq(y-WIDTH/2.0) - sq(z-WIDTH/2.0)>0 && level < MAX_LEVEL );
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++) {
face vector av = a;
Linear forcing term
We compute the average velocity and add the corresponding linear forcing term.
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);
}
boundary((scalar *){av});
}
Outputs
We log the evolution of viscous dissipation, kinetic energy, and microscale Reynolds number.
event logfile (i++) {
//~ printf("\n%d\t%g\t%g",i,t,dt);
//~ fflush(stdout);
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","w");
if (i == 0) {
//~ fprintf (ferr, "t dissipation energy Reynolds\n");
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));
fflush (fd);
}
We use adaptivity.
#if TREE
event adapt (i++) {
double uemax = 0.2*normf(u.x).avg;
// double tr1emax = 0.2*normf(c1).avg;
double femax = 1e-2;
adapt_wavelet ((scalar *){f, u},
(double[]){femax, uemax, uemax, uemax},MAX_LEVEL);
}
#endif
We output a full snapshot every time unit.
double endtimemark = 300;
event snapshot (t=0; t <=MAXTIME; t+=1)
{
char namedump[80];
sprintf(namedump,"./dat/dump_%g",t);
dump(file = "dump");
}
We also want to count the drops and bubbles in the flow.
int endflag = 0;
event countDropsBubble(t+=0.01)
{
scalar m1[]; //droplets
scalar m2[]; //bubbles
foreach(){
m1[] = f[] > 0.05; //i.e. set m true if f[] is close to unity (droplets)
m2[] = f[] < 0.99; //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]; //bubble
double ar1[n1];
double ar2[n2];
We initialize:
for (int j=0; j<n1; j++)
{
v1[j] = ar1[j] = b1[j].x = b1[j].y = b1[j].z = 0.0;
}
for (int j=0; j<n2; j++)
{
v2[j] = ar2[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
if (f[] > 1e-6 && f[] < 1. - 1e-6) {
coord n = interface_normal (point, f), p;
double alpha = plane_alpha (f[], n);
ar1[j] += pow(Delta, dimension - 1)*plane_area_center (n, alpha, &p);
}
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[]);
if (f[] > 1e-6 && f[] < 1. - 1e-6) {
coord n = interface_normal (point, f), p;
double alpha = plane_alpha (f[], n);
ar2[j] += pow(Delta, dimension - 1)*plane_area_center (n, alpha, &p);
}
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, ar1, 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, ar2, n2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce (MPI_IN_PLACE, b2, 3*n2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#endif
Output 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 %g\n", i, t,
j, v1[j], b1[j].x/v1[j], b1[j].y/v1[j],b1[j].z/v1[j],ar1[j]);
}
for (int j=0; j<n2; j++)
{
fprintf (fbubb, "%d %g %d %g %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],ar2[j]);
}
double temp;
// Simulation end criteria
if (n2 > 1){ //if more than one droplet exists
//sort droplet volume ascending
for (int k=0; k<n2; k++){
for (int j=0; j<(n2-k); j++){
if(v2[j] > v2[j+1]){
temp = *(v2 + j);
*(v2 + j) = *(v2 + j + 1);
*(v2 + j + 1) = temp;
}
}
}
double delta = L0/(pow(2,MAX_LEVEL));
// Thanks to Alienor who pointed out that looking the size of smallest droplet can result in judging ficticious breakup
// //hence we will now look for second largest droplet which should be more than 4 delta cube
if (endflag == 0 && v2[n2-2] > 4.*pow(delta,3)){ //if first breakup and smallest droplet has volume more than 4 times cell volume
endflag = 1; //declare fist breakup
endtimemark = t; //not the time for first breakup
dump(file = "breakup"); //save the breakup dump file
}
else if (endflag == 1){
if (t > (endtimemark + 30)){ // end simulation after 10tc
printf("End time after breakup");
dump(file = "end");
exit(0);
}
}
}
}
scalar omg[];
void vorticity3D (vector u, scalar omg) {
foreach() {
omg[] = 0.;
foreach_dimension()
omg[] += sq((u.y[1] - u.y[-1] -
u.x[0,1] + u.x[0,-1])/(2*Delta));
omg[] = sqrt(omg[]);
}
boundary ((scalar*){omg});
}
event movie (t+=0.1)
{
char namedump[200];
sprintf(namedump,"Images/temp.ppm");
vorticity3D(u, omg);
view (fov = 44, camera = "iso", ty = .2,
width = 1200, height = 1200, samples = 4);
clear();
//box();
draw_vof("f");
squares ("omg", spread = 5, linear = true, map = jet, alpha = 1, n = {1, 0 , 0});
squares ("omg", spread = 5, linear = true, map = jet, alpha = 1, n = {0, 1 , 0});
squares ("omg", spread = 5, linear = true, map = jet, alpha = 1, n = {0, 0 , 1});
save (namedump);
sprintf(namedump,"convert Images/temp.ppm Images/snap-%g.png",t*100);
system(namedump);
}
End the simulation.
References
[farsoiya2023role] |
Palas Kumar Farsoiya, Zehua Liu, Andreas Daiss, Rodney O Fox, and Luc Deike. Role of viscosity in turbulent drop break-up. Journal of Fluid Mechanics, 972:A11, 2023. |