src/test/shrinking.c

    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'
    Bubble radius as a function of time (script)

    Animation of the pressure and temperature fields.

    [saade2023]

    Youssef Saade, Detlef Lohse, and Daniel Fuster. A multigrid solver for the coupled pressure-temperature equations in an all-Mach solver with VoF. Journal of Computational Physics, 476:111865, 2023. [ DOI | http | .pdf ]

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