/** In this simulation we observe the bursting of a bubble on a water-air interface. There is several option in this simulation we can activate to study different things arround the bursting Bubble process First, this simulation is long enough to allow us to observe all the jet drops coming from a bursting bubble event ## Available option * ACOUSTIC * Gravity * BRUIT * AdaptOn * ROMEO * DIV * IntTracking * Bview * MEASURE * ENERGY ### ACOUSTIC This option will increase the domain size and it will add a pressure measurement far away from the bubble. This was a request from an acoustician friends, in order to compare the pressure he measured experimentaly, and the one observe by Basilisk, under the incompressible assumption */ /** ### Gravity This option will set the gravity forces to true. They are defined by the Bond number (Bo), with $Bo = \frac{\rho g R^2}{\gamma}$, with $\rho$ being the liquid density, g the acceleration of gravity, R the radius of the bubble, and $\gamma$ the surface tension. ### BRUIT This option (BRUIT is noise in french) will add an initial noise in the simulation. This noise is a random velocity in the liquid field. This velocity field is centerd on 0 and has an amplitude based on the velocity of the first drop. ### AdaptOn When set to true, this option will allow the use of the adaptiv mesh of Basilisk ### ROMEO Juliette Pierre asked me if it was possible to track some bubble taht could be trapped below the simulation. This option allow the simulation to track those kind of bubble in the liquid phase ### DIV This option will allow us to observe the divergence field in the simulation ### IntTracking This option will allow the simulation to output the interface of the simulation for the first steps of the simulation ### Bview This option allow the use of bview to output images and video ### MEASURE This option allow the use of the measure module. This module is a post-processing tool to measure the position of the interface allong the axis of symmetry, but also the size, the position and the velocity of all the droplets coming from the simulation ### ENERGY This option allow the use of the energy module (which is not working). This module is a post-processing tool to measure the energy in the simulation ## Publication linked to this simulation The results on all the drops were published [here](https://journals.aps.org/prfluids/abstract/10.1103/PhysRevFluids.5.033605) in Physical Review of Fluids. A detail study on the statistic is available [here](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GL092919) in Geophysical Research Letters. ## Basilisk option We want to perform an axi-symetric simulation, with the centered Navier-Stokes solver, on a 2 fluid simulations. */ #include #include "axi.h" #include "navier-stokes/centered.h" #include "two-phase.h" #include "tension.h" /** tag.h is a module use to tag the different phase of bubble or liquid */ #include "tag.h" /** distance.h is a module use to get the height function */ #include "distance.h" /** We add profiling option to track the perfs of the simulation. Those module are not mandatory for the simulation to work properly, but they provide usefull information*/ #include "navier-stokes/perfs.h" #include "profiling.h" /** Last option: the gravity*/ #define Gravity 1 #if Gravity #include "reduced.h" #endif /** ## Physical Parameters There are 2 parameters for this simulation. The first one is the Laplace number. The second one is the Bond number. The default value used for the (James) Bond number is 0.07. It's used to compute the initial shape of the bubble.*/ double La = 10000; double Bond = 0.07; /** The surface tension $\sigma$ is defined by the parameter Sigma. We define it equal to 1.*/ #define Sigma 1. /** The simulation should end at t = END. Iend is used in case of noise to mark the end of an output*/ #define END 1.2 #define Iend 8500 /** ## Options Those option have been describe at the begining of this file */ #define ACOUSTIC 0 #define BRUIT 0 double amplitude = 0.01; // the noise amplitude /** Alternatively, we can add some noise on the interface (only if there is no noise in the liquid phase). This option is disabled and may not work anymore*/ // #if !BRUIT // #define InterfaceNoise 1 // #endif #define AdaptOn 1 #define ROMEO 0 #define DIV 0 #define IntTracking 0 /** ## Numerical Parameters LEVEL stand for the maximum refinement of the grid L0 is for the size of the simulation domain. We can see here the impact of the ACOUSTIC option*/ int LEVEL = 10; #if ACOUSTIC #define L0 25. #else #define L0 10. #endif /** The MinLevel option define the minimum level of refinement for the adaptive mesh */ #define MinLevel 4 /** If we want to output a simulation when there is creation of new liquid phase, we should define DumpSimulation to 1. Otherwise, it should be 0.*/ #define DumpSimulation 1 /** We add the boundary condition to allow the liquid drops to escape the simulation domain once the reach to top of the domain */ u.n[right] = neumann(0.); p[right] = dirichlet(0.); /** intemax, uemax1 and uemax2 are numerical parameters that controll the condition when the mesh will increase or decrease its refinement level*/ double intemax = 0.005; double uemax1 = 0.2; double uemax2 = 0.15; /** ## The simulation setup*/ int main(int argc, char const *argv[]) { size (L0); // size of the domain #if ACOUSTIC origin (-20, 0.); init_grid (1 << (11)); // initial grid in ACOUSTIC: $2^11$ cells in each direction #else origin (-L0/2., 0.); init_grid (1 << (9)); // initial grid: $2^9$ cells in each direction #endif /** By default, the value of the Laplace number is 10 000. But we can change that with an initial input argument. */ if (argc >= 2){ La = atof (argv[1]); } /** By default, the value of the Bond number is 0.07. But we can change that with an initial input argument. */ if (argc >=3){ Bond = atof (argv[2]); } /** By default, the value of the maximum level refinement is 10. But we can change that with an initial input argument. */ if (argc >=4){ LEVEL = atoi (argv[3]); } /** By default the adaptiv criterion on the velocity is at 0.15, but we can change that with an initial input argument.*/ if (argc >=5){ uemax2 = atof (argv[4]); } /** By default the initial noise amplitude is at 0.01, but we can change that with an initial input argument.*/ if (argc >=6){ amplitude = atof (argv[5]); } /** We divide by $La$ to obtain the Oh number, so $La$ must be positive. Concerning the $Bo$ number, we prefer to avoid numerical problems, so we forced him to be positive. */ assert(La > 0.); assert(Bond > 0.); double Oh = sqrt(1./La); /** We then define the properties of the fluids. In this case, the fluid 1 corresponds to water, and the fluid 2 corresponds to air. */ rho1 = 1., mu1 = Oh; rho2 = 1./998., mu2 = Oh/55.; /** We have defined the density, the viscosity and the shape with dimensionless number. Therefore, we have to fix the surface tension to 1. */ f.sigma = Sigma; /** Last, if the gravity is defined, we set the gravity action to the Bond number value. */ #if Gravity G.x = -Bond; #endif /** In case of simulation with noise, we initialize the random seed to make sure that all the simulation with noise will be different*/ srand(time(NULL)); /** We print the characteristics of the fluid in the log file. */ fprintf (stderr, "props %f %f %f %f %f %i\n \n", rho1, mu1, La, Bond, Sigma, LEVEL); run(); } /** ## Some functions We compute the velocity of Ganan-Calvo (PRL 2017) (with the correction from Deike et al PRF 2018). This will be the velocity of reference for the noise.*/ double refVelo(double La, double Bo){ /* */ return sqrt(La)*19*pow(( (1+2.2*Bo) * La * (pow(500,-1./2.)-pow(La,-1./2.)) ),-3./4.); } /** If we want to add noise in the simulation, there is a function for that.*/ double velocityNoise(double amplitude) { return( (2*(rand()/(double)RAND_MAX)-1)*amplitude ); } /** ## Initialisation We initialize the simulation with the bubble shape. We include the library findBond in the code. It will solve the bubble shape equation and return a liste of coordinate. Then, we use a second function, geometry that will take the hollow circle and the bubble shape, to proceed to the union of the 2 distance fields.*/ #include "findBond.h" #include "view.h" //include here because we didn't need to perform output before event init (t = 0) { if (!restore (file = "restart") ){ /** Here, dataShape will contain the coordinate of the water/air interface. */ coord* dataShape; Circle* hollow = NULL; Circle fillet; hollow = &fillet; Circle* cap = NULL; Circle topCap; cap = &topCap; dataShape = shapeBond(Bond, hollow, cap); // fprintf(stderr, "x= %f y= %f r= %f\n",fillet.x, fillet.y, fillet.r ); // for debug purpose /** We are using the distance function, applied on the list of coordinate dataShape. This will generate a distance function. */ scalar d[]; distance (d, dataShape); #if BRUIT double Ca = refVelo(La, Bond); fprintf(stderr, "ref Velocity = %f, velo ref adim = %f \n", Ca, Ca*sqrt(1/La)); fprintf(stderr, "amplitude = %f\n", amplitude); foreach() foreach_dimension(){ if (d[]<0.){ u.x[] = velocityNoise(amplitude)*Ca; } } /** Since we changed the velocity field, we need to redifine the boundary condition*/ boundary ({u.x}); boundary ({u.y}); dump("noiseInit"); FILE * fp2 = fopen("noiseField", "w"); // output_field({u.x, u.y}, fp2, n = 1<<9); #endif #if AdaptOn while (adapt_wavelet ({d}, (double[]){intemax}, LEVEL,9).nf); #endif /** The distance function is defined at the center of each cell, we have to calculate the value of this function at each vertex. */ #if AdaptOn do { #endif vertex scalar phi[]; foreach_vertex(){ double a = (d[] + d[-1] + d[0,-1] + d[-1,-1])/4.; double b = -sq (fillet.r) + sq (x - fillet.y) + sq (y - fillet.x); double c = min(a,b); phi[] = -c; } /** We can now initialize the volume fraction of the domain. */ fractions (phi, f); #if AdaptOn } while(adapt_wavelet({f}, (double[]){intemax}, LEVEL, 9).nf); #endif /** We output the initial situation */ static FILE * fp = fopen("initial.png", "w"); output_ppm(f, fp, 400, min = 0, max = 1, box = {{-2.5,0},{2.5,3}}); // a png file, without bview dump("initial"); clear(); view (fov = 4.69417, quat = {0,0,-0.707108,0.707106}, tx = 0., ty = 0.0855609, bg = {1,1,1}, width = 1280, height = 720, samples = 4); draw_vof("f", filled = 1); mirror({0,1,0}, 0) draw_vof("f", filled = 1); save("initial.ppm"); #if BRUIT clear(); view (fov = 12, quat = {0,0,-0.707108,0.707106}, tx = 0., ty = 0.0855609, bg = {1,1,1}, width = 1280, height = 720, samples = 4); draw_vof("f"); squares("u.x", min = -0.1, max = 0.1, map = cool_warm); mirror({0,1,0}, 0){ draw_vof("f"); squares("u.y", min = -0.1, max = 0.1, map = cool_warm);} save("initialNoise.png"); #endif //BRUIT } else { // restart Case FILE * fp = fopen ("initialRestore.png", "w"); output_ppm(cm, fp, 400); } } /** ## Outputing noise in the simulation When we have noise, we output in mp4 file the first Iend steps. We show the velocity field and the vorticity field*/ #if (BRUIT || InterfaceNoise) event noiseMovie (i++; i < Iend) { output_field(); { static FILE * fpNoise = popen ("ppm2mp4 noise.mp4", "w"); clear(); view (fov = 12.699, quat = {0,0,0.707107,-0.707107}, ty = 0.00253494, bg = {1,1,1}, width = 1920, height = 1080, samples = 4); draw_vof("f"); squares("u.x", map = cool_warm, min = -0.5, max = 0.5); mirror({0,1,0}, 0) { draw_vof("f"); squares("u.y",map = cool_warm, min = -0.5, max = 0.5); } save(fp = fpNoise); } { scalar omega[]; vorticity (u, omega); static FILE * fpNoise = popen ("ppm2mp4 noiseOmega.mp4", "w"); clear(); view (fov = 12.699, quat = {0,0,0.707107,-0.707107}, ty = 0.00253494, bg = {1,1,1}, width = 1920, height = 1080, samples = 4); draw_vof("f"); squares("omega", min = -5, max = 5, linear = true); mirror({0,1,0}, 0) { draw_vof("f"); squares("omega", min = -5, max = 5, linear = true); } save(fp = fpNoise); } } #endif /** ## Adaptation At each step, we will adapt the mesh. We make sure that the mesh size will be at least $1/2^MinLevel$ */ #if AdaptOn event adapt (i++) { if (i<20) adapt_wavelet ({f,u}, (double []){intemax, uemax1, uemax1, uemax1}, maxlevel = LEVEL, minlevel = MinLevel); else adapt_wavelet ({f,u}, (double []){intemax, uemax2, uemax2, uemax2}, maxlevel = LEVEL, minlevel = MinLevel); } #endif /** ## General output We output a few files. One will display the interface to observe the evolution of the simulation. The other are video files, output when Bview is define to 1*/ int jj = 0; event outputInterface (i += 500) { // event outputFilm (t += 0.001) { // it's better to have an output based on the steps rather than the time /** We can output the interface in a txt file */ // { // static FILE * fp = fopen ("interface", "w"); // jj++; // char interfaceFile[80]; // sprintf(interfaceFile, "interface-%03d",jj); // static FILE * fp = fopen (interfaceFile, "w"); // FILE * fp = fopen (interfaceFile, "w"); // output_facets (f, fp); // fprintf(fp, "\n"); // } /** If we want to use bview to observe the evolution of the simulation, we define Bview to 1.*/ #define Bview 1 #if Bview scalar omega[]; vorticity (u, omega); /** We display the evolution of the interface. The fluid will be in black, the gas in white*/ { static FILE * fp1 = popen ("ppm2mp4 interface.mp4", "w"); clear(); view (fov = 12.699, quat = {0,0,0.707107,-0.707107}, ty = 0.00253494, bg = {1,1,1}, width = 1920, height = 1080, samples = 4); draw_vof("f", filled = 1); mirror({0,1,0}, 0) draw_vof("f", filled = 1); save(fp = fp1); } /* We show the vorticity field*/ { static FILE * fp2 = popen ("ppm2mp4 omega.mp4", "w"); // scalar omega[]; // vorticity (u, omega); // scalar l[]; // foreach() // l[] = level; clear(); view (fov = 12.699, quat = {0,0,0.707107,-0.707107}, ty = 0.00253494, bg = {1,1,1}, width = 1920, height = 1080, samples = 4); draw_vof("f"); squares("omega", map = jet, min = -50, max = 50, linear = true); // squares("omega", map = bwr, min = -50, max = 50, linear = true); mirror({0,1,0}, 0) { draw_vof("f"); squares("-omega", map = jet, min = -50, max = 50, linear = true); // squares("-omega", map = bwr, min = -50, max = 50, linear = true); } save(fp = fp2); } // The colormap bwr is not yet in the general version of basilisk. /** The velocity field.*/ // { // static FILE * fp2 = popen ("ppm2mp4 ux.mp4", "w"); // scalar omega[]; // vorticity (u, omega); // scalar l[]; // foreach() // l[] = level; // clear(); // view (fov = 12.699, quat = {0,0,0.707107,-0.707107}, ty = 0.00253494, bg = {1,1,1}, width = 1920, height = 1080, samples = 4); // draw_vof("f"); // squares("u.x", map = jet, min = -30, max = 30); // mirror({0,1,0}, 0) { // draw_vof("f"); // squares("u.x", map = jet, min = -30, max = 30); // } // save(fp = fp2); // } /** The vorticity field, with the cool_warm colormap. The liquid is in black, the vorticity is only showed in the gas*/ // { static FILE * fp2 = popen // ("ppm2mp4 omegaFilled.mp4", "w"); // scalar l[]; // foreach() // l[] = level; // clear(); // view (fov = 12.699, quat = {0,0,0.707107,-0.707107}, ty = 0.00253494, bg = {1,1,1}, width = 1920, height = 1080, samples = 4); // draw_vof("f", filled = 1); // squares("omega", map = cool_warm, min = -30, max = 30); // mirror({0,1,0}, 0) { // draw_vof("f", filled = 1); // squares("omega", map = cool_warm, min = -30, max = 30); // } // save(fp = fp2); // } /** The "level" field. The color is based on the level of resolution*/ { static FILE * fp3 = popen ("ppm2mp4 level.mp4", "w"); scalar l[]; foreach() l[] = level; clear(); view (fov = 12.699, quat = {0,0,0.707107,-0.707107}, ty = 0.00253494, bg = {1,1,1}, width = 1920, height = 1080, samples = 4); draw_vof("f", filled= 1); squares("l", min = 1, max = LEVEL); mirror({0,1,0}, 0) { draw_vof("f"); squares("l", min = 1, max = LEVEL); } save(fp = fp3); } /** The "level" field again, but with a zoom.*/ { static FILE * fplz = popen ("ppm2mp4 levelZoom.mp4", "w"); scalar l[]; foreach() l[] = level; clear(); view (fov = 4.42681, quat = {0,0,0,1}, tx = 0.0579169, bg = {1,1,1}, width = 1920, height = 1080, samples = 4); draw_vof("f"); squares("l", min = 1, max = LEVEL); cells(); mirror({0,1,0}, 0) { draw_vof("f"); squares("l", min = 1, max = LEVEL); cells(); } save(fp = fplz); } /** The vorticity field, with a zoom this time*/ { static FILE * fpz = popen("ppm2mp4 omegaZoom.mp4", "w"); clear(); view (fov = 6.69986, quat = {0,0,0,1}, tx = -0.05, ty = 0, bg = {1,1,1}, width = 1920, height = 1080, samples = 4); squares("omega", min = -200, max = 200, map = jet, linear = true); // squares("omega", min = -200, max = 200, map = bwr, linear = true); draw_vof("f"); mirror(n = {0,1,0}, alpha = 0){ squares("-omega", min = -200, max = 200, map = jet, linear = true); // squares("omega", min = -200, max = 200, map = bwr, linear = true); draw_vof("f"); } save(fp = fpz); } #if ROMEO /** The liquid/gas interface, below the bubble to observe a potential tiny bubble.*/ { static FILE * romeo = popen("ppm2mp4 zoom.mp4", "w"); scalar l[]; foreach() l[] = level; clear; view (fov = 0.647604, quat = {0,0,-0.707107,0.707107}, tx = -0.000410862, ty = 0.160506, bg = {1,1,1}, width = 1476, height = 870, samples = 4); draw_vof("f", filled = 1); squares("l", min = 1, max = LEVEL); mirror({0,1,0}, 0) { draw_vof("f"); squares("l", min = 1, max = LEVEL); } save(fp = romeo); } #endif #endif { char inter[80]; sprintf(inter, "interface-t-%g", t); FILE * fp = fopen(inter, "w"); output_facets(f, fp); fclose(fp); } dump ("dump"); // for a potential restart } /** This small function is used to detect a time step. This allow to have an output at a wanted time steps, without forcing the solver to reach it precisely*/ int timeDetect (double t, double tPrevious, double target) { if (tPrevious == -1) return -1; if ((tPrevious - target <0.) && (t - target > 0.)) return 1; else return -1; } double tAvant = -1; /** We output a movie that travel allong the jet The simulation evolution looks like:


