src/test/sessile.c

Sessile drop

A sessile drop is a drop of liquid at rest on a solid surface. In the absence of gravity, the shape of the drop is controlled by surface tension only. An important parameter is the “contact angle” θ between the solid surface and the interface. In the absence of gravity, the drop is hemispherical and it is easy to show that the relation between the radius of the drop R and its volume V is (for two-dimensional drops) V=R2(θsinθcosθ)

To test this relation, a drop is initialised as a half-disk (i.e. the initial contact angle is 90$^\circ$) and the contact angle is varied between 15$^\circ$ and 165$^\circ$. The drop oscillates and eventually relaxes to its equilibrium position. This equilibrium is exact to within machine accuracy. The curvature along the interface is constant.

Note that shallower angles are not accessible yet.

Equilibrium shapes for 15^\circ \leq \theta \leq 165^\circ

Equilibrium shapes for 15θ165

Note that in the figure above the jagged interface representation near the contact line is a visualisation artefact.

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

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

To set the contact angle, we allocate a height-function vector field and set the contact angle boundary condition on its tangential component.

vector h[];
double theta0 = 30;
h.t[bottom] = contact_angle (theta0*pi/180.);

int main()
{
  size (2);

We use a constant viscosity.

  const face vector muc[] = {.1,.1};
  μ = muc;

We must associate the height function field with the VOF tracer, so that it is used by the relevant functions (curvature calculation in particular).

  f.height = h;

We set the surface tension coefficient and run for the range of contact angles.

  f.σ = 1.;

  for (theta0 = 15; theta0 <= 165; theta0 += 15)
    run();
}

The initial drop is a quarter of a circle.

event init (t = 0)
{
  fraction (f, - (sq(x) + sq(y) - sq(0.5)));
}

#if 0
event logfile (i++)
{
  fprintf (fout, "%g %g\n", t, normf(u.x).max);
}

event snapshots (t += 1)
{
  p.nodump = false;
  dump();
}
#endif

At equilibrium (t = 10 seems sufficient), we output the interface shape and compute the (constant) curvature.

event end (t = 10)
{
  output_facets (f, stdout);
  
  scalar κ[];
  curvature (f, κ);
  stats s = statsf (κ);
  double R = s.volume/s.sum, V = 2.*statsf(f).sum;
  fprintf (ferr, "%d %g %.5g %.5g\n", N, theta0, R/sqrt(V/π), s.stddev);
}

We compare R/R0 to the analytical expression, with R0=V/π.

See also