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)

    Bubble radius as a function of time (script)

    Animation of the pressure and temperature fields.

    [saade2023]

    hal-03950917 not found or server error.

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