src/test/spurious.c

Circular droplet in equilibrium

This is the classical “spurious” or “parasitic currents” test case discussed in Popinet, 2009.

We use the Navier–Stokes solver with VOF interface tracking and surface tension.

#define JACOBI 1

#include "navier-stokes/centered.h"
#include "vof.h"
#include "tension.h"

The interface is represented by the volume fraction field c.

scalar c[], * interfaces = {c};

The diameter of the droplet is 0.8. The density is constant (equal to unity by default), and the viscosity is defined through the Laplace number \displaystyle La = \sigma\rho D/\mu^2 with \sigma set to one. The simulation time is set to the characteristic viscous damping timescale.

#define DIAMETER 0.8
#define MU sqrt(DIAMETER/LAPLACE)
#define TMAX (sq(DIAMETER)/MU)

We will vary the number of levels of refinement (to study the convergence), the Laplace number and DC a convergence parameter which measures the variation in volume fraction between successive timesteps (to evaluate whether we are close to a steady solution).

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.

  TOLERANCE = 1e-6 [*];
  stokes = true;
  c.sigma = 1;
  LEVEL = 5;
  N = 1 << LEVEL;
  for (LAPLACE = 120; LAPLACE <= 12000; LAPLACE *= 10)
    run();

We now fix the Laplace number and look for stationary solutions (i.e. the volume fraction field is converged to within 1e-10) for varying spatial resolutions.

  LAPLACE = 12000; DC = 1e-10;
  for (LEVEL = 3; LEVEL <= 7; LEVEL++) 
    if (LEVEL != 5) {
      N = 1 << LEVEL;
      run();
    }
}

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

scalar cn[];

event init (i = 0) {

We set the constant viscosity field…

  const face vector muc[] = {MU,MU};
  mu = muc;

… 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.

  fraction (c, sq(DIAMETER/2) - sq(x) - sq(y));
  foreach()
    cn[] = c[];
}

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

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

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

And we output the evolution of the maximum velocity.

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

event error (t = end) {

At the end of the simulation, we compute the equivalent radius of the droplet.

  double vol = statsf(c).sum;
  double radius = sqrt(4.*vol/pi);

We recompute the reference solution.

  scalar cref[];
  fraction (cref, sq(DIAMETER/2) - sq(x) - sq(y));

And compute the maximum error on the curvature ekmax, the norm of the velocity un and the shape error ec.

  double ekmax = 0.;
  scalar un[], ec[], kappa[];
  curvature (c, kappa);
  foreach() {
    un[] = norm(u);
    ec[] = c[] - cref[];
    if (kappa[] != nodata) {
      double ek = fabs (kappa[] - (/*AXI*/ + 1.)/radius);
      if (ek > ekmax)
	ekmax = ek;
    }
  }

We output these on standard error (i.e. the log file).

  norm ne = normf (ec);
  fprintf (stderr, "%d %g %g %g %g %g %g\n", 
	   LEVEL, LAPLACE, 
	   normf(un).max*sqrt(DIAMETER), 
	   ne.avg, ne.rms, ne.max,
	   ekmax);
}

#if 0
event gfsview (i += 10) {
  static FILE * fp = popen ("gfsview2D spurious.gfv", "w");
  output_gfs (fp);
}
#endif

We use an adaptive mesh with a constant (maximum) resolution along the interface.

#if TREE
event adapt (i <= 10; i++) {
  adapt_wavelet ({c}, (double[]){0}, maxlevel = LEVEL, minlevel = 0);
}
#endif

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"

{u}

The equilibrium shape and curvature converge toward the exact shape and curvature at close to second-order rate.

set xlabel 'D'
set ylabel 'Shape error'
set logscale x
set xtics 2
set pointsize 1
plot [5:120]'< sort -n -k1,2 log' u (0.8*2**$1):5 w lp t "RMS", \
            '< sort -n -k1,2 log' u (0.8*2**$1):6 w lp t "Max", \
             0.2/(x*x) t "Second order"

Convergence of the error on the equilibrium shape of the droplet with resolution. The diameter is given in number of grid points. (script)

set ylabel 'Relative curvature error'
plot [5:120]'< sort -n -k1,2 log' u (0.8*2**$1):($7/2.5) w lp t "Max", \
             0.6/(x*x) t "Second order"

Convergence of the relative error on the equilibrium curvature value with resolution. The diameter is given in number of grid points. (script)

See also

• Same test with Gerris