This is a standard file to reproduce the simulations presented in:

    Rivière, Mostert, Perrard and Deike, Sub-Hinze scale bubble production in turbulent bubble break up. Journal of Fluid Mechanics, 917:A40, 2021.

    Perrard, Rivière, Mostert and Deike, Bubble deformation by a turbulent flow. Journal of Fluid Mechanics, 920:A15, 2021.

    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.

    We use the \lambda_2 criterion and Basilisk View for visualisation of vortices.

    #include "lambda2.h"
    #include "view.h"

    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.0

    Define domain size.

    #define WIDTH 120.0

    We 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);
        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);
      N = 1 << maxlevel;

    Set tolerance on flow divergence.

      TOLERANCE = 1e-4;

    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.

    event init (i = 0) {

    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.

      else if(restore (file = "dump")){

    If there is a dump file named “restart”, then insert the bubble.

          fraction (f, sq(x-WIDTH/2.0) + sq(y-WIDTH/2.0) + sq(z-WIDTH/2.0) - sq(R0));
              	u.x[] = f[]*u.x[];

    Include accelerations:

    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;
          av.x[] += f[]*amp_force*((u.x[] + u.x[-1])/2. - ubar.x);


    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
        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];  //bubble

    We 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};
            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};
              b2[j].x += dv()*(1.0-f[])*p.x;

    Reduce for MPI.

    #if _MPI
      MPI_Allreduce (MPI_IN_PLACE, b1, 3*n1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
      MPI_Allreduce (MPI_IN_PLACE, b2, 3*n2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

    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\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);
      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);
          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");
      return 0;

    We refine on the velocity field in the precursor, and on the VOF field in the full simulation.

    #if TREE
    event adapt (i++) {
      double uemax = 0.2*normf(u.x).avg;
    unused variable ‘femax’ [-Wunused-variable]
      double femax = 0.01;
      adapt_wavelet ((scalar *){u}, (double[]){uemax,uemax,uemax}, maxlevel);
      //adapt_wavelet ((scalar *){f, u}, (double[]){femax, uemax, uemax, uemax},maxlevel);

    We output a full snapshot every time unit.

    int j = 0;
    event snapshot (t=1; t <=MAXTIME; i+=10)
      char namedump[80];
      dump(file = namedump);
    if (j==2)
    event allsnapshots (t=60; t<=MAXTIME; t +=1)
    char step[100];
    p.nodump = false;
    sprintf(step, "time=%g",t);

    End the simulation after 300 time units.

    event end (t=MAXTIME) {
      dump ("end");