PL. 115 · FIELDS / ODE / POINT VORTICES
Lamb–Oseen Vortices
Hermann von Helmholtz, 1858 & Gustav Kirchhoff, 1876 · Carl Wilhelm Oseen, 1912 & Horace Lamb, 1932
OPEN THE LIVE PLATE ▸DEFINITION
u(x) = Σⱼ (Γⱼ/2π) · ẑ × (x − xⱼ)/|x − xⱼ|² · [1 − exp(−|x − xⱼ|²/a²)] a²(t) = a₀² + 4νt (the Lamb–Oseen core, spreading by diffusion) ẋᵢ = u(xᵢ) with the i-th term omitted — a vortex does not move itself
NOTES
Helmholtz showed in 1858 that a vortex in an ideal fluid is permanent: it travels with the flow and carries its circulation with it forever, so a handful of them becomes a closed system of ordinary differential equations, each one swept along by the field the others make. Kirchhoff wrote that system down in Hamiltonian form. The flaw is the word point — the velocity a point vortex induces runs to infinity at its own centre, so two that pass close enough exchange an unbounded kick. Oseen removed the singularity in 1912 by solving the viscous problem instead: released into a real fluid, a point vortex diffuses, its vorticity spreading as a Gaussian whose width grows as the square root of four times the viscosity times the elapsed time, and the velocity it induces now climbs smoothly from zero at the centre to the ideal circulation over two pi r outside the core. Lamb gave the same solution in Hydrodynamics, which is why it carries both names. That is the field these tracers are dropped into — they have no circulation of their own and only report what the vortices do. Opposite-signed pairs translate in a straight line, same-signed pairs orbit each other, and three or more are chaotic. The rings are the cores at their present radius and they only ever grow; when one passes about a quarter of the way to its nearest neighbour, which is the ratio at which two real co-rotating vortices of this kind are observed to merge, the point-vortex idealisation has stopped describing anything and the plate says so and starts over. Plates 23 and 64 read their currents out of a noise field; this one integrates a field whose sources are themselves moving.
PROVENANCE
- Origin — the dynamics
- H. Helmholtz, "Über Integrale der hydrodynamischen Gleichungen, welche den Wirbelbewegungen entsprechen", Journal für die reine und angewandte Mathematik (Crelle) 55, 1858, 25–55. The N-vortex system as a Hamiltonian one is G. Kirchhoff, Vorlesungen über mathematische Physik: Mechanik, Teubner, Leipzig, 1876
- Origin — the viscous core
- C. W. Oseen, "Über die Wirbelbewegung in einer reibenden Flüssigkeit", Arkiv för matematik, astronomi och fysik 7, 1912. Secondary sources disagree about the page range — 14–21 and 14–26 both circulate — so none is given here rather than one that might be wrong. The same solution is Art. 334a of H. Lamb, Hydrodynamics, 6th edition, Cambridge University Press, 1932
- Standing
- Public domain — a closed-form solution of the Navier–Stokes equations published in 1912, and an ODE system published in 1858. No patent is possible on any part of it
- No source link
- Oseen 1912 has no DOI, which is expected of a 1912 Swedish academy series. The Crossref record for the Helmholtz paper matches on title, journal, volume, issue, pages and year but carries no author field at all, and this atlas requires all four of title, author, year and venue before it will print a link. A citation that is nearly right is worse here than none
- Where the life cycle number comes from
- The core-to-separation ratio at which two co-rotating Gaussian vortices are observed to merge is 0.22–0.24, measured experimentally and numerically: P. Meunier, S. Le Dizès & T. Leweke, "Physics of vortex merging", Comptes Rendus Physique 6, 2005, 431–450 (doi:10.1016/j.crhy.2005.06.003). This model cannot merge — point vortices with regularised cores never do — so reaching that ratio is the plate announcing that it has run past its own validity, not depicting a merger
- Checked, not just plotted
- Both two-vortex behaviours have closed-form solutions, so they were driven through this plate own integrator and compared with them rather than trusted. An opposite-signed pair at separation 0.24 translates in a straight line at 0.145892 frame-widths per second against a predicted Γ/2πd times the core factor of 0.145892, with zero lateral drift and zero change of separation; a same-signed pair holds its centroid exactly and turns 11.609719 revolutions in sixty seconds against a predicted 11.609719, a period error of six parts in a billion. The core factor was checked separately at four core ratios: at a separation of half a core radius the pair moves 4.52 times slower than an ideal point vortex would carry it, which is the 1912 solution doing the thing it was written to do
- Constants
- Circulation sets how fast everything moves, viscosity how fast the cores spread and therefore how long an epoch lasts, and the initial core radius how tight the swirls start. Vortex count and tracer count are locked: both are cost, and the configuration regenerate would otherwise reach for is already derived from the seed
HOUDINI · VEX
The same published mathematics as a Detail Wrangle body. Paste it into a Wrangle with Run Over set to Detail; every constant is the published value plus a tweak channel, so Create Spare Parameters gives a slider that starts where the paper does.
// FORMA — PL. 115 · LAMB–OSEEN VORTICES — Hermann von Helmholtz, 1858 & Gustav Kirchhoff, 1876 · Carl Wilhelm Oseen, 1912 & Horace Lamb, 1932
// u(x) = Σⱼ (Γⱼ/2π) · ẑ × (x − xⱼ)/|x − xⱼ|² · [1 − exp(−|x − xⱼ|²/a²)]
// a²(t) = a₀² + 4νt (the Lamb–Oseen core, spreading by diffusion)
// ẋᵢ = u(xᵢ) with the i-th term omitted — a vortex does not move itself
// Paste into a Detail Wrangle (Run Over: Detail), no inputs needed.
// Written from the published mathematics, not adapted from any code.
// Constants arrive at their published values. Press the node's Create
// Spare Parameters button and every tweak becomes a slider — starting
// at 0, the published figure, and moving in the constant's own units.
// https://forma-gen.com/#plate=lamboseen
float p_vortices = 4 + chf('vortices_tweak'); // vortices · live 2 .. 8
float p_circ = 0.22 + chf('circ_tweak'); // circulation Γ · live 0.08 .. 0.6
float p_visc = 7 + chf('visc_tweak'); // viscosity ν × 10⁻⁵ · live 2 .. 10
float p_core = 0.022 + chf('core_tweak'); // initial core radius · live 0.01 .. 0.034
float p_tracers = 900 + chf('tracers_tweak'); // tracer particles · live 300 .. 2000
// The plate's own colour: FORMA's FIELDS accent as a cosine ramp,
// brightest near t = 0 and t = 1, near-black around t = 0.5.
vector forma_ramp(float t){
return set(
0.46 + 0.5 * cos(6.28318530718 * (t + 0)),
0.3031 + 0.3294 * cos(6.28318530718 * (t + 0.05)),
0.11 + 0.1196 * cos(6.28318530718 * (t + 0.1)));
}
// N point vortices advecting one another through the Biot-Savart velocity
// each of them induces, with Oseen's viscous core on every term so nothing is
// singular: v_theta = G/(2 pi r) * (1 - exp(-r^2 / a^2)), and a^2 = a0^2 + 4 nu t,
// the Gaussian spreading by diffusion exactly as the 1912 solution says. RK4
// on the vortices, midpoint on the tracers — the plate's own split, for the
// plate's own reason: the vortex system is the specimen and there are at most
// eight of them, while the tracers are witnesses to the field and are where
// all the cost is. Opposite-signed pairs translate, same-signed pairs
// co-rotate. Each tracer emits its recent path as a polyline; each vortex
// emits its core at the final radius as a closed ring, white for positive
// circulation and the order accent for negative, which is the plate's own
// HILITE and TRACE. Deterministic: every draw is random(counted seed).
//
// THREE THINGS ARE DELIBERATELY NOT THE PAGE, and all three are recorded
// here rather than left to be discovered:
//
// 1. Y IS NOT NEGATED. Almost every 2-D port here flips y because canvas y
// runs down and Houdini y runs up. This plate is the exception: a vortex
// system is chiral, so the plate keeps its own domain y-UP and negates
// once, at the blit. The domain therefore already matches Houdini and
// passing it straight through is what makes the picture agree — flipping
// again would mirror it and reverse every circulation.
// 2. NO REJECTION LOOP. The plate redraws its layout until no two vortices
// start closer than SEP_LO, because it runs for thousands of frames and
// a crowded start ends its life cycle early. A cook runs a fixed number
// of steps and stops, so there is no life cycle to protect and the first
// layout is taken.
// 3. THE HASH. mulberry32 is exact-int32 bit mixing and VEX has no shifts,
// so the seeded stream here is VEX's own random(). The configurations
// differ from the page's in which pairs land where, not in what they are
// or how they behave.
//
// Helpers carry this port's own prefix so it can share a Houdini scene with
// flowfield and curl, whose fields it will be compared against.
float SCALE = 560.0; // the world size flowfield and curl use
float DT = 1.0 / 60.0; // the plate's own step
int forma_warm = 300; // steps run before any trail is recorded
int forma_rec = 130; // steps recorded — the span the plate's own
// fading buffer still shows at alpha 0.19/6
float HX = 0.5, HY = 0.5; // half-extents of a square frame, in the
// plate's units of one frame side
// The velocity every vortex feels at the given positions, which is what one
// RK4 stage needs: all of them move together. A vortex is omitted from its
// own sum, because a circularly symmetric vortex induces nothing on itself.
function float[] forma_lo_stage(float px[]; float py[]; float gg[]; int nv; float rc2){
float out[];
resize(out, nv * 2);
for (int i = 0; i < nv; i++){
float ux = 0.0, uy = 0.0;
for (int j = 0; j < nv; j++){
if (j == i) continue;
float dx = px[i] - px[j], dy = py[i] - py[j];
float r2 = dx * dx + dy * dy;
if (r2 < 1e-12) continue;
float f = gg[j] * (1.0 - exp(-r2 / rc2)) / (r2 * 6.283185307);
ux -= dy * f; // z_hat x (dx, dy) = (-dy, dx)
uy += dx * f;
}
out[i * 2] = ux; out[i * 2 + 1] = uy;
}
return out;
}
// The same field at one arbitrary point — a tracer, which carries no
// circulation of its own and only reports what the vortices do.
function vector forma_lo_vel(float qx, qy; float px[]; float py[]; float gg[]; int nv; float rc2){
float ux = 0.0, uy = 0.0;
for (int j = 0; j < nv; j++){
float dx = qx - px[j], dy = qy - py[j];
float r2 = dx * dx + dy * dy;
if (r2 < 1e-12) continue;
float f = gg[j] * (1.0 - exp(-r2 / rc2)) / (r2 * 6.283185307);
ux -= dy * f;
uy += dx * f;
}
return set(ux, uy, 0.0);
}
int NV = int(rint(p_vortices));
int NT = int(rint(p_tracers));
int seed = 1912; // the plate's own seed, counted upward
// The layout: vortices in pairs, because that is where the two textbook
// behaviours live. Pair centres go on a ring wide enough that neighbouring
// pairs stay about one separation apart while the outermost vortex stays in
// frame; scattering them independently lets two pairs overlap.
int pairs = (NV + 1) / 2;
float SEP_LO = 0.19, SEP_HI = 0.27;
float ring = 0.0;
if (pairs > 1) ring = min(0.44 - SEP_HI * 0.5, 0.95 * SEP_HI / sin(M_PI / float(pairs)));
float px[], py[], gg[];
resize(px, NV); resize(py, NV); resize(gg, NV);
for (int i = 0; i < NV; i += 2){
int k = i / 2;
float ang = (float(k) + random(seed) * 0.55) / float(pairs) * 2.0 * M_PI; seed++;
float cx = cos(ang) * ring, cy = sin(ang) * ring;
if (i + 1 < NV){
float sep = SEP_LO + random(seed) * (SEP_HI - SEP_LO); seed++;
float g = p_circ * (0.7 + random(seed) * 0.6); seed++;
// about a third of pairs are opposite-signed: a dipole is the more
// striking behaviour and also the shorter-lived, since it propagates
int opposite = 0;
if (random(seed) < 0.35) opposite = 1;
seed++;
// a dipole flies perpendicular to its own axis, so choosing the axis
// chooses the flight; it is aimed thirty to sixty degrees off the
// frame centre, because aiming exactly inward fires it at its
// neighbour and aiming at random fires half of them off the edge
float th;
if (opposite == 1){
float sgn = 1.0;
if (random(seed) < 0.5) sgn = -1.0;
seed++;
float off = sgn * (0.55 + random(seed) * 0.5); seed++;
th = atan2(-cy, -cx) + off + M_PI * 0.5;
} else {
th = random(seed) * 2.0 * M_PI; seed++;
}
px[i] = cx + cos(th) * sep * 0.5;
py[i] = cy + sin(th) * sep * 0.5;
px[i + 1] = cx - cos(th) * sep * 0.5;
py[i + 1] = cy - sin(th) * sep * 0.5;
gg[i] = g;
gg[i + 1] = g;
if (opposite == 1) gg[i + 1] = -g;
} else {
// an odd count leaves one vortex on its own, orbiting the rest
px[i] = cx; py[i] = cy;
float sgn = 1.0;
if (random(seed) < 0.5) sgn = -1.0;
seed++;
gg[i] = p_circ * sgn;
}
}
float tx[], ty[], hue[];
resize(tx, NT); resize(ty, NT); resize(hue, NT);
for (int i = 0; i < NT; i++){
tx[i] = (random(seed) * 2.0 - 1.0) * HX; seed++;
ty[i] = (random(seed) * 2.0 - 1.0) * HY; seed++;
hue[i] = random(seed); seed++;
}
// One polyline per tracer, opened once recording starts and reopened wherever
// a tracer leaves the frame and respawns. The first point of a run is HELD
// rather than committed: a tracer that respawns on the very step after it
// opened would otherwise leave a one-vertex polyline behind, which is a prim
// that draws nothing. A single point is not a path.
int prim[]; resize(prim, NT);
int open[]; resize(open, NT);
float heldx[], heldy[]; resize(heldx, NT); resize(heldy, NT);
int held[]; resize(held, NT);
for (int i = 0; i < NT; i++){ prim[i] = -1; open[i] = 0; held[i] = 0; }
float tau = 0.0;
float acore = p_core;
float ax[], ay[];
resize(ax, NV); resize(ay, NV);
for (int s = 0; s < forma_warm + forma_rec; s++){
// the core spreads first, so the field everything sees this step is the
// field at this time: a^2 = a0^2 + 4 nu t, Oseen's solution begun from a
// finite core rather than a singular one
tau += DT;
acore = sqrt(p_core * p_core + 4.0 * p_visc * 1.0e-5 * tau);
float rc2 = acore * acore;
float k1[] = forma_lo_stage(px, py, gg, NV, rc2);
for (int i = 0; i < NV; i++){
ax[i] = px[i] + k1[i * 2] * DT * 0.5;
ay[i] = py[i] + k1[i * 2 + 1] * DT * 0.5;
}
float k2[] = forma_lo_stage(ax, ay, gg, NV, rc2);
for (int i = 0; i < NV; i++){
ax[i] = px[i] + k2[i * 2] * DT * 0.5;
ay[i] = py[i] + k2[i * 2 + 1] * DT * 0.5;
}
float k3[] = forma_lo_stage(ax, ay, gg, NV, rc2);
for (int i = 0; i < NV; i++){
ax[i] = px[i] + k3[i * 2] * DT;
ay[i] = py[i] + k3[i * 2 + 1] * DT;
}
float k4[] = forma_lo_stage(ax, ay, gg, NV, rc2);
for (int i = 0; i < NV; i++){
px[i] += DT / 6.0 * (k1[i * 2] + 2.0 * k2[i * 2] + 2.0 * k3[i * 2] + k4[i * 2]);
py[i] += DT / 6.0 * (k1[i * 2 + 1] + 2.0 * k2[i * 2 + 1] + 2.0 * k3[i * 2 + 1] + k4[i * 2 + 1]);
}
for (int i = 0; i < NT; i++){
// midpoint: evaluate where the tracer would be half a step on, and
// use that velocity for the whole step. Forward Euler spirals a
// tracer outward near a core and draws a decay the equations do not
// contain.
vector v1 = forma_lo_vel(tx[i], ty[i], px, py, gg, NV, rc2);
vector v2 = forma_lo_vel(tx[i] + v1.x * DT * 0.5, ty[i] + v1.y * DT * 0.5,
px, py, gg, NV, rc2);
if (s >= forma_warm){
vector col = forma_ramp(0.84 + 0.22 * hue[i]);
if (open[i] == 0){
// no committed run yet: hold this point, or commit the held
// one alongside it now that the run has two
if (held[i] == 0){
heldx[i] = tx[i]; heldy[i] = ty[i]; held[i] = 1;
} else {
prim[i] = addprim(0, "polyline");
open[i] = 1; held[i] = 0;
// no y flip: the plate holds this domain y-up already
int p0 = addpoint(0, set(heldx[i] * SCALE, heldy[i] * SCALE, 0.0));
setpointattrib(0, "Cd", p0, col);
setpointattrib(0, "Alpha", p0, 0.34);
addvertex(0, prim[i], p0);
int p1 = addpoint(0, set(tx[i] * SCALE, ty[i] * SCALE, 0.0));
setpointattrib(0, "Cd", p1, col);
setpointattrib(0, "Alpha", p1, 0.34);
addvertex(0, prim[i], p1);
}
} else {
int pt = addpoint(0, set(tx[i] * SCALE, ty[i] * SCALE, 0.0));
setpointattrib(0, "Cd", pt, col);
setpointattrib(0, "Alpha", pt, 0.34);
addvertex(0, prim[i], pt);
}
}
tx[i] += v2.x * DT;
ty[i] += v2.y * DT;
if (tx[i] < -HX || tx[i] > HX || ty[i] < -HY || ty[i] > HY){
tx[i] = (random(seed) * 2.0 - 1.0) * HX; seed++;
ty[i] = (random(seed) * 2.0 - 1.0) * HY; seed++;
open[i] = 0; held[i] = 0; // a respawn is not a path it flew
}
}
}
// the cores at the radius they have reached, which is the exhibit: they only
// ever grow. White is positive circulation and the order accent negative, so
// a dipole can be told from a co-rotating pair by eye — the plate's own
// HILITE and TRACE, since a wrangle has no access to either.
int RING_PTS = 64;
for (int i = 0; i < NV; i++){
// a vortex that has propagated out of frame is clipped by the canvas on
// the page, so it is not drawn here either: one ring far outside the
// field would set the bounds of the whole cook
if (abs(px[i]) > HX * 1.15 || abs(py[i]) > HY * 1.15) continue;
vector col = set(1.0, 1.0, 1.0);
if (gg[i] < 0.0) col = forma_ramp(1.0);
int cp = addprim(0, "polyline");
for (int a = 0; a <= RING_PTS; a++){
float th = float(a) / float(RING_PTS) * 2.0 * M_PI;
int pt = addpoint(0, set((px[i] + cos(th) * acore) * SCALE,
(py[i] + sin(th) * acore) * SCALE, 0.0));
setpointattrib(0, "Cd", pt, col);
setpointattrib(0, "Alpha", pt, 0.5);
addvertex(0, cp, pt);
}
}