sandbox/aberny/bursting.c

    Bursting bubble

    With this code, we want to simulate a symplified version of a bursting bubble in a water bath. We based this simulation on the work of Duchemin et al. published in 2002.

    #include "axi.h"
#include "navier-stokes/centered.h"
// #define mu(f) ((mu1*mu2)/(clamp(f,0,1)*(mu2-mu1)+mu1))
#include "two-phase.h"
#include "tension.h"
#include "tag.h"

    Definition of the geometrical parameters. The 2 radii correspond respectively to the radius of the bubble and to the radius of a smooth interface.

    #define LEVEL 14
#define L0 10
#define R1 1.

    The parameters L1, L2 and L3 are used to define the hole between the upper interface and the spherical droplet. Once the hole is done, we smooth the interface by adding a circle between the upper interface and the drop. This circle as to be continuous with the upper interface but also with the bubble interface. 2\times L2 defines the distance between the top of the bubble and the interface. L1 defines the radius of the smoothing circle, but also the distance between the center of this cylinder with the upper interface.

    #define L2 0.05
#define L3 0.15
#define L1 (-(sq(R1+2.*L2)-sq(R1)+sq(L3))/(4.*(L2+R1)))

#define Sigma 1.
    invalid storage class for function ‘_boundary9_homogeneous’
    invalid storage class for function ‘_boundary9’
    invalid storage class for function ‘_boundary10_homogeneous’ 
    uf.n[left] = 0.;
    invalid storage class for function ‘_boundary10’
    invalid storage class for function ‘_boundary11_homogeneous’ 
    uf.n[right] = 0.;
    invalid storage class for function ‘_boundary11’
    invalid storage class for function ‘_boundary12_homogeneous’ 
    uf.n[top] = 0.;
    invalid storage class for function ‘_boundary12’ 
    uf.n[bottom] = 0.;
    invalid storage class for function ‘metric’
    invalid storage class for function ‘metric_expr0’ 
    double uemax = 0.3;
    invalid storage class for function ‘defaults’
    invalid storage class for function ‘defaults_expr0’ 
    int main(int argc, char *argv[]) {
    invalid storage class for function ‘init’ 
      size (L0);
    invalid storage class for function ‘init_expr0’ 
      origin (-L0/2., 0.);
    invalid storage class for function ‘set_dtmax’ 
      init_grid (1 << (9));
    invalid storage class for function ‘set_dtmax_expr0’
    invalid storage class for function ‘stability’
    invalid storage class for function ‘stability_expr0’
    By default, the value of the Laplace number is 10 000. But, we can
    invalid storage class for function ‘vof’
    change that with an initial input argument.
    invalid storage class for function ‘vof_expr0’ 
      double La = 10000.;
    invalid storage class for function ‘tracer_advection’ 
      if (argc >= 2){
    invalid storage class for function ‘tracer_advection_expr0’ 
          La = atof (argv[1]);
    invalid storage class for function ‘tracer_diffusion’ 
        }
    invalid storage class for function ‘tracer_diffusion_expr0’
    invalid storage class for function ‘properties’
    invalid storage class for function ‘properties_expr0’
    The Ohnesorge number is defined such that: $Oh =
    invalid storage class for function ‘advection_term’
    invalid storage class for function ‘advection_term_expr0’ 
      double Oh = sqrt(1./La);
    invalid storage class for function ‘viscous_term’
    invalid storage class for function ‘viscous_term_expr0’
    invalid storage class for function ‘acceleration’
    We then define the properties of the fluids. In this case, the fluid
    invalid storage class for function ‘acceleration_expr0’
    1 correspond to water, and the fluid 2 correspond to air.
    invalid storage class for function ‘projection’
    invalid storage class for function ‘projection_expr0’ 
      rho1 = 1., mu1 = Oh;
    invalid storage class for function ‘end_timestep’ 
      rho2 = 1./998., mu2 = Oh/55;
    invalid storage class for function ‘end_timestep_expr0’
    invalid storage class for function ‘adapt’
    invalid storage class for function ‘adapt_expr0’
    The surface tension is defined with the Ohnesorge number
    invalid storage class for function ‘defaults_0’
    invalid storage class for function ‘defaults_0_expr0’ 
      f.sigma = Sigma;
    invalid storage class for function ‘stability_0’
    invalid storage class for function ‘stability_0_expr0’
    invalid storage class for function ‘vof_0’
    We print the characteristics of the fluid in the logfile.
    invalid storage class for function ‘vof_0_expr0’ 
      fprintf (stderr, "props %f %f %f %f %f\n", rho1, mu1, La, R1, Sigma);
    invalid storage class for function ‘defaults_1’
    invalid storage class for function ‘defaults_1_expr0’ 
      run();
    invalid storage class for function ‘properties_0’ 
    }
    invalid storage class for function ‘properties_0_expr0’
    invalid storage class for function ‘defaults_2’ 
    double geometry(double x, double y) {
    invalid storage class for function ‘defaults_2_expr0’
    invalid storage class for function ‘acceleration_0’
    invalid storage class for function ‘acceleration_0_expr0’
    The geometry need the definition of 2 circles (C*) and 3 lines
    invalid storage class for function ‘stability_1’
    (D*). Then we need to proceed correctly to the intersection and the
    invalid storage class for function ‘stability_1_expr0’
    union of those parameters.
    invalid storage class for function ‘acceleration_1’
    invalid storage class for function ‘acceleration_1_expr0’ 
      double C1 = sq(x + R1 + 2*L2) + sq(y) - sq(R1);
    invalid storage class for function ‘init_0’ 
      double C2 = sq(x - L1) + sq(y - L3) - sq(L1);
    invalid storage class for function ‘init_0_expr0’ 
      double D1 = - x - 1e-8;
    invalid storage class for function ‘adapt_0’ 
      double D2 = y - L3;
    invalid storage class for function ‘adapt_0_expr0’ 
      double D3 = - x + (2*L1);
    invalid storage class for function ‘droplets’
    invalid storage class for function ‘droplets_expr0’
    invalid storage class for function ‘droplets_expr1’
    Our geometry is:
    invalid storage class for function ‘interface’
    \displaystyle geom = C1 \cup (D3 \cap (D1 \cup D2))
    invalid storage class for function ‘interface_expr0’
    invalid storage class for function ‘logfile’ 
      double D1D2 = min(D1, D2);
    invalid storage class for function ‘logfile_expr0’ 
      double D1D2D3 = max(D1D2, D3);
    invalid storage class for function ‘finalState’ 
      double D1D2D3C1 = min(D1D2D3, C1);
    invalid storage class for function ‘finalState_expr0’ 
      double D1D2D3C1C2 = min(-D1D2D3C1, C2);
    invalid storage class for function ‘_set_boundary4’ 
      return -D1D2D3C1C2;
    invalid storage class for function ‘_set_boundary5’ 
    }
    invalid storage class for function ‘_set_boundary6’
    invalid storage class for function ‘_set_boundary7’
    invalid storage class for function ‘_set_boundary8’
    We initialise the geometry of the interface.
    invalid storage class for function ‘_set_boundary9’
    invalid storage class for function ‘_set_boundary10’
    invalid storage class for function ‘_set_boundary11’
    Initial interface of the simulation. The length of the rim is half of the size of the hole
    invalid storage class for function ‘_set_boundary12’
    ISO C forbids nested functions [-Wpedantic]
    ‘_init_solver’ defined but not used [-Wunused-function] 
    double intemax = 0.005;

