src/test/bubble.h

    A generic compressible gas bubble in a liquid

    #include "compressible/two-phase.h"
#include "compressible/Mie-Gruneisen.h"
#include "tension.h"
#include "compressible/tension.h"
#include "rayleigh-plesset.h"

    Static refinement is used with a level of refinement varying from MINLEVEL far from the bubble to MAXLEVEL close to the bubble.

    int MINLEVEL = 4, MAXLEVEL;

    The liquid and gas densities as well as the initial gas pressure pg0 and the pressure “at infinity” pinf.

    double rhoL = 1., rhoG = 0.001;
double pg0 = 100., pinf = 100.;

    The final time and the bubble radius.

    double tend = 1.;
double R0 = 1.;

    Boundary conditions: imposed pressure and inflow “at infinity”.

    p[right]   = dirichlet (pinf);
q.n[right] = neumann (0.);

#if dimension > 1
p[top]     = dirichlet (pinf);
q.n[top]   = neumann (0.);
#endif

event init (i = 0)
{

    We compute the reference Rayleigh-Plesset and Keller-Miksis solutions.

      static bool first = true;
  if (first) {
    struct RPdata RPd = {
      .rhol  = rhoL,
      .pliq = pinf,
      .p0 = pg0,
      .sigma = f.sigma,
      .gamma = gamma2,
      .R0 = R0,
      .visc = mu1,
      .cson = sqrt(gamma1*(pinf + PI1)/rhoL)
    };
    
    FILE * fp = fopen("RP.dat", "w");
    Integrate_RP (fp, 0., tend, &RPd);
    fclose (fp);
  
    fp = fopen("RPinviscid.dat", "w");
    RPd.visc = 0;
    Integrate_RP (fp, 0., tend, &RPd);
    fclose (fp);

    first = false;
  }

    The static mesh refinement.

    #if TREE
  for (int l = MINLEVEL ; l <= MAXLEVEL; l++)
    refine (level < l && sqrt(sq(x) + sq(y)) < (2.5*R0 + 4.*sqrt(2.)*L0/(1 << (l - 1))));
#endif

    The initial volume fraction, densities, energies and pressure.

      fraction (f, sq(x) + sq(y) - sq(R0));
  
  foreach() {
    frho1[] = f[]*rhoL;
    frho2[] = (1. - f[])*rhoG;

    The initial liquid pressure field is approximated from the solution in the incompressible limit.

        double r = sqrt(sq(x) + sq(y));
    double pL = pinf*(1. - R0/r) + (pg0 - 2*f.sigma/R0)*R0/r;
	
    fE1[] = f[]*(pL + PI1*gamma1)/(gamma1 - 1.);
    fE2[] = (1. - f[])*pg0/(gamma2 - 1.);

    p[] = average_pressure (point);
  }
}

    Some statistics.

    event logfile (i++)
{
  scalar pgas[], keliq[];
  double volume = 0.;
  
  foreach (reduction(+:volume)) {
   
    double Ek = 0.;
    foreach_dimension()
      Ek += sq(q.x[]);
    Ek /= 2.*(frho1[] + frho2[]);

    The gas pressure is recovered from the energy.

        pgas[] = average_pressure (point)*(1. - f[]);
    keliq[] = Ek*f[];
    volume += dv()*(1. - f[]);
  }

  if (i == 0)
    fprintf (stderr, "t volume statsf(keliq).sum statsf(pgas).sum/volume\n");
  fprintf (stderr, "%10.6f %10.6f %10.8f %10.6f\n",
	   t, volume, statsf(keliq).sum, statsf(pgas).sum/volume);
}

    The bubble shape as a function of time.

    event profiles (t += tend/10; t <= tend) {
  output_facets (f);
}