sandbox/fuster/Allmach3.0/oscillation.c

Shape oscillation of an inviscid droplet

This test case, discussed in Popinet, 2009 and included in the incompressible version , it is used here to test the compressible solver

A two-dimensional elliptical droplet (density ratio 1/1000) 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 the compressible solver with VOF interface tracking and surface tension.

#include "grid/multigrid.h"
#include "two-phase-compressible.h"
#include "compressible-tension.h"

The diameter of the droplet is 0.2.

#define D 0.2

We will vary the level of refinement to study convergence.

FILE * fp = NULL;
int LEVEL;
double rho1, rho2;
double p0L = 1.;
double CFLac = 5;

event stability (i++) {
  
  double cson = fmax(sqrt(gamma1*(p0L + PI1)/rho1), sqrt(gamma2*p0L/rho2));

  foreach () {
      DT    = CFLac/cson/pow(2,LEVEL);
      dtmax = CFLac/cson/pow(2,LEVEL);
  }

}

int main() {

The density is variable.

  rho1 = 1, rho2 = 1e-3;

  size(2.);
  X0 = Y0 = -L0/2.;

The surface tension is unity. Decreasing the tolerance on the Poisson solve improves the results. We cleanup existing files and vary the level of refinement.

  f.sigma = 1.;
  
  PI1 = 300.;
  gamma1 = 7.14;
  gamma2 = 1.4;

  TOLERANCE = 1e-4;
  remove ("error");
  remove ("laplace");
  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");
}

event init (i = 0) {

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

  fraction (f, D/2.*(1. + 0.05*cos(2.*atan2(y,x))) - sqrt(sq(x) + sq(y)));
  
  double p0 = p0L + f.sigma/D*2;
  foreach() {
    frho1[]  = rho1*f[];
    frho2[]  = rho2*(1. - f[]);
    double pL = p0L + f.sigma;
    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.;
  }

}

At each timestep we output the kinetic energy.

event logfile (i++; t <= 0.5) {
  double ke = 0.;
  foreach (reduction(+:ke))
    ke += sq(Delta)*(sq(q.x[]) + sq(q.y[]))/(frho1[]+frho2[]);
  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.

We also compute the relative error on the pulsation, using the theoretical value \omega_0 as reference.

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./1000.\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 sprintf (\"fit %%g %%.6f %%.2f %%.0f\\n\", 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);
}

Results

set xlabel 'Time'
set ylabel 'Kinetic energy'
set logscale y
set yrange[0.0001:0.01]
plot [0:0.5][8e-5:]'k-8' t "25.6" w l, 'k-7' t "12.8" w l,               \
  'k-6' t "6.4" w l, 'k-5' t "3.2" 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)

set xlabel 'Diameter (grid points)'
set ylabel 'Frequency error (%)'
set logscale x 2
set au y
unset grid
set xzeroaxis
set key spacing 1.5 top right
ftitle(a,b,c) = sprintf("%.0f/x^{%4.2f} (%s)", exp(a), -b, c)
f(x)=a+b*x
fit f(x) 'error' u (log($1)):(log(abs($2)*100.)) via a,b
plot 'error' u ($1):(abs($2)*100.) t "" w p pt 5 ps 2, \
      exp(f(log(x))) t ftitle(a,b,"standard")

Relative error in the oscillation frequency as a function of resolution. (script)

The amount of numerical damping can be estimated by computing an equivalent viscosity. With viscosity, kinetic energy is expected to decrease as: \displaystyle \exp(-C\nu/D^2t) where C is a constant, \nu the viscosity and D the droplet diameter. Using curve fitting the damping coefficient b=C\nu/D^2 can be estimated (black curves on Figure ). 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 xlabel 'Diameter (grid points)'
set ylabel 'Equivalent Laplace number'
set grid
set key bottom right
plot 'laplace' t "standard" w p pt 5 ps 2

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

See also

• Same test with the incompressible version
• Same test with Gerris