fast-oscillation.c

Shape oscillation of an inviscid droplet with density ratio of 1/10
- Results

Shape oscillation of an inviscid droplet with density ratio of 1/10

A two-dimensional elliptical droplet (density ratio 1/10) is released in a large domain. Under the effect of surface-tension the shape of the droplet oscillates around its (circular) equilibrium shape. The fluids inside and outside the droplet are inviscid so ideally no damping of the oscillations should occur. As illustrated on the figures some damping occurs in the simulation due to numerical dissipation.

This simulation is also a stringent test case of the accuracy of the surface tension representation as no explicit viscosity can damp eventual parasitic currents.

We use a the Navier–Stokes solver with VOF interface tracking and surface tension using FastP* method for solving poisson equation.

#include "fast-navierstokes.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.2. The density inside the droplet is one and outside 10^-1.

#define D 0.2
#define rho(c) (clamp(c,0,1)*(1. - 1e-1) + 1e-1)

FILE * fp = NULL;

int LEVEL;

int main() {

The surface tension is unity. We cleanup existing files and vary the level of refinement.

  c.sigma = 1;
  remove ("error");
  remove ("laplace");

We choose FastP* method for solving poisson equation

  poisson_solver=2;

We set rho0 parameter as smaller density

  rho0=1e-1;

We reducing the TOLERANCE to imporve resutls

TOLERANCE=1e-6;

We will vary the level of refinement to study convergence.

  for (LEVEL = 5; LEVEL <= 8; LEVEL++) {
    N = 1 << LEVEL;

We open a file indexed by the level to store the time evolution of the kinetic energy.

    char name[80];
    sprintf (name, "k-%d", LEVEL);
    fp = fopen (name, "w");
    run();
    fclose (fp);
  }

We use grep to filter the lines generated by gnuplot containing the results of the fits (see below).

  system ("grep ^fit out >> log");
}

We have the option of using some “smearing” of the density jump.

#if 0
#define cf c
#else
scalar cf[];
#endif

The density is variable. We allocate a new field to store its inverse.

face vector alphav[];

event init (i = 0) {

We initialise the shape of the interface, a slightly elliptic droplet.

  fraction (c, D/2.*(1. + 0.05*cos(2.*atan2(y,x))) - sqrt(sq(x) + sq(y)));
#ifndef cf
  foreach()
    cf[] = c[];
#endif

The density is variable, and defined by the function below.

  alpha = alphav;
}

The density is defined at each timestep via the properties() event declared by the Navier–Stokes solver.

event properties (i++) {

When using smearing of the density jump, we initialise cf with the vertex-average of c.

#ifndef cf
  foreach()
    cf[] = (4.*c[] + 
	    2.*(c[0,1] + c[0,-1] + c[1,0] + c[-1,0]) +
	    c[-1,-1] + c[1,-1] + c[1,1] + c[-1,1])/16.;
  boundary ({cf});
#endif

The inverse of the density \alpha is then given by the face-averaged value of cf and the arithmetic average of density defined by rho().

  foreach_face() {
    double cm = (cf[] + cf[-1])/2.;
    alphav.x[] = 1./rho(cm);
  }
  boundary ((scalar *){alphav});  
}

By default tracers are defined at t-\Delta t/2. We use the first keyword to move VOF advection before the logfile output i.e. at t+\Delta/2. This improves the results.

event vof (i++, first);

At each timestep we output the kinetic energy.

event logfile (i++; t <= 1) {
  double ke = 0.;
  foreach (reduction(+:ke))
    ke += sq(Delta)*(sq(u.x[]) + sq(u.y[]))*rho(cf[]);
  fprintf (fp, "%g %g %d\n", t, ke, mgp.i);
  fflush (fp);
}

At the end of the simulation, we use gnuplot to fit a function of the form \displaystyle k(t) = ae^{-bt}(1-\cos(ct)) to the kinetic energy. This gives estimates of the oscillation pulsation c and of the damping b.

event fit (t = end) {
  FILE * fp = popen ("gnuplot 2>&1", "w");
  fprintf (fp, 
           "k(t)=a*exp(-b*t)*(1.-cos(c*t))\n"
           "a = 3e-4\n"
           "b = 1.5\n"
	   "\n"
           "D = %g\n"
           "n = 2.\n"
           "sigma = 1.\n"
           "rhol = 1.\n"
           "rhog = 1./10.\n"
           "r0 = D/2.\n"
           "omega0 = sqrt((n**3-n)*sigma/((rhol+rhog)*r0**3))\n"
	   "\n"
           "c = 2.*omega0\n"
           "fit k(x) 'k-%d' via a,b,c\n"
	   "level = %d\n"
	   "res = D/%g*2.**level\n"
           "print \"fit \", res, a, b, c, D\n"
	   "\n"
	   "set table 'fit-%d'\n"
	   "plot [0:1] 2.*a*exp(-b*x)\n"
	   "unset table\n"
	   "\n"
	   "set print 'error' append\n"
	   "print res, c/2./omega0-1., D\n"
	   "\n"
	   "set print 'laplace' append\n"
	   "empirical_constant = 30.\n"
	   "print res, (1./(b**2.*D**3.))*empirical_constant**2, D\n"
	   "\n",
	   D, LEVEL, LEVEL, L0, LEVEL);
  pclose (fp);
}


#if TREE
event adapt (i++) {
  adapt_wavelet ({c,u}, (double[]){5e-3,1e-3,1e-3}, LEVEL);
  event ("properties");
}
#endif

Results

set output 'k.png'
set xlabel 'Time'
set ylabel 'Kinetic energy'
set logscale y
plot [0:1][8e-5:1e-3]'k-8' t "51.2" w l, 'k-7' t "25.6" w l,               \
  'k-6' t "12.8" w l, 'k-5' t "6.4" w l,			       \
  'fit-8' t "" w l lt 7, 'fit-7' t "" w l lt 7, 'fit-6' t "" w l lt 7, \
  'fit-5' t "" w l lt 7

Evolution of the kinetic energy as a function of time for the spatial resolutions (number of grid points per diameter) indicated in the legend. The black lines are fitted decreasing exponential functions. (script)

An equivalent Laplace number can then be computed as: \displaystyle La=\frac{\sigma D}{\rho\nu^2}=\frac{\sigma C^2}{\rho b^2 D^3} The equivalent Laplace number depends on spatial resolution as illustrated below.

set output 'laplace.png'
set xlabel 'Diameter (grid points)'
set ylabel 'Equivalent Laplace number'
set grid
plot 'laplace' t "" w p pt 5 ps 2

Equivalent Laplace number estimated from the numerical damping of kinetic energy. (script)