3D Sessile drop

This is the 3D equivalent of the 2D test case.

The volume of a spherical cap of radius $R$ and (contact) angle $\theta$ is $V=\frac{\pi }{3}{R}^{3}\left(2+\mathrm{cos}\theta \right)\left(1-\mathrm{cos}\theta {\right)}^{2}$ or equivalently $\frac{R}{{R}_{0}}={\left(\frac{1}{4}\left(2+\mathrm{cos}\theta \right)\left(1-\mathrm{cos}\theta {\right)}^{2}\right)}^{-1/3}$ with ${R}_{0}$ the equivalent radius of the droplet ${R}_{0}={\left(\frac{3V}{4\pi }\right)}^{1/3}$ To test this relation, a drop is initialised as a half-sphere (i.e. the initial contact angle is 90\$^\circ\$) and the contact angle is varied between 30\$^\circ\$ and 150\$^\circ\$. The drop oscillates and eventually relaxes to its equilibrium position. The curvature along the interface is close to constant.

Relaxation toward a ${30}^{\circ }$ contact angle.

Note that shallower angles are not accessible yet.

#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "contact.h"
#include "vof.h"

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

#include "tension.h"
#include "view.h"

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

double theta0 = 45;

vector h[];
h.t[back] = contact_angle (theta0*pi/180.);
h.r[back] = contact_angle (theta0*pi/180.);

#define MAXLEVEL 5

int main()
{

We use a constant viscosity.

const face vector muc[] = {.1,.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.;

N = 1 << MAXLEVEL;
for (theta0 = 30; theta0 <= 150; theta0 += 30)
run();
}

The initial drop is a quarter of a sphere.

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

We log statistics on the maximum velocity, curvature and volume. If the standard deviation of curvature falls below ${10}^{-2}$, we assume that the steady shape is reached and we stop the calculation.

event logfile (i += 10; t <= 10)
{
scalar κ[];
cstats cs = curvature (f, κ);
foreach()
if (f[] <= 1e-3 || f[] >= 1. - 1e-3)
κ[] = nodata;
stats s = statsf (κ);
fprintf (fout, "%g %g %g %g %g %g %g %d %d %d %d\n", t, normf(u.x).max,
s.min, s.sum/s.volume, s.max, s.stddev, statsf(f).sum,
cs.h, cs.f, cs.a, cs.c);
fflush (fout);
if (s.stddev < 1e-2)
return 1; // stops
}

#if 0
event snapshots (i += 10)
{
scalar κ[];
curvature (f, κ);
p.nodump = false;
dump (buffered = true);
}
#endif

At equilibrium we output the (almost constant) radius, volume, maximum velocity and time.

event end (t = end)
{
scalar κ[];
curvature (f, κ);
stats s = statsf (κ);
double R = s.volume/s.sum, V = 4.*statsf(f).sum;
fprintf (ferr, "%d %g %.5g %.3g %.5g %.3g %.5g\n",
N, theta0, R, s.stddev, V, normf(u.x).max, t);
}

We make a movie of the relaxing interface for $\theta ={30}^{\circ }$. We use symmetries since only a quarter of the drop is simulated.

event movie (i += 5; t <= 3)
{
if (theta0 == 30.) {
view (fov = 26.6776, quat = {0.474458,0.144142,0.234923,0.836017},
tx = -0.0137556, ty = -0.00718937, bg = {1,1,1},
width = 758, height = 552);
draw_vof ("f", edges = true);
cells (lc = {1,0,0});
mirror (n = {1,0,0}) {
draw_vof ("f", edges = true);
cells (lc = {1,0,0});
}
mirror (n = {0,1,0}) {
draw_vof ("f", edges = true);
cells (lc = {1,0,0});
mirror (n = {1,0,0}) {
draw_vof ("f", edges = true);
cells (lc = {1,0,0});
}
}
save ("movie.mp4");
}
}

We use refinement based on a smooth version of the volume fraction. This guarantees constant refinement around the interface, which seems to be necessary to reach balance (this should be improved). Note however that this is specific to this test case and should not generally be used in “production” runs, which should work fine with the default criterion.

#if TREE
#if 1
scalar f1[];
foreach()
f1[] = f[];
boundary ({f1});
adapt_wavelet ({f1}, (double[]){1e-3}, minlevel = 3, maxlevel = MAXLEVEL);
#else
adapt_wavelet ({f}, (double[]){1e-4}, minlevel = 3, maxlevel = MAXLEVEL);
#endif
}
#endif

We compare $R/{R}_{0}$ to the analytical expression.