sandbox/fuster/Allmach3.0/spurioustest.c

    Circular droplet in equilibrium

    This is the classical “spurious” or “parasitic currents” test case discussed in Popinet, 2009 to test the compressible Navier–Stokes solver with VOF interface tracking and surface tension.

    #define JACOBI 1

#include "two-phase-compressible.h"
#include "compressible-tension.h"

#define DIAMETER 0.8
#define MU sqrt(DIAMETER/LAPLACE)
/*Change 0.1 to 5 to obtain the results plotted at the end of this page */
#define TMAX (5.*sq(DIAMETER)/MU)

int LEVEL;
double LAPLACE;
double DC = 0.;
FILE * fp = NULL;

int main() {

    We neglect the advection terms and vary the Laplace, for a constant resolution of 5 levels.

      PI1 = 300.;
  gamma1 = 7.14;
  gamma2 = 1.4;

  TOLERANCE = 1e-6 [*];
  L0 = 1. [0];
  DT = HUGE [0];
  
  f.sigma = 1;
  f.gradient = zero;
  LEVEL = 5;
  N = 1 << LEVEL;
  for (LAPLACE = 120; LAPLACE <= 12000; LAPLACE *= 10) {

    We set the constant viscosity field…

        mu1 = mu2 = MU;
    run();
  }
}

    We allocate a field to store the previous volume fraction field (to check for stationary solutions).

    scalar cn[];

event init (i = 0) {

    … open a new file to store the evolution of the amplitude of spurious currents for the various LAPLACE, LEVEL combinations…

      char name[80];
  sprintf (name, "La-%g-%d", LAPLACE, LEVEL);
  if (fp)
    fclose (fp);
  fp = fopen (name, "w");

    … and initialise the shape of the interface and the initial volume fraction field.

      double p0L = 1.;
  double p0 = p0L + f.sigma/DIAMETER*2;
  fraction (f, sq(DIAMETER/2) - sq(x) - sq(y));
  foreach() {
    cn[] = f[];
    frho1[]  = f[];
    frho2[]  = (1. - f[]);
    double pL = p0L;
    p[]   = pL*f[] + p0*(1.-f[]);
    fE1[]   = f[]*(pL/(gamma1 - 1.) + PI1*gamma1/(gamma1 - 1.));
    fE2[]   = (1.-f[])*p0/(gamma2 - 1.);
    q.x[] = 0.;
    q.y[] = 0.;
  }
  boundary ({cn});
}

event logfile (i++; t <= TMAX)
{

    At every timestep, we check whether the volume fraction field has converged.

      double dc = change (f, cn);
  if (i > 1 && dc < DC)
    return 1; /* stop */

    And we output the evolution of the maximum velocity.

      scalar un[];
  foreach()
    un[] = norm(q);
  fprintf (fp, "%g %g %g\n",
     MU*t/sq(DIAMETER), normf(un).max*sqrt(DIAMETER), dc);
}

    Results

    The maximum velocity converges toward machine zero for a wide range of Laplace numbers on a timescale comparable to the viscous dissipation timescale, as expected.

    set xlabel 't{/Symbol m}/D^2'
set ylabel 'U(D/{/Symbol s})^{1/2}'
set logscale y
plot 'La-120-5' w l t "La=120", 'La-1200-5' w l t "La=1200", 'La-12000-5' w l t "La=12000"

    Evolution of the amplitude of the capillary currents $(|(script){u}

    The high frequency oscillations scale with the bubble resonance frequency \omega_0 = \frac{1}{R_0}\sqrt{\frac{3 \gamma p_0}{\rho_l}}

    set xlabel 't f_0'
set ylabel 'U(D/{/Symbol s})^{1/2}'
set logscale y
DIAMETER=0.8
plot [0:10] 'La-120-5' u ($1*sqrt(120*3*1.4*(1. + 2./DIAMETER))/(pi*DIAMETER)):2 w l t "La=120", \
     'La-1200-5' u ($1*sqrt(1200*3*1.4*(1. + 2./DIAMETER))/(pi*DIAMETER)):2 w l t "La=1200", \
     'La-12000-5' u ($1*sqrt(12000*3*1.4*(1 + 2./DIAMETER))/(pi*DIAMETER)):2 w l t "La=12000"
    In the unit system used this implies \tau'=t \omega_0/(2 \pi) (script)

    In the unit system used this implies \tau'=t \omega_0/(2 \pi) (script)

    See also