Fluid Flow Simulator

Put a finger in a slow stream and the water slides around it and closes up behind, almost as if nothing were there. Speed the stream up and the picture changes. The water no longer closes up. Eddies form behind the finger, let go, and drift away, first from one side and then from the other, in a rhythm steady enough to feel as a tremor in the finger.

What decides which of these flows you get is not the speed alone, nor the size of the finger, nor the thickness of the liquid. It is a single combination of all three. That combination is the Reynolds number, and much of what is beautiful about the wake behind an obstacle can be read off from it.

One ratio for every stream

The Reynolds number is the flow speed times a length, divided by the kinematic viscosity of the fluid. It compares two tendencies. Inertia keeps a moving parcel of fluid going; viscosity drags neighbouring parcels toward a common speed and smooths differences out. At low values the smoothing wins and the flow is calm and laminar. At high values inertia wins, and the flow breaks into eddies and, eventually, turbulence. George Stokes introduced the quantity in 1851, and Arnold Sommerfeld named it in 1908 after Osborne Reynolds, whose 1883 experiments with dye in a glass pipe made it famous.

The power of the ratio is that it carries no units. A wire in a breeze and a bridge pier in a river behave alike when their Reynolds numbers match, even though one is a millimetre across and the other several metres. Engineers lean on this when they test a small model in a wind tunnel, and a simulation leans on it even harder, because a simulation has no metres at all.

For a circular cylinder the length is the diameter, and the cylinder’s wake is the classic case. Below a Reynolds number of about 47 the flow behind it is steady: two small eddies sit side by side in its lee and stay put. Above that, the pair can no longer hold still. One eddy grows, breaks away and is carried off; then its partner on the other side does the same. The result is the double row of alternating swirls known as a Karman vortex street.

Counting particles instead of solving equations

The usual way to compute a flow is to write down the Navier-Stokes equations and solve them on a grid. The lattice Boltzmann method takes a stranger route. It never solves those equations directly. Instead it tracks how much fluid is moving in each of a handful of fixed directions at every point of a grid, and lets those amounts hop and collide.

The idea grew out of lattice gas automata. In papers from 1973 and 1976, Jean Hardy, Yves Pomeau and Olivier de Pazzis moved true-or-false particles around a square grid; in 1986 Uriel Frisch, Brosl Hasslacher and Pomeau switched to a hexagonal grid, which behaved much more like a real fluid. Those models were noisy, because each site either held a particle or did not. The lattice Boltzmann method replaced each yes-or-no particle with its average, a smooth density for every direction, and the noise went away.

In two dimensions the common lattice is called D2Q9: nine populations per site, one at rest, four moving along the axes and four along the diagonals. Each step has two parts. In streaming, every population moves one site in its own direction. In collision, the nine populations at a site relax part of the way toward an equilibrium set by the local density and velocity. The equilibrium weights are 4/9 for the population at rest, 1/9 for each axis direction and 1/36 for each diagonal. The simplest collision rule, known after Bhatnagar, Gross and Krook as BGK, relaxes by a fixed fraction each step, set by a relaxation time tau.

From such bookkeeping, a careful expansion recovers the equations of a viscous fluid, and the viscosity falls out as (tau minus one half) divided by three in lattice units. That relation is the hinge of the whole method. Choose a speed and a length in cells, choose the Reynolds number you want, and tau follows. The same expansion also warns that the equilibrium is only right for flows well below the lattice speed of sound, which is why lattice speeds are kept small.

Testing a toy against exact answers

A simulation earns trust by matching cases whose answers are known. Two are standard. Flow pushed by a pressure difference between two flat walls settles into a parabola, plane Poiseuille flow. A grid of counter-rotating vortices in a periodic box, the Taylor-Green vortex, keeps its shape while its speed decays exactly as the exponential of minus twice the viscosity times the wavenumber squared times time. The solver behind this page matches the parabola to within 0.06 percent of the peak speed on a 32-cell channel, and the decay rate to within 0.03 percent; a periodic box loses less than a trillionth of its mass in ten thousand steps.

