spurious.c

## 3D Circular droplet in equilibrium

This case is almost exactly the same as the one proposed by S. Popinet in 3D. All credit to him !

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

#define JACOBI 1

#include "grid/multigrid3D.h"
#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;
      init_grid(N);
      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) -sq(z));
  foreach()
    cn[] = c[];
  boundary ({cn});
}

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) - sq(z));

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

  scalar le[];
  foreach(){
    le[] = level;
  }
  char name_grid[100];
  sprintf (name_grid, "grid-%g-%d.png", LAPLACE, LEVEL);
  FILE* fgrid = fopen (name_grid, "w");
  output_ppm (le, file = name_grid,min=0,max=LEVEL,n = 100,box = {{0.,0.},{1.,1.}});
  fclose (fgrid);

}

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

## Results

The maximum velocity converges toward machine zero for a wide range of Laplace numbers (except for La=120)…