event init (t = 0) {
  double iteration = 0;
  do {
    iteration++;
    fraction(f, geometry(x,y))
    expected ‘;’ before ‘}’ token 
      } while( adapt_wavelet({f,u}, (double []){intemax,uemax,uemax,uemax},
     maxlevel = LEVEL, 5).nf != 0 && iteration <= 10);

  static FILE * fp = fopen("initial.png", "w");
  output_ppm(f, fp, 400, min = 0, max = 1, box = {{-2.5,0},{0.5,3}});
  static FILE * fp2 = fopen("initial.gfs", "w");
  output_gfs(fp2);
    expected ‘while’ before ‘return’ 
    }

    We use an adaptive mesh. Our refinement criteria are based on the interface and the velocity.

    invalid storage class for function ‘adapt_0_expr0’
    invalid storage class for function ‘adapt_0’ 
    event adapt (i++) {
  adapt_wavelet ({f,u}, (double []){intemax,uemax,uemax,uemax},
		 maxlevel = LEVEL, 5);
}

    We output the shape of the interface at regular intervals.

    This gives the following image

    unset key
set size ratio -1
plot [-2:2]'interface' u 2:1 w l, 'interface' u (-$2):1 w l lt 1
    Evolution of the interface (script)

    Evolution of the interface (script)

    We output the position of the interface of the fluid allong the axis of symmetry at every timestep. We track this position with the height function. We also compute the corresponding velocity.

    variable ‘xDropPrev’ set but not used [-Wunused-but-set-variable] 
    static double tprev = -1, xprev = -1,xDropPrev = -1;
    unused variable ‘toto’ [-Wunused-variable]
    ‘toto’ defined but not used [-Wunused-variable] 
    static double compt = 0, toto = 0;
    invalid storage class for function ‘droplets_expr0’
    invalid storage class for function ‘droplets_expr1’
    invalid storage class for function ‘droplets’ 
    event droplets (i ++; t<=2) {
  scalar m[];
  foreach()
    m[] = f[] > 1.5e-1;
  int n = tag(m);

    We want to observe the formation of the first drop. We detect when the tag function return more than 1 phase of the internal fluid. We then output 15 steps in gfs file to observe the behaviour of the fluid at this moment.

      if (n>1 && t>0.15) {
    compt++;
    if (compt>=1 && compt<=15) {
      char rupture[80];
    format ‘%ld’ expects argument of type ‘long int’, but argument 3 has type ‘int’ [-Wformat=] 
          sprintf(rupture, "rupture-%07ld.gfs",i);
      FILE* fprup = fopen (rupture, "w");
      output_gfs(fprup);
    }
  }

    We tag the droplet, when there is some.

      double v[n];
  coord b[n];
  double dropVelo[n];

  for (int j = 0; j < n; j++)
    dropVelo[j] = v[j] = b[j].x = b[j].y = b[j].z = 0.;

  foreach_leaf()
    if (m[] > 0) {
      int j = m[] - 1;
      v[j] += dv()*f[];
      dropVelo[j] += dv()*f[]*u.x[];
      coord p = {x,y,z};
      foreach_dimension()
      b[j].x += dv()*f[]*p.x;
    }
    We compute the droplets velocity. We need at least 2 steps with a droplet to have a correct velocity. We add for that an extra variable, which is counting the number of steps with the presence of a drop. We also compute the volume of the first drop.
    expected declaration or statement at end of input 
      double xDrop, veloDrop, volume, radius;

  for (int j = 1; j < n; j++) {

    Computation of the drop velocity, and the drop position.

        xDrop = b[j].x/v[j];
    veloDrop = dropVelo[j]/v[j];

    Computation of the drop radius.

        volume = v[j];
    radius = sqrt(2*volume/M_PI);

    fprintf (fout, "%d %g %d %g %g %g %g %g %d\n", i, t,
       j, v[j], xDrop, veloDrop, volume, radius, n);
  }
  fflush (fout);

    We also compute the velocity of the top of the jet. We use for that the height function.

      vector h[];
  double xMax = -HUGE;;
  heights (f, h);
  foreach()
    if (y <= Delta && h.x[] != nodata && m[]==0) {
      double xi = x + height(h.x[])*Delta;
      if (xi > xMax)
        xMax = xi;
    }
  double velo = tprev >= 0 ? (xMax - xprev)/(t - tprev) : 0.;
  double deltaT = t - tprev;
  printf ("%d %g %g %g %g\n", i, t, deltaT, xMax, velo);
  fflush (stdout);

  if (xMax > 0. && xprev < 0.) {
    fprintf (stderr, "end %g %g\n", t, velo);
    char name[80];
    sprintf(name, "jetX-0.gfs");
    FILE* fp = fopen (name, "w");
    output_gfs (fp);
  }

    The stop condition of the simulation is based on 2 conditions. First, we need to overpass the axi x=0. Secondly, we need to observe at least 2 steps consecutive steps with a drops.

      if (xMax>=0. && n>1) {
    fprintf(stderr, "firstDrop %g %g %g %g\n", t, veloDrop, volume, radius);
    return 1;
  }

  tprev = t, xprev = xMax, xDropPrev = xDrop;
}

    This gives the following figure

    set xlabel 'Time'
