# 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;
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;
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;
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.;
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;
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;
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;
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;
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;
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;
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() {
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));
}
char fname;
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;
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.

The marvelous design of this website is taken from Suckless.org