sandbox/piersonj/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).
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)…