set ylabel 'Velocity'
set grid
plot [0.01:]'out' u 1:4 w l t ''
    Evolution of the velocity of the centerline (script)

    Evolution of the velocity of the centerline (script)

    invalid storage class for function ‘interface_expr0’
    invalid storage class for function ‘interface’ 
    event interface(i+= 10) {
  static FILE * fp = popen ("ppm2mpeg > tracer.mpeg", "w");
  static FILE * fp2 = popen ("ppm2mpeg > grid.mpeg", "w");
  output_ppm(f, fp, 512, min=0, max=1, box = {{-2.5,0},{0.5,3}});
  scalar l[];
  foreach()
    l[] = level;
  output_ppm(l,fp2, min=5, max=LEVEL, box = {{-2.5,0},{0.5,3}});
}

    We can follow the evolution of the simulation with a video

    The bursting bubble

    #if 0
event gfsview (i += 10) {
  static FILE * fp = popen ("gfsview2D bursting.gfv", "w");
  output_gfs (fp);
}
#endif
    invalid storage class for function ‘logfile_expr0’
    invalid storage class for function ‘logfile’ 
    event logfile (i += 10) {
  fprintf (stderr, "%d %g %d %d %d\n", i, t, mgp.i, mgpf.i, mgu.i);
}

    We output the final state of the simulation. In the final state, we include the tag function, to see where basilisk identify the different phase.

    invalid storage class for function ‘finalState_expr0’
    invalid storage class for function ‘finalState’ 
    event finalState (t = end) {
  scalar m[];
  foreach()
    m[] = f[] > 1.5e-1;
    unused variable ‘n’ [-Wunused-variable] 
      int n = tag(m);
  output_gfs (file = "final.gfs");
  static FILE * fp = fopen("final.ppm", "w");
  output_ppm(f, fp, 512, min=0, max=1, box = {{-2.5,0},{0.5,3}});
  dump (file = "dump");
}

    We can control the refinement parameter by loocking at the evolution of the mesh during the simulation.

    Evolution of the mesh

    Evolution of the mesh