# A bubble shrinking due to thermal effects

This reproduces Figure 3.a of Saade et al, 2023 where more detailed explanations can be found (section “4.1 Epstein-Plesset like problem for temperature”).

set xlabel 't/tau'
set ylabel 'R/R_0'
plot "log" u 1:2 w l t 'Spherical', "../shrinking-axi/log" u 1:2 w l t 'Axisymmetric'

Animation of the pressure and temperature fields.

We run both the axisymmetric and the spherically-symmetric versions.

#if AXIS
# include "axi.h"
# include "compressible/thermal.h"
# include "view.h"
# define LEVEL 8
#else
# include "spherisym.h"
# include "compressible/thermal.h"
# define LEVEL 9
#endif
#include "compressible/NASG.h"

The initial density of the gas is chosen such that the initial temperature inside the bubble is twice the far-field temperature T_\infty.

double rhoL = 1., rhoR = 0.02011771644;
double p0L = 1.;
double p0 = 1.;
double tend = 1.;
double R0 = 1.;
double tau;

The problem is rendered dimensionless using the ambient pressure, the liquid density, the far-field temperature and the bubble initial radius. The values employed for this simulation are respectively listed.

double pdim = 5e6;
double rhodim = 975.91;
double Tdim = 350;
double Rdim = 1e-4;

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

#if AXIS
p[left]    = dirichlet(p0L);
q.n[left]  = neumann(0.);

p[top]    = dirichlet(p0L);
q.n[top]  = neumann(0.);
#endif

Although the thermal solver is implicit and unconditionally stable, a diffusive CFL condition is employed for better accuracy.

event stability (i++) {
dtmax = rhoR*cp2*sq(L0/pow(2,LEVEL))/kappa2/2.;
}

int main()
{
L0 = 8.;
#if AXIS
X0 = -L0/2.;
#endif

Liquid water parameters in the Noble-Abel Stiffened Gas (NASG) equation of state.

  gamma1 = 1.187;
PI1 = 7028e5/pdim;
b1 = 6.61e-4*rhodim;
q1 = -1177788*rhodim/pdim;

Specific heats and thermal conductivity of the fluids.

  cv1 = 3610*rhodim*Tdim/pdim; cv2 = 729.1*rhodim*Tdim/pdim;
cp1 = 4285*rhodim*Tdim/pdim; cp2 = 1063*rhodim*Tdim/pdim;

kappa1 = 0.6705/(Rdim/Tdim*sqrt(cube(pdim)/rhodim));
kappa2 = 0.03153/(Rdim/Tdim*sqrt(cube(pdim)/rhodim));

mu1 = 3.7e-4/(Rdim*sqrt(rhodim*pdim));
mu2 = mu1*1e-2;

The diffusive time scale \tau based on the gas properties.

  tau = rhoR*cp2/kappa2;

tend *= tau;

#if TREE
N = 1 << 4;
#else
N = 1 << LEVEL;
#endif

run();
}

event init (t = 0)
{
if (!restore (file = "restart")) {

The static mesh refinement.

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

Initialization of a bubble with initial radius R0.

    fraction (f, - (sq(R0) - sq(x) - sq(y)));

foreach() {
frho1[]  = f[]*rhoL;
frho2[]  = (1. - f[])*rhoR;

double pL = p0L;
p[] = pL*f[] + p0*(1. - f[]);
T[] = average_temperature (point, p[]);

fE1[] = (pL + gamma1*PI1)/(gamma1 - 1.)*(f[] - frho1[]*b1) + frho1[]*q1;
fE2[] = (1. - f[])*(p0/(gamma2 - 1.));
}
}
}

We log the evolution of the bubble radius.

event centroid (i += 20)
{
double volume = 0.;
foreach(reduction(+:volume))
volume += dv()*(1. - f[]);
#if AXIS
volume /= 2.;
#endif
fprintf (stderr ,"%g %g\n", t/tau, pow(3.*volume,1./3.));
}

Output of some statistics about the fields.

event logfile (i++)
{
stats sp = statsf (p), su = statsf (q.x), sT = statsf (T);
if (i == 0)
fprintf (stdout, "t dt max(p) max(T) max(u)\n");
fprintf (stdout, "%g %g %g %g %g\n",
t/tau, dt/tau, sp.max, sT.max, su.max);
}

On the fly movie generation.

#if AXIS
event movie (t += 0.01*4737.81)
{
view (fov = 12.5, quat = {0,0,-cos(M_PI/4.),cos(M_PI/4.)}, width = 640, height = 990);
draw_vof ("f");
squares ("p", min = 1., map = cool_warm);
char s[80];
sprintf (s, "t = %.2f", t/4737.81);
mirror({0,1}) {
draw_vof ("f");
squares ("T", min = 1.0896, max = 2.1792, map = cool_warm);
draw_string (s, pos = 2, size = 16, lc = {255,255,255}, lw = 4);
draw_string ("Temperature", pos = 3, size = 25, lc = {255,255,255}, lw = 4);
draw_string ("Pressure", size = 25, lc = {255,255,255}, lw = 4);
}
save ("T2.mp4");

Saving dump files for post-processing. (Uncomment)

#if 0
char name[80];
sprintf (name,"dump-%g",t/4737.81);
dump (name, list = (scalar *){f,p,T});
#endif
}
#endif // AXIS

event ending (t = tend);