The weights matter more than they look. Change the diagonal weight from 1/36 to 1/9 and the populations at a site no longer sum to the right density. In the channel test that single change does not produce a slightly wrong parabola. The run blows up.

Boundaries need equal care. At the obstacle, a population that would enter a solid cell is sent straight back the way it came, a rule called bounce-back, which puts the wall halfway between cells. Counting the momentum those reflected populations carry gives the force on the obstacle, the drag along the stream and the lift across it. At the inlet, Qisu Zou and Xiaoyi He showed how to rebuild the unknown populations so that a chosen velocity is met exactly, by assuming that the non-equilibrium part bounces back too.

The method has a fragile side. James Sterling and Shiyi Chen found that runs become unstable when tau creeps too close to one half, which is exactly where the Reynolds number gets large. Here a margin of 0.54 is kept, chosen from test runs in which 0.52 and 0.53 blew up within a few thousand steps. When a requested Reynolds number would push tau below it, the grid has to get finer, or the request is refused.

Reading the rhythm of a wake

The rate at which a wake sheds vortices has its own dimensionless number. Multiply the shedding frequency by the diameter and divide by the flow speed, and you get the Strouhal number, named for Vincenc Strouhal, who studied wires singing in the wind in 1878. Those wires sing because each shed vortex gives the wire a small sideways push, and the pushes alternate. The same alternating force makes power lines hum and can shake tall buildings and masts when it matches their natural frequency.

To measure it, the solver records the lift on the obstacle every ten steps, drops the start-up, and runs a fast Fourier transform, an algorithm usually credited to James Cooley and John Tukey in 1965. The tallest peak in the spectrum is the shedding frequency. A peak-fitting step refines it between the frequency bins, and counting zero crossings of the lift gives an independent check.

Here is a worked example you can repeat. Set the Reynolds number to 100 and choose the wide grid of 640 by 320 cells, with a cylinder 16 cells across, so it blocks 5 percent of the channel. The inlet speed is 0.1 lattice units, so the viscosity must be 0.1 times 16 divided by 100, which is 0.016, and tau is 0.548.

After 25,000 steps the spectrum peaks at a Strouhal number of 0.166, the zero crossings agree, and the mean drag coefficient is 1.46. Four published simulations gathered by Cui and colleagues give values from 0.163 to 0.17 at this Reynolds number, with drag from 1.33 to 1.46. Another group’s grid study of its own solver gives 0.161 to 0.165.

Now take the standard grid, where the same cylinder fills 10 percent of a narrower channel. The Strouhal number rises to 0.175 and the drag to 1.57. Nothing has gone wrong. The walls squeeze the flow past the cylinder, it speeds up, and the wake beats faster. A channel is not open water, and the blockage belongs in any honest comparison.

The ratio also travels. Water at 25 degrees Celsius has a density of about 997 kilograms per cubic metre and a viscosity of about 0.89 millipascal seconds. A wire 2 millimetres across in a current of 4.5 centimetres per second therefore sits at a Reynolds number of about 101. With a Strouhal number near 0.166, it should shed about 3.7 vortex pairs every second.

What a flat channel cannot show

The biggest limit is the second dimension. A real cylinder has length, and its wake is free to twist along it. Simulations of the cylinder wake find a first three-dimensional instability, called mode A, near a Reynolds number of 190, and a second, mode B, near 240. Past those values a two-dimensional model still draws a tidy vortex street, but the real wake no longer looks like one, and its forces and shedding rate drift away from the flat answer.

There are smaller limits too. Obstacles are built from square cells, so a cylinder is a staircase at its edge, and its effective size depends on the grid. The inlet and outlet sit a few diameters away, and the side walls are close enough to matter. The flow is slightly compressible, because the method works at a modest Mach number on purpose. None of this is hidden: the blockage, tau and Mach number are shown beside every run.

Within those limits, the lattice does something remarkable. Nine numbers per site, streamed and relaxed with simple arithmetic, produce eddies that let go at a rate within a few percent of careful published work. The street that Theodore von Karman described more than a century ago, and that Arnulph Mallock had studied before him, appears from nothing more than counting.