Contents

Fluid Motion Simulations and Artwork
Water Ripple Simulation  calculated in Mathematica 4.2, 12/2/04; rendered in POVRay 3.6.1, 7/29/05
Here is a simple water ripple simulation showing single slit wave diffraction. The following Mathematica code solves the wave equation with damping using the finite difference method. You can read more about this algorithm on Hugo Elias’ website. (Note: technically this simulation should use Neumann boundary conditions but I decided it was simpler to demonstrate using Dirichlet boundary conditions). See also the Shallow Water Equation.
(* runtime: 18 seconds, c is the wave speed and b is a damping factor *)
n = 64; c = 1; b = 5; dx = 1.0/(n  1); Courant = Sqrt[2.0]/2;dt = Courant dx/c; z1 = z2 = Table[0.0, {n}, {n}];
Do[{z1, z2} = {z2, z1}; z1[[n/2, n/4]] = Sin[16Pi t]; Do[If[0.45 < i/n < 0.55  ! (0.48 < j/n < 0.52), z2[[i, j]] = z1[[i, j]] + (z1[[i, j]]  z2[[i, j]] + Courant^2 (z1[[i  1, j]] + z1[[i + 1, j]] + z1[[i, j  1]] + z1[[i, j + 1]]  4 z1[[i, j]]))/(1 + b dt)], {i, 2, n  1}, {j, 2, n  1}]; ListPlot3D[z1, Mesh > False, PlotRange > {1, 1}], {t, 0, 1, dt}];
Links
LSSM  POVRay animations using this same technique, by Tim Nikias Wenclawiak
POVRay Ripples  another one by Rune Johansen
Water Ripples  Delphi program by Jan Horn
Waves  C++ program by Maciej Matyka

