src/test/meniscus3D.c

3D meniscus

An initial perturbed free-surface is contained within a cube. The contact angles on the faces of the cube are set by imposing a Neumann boundary condition on \eta as \displaystyle \frac{\partial \eta}{\partial n} = \frac{1}{\sqrt{R^2 - x^2 - y^2}}, with n the direction normal to the face and R = \sqrt{3}. These boundary conditions correspond to the static equilibrium solution for the intersection of a sphere of radius R with the cube. The contact angle varies between \text{atan} \left( 1 / \sqrt{2} \right) \approx 35 degrees on the axis of symmetry and 45 degrees (\partial_n \eta = 1) in the corners. A constant surface tension is imposed and the viscous fluid relaxes toward its equilibrium position.

Time-implicit timestepping is used and greatly decreases the number of timesteps required to reach convergence.

Both the maximum velocity and standard deviation on the curvature decrease exponentially.

set xlabel 'Time'
set logscale y
set format y "10^{%L}"
set xrange [0:5]
plot 'log' u 1:3 w l t 'Maximum velocity', \
     'log' u 1:5 w l t 'Standard deviation of curvature'

Convergence of the maximum velocity and standard deviation on curvature (script)

The final value of the mean curvature is close to the theoretical value 1/R = 1/\sqrt{3}.

set ylabel 'Mean curvature - 1/R'
plot 'log' u 1:($4/2. - 1./sqrt(3.)) w l t ''

Convergence of the mean curvature toward its theoretical value (script)

#include "grid/multigrid.h"
#include "layered/hydro-tension.h"
#include "layered/implicit.h"
#include "view.h"

int main()
{
  L0 = 2.;
  origin (-L0/2., -L0/2.);

  G = 0.;
  nu = 1.;

We use a large timestep and decrease the tolerance to get clear convergence of the maximum velocity and standard deviation of curvature.

  TOLERANCE = 1e-8;
  CFL_H = 40;
  
  N = 64;
  nl = 1;

  run();
}

Contact angle boundary conditions

This is necessary to balance the surface tension acceleration on the side walls. This needs to match the contact angle boundary conditions below.

event half_advection (i++)
{
  foreach_dimension() {
    ha.n[left] = 0.;
    ha.n[right] = 0.;
    hu.n[left] = 0.;
    hu.n[right] = 0.;
    hf.n[left] = 0.;
    hf.n[right] = 0.;
  }
}

event init (i = 0)
{
  foreach_dimension() {
    eta[left]  = neumann (1./sqrt(3. - x*x - y*y));
    eta[right] = neumann (1./sqrt(3. - x*x - y*y));
  }

Initial conditions

The initial free-surface shape.

  foreach() {
    double H = 0.5 + 0.1*cos(2.*pi*x)*cos(2.*pi*y);
    zb[] = - 0.5;
    foreach_layer()
      h[] = H/nl;
  }
}

Outputs

We compute the interface curvature (times \sigma, which is unity). This needs to be done in the pressure() event since the \sigma fields used by the sigma_kappa() macro are temporary.

scalar kappa[];

event pressure (i++)
{
  foreach()
    kappa[] = sigma_kappa (eta, 0);
}
  
event logfile (i++)
{
  stats s = statsf (kappa);
  fprintf (stderr, "%g %g %g %g %g %d\n", t, dt, normf(u.x).max,
	   s.sum/s.volume, s.stddev, mgH.i);
}

event movie (i++; t <= 5)
{
  view (fov = 22, quat = {0.475152,0.161235,0.235565,0.832313},
	tx = 0.02, ty = 0., width = 800, height = 600);
  char s[80];
  sprintf (s, "t = %.2f", t);
  draw_string (s, size = 80);
  squares ("eta", linear = true, z = "eta", min = -0.3, max = 0.3);
  save ("movie.mp4");
}