The bubble bursting

We output a file at regular time step until the end of the simulation */ event simuInfo (i += 500; t<= END) { // event simuInfo (i += 100; t<= END) { // event simuInfo(i+=2500; i<=Iend){ /** It could be usefull to observe the vorticity of the flow. That's why we compute it.*/ vector h[]; heights (f, h); scalar omega[]; vorticity (u, omega); char dumpFile[80]; sprintf (dumpFile, "simulation-%07d", i); dump (file = dumpFile); } /** A general evolution pictures can be: ~~~gnuplot Evolution of the interface unset key set size ratio -1 plot [-2:2]'interface' every :10 u 2:1 w l, 'interface' u (-$2):1 w l lt 1 ~~~ ## Measure measure.h includes all we need for the droplet measurement. It provides a new structure to store the droplet data, it's initialisation and a comparison function for the qsort function (comming from stdlib.h). This library also contains an event for the physical measure. If we just want the video output, we just have to set MEASURE to 0*/ #define MEASURE 1 #if MEASURE #include "measure.h" #endif #define ENERGY 0 #if ENERGY #include "energy.h" #endif #if ACOUSTIC #include "acoustic.h" #endif event logfile (i++) { fprintf(stderr, "%g %i %i %i\n", t, i, mgpf.i, mgpf.nrelax); } /** ## End As verification, we output the final state. */ // event finalState (i = Iend) { // dump(file = "final"); // } event finalState (t = end) { dump(file = "final"); } /** ## Post-Processing We are interested in the velocity of the interface along the axis of symmetry. ~~~gnuplot Velocity along the axis of symmetry set ratio -1 set xlabel 'Time' set ylabel 'Velocity' set yrange [0:0.4] set grid plot 'out' u 2:6 ~~~ We can also follow the evolution of the position of the different elements (drops and interface). ~~~gnuplot Position of the elements along the axis of symmetry set xlabel 'Time' set ylabel 'Position on the axis' plot 'out' u 2:4 w l t "jet position", 'drop' u 2:4 every ::0::0 t "first drop", 'drop' u 2:4 every ::1::1 t "second drop", 'drop' u 2:4 every ::2::2 t "third drop" ~~~ Max number of cells: ~~~bash max=`awk '{max_val=($6