Unstable Vortex Filaments  C#, POVRay 3.6.1, 9/6/12
Here is a simulation of two unstable leapfrogging vortex filaments in 3D, assuming inviscid incompressible potential flow. The vortex filaments were simulated using the BiotSavart law for vortex segments. Here is some Mathematica code:
(* runtime: 38 seconds *)
abs[x_] := Sqrt[x.x]; Normalize[x_] := x/abs[x]; Mean[x_] := Plus @@ x/Length[x];
v[p_] := Module[{v = 0.0}, Scan[({p1,p2} = #; r12 = p2  p1; zhat = Normalize[r12]; r0 = p1  p + (p  p1).zhat zhat; r = abs[r0]; If[r > 0.001, gamma = 1.0/abs[r12]; rcore = 0.15Sqrt[gamma]; rhat = Normalize[r0]; cos1 = rhat.Normalize[p1  p]; cos2 = rhat.Normalize[p2  p]; thetahat = Cross[zhat, rhat]; v += (gamma/r)(1  Exp[(r/rcore)^2])(cos1 + cos2)thetahat]) &, Partition[filament, 2, 1, 1]]; v];
dt = 0.001; filament = Table[{Cos[theta], Sin[theta], 0.0}, {theta, 0.0, 2Pi  0.001, Pi/18}];
Do[filament = Map[(# + v[#] dt) &, filament]; If[t > 190, p0 = Mean[filament]; Show[Graphics3D[Line[Map[#  p0 &, filament]], PlotRange > 2{{1, 1}, {1, 1}, {1, 1}}]]], {t, 1, 230}]
Links
Vortex Movies  by Tee Tai Lim, I especially like the movies of leapfrogging and vortex ring collisions
Turbulence Infinite  a single incredibly knotted vortex filament, by Mark Stock
Fuzzy Smoke  vortex filament simulation, by Daniel Piker

LeapFrogging Bubble Rings  Mathematica 4.2, POVRay 3.6.1, 12/31/05
LeapFrogging Vortices  Mathematica 4.2 version: 12/31/05; C++ version: 10/22/07
Here is a simple technique for modeling 2D vortices assuming inviscid incompressible potential flow (irrotational). This technique was inspired from Kerry Mitchell’s paper based on Ultra Fractal. Here is some Mathematica code to plot the entrained fluid:
(* runtime: 25 seconds, increase n for better resolution *)
tmax = 0.85; rcore = 0.1; Klist = {1, 1, 1, 1}; zlist0 = {1  0.5I, 1 + 0.5I, 0.5  0.5I, 0.5 + 0.5I}; m = Length[zlist0];
v[K_, z_, z0_] := Module[{r2 = Abs[z  z0]^2}, I K(z0  z)/r2(1  Exp[r2/rcore^2])];
zlist = NDSolve[Flatten[Table[{Subscript[z,i]'[t] == Sum[If[i == j, 0, v[Klist[[j]], Subscript[z, i][t], Subscript[z, j][t]]], {j, 1, m}], Subscript[z, i][0] ==zlist0[[i]]}, {i, 1, m}]], Table[Subscript[z, i][t], {i, 1, m}], {t, 0, tmax}][[1, All, 2]];
n = 23; image = Table[NDSolve[{z'[t] == Sum[v[Klist[[i]], z[t], zlist[[i]]], {i, 1, m}], z[tmax] == x + I y}, z[t], {t, 0, tmax}, MaxSteps > 5000][[1, 1, 2]] /. t > 0, {y, 1.125, 1.125, 2.25/n}, {x, 0.35, 1.9, 2.25/n}];
ListDensityPlot[Map[Sign[Im[#]]Arg[#] &, image, {2}], Mesh > False, Frame > False, AspectRatio > Automatic]
Here is some Mathematica code to numerically solve this using the 4th order RungeKutta method:
(* runtime: 21 seconds, increase n for better resolution *)
n = 23; tmax = 0.85; dt = tmax/100; rcore = 0.1; Klist = {1, 1, 1, 1}; zlist = {1  0.5I, 1 + 0.5I, 0.5  0.5I, 0.5 + 0.5I}; m = Length[zlist];
v[K_, z_, z0_] := Module[{r2 = Abs[z  z0]^2}, I K(z0  z)/r2(1  Exp[r2/rcore^2])];
Runge[z_] := Module[{}, k1 = dt vtot[z]; k2 = dt vtot[z + k1/2]; k3 = dt vtot[z + k2/2]; k4 = dt vtot[z + k3]; (k1 + 2 k2 + 2 k3 + k4)/6];
vtot[zlist_] := Table[Sum[If[i == j, 0, v[Klist[[j]], zlist[[i]], zlist[[j]]]], {j, 1, m}], {i, 1, m}]; zlists = Transpose[Table[zlist += Runge[zlist], {t, 0, tmax, dt}]];
vtot[z_] := Plus @@ v[Klist, z, zlists[[All, t]]]; image = Table[z = x + I y; Do[z = Runge[z], {t, 100, 1, 1}]; z, {y, 1.125, 1.125, 2.25/n}, {x, 0.35, 1.9, 2.25/n}];
ListDensityPlot[Map[Sign[Im[#]]Arg[#] &, image, {2}], Mesh > False, Frame > False, AspectRatio > Automatic]
Links
YouTube animation  a very similar looking animation
My Educational Project  some basic potential flows using Matlab and Mathematica
Nuclear Bombs  mushroom clouds are giant vortex rings
Birth of the Jellyfish  nice animation by Kerry Mitchell

Driven Cavity  calculated in Fortran 90, plotted in Mathematica 4.2, 11/2/05
Here is the flow inside a square box where the flow across the top of the box is moving to the right at Reynolds number 1000. This program uses the finite volume method to solve the NavierStokes equations assuming steady, incompressible, viscous, laminar flow. This was calculated on a 300×300 nonuniform grid and took 3 hours to run on my laptop. The animation shows the motion along the streamlines/pathlines. The following Mathematica code is not completely accurate at the boundaries, but it gives the basic idea:
(* runtime: 20 minutes *)
n = 50; dx = 1.0/(n  1); RE = 1000; relax = 0.4; residu = residv = residp = 1;
u = Table[If[j == 1  i == 2  i == n, 0, 1], {i, 1, n}, {j, 1, n}];
v = p = pstar = ap = an = as = aw = ae = sc = sp = du = dv = Table[0, {n}, {n}];
lisolv[i0_, j0_, phi0_] := Module[{phi = phi0, p = q =Table[0, {n}]}, Do[q[[j0  1]] = phi[[i, j0  1]]; Do[p[[j]] = an[[i, j]]/(ap[[i, j]] as[[i, j]]p[[j  1]]); q[[j]] = (as[[i, j]]q[[j  1]] + ae[[i, j]] phi[[i + 1, j]] + aw[[i, j]]phi[[i  1, j]] + sc[[i, j]])/(ap[[i, j]]  as[[i, j]]p[[j  1]]), {j, j0, n  1}]; Do[phi[[i,j]] = p[[j]]phi[[i, j + 1]] + q[[j]], {j, n  1, j0, 1}], {i, i0, n  1}]; phi];
calc[i0_, j0_] := Module[{}, {cn, cs, ce, cw} = 0.5 {v[[i,j + 1]] + v[[i0, j0 + 1]], v[[i, j]] + v[[i0, j0]], u[[i + 1, j]] + u[[i0 + 1, j0]], u[[i, j]] + u[[i0, j0]]}dx; cp = Max[0, cn  cs + ce  cw]; {an[[i, j]], as[[i, j]], ae[[i, j]], aw[[i, j]]} = Map[Max[0, 1  0.5 Abs[#]] &, RE{cn, cs, ce, cw}]/RE + Map[Max[0, #] &, {cn, cs, ce, cw}]; sp[[i, j]] = cp; ap[[i, j]] = an[[i, j]] + as[[i, j]] + ae[[i, j]] + aw[[i, j]]  sp[[i, j]]];
While[Max[residu, residv, residp] > 0.001, residu = residv = residp = 0; Do[calc[i  1, j]; sc[[i, j]] = cp u[[i, j]] + dx(p[[i  1, j]]  p[[i, j]]); du[[i, j]] = relax dx/ap[[i, j]]; residu += Abs[an[[i, j]]u[[i, j + 1]] + as[[i, j]]u[[i, j  1]] + ae[[i, j]]u[[i + 1, j]] + aw[[i, j]]u[[i  1, j]]  ap[[i, j]]u[[i, j]] + sc[[i, j]]], {i, 3, n  1}, {j, 2, n  1}]; ap /= relax; sc += (1  relax) ap u; u = lisolv[3, 2, u]; Do[calc[i, j  1]; sc[[i, j]] = cp v[[i, j]] + dx(p[[i, j  1]]  p[[i, j]]); dv[[i, j]] = relax dx/ap[[i, j]];residv += Abs[an[[i, j]]v[[i, j + 1]] + as[[i, j]]v[[i, j  1]] + ae[[i, j]]v[[i + 1, j]] + aw[[i, j]]v[[i  1, j]]  ap[[i, j]]v[[i, j]] + sc[[i, j]]], {i, 2, n  1}, {j, 3, n  1}]; ap /= relax; sc += (1  relax)ap v; v = lisolv[2, 3, v];Do[an[[i, j]] = dx dv[[i, j + 1]]; as[[i, j]] = dx dv[[i, j]]; ae[[i, j]] = dx du[[i + 1, j]]; aw[[i, j]] = dx du[[i, j]]; sp[[i, j]] = 0; sc[[i, j]] = ((v[[i, j + 1]]  v[[i, j]])dx + (u[[i + 1, j]]  u[[i, j]]) dx); residp += Abs[sc[[i, j]]]; ap[[i, j]] = an[[i, j]] + as[[i, j]] + ae[[i, j]] + aw[[i, j]]  sp[[i, j]]; pstar[[i, j]] = 0, {i, 2, n  1}, {j, 2, n  1}]; Do[pstar = lisolv[2, 2, pstar], {2}]; Do[If[i != 2, u[[i, j]] += du[[i, j]](pstar[[i  1, j]]  pstar[[i, j]])]; If[j != 2, v[[i, j]] += dv[[i, j]](pstar[[i, j  1]]  pstar[[i, j]])]; p[[i, j]] += 0.3(pstar[[i, j]]  pstar[[n  1, n  1]]), {i, 2, n  1}, {j, 2, n  1}]];

Here is some Mathematica code to plot streamlines. The blue lines are rotating clockwise and the red lines are rotating counterclockwise. The plot on the left agrees fairly well with Ercan Erturk’s results.
psi = Table[0, {n  1}, {n  1}]; Do[psi[[i, j]] = psi[[i, j  1]] + dx(u[[i + 1, j  1]] + u[[i + 1, j]])/2, {j, 2, n  1}, {i, 1, n  1}];
psi = ListInterpolation[psi, {{0, 1}, {0, 1}}];
ContourPlot[psi[x, y], {x, 0, 1}, {y, 0, 1}, PlotPoints > 50, PlotRange > All, ContourShading > False,Contours > {0.08, 0.077, 0.07, 0.06, 0.045, 0.025, 0.01, 0.0025, 0, 8*^6, 2*^6, 3*^5, 8*^5, 1*^4, 3*^4, 6*^4, 9*^4}, ContourStyle > Table[{Hue[2(1  x)/3]}, {x, 0, 1, 1/16}]]

Here is some Mathematica code to plot vorticity contours. The plot on the left agrees fairly well with Ghia’s results.
omega = ListInterpolation[Table[(v[[i, j]]  v[[i  1, j]])/dx  (u[[i, j]]  u[[i, j  1]])/dx, {i, 2, n  1}, {j, 2, n  1}], {{0, 1}, {0, 1}}];
ContourPlot[omega[x, y], {x, 0, 1}, {y, 0, 1},PlotPoints > 50, PlotRange > All, ContourShading > False, Contours > Range[14, 7],ContourStyle > Table[{Hue[2(1  x)/3]}, {x, 0, 1, 1/21}]]
Computational Fluid Dynamics (CFD) Links
CFD Gallery  by Arthur Veldman, I especially like the simulations of turbulence
Turbulent Boundary Layer  large simulation by Said Elghobashi
Colorful Fluid Mixing Gallery  by André Bakker
CFD code  simple Fortran codes by Abdusamad Salih
FlowLab  educational resource by Fluent

Driven Cavity  Mathematica 4.2, 11/16/06
Here is the same driven cavity using the finite difference method. This code uses the Alternating Direction Implicit (ADI) method to solve the vorticity transport equation assuming unsteady, incompressible, viscous flow. This method seems to be much faster and easier to implement in this case than the finite volume method.
(* runtime: 1 minute *)
<< LinearAlgebra`Tridiagonal`; n = 20; RE = 1000; dx = 1.0/(n  1); dt = 0.25; cxy = dt/(4dx); sxy = dt/(2RE dx^2); lambda = 1.5; free = Range[2, n  1];
omega = v = psi = Table[0, {n}, {n}]; u = Table[If[j == n, 1, 0], {n}, {j, 1, n}];
Do[omega[[All, 1]] = 2psi[[All, 2]]/dx^2; omega[[All, n]] = 2 (psi[[All, n  1]] + dx)/dx^2; omega[[1, All]] = 2 psi[[2, All]]/dx^2; omega[[n, All]] = 2 psi[[n  1, All]]/dx^2; omega[[free, free]] = Transpose[Partition[TridiagonalSolve[Drop[Flatten[Table[If[i == n  1, 0, (cxy u[[i, j]] + sxy)], {j, 2, n  1}, {i, 2, n  1}]], 1], Flatten[Table[1 + 2sxy, {j, 2, n  1}, {i, 2, n  1}]], Drop[Flatten[Table[If[i == n  1,0, cxy u[[i + 1, j]]  sxy], {j, 2,n  1}, {i, 2, n  1}]], 1], Flatten[Table[(cxy v[[i, j  1]] + sxy)omega[[i, j  1]] + (1  2sxy)omega[[i, j]] + (cxy v[[i, j + 1]] + sxy)omega[[i, j + 1]] + If[i == 2, (cxy u[[i  1, j]] + sxy)omega[[i  1, j]],0]  If[i == n  1, (cxy u[[i + 1, j]]  sxy) omega[[i + 1, j]], 0], {j, 2, n  1}, {i, 2, n  1}]]], n  2]]; omega[[free, free]] = Partition[TridiagonalSolve[Drop[Flatten[Table[If[j == n  1, 0, (cxy v[[i, j]] + sxy)], {i, 2, n  1}, {j, 2, n  1}]], 1], Flatten[Table[1 + 2sxy, {i, 2, n  1}, {j, 2, n  1}]], Drop[Flatten[Table[If[j == n  1, 0, cxy v[[i, j + 1]]  sxy], {i, 2, n  1}, {j, 2, n  1}]], 1], Flatten[Table[(cxy u[[i  1, j]] + sxy)omega[[i  1, j]] + (1  2sxy)omega[[i, j]] + (cxy u[[i + 1, j]] +sxy)omega[[i + 1, j]] + If[j == 2, (cxy v[[i, j  1]] + sxy)omega[[i, j  1]], 0]  If[j == n  1, (cxy v[[i, j + 1]]  sxy) omega[[i, j + 1]], 0], {i, 2, n  1}, {j, 2, n  1}]]], n  2]; resid = 1; While[resid > 0.001, resid = 0; Do[resid += Abs[psi[[i, j]]  (psi[[i,j]] = lambda(psi[[i, j  1]] + psi[[i  1, j]] + psi[[i + 1, j]] + psi[[i, j + 1]]  dx^2omega[[i, j]])/4 + (1  lambda)psi[[i, j]])], {i, 2, n  1}, {j, 2, n  1}]]; Do[u[[i, j]] = (psi[[i, j + 1]]  psi[[i, j  1]])/(2dx); v[[i, j]] = (psi[[i + 1, j]]  psi[[i  1, j]])/(2dx), {i, 2, n  1}, {j, 2, n  1}]; ListContourPlot[Transpose[psi], PlotRange > All, ContourShading > False], {100}];
Link: Driven Cavity  simple Fortran program by Zheming Zheng that uses this method

2D Smoothed Particle Hydrodynamics (SPH)  calculated in C++, rendered in POVRay 3.6.1, 12/10/08
The following Mathematica code for 2D SPH uses the same kernal functions as the 3D code:
(* runtime: 6.4 minutes *)
Norm[x_] := x.x; PlotColor[x_] := Hue[2(1  Min[1, Max[0, x]])/3];
dt = 0.004; g = {0, 9.8}; mu = 0.2; rho0 = 600; m = 0.00020543; r0 = 0.004; rsmooth = 0.01; rsmooth2 = rsmooth^2; vmax = 200; vmax2 = vmax^2;
WPoly6Kernal = 315/(64Pi rsmooth^9); SpikyKernal = 45/(Pi rsmooth^6); LaplacianKernal = 45/(Pi rsmooth^6); IntStiff = 0.5; ExtStiff = 20000; damp = 256;
p = 1; v = 2; v2 = 3; F = 4; P = 5; nu = 6; xmax = 0.1;
dx = 0.87(m/rho0)^(1/3); particles = Flatten[Table[{{x, y}, {0, 0}, {0, 0}, {0, 0}, 0, 0}, {y, xmax, 0,dx}, {x, 0, xmax, dx}], 1]; n = Length[particles];
Do[Do[rho = WPoly6Kernal m Sum[If[i != j, Max[0, (rsmooth2  Norm[particles[[i, p]]  particles[[j, p]]])^3], 0], {j, 1, n}]; particles[[i, P]] = (rho  rho0)IntStiff; particles[[i, nu]] = If[rho > 0, 1/rho, 0], {i, 1, n}];Do[A = particles[[i]]; particles[[i, F]] = Sum[If[i != j, B = particles[[j]]; dx = A[[p]]  B[[p]]; r2 = Norm[dx];If[r2 < rsmooth2, r = Sqrt[r2]; c = rsmooth  r; c A[[nu]] B[[nu]](0.5 c SpikyKernal(A[[P]] + B[[P]])dx/r + LaplacianKernal mu(B[[v2]]  A[[v2]])), 0], 0], {j, 1, n}], {i, 1, n}];Do[A = particles[[i]]; a = A[[F]]m; If[Norm[a] > vmax2, a = vmax a/Sqrt[Norm[a]]]; Do[d = xmax  sign A[[p, dim]]; If[d < 2r0  0.00001, normal = {0, 0}; normal[[dim]] = sign;a += (ExtStiff(2r0  d)  damp normal.A[[v2]]) normal], {sign, 1, 1, 2}, {dim, 1, 2}];a += g; particles[[i, v2]] = A[[v]] + a dt/2; particles[[i, v]] = A[[v]] + a dt; particles[[i, p]] += particles[[i, v]] dt, {i, 1, n}]; Show[Graphics[Table[{PlotColor[particles[[i, P]]/1500], Point[particles[[i, p]]]}, {i, 1,n}], AspectRatio > Automatic, PlotRange > xmax{{1, 1}, {1, 1}}]], {100}]
Links
Fluid v1.0  C++ program for fluid surfaces by Maciej Matyka, uses Marker And Cell (MAC) method
Splash  analytical Mathematica plot by Dr. Jim Swift
Particle Fluid  Java applet using a simple technique, by Chris Laurel
Lobe Compressor  SPH simulation by Simerics

3D Smoothed Particle Hydrodynamics (SPH)  calculated in C++, rendered in POVRay 3.6.1, 12/10/08
When a droplet falls into shallow water, it creates a crown or "coronet". This droplet simulation was calculated using Smoothed Particle Hydrodynamics (SPH). SPH is one of the most impressivelooking fluid simulation techniques.
Droplet Links
Liquid Sculpture  beautiful high speed photographs, by Martin Waugh, see also this video
Water Figures  beautiful highspeed camera splashes by Fotoopa
Other Links
Fluids v.3  fast SPH C++ program by Rama Hoetzlein
Physics Demos  fluid Java applets by Grant Kot
Fluid Animations  amazing animations by Ron Fedkiw, with Eran Guendelman, Andrew Selle, Frank Losasso, et al.
POVFlow  by Fidos
GPU Fluid Solver  by Keenan Crane, uses Graphics Processing Unit (GPU) for fast calculations
Free Surface Fluid Simulations  impressive simulations using the LatticeBoltzmann method with level sets, by Nils Thuerey, author of Blender’s fluid package
ball in water  nice 3D POVRay animation by Fidos
Harmonic Fluids  simulation of fluid sounds
Smoke Animations  by Jos Stam, Henrik Jensen, Ron Fedkiw, et al.
Fire Animations  by Duc Nguyen, Henrik Jensen, Ron Fedkiw
Bouncing liquid jets  University of Texas
Mark Stock  impressive fluid simulation artwork
Rune’s Particle System  nice POVRay program by Rune Johansen

Von Kármán Vortex Street  C++, 1/11/05
When fluid passes an object, it can leave a trail of vortices called a Von Kármán Vortex Street. This animation shows the vorticity where blue is clockwise and red is counterclockwise. This simulation assumes unsteady, incompressible, viscid, laminar flow. For solenoidal flows, mass conservation can be achieved by taking the Fast Fourier Transform (FFT) of the velocity and then removing the radial component of the wave number vectors. The code was adapted from Jos Stam’s C Source Code and paper which is patented by Alias. Thanks to xzczd for debugging this code. Click here to see an optimized implementation:
(* runtime: 4 minutes *)
n = 64; dt = 0.3; mu = 0.001; v = Table[{0, 0}, {n}, {n}];
ListAnimate@Table[Do[If[i < 5, v[[i, j]] = {0.1, 0}]; If[(i  n/4)^2 + (j  n/2)^2 < 4^2, v[[i, j]] = {0, 0}], {i, 1, n}, {j, 1, n}]; ui = ListInterpolation[v[[All, All, 1]]]; vi = ListInterpolation[v[[All, All, 2]]]; v = Table[{i2, j2} = {i, j}  n dt v[[i, j]]; {ui[i2, j2], vi[i2, j2]}, {i, 1, n}, {j, 1, n}]; v = Transpose[Map[Fourier[v[[All, All, #]]] &, {1, 2}], {3, 1, 2}]; v = Table[x = Mod[i + n/2  1, n]  n/2; y = Mod[j + n/2  1, n]  n/2; k = x^2 + y^2; Exp[k dt mu] If[k > 0, (v[[i, j]] . {y, x}/k) {y, x}, v[[i, j]]], {i, 1, n}, {j, 1, n}]; v = Transpose[Map[Re[InverseFourier[v[[All, All, #]]]] &, {1, 2}], {3, 1, 2}]; ListDensityPlot[Table[((v[[i + 1, j, 2]]  v[[i  1, j, 2]])  (v[[i, j + 1, 1]]  v[[i, j  1, 1]]))/2, {j, 2, n  1}, {i, 2, n  1}], Mesh > False, Frame > False, ColorFunction > (Hue[2 #/3] &), PlotRange > All], {t, 1, 350}]
Links
smoke program  uses this method, by Gustav Taxén
fire program  uses this method, Hugo Elias
Fluid Simulations  nice vortex street simulations by Maciej Matyka
Vortex street in the clouds  photo of a vortex street created by an island
Airplane Vortices  nice photo of trailing airplane vortices
Harmonic Balance CFD

Flapping Wing  Fortran 90, rendered in POVRay 3.6.1, 4/28/06
Joukowski Airfoil  C++, 11/8/07
These animations were created using a conformal mapping technique called the Joukowski Transformation. A Joukowski airfoil can be thought of as a modified Rankine oval. It assumes inviscid incompressible potential flow (irrotational). Potential flow can account for lift on the airfoil but it cannot account for drag because it does not account for the viscous boundary layer (D'Alembert's paradox). In these animations, red represents regions of low pressure. The left animation shows what the surrounding fluid looks like when the Kutta condition is applied. Notice that the fluid separates smoothly at the trailing edge of the airfoil and a low pressure region is produced on the upper surface of the wing, resulting in lift. The lift is proportional to the circulation around the airfoil. The right animation shows what the surrounding fluid looks like when there is no circulation around the airfoil (no lift). Notice the sharp singularity at the trailing edge of the airfoil. This singularity would not happen on a real airfoil. Instead the flow would separate and the airfoil would stall.

Joukowski Airfoil  Mathematica 4.2, 11/2/05
Here is an animation that shows how the streamlines change when you increase the circulation around the airfoil. (Please note: The background fluid motion in this animation is just for effect and is not accurate!) Here is some Mathematica code to plot the streamlines and pressure using Bernoulli’s equation:
(* runtime: 13 seconds *)
Clear[z]; U = rho = 1; chord = 4; thk = 0.5; alpha = Pi/9; y0 = 0.2; x0 = thk/5.2; L = chord/4; a = Sqrt[y0^2 + L^2]; gamma = 4Pi a U Sin[alpha + ArcCos[L/a]];
w[z_, sign_] := Module[{zeta = (z + sign Sqrt[z^2  4 L^2])/2}, zeta = (zeta  x0  I y0)Exp[I alpha]/Sqrt[(1  x0)^2 + y0^2]; U(zeta + a^2/zeta) + I gamma Log[zeta]/(2Pi)];
sign[z_] := Sign[Re[z]]If[Abs[Re[z]] < chord/2 && 0 < Im[z] < 2y0(1  (2Re[z]/chord)^2), 1, 1];w[z_] := w[z, sign[z]]; V[z_] = D[w[z, sign], z] /. sign > sign[z];
ContourPlot[Im[w[(x + I y)Exp[I alpha]]], {x, 3, 3}, {y, 3, 3}, Contours > Table[x^3 + 0.0208, {x, 2, 2, 0.1}], PlotPoints > 100, ContourShading > False, AspectRatio > Automatic]
DensityPlot[0.5rho Abs[V[(x + I y)Exp[I alpha]]]^2, {x, 3, 3}, {y, 3, 3}, PlotPoints > 275, Mesh > False, Frame > False, ColorFunction > (If[# == 1, Hue[0, 0, 0], Hue[(5#  1)/6]] &), AspectRatio > Automatic]
Links
Joukowski Animation  nice animation showing how the fluid moves
Joukowski Airfoil  nice Java applet by NASA
Nikolai Joukowski  used this technique to find the lift on an airfoil in 1906, long before modern computers

NACA Airfoil  Mathematica 4.2, 12/16/05
Here is a 9415 NACA airfoil. It looks similar to the Joukowski airfoil, but the shape is slightly better. Here is some Mathematica code:
(* runtime: 0.01 second *)
Clear[x]; n = 40; c = 0.09; p = 0.4; t = 0.15;
camber := c If[x < p, (2p x  x^2)/p^2, ((1  2p) + 2p x  x^2)/(1  p)^2]; theta = ArcTan[D[camber, x]];
p = Table[x = 0.5(1  Cos[Pi s]); x1 = 1.00893x; thk = 5t(0.2969Sqrt[x1]  0.126x1  0.3516x1^2 + 0.2843x1^3  0.1015x1^4); {x, camber} + Sign[s]thk {Sin[theta],Cos[theta]}, {s, 1, 1, 2/(n  1)}];
ListPlot[p, PlotJoined > True, AspectRatio > Automatic]
We can approximate the pressure profile using the vortex panel method assuming inviscid incompressible potential flow (irrotational). The following Mathematica code was adapted from my Fortran program for my 4/25/99 research project on optimal airfoil design:
(* runtime: 0.25 second *)
alpha = Pi/9; pc = Table[(p[[i]] + p[[i + 1]])/2, {i, 1, n  1}]; s = Table[v = p[[i + 1]]  p[[i]];Sqrt[v.v], {i, 1, n  1}]; theta = Table[v = p[[i + 1]]  p[[i]]; ArcTan[v[[1]], v[[2]]], {i, 1, n  1}]; sin = Sin[theta]; cos = Cos[theta]; Cn1 = Cn2 = Ct1 = Ct2 =Table[0, {n  1}, {n  1}];
Do[If[i == j, Cn1[[i, j]] = 1; Cn2[[i, j]] = 1; Ct1[[i, j]] = Ct2[[i, j]] = Pi/2,v = pc[[i]]  p[[j]]; a = v.{cos[[j]], sin[[j]]}; b = v.v;t = theta[[i]]  theta[[j]];c = Sin[t]; d = Cos[t]; e = v.{sin[[j]], cos[[j]]}; f =Log[1 + s[[j]](s[[j]] + 2a)/b]; g = ArcTan[b + a s[[j]], e s[[j]]];t = theta[[i]]  2theta[[j]]; p1 = v.{Sin[t], Cos[t]}; q = v.{Cos[t], Sin[t]}; Cn2[[i, j]] = d + (0.5q f  (a c + d e)g)/s[[j]]; Cn1[[i, j]] = 0.5d f + c g  Cn2[[i, j]];Ct2[[i, j]] = c + 0.5p1 f/s[[j]] + (a d  c e)g/s[[j]]; Ct1[[i, j]] = 0.5c f  d g  Ct2[[i, j]]], {i, 1, n  1}, {j, 1, n  1}];
gamma = LinearSolve[Table[If[i == n, If[j == 1  j == n, 1, 0], If[j == n, 0, Cn1[[i, j]]] + If[j == 1, 0, Cn2[[i, j  1]]]], {i, 1, n}, {j, 1, n}], Table[If[i ==n, 0, Sin[theta[[i]]  alpha]], {i, 1, n}]];
ListPlot[Table[q = Cos[theta[[i]]  alpha] + Sum[(If[j == n, 0, Ct1[[i, j]]] + If[j == 1, 0, Ct2[[i, j  1]]])gamma[[j]], {j, 1, n}]; {pc[[i, 1]], 1  q^2}, {i, 1, n  1}], PlotJoined > True]
Here is some Mathematica code for a pretty pressure plot:
(* runtime: 33 seconds *)
w[z_] := z Exp[I alpha] + I Sum[s[[j]]gamma[[j]]Log[z  pc[[j, 1]]  I pc[[j, 2]]], {j, 1, n  1}]; V[z_] = D[w[z], z];
DensityPlot[Abs[V[(x + I y)Exp[I alpha]]]^2, {x, 0.25, 1.25}, {y, 0.75, 0.75}, PlotPoints > 275, Mesh > False, Frame > False, ColorFunction > (Hue[(5#  1)/6] &), AspectRatio > Automatic]
Links
Panel Code  by Kevin Jones
Vortex Panel Java applet  by William Devenport

KelvinHelmholtz Instability Waves  Mathematica 4.2, C++, POVRay 3.6.1, 11/16/06
KelvinHelmholtz instabilities are formed when one fluid layer is on top of another fluid layer moving with a different velocity. These instabilities take the form of small waves that can eventually grow into vortex rollers. This is a purely potential flow (irrotational) phenomenon. The following Mathematica code was adapted from Zheming Zheng’s Fortran program. It uses vortex blobs to simulate smooth rollers and periodic point vortices to simulate the far field. A simple predictorcorrector scheme is used to integrate the differential equations:
(* runtime: 9 minutes *)
n = 40; dt = 0.05; L = 2Pi; gamma = L/n; rcore = 0.5; plist = Table[{x, 0.01 Sin[2Pi x/L]}, {x, 0, L(1  1/n), L/n}];
vcalc[plist_] := Table[Sum[{x, y} = plist[[j]]  plist[[i]]; gamma/(2Pi) Sum[If[i == j && k == 0, 0, r2 = (x + k L)^2 +y^2; {y, x + k L}/r2(1  Exp[r2/rcore^2]  1)], {k, 3, 3}] + If[i ==j, 0, gamma/(2L) {Sinh[2Pi y/L], Sin[2Pi x/L]}/(Cosh[2Pi y/L]  Cos[2Pi x/L])], {j, 1, n}], {i, 1, n}];
Do[vlist = vcalc[plist]; vlist0 =vlist; vlist = vcalc[plist + dt vlist]; plist += dt(vlist0 + vlist)/2; ListPlot[plist, PlotJoined > True, PlotRange > {2, 2}, AspectRatio > Automatic], {i, 1, 400}];

Here is some Mathematica code to plot streamlines assuming periodic point vortices:
(* runtime: 2 seconds *)
Clear[psi]; psi[x_, y_] := Plus @@ Map[gamma/(4 Pi) Re[Log[Cos[2Pi(x  #[[1]])/L]  Cosh[2Pi (y  #[[2]])/L]]] &, plist];
ContourPlot[psi[x, y], {x, 0, L}, {y, L/3, L/3}, PlotRange > All, ContourShading > False, PlotPoints > 20, AspectRatio > Automatic]
Link: KelvinHelmholtz waves in the clouds  Birmingham, Alabama

KelvinHelmholtz Instability Waves  Mathematica 4.2, C++, 11/16/06
Magnus Effect  C++, Mathematica 4.2, POVRay 3.6.1, 10/24/07
A Curveball is a spinning ball that accelerates sideways due to the Magnus Effect. Here is some Mathematica code to plot the streamlines around a lifting cylinder, assuming inviscid incompressible potential flow (irrotational):
(* runtime: 0.3 second *)
U = 1; gamma = 4; a = 1; w[z_] := U(z + a^2/z) + I gamma Log[z]/(2Pi);
ContourPlot[Im[w[x + I y]], {x, 2, 2}, {y, 2, 2}, Contours > Table[psi, {psi, 3, 3, 0.25}], PlotPoints > 100, ContourShading > False, AspectRatio > Automatic]
Here is some Mathematica code to plot the entrained fluid using the 4th order RungeKutta method:
(* runtime: 17 seconds, increase n for better resolution *)
Clear[v]; n = 40; U = 1; gamma = 4; a = 1; tmax = 5.0; dt = tmax/100; v[z_] := If[Abs[z] < 1, 0, Conjugate[U (1  (a/z)^2) + I gamma/(2Pi z)]];
image = Table[z = x + I y; Do[k1 = dt v[z]; k2 = dt v[z + k1/2]; k3 = dt v[z + k2/2]; k4 =dt v[z + k3]; z = (k1 + 2 k2 + 2 k3 + k4)/6, {t, tmax, 0, dt}]; z, {y, 2, 2, 4.0/n}, {x, 2, 2, 4.0/n}];
ListDensityPlot[Map[Abs[Mod[#, 1]] &, image, {2}], Mesh > False, Frame > False, AspectRatio > 1]
Link: My Educational Project  some basic potential flows using Matlab and Mathematica

1D Shock Tube  Fortran 90, Mathematica 4.2, POVRay 3.6.1, 3/15/07
A shock tube is a tube containing high and low pressure gas separated by a thin diaphragm. A shock wave is produced when the diaphragm is quickly removed. The color in the upper plot shows the pressure. The lower plot shows the density. The following Mathematica code solves Euler’s equations in 1D using the finite volume method with the JamesonSchmidtTurkel (JST) scheme and RungeKutta time stepping.
(* runtime: 5 seconds *)
gamma = 1.4;
R[W_] := Module[{}, rho = W[[All, 1]]; u = W[[All, 2]]/rho; p = (gamma  1)(W[[All, 3]]  rho u^2/2); F = u W + Transpose[{Table[0, {n}], p, u p}]; h = Table[(F[[Min[n, i + 1]]] + F[[i]])/2, {i,1, n}]; Q = Table[h[[Max[i, 2]]]  h[[Max[i, 2]  1]], {i, 1, n}]; nu = Table[i = Max[2, Min[n  1, i]]; Abs[(p[[i + 1]]  2p[[i]] + p[[i  1]])/(p[[i + 1]] + 2p[[i]] + p[[i  1]])], {i, 1, n}]; S = Table[Max[nu[[Min[n, i + 1]]], nu[[i]]], {i, 1, n}]; alpha1 = 1/2; beta1 = 1/4;alpha2 = alpha1; beta2 = beta1; epsilon2 = Map[Min[alpha1, alpha2#] &, S];epsilon4 = Map[Max[0, beta1  beta2#] &, epsilon2];dW = Table[W[[Min[n  1, i] + 1]]  W[[Min[n  1, i]]], {i, 1, n}];dW3 = Table[i = Max[2, Min[n  2, i]]; W[[i  1]] + 3W[[i]]  3W[[i + 1]] + W[[i + 2]], {i, 1, n}];d = (epsilon2 dW  epsilon4 dW3)(Abs[u] + a); Dflux = Table[d[[Max[2, i]]]  d[[Max[2, i]  1]], {i, 1, n}]; (Q  Dflux)/dx];
n = 50; dx = 1.0/n; a = 1.0; dt = dx/(1.0 + a); W = Table[{If[i > n/2, 0.125, 1], 0, If[i > n/2, 0.1, 1]/(gamma  1)}, {i, 1, n}];
Do[W = dt R[W  dt R[W  dt R[W]/4]/3]/2;ListPlot[W[[All, 1]], PlotJoined > True,PlotRange > {0, 1}, AxesLabel > {"i", "rho"}], {t, 0, 100dt, dt}];

Jet  C++, 6/6/07
This simulation is based on a simple method presented in this paper and sample C++ code by Jos Stam.
(* runtime: * seconds *)

SchwarzChristoffel Flow  C++, 11/*/07
Here is a potential flow field of a vortex near a flat wall mapped to stair steps using a SchwarzChristoffel transformation. I got the inspiration for this stair step shape from Matthias Weber's website. Here is some Mathematica code:
(* runtime: 0.7 second *)
Clear[z]; gamma = 1; z0 = 1 + I; z1 = I;z2 = 0; z3 = 1; zeta1 = 1; zeta2 = 0; zeta3 = 1; theta1 = 3Pi/2; theta2 = Pi/2; theta3 = 3Pi/2;
z = z[zeta] /. DSolve[{z'[zeta] == K(zeta  zeta1)^(theta1/Pi  1)(zeta  zeta2)^(theta2/Pi  1)(zeta zeta3)^(theta3/Pi  1), z[zeta1] == z1}, z[zeta], zeta][[1]];
z = z /. K > (K /. Solve[(z /. zeta > zeta2) == z2, K][[1]]); zeta0 = zeta /. FindRoot[z == z0, {zeta, zeta2 + I}];
Z[w_] = z /. zeta > (zeta /. Solve[w == (I gamma/(2Pi))(Log[zeta  zeta0]  Log[zeta  Conjugate[zeta0]]), zeta][[1]]);
ToPoint[z_] := {Re[z], Im[z]}; grid = Table[ToPoint[Z[phi + I psi]], {phi, 1, 0, 0.05}, {psi, 1*^6, 1, 0.02}];
Show[Graphics[{Map[Line, grid], Map[Line, Transpose[grid]]}, AspectRatio > Automatic, PlotRange > {{1, 2}, {1, 2}}]]
Links
My Educational Project  some basic potential flows using Matlab and Mathematica
SchwarzChristoffel Maps  nice plots by Matthias Weber used for minimal surfaces
SchwarzChristoffel Transformations  many examples by John Mathews

Periodic Steady Vortices in a Stagnation Point Flow  C++, 11/17/07
Vortices  Mathematica 4.2, C++, 10/31/06
Here is another code to solve Euler’s equations (inviscid flow) and plot the pathlines. This method was adapted from Stephen MontgomerySmith’s Euler2D program. The vortex effects are found using the BiotSavart law and the differential equations are solved using the Adams Bashforth method. Click here to download some Matlab code for this. Shown below is some Mathematica code:
(* runtime: 19 seconds *)
ToPoint[z_] := {Re[z], Im[z]};dt = 0.01; rcore = 0.1; dz0 = 0; Klist = {1, 1, 1, 1};
zlist = Join[{1  0.5I, 1 + 0.5I, 0.5  0.5I, 0.5 + 0.5I}, Flatten[Table[x + I y, {x, 1, 1, 0.1}, {y, 1, 1, 0.1}], 1]];
paths = Map[# Table[1, {10}] &, zlist];
Do[dz = Map[Sum[zj = zlist[[j]]; If[# != zj, r2 = Abs[#  zj]^2; I Klist[[j]](zj  #)/r2 (1 Exp[r2/rcore^2]), 0], {j, 1, Length[Klist]}] &, zlist]; zlist += dt (1.5dz  0.5dz0); dz0 = dz; paths = Table[Prepend[Delete[paths[[i]], 1], zlist[[i]]], {i, 1, Length[zlist]}];Show[Graphics[Map[Line[Map[ToPoint, #]] &, paths], PlotRange > {{1, 1}, {1, 1}}, AspectRatio > Automatic]], {85}];
Link: math paper by Stephen MontgomerySmith, seems fairly advanced to me

Ferrofluid  4/6/05
Whirlpools  3/16/08, 3/31/12
Smoke Ring  10/13/07
Supersonic Flow  CFL3D (Computational Fluids Laboratory 3D), Tecplot 360, 2/1/07
Here are some weird tests of CFL3D. The right image was my attempt to simulate shock diamonds. The color represents the Mach number.

Channel Flow  Fluent 5.5.14, Gambit, Mathematica 4.2, 5/14/05
Fluid Motion Links
Gallery of Fluid Dynamics  excellent web site on fluid dynamics by Mark Cramer
Gallery of Fluid Motion  from the American Institute of Physics
Best Fluid Motion Pictures Named  National Geographic article
Oobleck  amazing nonNewtonian fluid that behaves like a solid when you move fast
Hydraulophone video  a hydraulophone works just like a flute except it uses water instead of air, the world's largest hydraulophone is at the Ontario Science Centre in Canada
