Two-phase flow around RV Tangaroa

    This is an improved version of the famous Gerris example, illustrating the combination of complex solid boundaries, air-water turbulent flows and reduced gravity approach.

    We use the centered Navier–Stokes solver, two-phase flow and the momentum-conserving option. Note that the momentum-conserving option is crucial to obtain stable solutions for this air-water density ratio configuration.

    We also need to compute distance functions (to describe the ship geometry), use reduced gravity and visualisation functions.

    #include "distance.h"
    #include "reduced.h"
    #include "view.h"
    #include "lambda2.h"

    On supercomputers we need to control the maximum runtime and we check performances.

    #include "maxruntime.h"
    #include "navier-stokes/perfs.h"

    Importing the geometry

    This function computes the solid fraction given a pointer to an STL file, a tolerance (maximum relative error on distance) and a maximum level.

    void fraction_from_stl (scalar f, FILE * fp, double eps, int maxlevel)

    We read the STL file and compute the bounding box of the model.

      coord * p = input_stl (fp);
      coord min, max;
      bounding_box (p, &min, &max);
      double maxl = -HUGE;
        if (max.x - min.x > maxl)
          maxl = max.x - min.x;

    We initialize the distance field on the coarse initial mesh and refine it adaptively until the threshold error (on distance) is reached.

      scalar d[];
      distance (d, p);
      while (adapt_wavelet ({d}, (double[]){eps*maxl}, maxlevel, 5).nf);

    We also compute the volume fraction from the distance field. We first construct a vertex field interpolated from the centered field and then call the appropriate VOF functions.

      vertex scalar phi[];
        phi[] = (d[] + d[-1] + d[0,-1] + d[-1,-1] +
    	     d[0,0,-1] + d[-1,0,-1] + d[0,-1,-1] + d[-1,-1,-1])/8.;
      fractions (phi, f);

    Main function

    We can change both the maximum level of refinement and the Froude number at runtime.

    RV Tangaroa is 70 metres long. If we assume that it moves at 20 knots (twice its actual cruise speed), this gives a Froude number of approx 0.4.

    int LEVEL = 9;
    double FROUDE = 0.4;

    We need additional (fraction) fields for the ship geometry and for the (inflow) boundary condition.

    scalar tangaroa[], f0[];
    int main (int argc, char * argv[])
      maxruntime (&argc, argv);
      if (argc > 1)
        LEVEL = atoi(argv[1]);
      if (argc > 2)
        FROUDE = atof(argv[2]);
      init_grid (32);
      rho1 = 1.; // water
      rho2 = 1./815.; // air

    The length of the ship is unity and the domain is five times larger. We change the origin so that the ship is not too close to the inflow.

      size (5.);
      origin (-L0/2.,-L0/3.,-L0/2.);

    We need to tell the code that both tangaroa and f0 are volume fraction fields.

      for (scalar s in {tangaroa,f0})
        s.refine = s.prolongation = fraction_refine;

    Since the ship length is one and the velocity one, the acceleration of gravity is…

      G.z = - 1./sq(FROUDE);

    Boundary conditions

    The inflow condition fixes the velocity (unity) and the water level (using f0).

    u.n[bottom] = dirichlet(1);
    p[bottom]   = neumann(0.);
    pf[bottom]  = neumann(0.);
    f[bottom]   = f0[];

    Outflow uses standard Neumann/Dirichlet conditions.

    u.n[top]  = neumann(0.);
    p[top]    = dirichlet(0.);
    pf[top]   = dirichlet(0.);

    Boundary conditions for the solid and fraction tracers.

    tangaroa[back] = 0;
    f[back] = 1;

    Not sure whether this is really useful.

    uf.n[left] = 0.;
    uf.n[right] = 0.;

    Initial conditions

    We can optionally restart, otherwise we open the STL file and initialize the corresponding fraction. We also initialize the f0 field used for the inflow condition and set the initial water level and velocity field.

    event init (t = 0) {
      if (!restore (file = "restart")) {
        FILE * fp = fopen ("tangaroa.stl", "r");
        fraction_from_stl (tangaroa, fp, 5e-4, LEVEL);
        fclose (fp);
        fraction (f0, - z);
        foreach() {
          f[] = f0[];
          u.y[] = 1.;

    We setup the default display.

      display ("view (quat = {0.542,0.150,0.209,0.799},"
    	   "      fov = 30, near = 0.01, far = 1000,"
    	   "      tx = -0.0074, ty = -0.0138, tz = -0.3211);"
    	   "draw_vof (c = 'f', color = 'z', min = -0.1, max = 0.1, "
    	   "          linear = true, fc = {0.45,0.78,0.92});"
    	   "draw_vof (c = 'tangaroa');", true);

    Boundary conditions on the ship

    We use a simple (but crude) imposition of u=0 in the solid.

    event velocity (i++) {
          u.x[] = (1. - tangaroa[])*u.x[];


    We generate animations of the ship surface (as represented by the solid fraction) and of the air-water interface, colored with the height.

    Several classical features of ship wakes are recognisable: breaking bow wave, breaking stern divergent wave, turbulent boundary layer, Kelvin waves etc…

    Evolution of the air-water interface

    We also use the \lambda_2 criterion to display the turbulent vortical structures generated by the airflow. The air recirculation at the top of the steep primary Kelvin waves is particularly noticeable.

    Turbulent vortical structures

    The computations above were done on the Irene supercomputer using 12 levels of refinement.

    event movie (t += 0.01; t <= 10)
      view (fov = 5.86528,
    	quat = {0.515965,0.140691,0.245247,0.808605},
    	tx = -0.07438, ty = -0.0612925,
    	width = 1024, height = 768);
      draw_vof ("tangaroa");
      scalar Z[];
      Z[back] = dirichlet (z);
        Z[] = z;
      draw_vof ("f", color = "Z", min = -0.1, max = 0.1, linear = true);
      save ("movie.mp4");
      draw_vof ("tangaroa", fc = {0.5,0.5,0.5});
      draw_vof ("f", color = "Z", min = -0.1, max = 0.1, linear = true);
      scalar l2[];
      lambda2 (u, l2);
      isosurface ("l2", -100);
      save ("l2.mp4");
    #if DUMP
    event snapshot (i += 100)
      char name[80];
      sprintf (name, "dump-%d", i);
      scalar l2[];
      lambda2 (u, l2);
      dump (file = name);

    Mesh adaptation

    This computation is only feasible thanks to mesh adaptation, based both on volume fraction and velocity accuracy.

    event adapt (i++) {
      double uemax = 0.1;
      adapt_wavelet ({f,tangaroa,u},
    		 (double[]){0.01,0.01,uemax,uemax,uemax}, LEVEL, 5);

    Running in parallel on Irene

    Running with MPI-parallelism is a bit more complicated than usual since the distance() function is not parallelised yet. A reasonably simple workaround is to first generate a restart/dump file on the local machine, without MPI, using something like:

    CFLAGS=-DDUMP=1 make tangaroa.tst

    (and also adjust the maximum level), then kill the running code (using Ctrl-C) and do:

    qcc -source -D_MPI=1 tangaroa.c
    scp _tangaroa.c
    scp tangaroa/dump-0

    then on irene (to run on 480 cores for 10 hours, with 12 levels of refinement):

    ccc_mpp -u popinets
    sed 's/WALLTIME/36000/g' | ccc_msub -n 480

    with the following script

    #MSUB -r tangaroa
    #MSUB -@,end
    #MSUB -o basilisk_%I.out
    #MSUB -e basilisk_%I.log
    #MSUB -q skylake
    #MSUB -A gen7760
    #MSUB -m scratch
    ##MSUB -w
    set -x
    mpicc -Wall -std=c99 -O2 _tangaroa.c -o tangaroa -I$HOME -L$HOME/gl -lglutils -lfb_osmesa -lOSMesa -lGLU -lm
    ccc_mprun -n ${BRIDGE_MSUB_NPROC} ./tangaroa -m WALLTIME 12 0.4 \
        2>> log-${BRIDGE_MSUB_NPROC} >> out-${BRIDGE_MSUB_NPROC}

    Generating MP4 movies

    Note that when running on Irene the ffmpeg MP4 encoder is not available and the logfile will contain the warning:

    open_image(): cannot find 'ppm2mp4' or 'ffmpeg'/'avconv'
      falling back to raw PPM outputs

    The MP4 files defined above will be renamed movie.mp4.ppm and l2.mp4.ppm. As the extension indicates, these (large) files are now raw (uncompressed) PPM images. To convert them to compressed MP4, you will need to copy them to a machine where ffmpeg (and Basilisk) are installed (i.e. your local machine) and do:

    ppm2mp4 movie.mp4 < movie.mp4.ppm
    ppm2mp4 l2.mp4 < l2.mp4.ppm