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);
}