12. A gallery of fluid motion

On this page we solve some example flow problems, present their case file (.c) contents and the resulting movies. For a dummy file case.c the workflow is.

$ FILE=case
$ gcc -O3 -DN=250 $FILE.c -lm 
$ ./a.out
...
$ convert $FILE-*.ppm $FILE.mp4 
$ rm $FILE-*.ppm

A Kelvin-Helmholtz instability

A shear layer can roll up when it is perturbed. For our case, we consider a jet in the x direction, within which we initialize a tracer field.

#include "ns.h"
#include "visual.h"
// The jet is two-thirds of LO wide
#define JET (fabs(y) < L0/6.)

void output_frame() {
  double omega[N][N];
  vorticity (omega);
  char fname[99];
  sprintf (fname, "kh-%04d.ppm", iter);
  output_ppm (tracer, fname, -1, 1, false);
  // Log progress
  printf ("%d %g\n", iter, t);
}

int main() {
  L0 = 5;
  Y0 = -L0/2.;
  nu = 2e-4;
  t_end = 20;
  t_interval = .2;
  // Overload function
  every_t_interval = output_frame;
  // Initialize
  foreach() {
    val (ux, 0, 0) =  JET ? 0.5 : -0.5;
    val (uy, 0, 0) += 1e-3*noise();
    val (tracer, 0, 0) = JET ? 1. : 0.;
  }
  // call solver..
  run();
}

Decaying two-dimensional turbulence

Two-dimensional turbulence is characterized by a cascade of small vortices towards larger ones over time. In our case, we initialize a Taylor-Green solution and perturb it such that we can see the unstable deformation and the resulting chaos.

#include "ns.h"
#include "visual.h"

void output_frame() {
  double omega[N][N];
  vorticity (omega);
  char fname[99];
  sprintf (fname, "turbulence-%04d.ppm", iter);
  output_ppm (omega, fname, -.2, .2, false);
  // Log progress
  printf ("%d %g\n", iter, t);
}

int main() {
  L0 = 2.*pi;
  t_end = 500;
  t_interval = 2.;
  // Overload function
  every_t_interval = output_frame;
  // Initialize
  foreach() {
    val (ux, 0, 0) =  sin(5*x)*cos(5*y);
    val (uy, 0, 0) = -sin(5*y)*cos(5*x);
    val (uy, 0, 0) += 0.001*noise();
  }
  // call solver..
  run();
}

Colliding dipolar vortices

Two vortex dipoles can collide and exchange vortices. Instead of following the previous Lamb-Chaplygin model, we project two localized Gaussian jets, targeted for a head-on collision

#include "ns.h"
#include "visual.h"

void output_frame() {
  double omega[N][N];
  vorticity (omega);
  char fname[99];
  sprintf (fname, "dip-%04d.ppm", iter);
  output_ppm (omega, fname, -.5, .5, false);
  // Log progress
  printf ("%d %g\n", iter, t);
}

int main() {
  L0 = 15;
  X0 = Y0 = -L0/2.;
  t_end = 50;
  t_interval = .5;
  // Overload function
  every_t_interval = output_frame;
  // Initialize
  foreach() 
    val (ux, 0, 0) =  exp(-sq(x + 3) - sq(y)) - exp(-sq(x - 3) - sq(y));
  projection (ux, uy, p, 1);
  // reset p
  foreach()
    val(p, 0, 0) = 0;
  // call solver..
  run();
}

Vortex rebound from a wall

A dipole may also find a solid obstacle on its path. A crude way to implement such a condtion is to prescribe the solution of the stationary wall every timestep.

#include "ns.h"
#include "visual.h"
// The wall is located at the lhs ofthe domain
// such that it appears at the rhs
#define WALL (x < (X0 + L0/10.))

void output_frame() {
  double omega[N][N];
  vorticity (omega);
  char fname[99];
  sprintf (fname, "dip-wall-%04d.ppm", iter);
  output_ppm (omega, fname, -.1, .1, false);
  // Log progress
  printf ("%d %g\n", iter, t);
}

void boundary() {
  foreach() {
    if (WALL) {
      val (ux, 0, 0) = 0;
      val (uy, 0, 0) = 0;
    }
  }
}

int main() {
  L0 = 20;
  X0 = Y0 = -L0/2.;
  t_end = 250;
  t_interval = .5;
  nu = 1e-5;
  // Overload functions
  every_t_interval = output_frame;
  every_iter = boundary;
  // Initialize
  foreach() 
    val (ux, 0, 0) =  exp(-sq(x - 2) - sq(y));
  projection (ux, uy, p, 1);
  // reset p
  foreach()
    val(p, 0, 0) = 0;
  // call solver..
  run();
}

Rising warm plume

We can couple a tracer to the acceleration field to describe a buoyancy force with the Boussinesq approximation. In this case, we initialize a warm (buoyant) Gaussian plume

#include "ns.h"
#include "visual.h"

void output_frame() {
  char fname[99];
  sprintf (fname, "plume-%04d.ppm", iter);
  output_ppm (tracer, fname, -.1, .1, false);
  // Log progress
  printf ("%d %g\n", iter, t);
}

void buoyancy() {
  foreach() 
    val(acceleration_y, 0, 0) = val(tracer, 0, 0);
}

int main() {
  L0 = 20;
  X0 = Y0 = -L0/2.;
  t_end = 20;
  t_interval = .1;
  nu = 1e-5;
  // Overload functions
  every_t_interval = output_frame;
  every_iter = buoyancy;
  // Initialize
  foreach() 
    val (tracer, 0, 0) =  exp(-sq(x) - sq(y + 5));
  // call solver..
  run();
}

A Rayleigh-Taylor instability

A heavy fluid resting on top of a lighter fluid creates an unstable stratification. Using the Boussinesq approximation, we initialize a sharp interface between a heavier and a lighter region.

#include "ns.h"
#include "visual.h"

void output_frame() {
  double grad[N][N];
  foreach() {
    val(grad, 0, 0) = sq((val(tracer, 1, 0) - val(tracer, -1, 0))/(2*Delta));
    val(grad, 0, 0) += sq((val(tracer, 0, 1) - val(tracer, 0, -1))/(2*Delta));
    val(grad, 0, 0) = -(val (grad, 0, 0));
  }
  char fname[99];
  sprintf (fname, "rt-%04d.ppm", iter);
  output_ppm (grad, fname, -100, 100, false);
  // Log progress
  printf ("%d %g\n", iter, t);
}

void buoyancy() {
  foreach() 
    val (acceleration_y, 0, 0) = val (tracer, 0, 0);
}

int main() {
  L0 = 1;
  Y0 = -L0/2.;
  t_end = 5;
  t_interval = .05;
  nu = kappa = 1e-5;
  // Overload functions
  every_t_interval = output_frame;
  every_iter = buoyancy;
  // Initialize
  double sqN = 1;
  foreach() {
    val (tracer, 0, 0) =  (y < 0 ? L0/2*sqN : - L0/2.*sqN) + sqN*y;
    val (uy, 0, 0) = 0.01*noise();
  }
  // Finding the hydropstatic pressure field requires some addional sweeping
  max_sweeps = 1e5;
  // call solver..
  run();
}

Because I cant help myself. I also discuss the good, the bad and the Ugly of our resulting solver here.