sandbox/crobert/2_Implicit/1_Surface_tension/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 to 45 degrees. A constant surface tension is imposed and the viscous fluid relaxes toward its equilibrium position: a 3D meniscus i.e. the intersection of a sphere and a cube with the given volume and contact angles.

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 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)

#include "grid/multigrid.h"
#include "crobert/2_Implicit/hydro-tension.h"
#include "layered/implicit.h"
#include "layered/remap.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.;
  }
  boundary ((scalar *){ha, hu, hf});  
}

event init (i = 0)
{
  foreach_dimension() {
    eta[left]    = neumann (1./tan(45.*pi/180.));
    eta[right]   = neumann (1./tan(45.*pi/180.));
  }

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");
}