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.

    Animation of the free-surface relaxation

    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)

    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;

    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;
          h[] = H/nl;


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