Viscosity Is Not What Kills a Vortex — Newton's 1/r and the Lamb–Oseen Core

A free vortex is a standard initial condition for checking a viscous solver. You seed $v_\theta = \Gamma/2\pi r$ and march. Viscosity is on, so the vortex should slowly weaken. Instead the velocity at radius 0.5 is unchanged to four decimal places after fifty times the initial diffusion time. The code is not broken. This post explains why, and what actually does kill a vortex. Half the answer was already in print in 1687.

What Book II of the Principia was aiming at#

The Principia runs to three books. Book I is the laws of force, Book III is universal gravitation. Book II, the one physicists rarely cite, is fluid mechanics. Its target was explicit: Descartes' vortex cosmology.

Descartes explained planetary orbits as bodies carried around by a fluid vortex. Space is packed with fine matter, and that matter turns in a vast circulating flow. Newton attacked the picture on fluid-mechanical grounds. Fluid motion always carries resistance, so without an external force the vortex must eventually die.

To make that argument he had to quantify resistance. So in Book II he posits that the resistance of a fluid is linear in the rate of shear.

Here $\tau$ is shear stress, $du/dy$ is the velocity gradient, and the constant of proportionality $\mu$ is the viscosity. This is the first place fluid viscosity is defined mathematically. Later it turned out that far more fluids violate the linear relation than obey it, which is why the obedient ones are called Newtonian.

The 1/r that falls out of a single cylinder#

Newton's derivation still reads as modern. An infinitely long cylinder spins at constant angular velocity inside a viscous fluid. At steady state the torque transmitted outward through any cylindrical surface of radius $r$ must be the same everywhere. Otherwise angular momentum piles up in the layer between two such surfaces.

For pure circular flow in cylindrical coordinates the shear stress is $\tau_{r\theta} = \mu\, r\, \frac{d}{dr}(v_\theta/r)$. It is not simply $\mu\, dv_\theta/dr$. Rigid rotation ($v_\theta \propto r$) must carry no shear, and this form satisfies that condition automatically. Torque per unit length is the stress times the moment arm $r$ times the circumference $2\pi r$.

Solving that ordinary differential equation gives two terms.

The first is rigid rotation, the second is a free vortex. If the outer boundary is infinitely far away and the fluid is at rest there, then $A=0$. What survives is $v_\theta \propto 1/r$. That is Newton's Book II result, and it is the same expression taught today as the steady solution of Taylor–Couette flow.

Kepler demands an exponent of 1/2#

Now the refutation closes. If Descartes' vortex carries the planets, its velocity profile has to reproduce Kepler's third law. The orbital period scales as $T \propto r^{3/2}$, so the speed scales as $v = 2\pi r / T \propto r^{-1/2}$.

Newton's fluid vortex gives $r^{-1}$. Kepler demands $r^{-1/2}$. The exponents disagree. A steady vortex in a viscous fluid cannot produce planetary orbits.

Move the exponent yourself in the dial below.

At n = −1 the spoke stays straight: that is rigid rotation. Push to n = 1 and the spoke winds up hard, yet the "viscous force" bar at lower right drops to green. At Kepler's n = 0.5 both bars are red. That is the observation point — holding that profile requires something to keep pushing.

Stress without force#

This resolves the opening puzzle. Substitute $v_\theta = C\,r^{-n}$ into the viscous term of the swirl equation.

That last $-v_\theta/r^2$ term is a curvature term. The Cartesian Laplacian has no counterpart to it. Collecting the bracket leaves $n^2-1$.

That coefficient is exactly zero at $n = \pm 1$. The case $n = -1$ is rigid rotation. The case $n = +1$ is the free vortex. So the free vortex feels no viscous force at all.

Not because the shear stress vanishes. The stress itself is $\tau_{r\theta} = -\mu\Gamma/\pi r^{2}$, which is enormous near the core. What happens is that the torque on the inner face of a fluid element exactly cancels the torque on its outer face. Force is not stress, it is the divergence of stress. The free vortex is the special profile whose stress divergence is zero.

The same statement works through vorticity. With $\omega = \frac{1}{r}\frac{d(r v_\theta)}{dr}$ and $v_\theta = \Gamma/2\pi r$, the product $r v_\theta$ is constant, so $\omega = 0$ everywhere for $r>0$. The viscous force in incompressible flow can be written $-\nu\,\nabla\times\boldsymbol{\omega}$. No vorticity, no viscous force.

The core sets the clock#

Vorticity is absent everywhere except one place: the origin. The circulation $\Gamma$ has to live somewhere, and in the idealized free vortex it is concentrated in a delta function at the axis. That is where viscosity actually goes to work.

Concentrate circulation $\Gamma$ at the origin and solve the diffusion problem, and you get the Lamb–Oseen vortex.

The exponential builds the core. For $r \gg \sqrt{4\nu t}$ the bracket goes to 1 and the free vortex is recovered. The core radius grows as $r_c = 2.2418\sqrt{\nu t}$, and the peak velocity falls as $t^{-1/2}$.

The circulation is the key quantity. As $r \to \infty$ it is always exactly $\Gamma$. Viscosity spreads vorticity; it does not remove it. Toggle the wall switch in the experiment below.

With the wall off, the core swells and the peak drops, but the circulation bar at lower right never leaves 1.000. Switch the wall on and that same bar starts to fall. Operate $\nu$ and the wall; watch the core radius and the gap between the two circulation histories.

Decay rates by radius, counted in Python#

Time to put numbers on the prose. March the axisymmetric swirl equation $\partial_t v_\theta = \nu(\partial_{rr} v_\theta + \partial_r v_\theta / r - v_\theta / r^2)$ with forward Euler on 400 cell-centered radial cells.

import numpy as np
 
NU, GAMMA, R, N = 1.0e-3, 1.0, 1.0, 400
dr = R / N
r = (np.arange(N) + 0.5) * dr
 
def lamb_oseen(rr, t):
    """Tangential velocity of a Lamb-Oseen vortex; small t is the bare 1/r vortex."""
    return GAMMA / (2 * np.pi * rr) * (1 - np.exp(-rr * rr / (4 * NU * t)))
 
def swirl_terms(v):
    """The three pieces of the swirl Laplacian, kept apart so the cancellation shows."""
    ghost = np.concatenate(([-v[0]], v, [v[-1] * r[-1] / (r[-1] + dr)]))
    return ((ghost[2:] - 2 * ghost[1:-1] + ghost[:-2]) / dr**2,
            (ghost[2:] - ghost[:-2]) / (2 * dr) / r,
            -v / r**2)
 
def swirl_operator(v, outer):
    """nu * (v_rr + v_r/r - v/r^2): all that is left of viscosity in pure swirl."""
    g_out = -v[-1] if outer == 'wall' else v[-1] * r[-1] / (r[-1] + dr)
    w = np.concatenate(([-v[0]], v, [g_out]))
    v_rr = (w[2:] - 2 * w[1:-1] + w[:-2]) / dr**2
    v_r = (w[2:] - w[:-2]) / (2 * dr)
    return NU * (v_rr + v_r / r - v / r**2)
 
def march_swirl(v, t_end, outer):
    dt = 0.2 * dr * dr / NU
    for _ in range(int(round(t_end / dt))):
        v = v + dt * swirl_operator(v, outer)
    return v
 
def circulation_at(v, radius):
    return 2 * np.pi * radius * np.interp(radius, r, v)
 
# 1. the bare 1/r vortex: three large terms that cancel
a, b, c = (np.interp(0.10, r, x) for x in swirl_terms(GAMMA / (2 * np.pi * r)))
print(f"free vortex at r=0.10:  v_rr={a:+8.2f}  v_r/r={b:+8.2f}  -v/r^2={c:+8.2f}  sum={a+b+c:+.2e}")
print(f"                        shear stress tau_rtheta = {-NU * GAMMA / (np.pi * 0.10**2):+.4f}")
 
# 2. march a real vortex and check it against the analytic solution
t0, t1 = 0.02, 1.0
v0 = lamb_oseen(r, t0)
v1 = march_swirl(v0, t1 - t0, 'free')
print(f"\nmarched {t0} -> {t1} s   L-inf vs Lamb-Oseen = {np.max(np.abs(v1 - lamb_oseen(r, t1))):.1e}")
print("    r    v(0.02)   v(1.00)    change")
for x in (0.02, 0.05, 0.15, 0.50, 0.95):
    p, q = np.interp(x, r, v0), np.interp(x, r, v1)
    print(f"{x:5.2f} {p:9.4f} {q:9.4f} {100 * (q - p) / p:+8.1f}%")
 
# 3. core spreads as sqrt(nu t), peak falls as 1/sqrt(t), circulation holds
print("\n    t   r_core   v_peak  v_peak*sqrt(t)  Gamma(0.9)")
v, tc = v0.copy(), t0
for t in (0.05, 0.20, 0.50, 1.00):
    v = march_swirl(v, t - tc, 'free'); tc = t
    i = int(np.argmax(v))
    print(f"{t:5.2f} {r[i]:8.4f} {v[i]:8.4f} {v[i] * np.sqrt(t):13.4f} {circulation_at(v, 0.9):11.4f}")
 
# 4. add a wall at r = R and the same vortex dies
print("\n    t   Gamma(0.9)  unbounded    walled")
vw, tc = v0.copy(), t0
for t in (1.0, 20.0, 70.0, 200.0):
    vw = march_swirl(vw, t - tc, 'wall'); tc = t
    print(f"{t:8.1f}     {circulation_at(lamb_oseen(r, t), 0.9):9.4f} {circulation_at(vw, 0.9):9.4f}")
print(f"\nslowest walled mode:  R^2/(nu*j11^2) = {R**2 / (NU * 3.8317**2):.1f} s")

free vortex at r=0.10:  v_rr= +318.81  v_r/r= -159.40  -v/r^2= -159.30  sum=+9.98e-02
                        shear stress tau_rtheta = -0.0318
 
marched 0.02 -> 1.0 s   L-inf vs Lamb-Oseen = 7.7e-04
    r    v(0.02)   v(1.00)    change
 0.02    7.9233    0.7562    -90.5%
 0.05    3.1851    1.4782    -53.6%
 0.15    1.0611    1.0573     -0.4%
 0.50    0.3183    0.3183     +0.0%
 0.95    0.1675    0.1675     +0.0%
 
    t   r_core   v_peak  v_peak*sqrt(t)  Gamma(0.9)
 0.05   0.0163   7.1699        1.6032      1.0000
 0.20   0.0312   3.5880        1.6046      1.0000
 0.50   0.0513   2.2697        1.6049      1.0000
 1.00   0.0713   1.6057        1.6057      1.0000
 
    t   Gamma(0.9)  unbounded    walled
     1.0        1.0000    0.9773
    20.0        1.0000    0.4182
    70.0        0.9446    0.2186
   200.0        0.6367    0.0342
 
slowest walled mode:  R^2/(nu*j11^2) = 68.1 s

The first block shows the cancellation. Three terms of order 300 sum to 0.1. That is $3\times10^{-4}$ in relative terms, and it is the discretization error of the second difference. Analytically the sum is exactly zero. The shear stress at the same location is not.

The second block answers the opening puzzle. At $r = 0.5$ and $r = 0.95$ the change is 0.0%, after fifty times the initial time. Only $r = 0.02$ lost 90%. Viscosity worked in the core and nowhere else.

The last two columns of the third block confirm the scaling. v_peak*sqrt(t) holds between 1.603 and 1.606, a 0.2% spread. The core radius matches $2.2418\sqrt{\nu t} = 0.0709$ against a measured 0.0713. Through all of it the circulation reads 1.0000.

Put a wall around it and the same vortex dies#

The fourth block is the conclusion. In the unbounded case the circulation at $r=0.9$ still reads 0.6367 at $t = 200$. That is not decay; it is the core having swollen out to that radius. Measure at a larger radius and it is still 1.

The walled case reads 0.0342. More than 95% is gone. What the wall does is remove angular momentum from the system. The decay time is set by the slowest mode. The mode satisfying $v_\theta = 0$ at the wall and regularity on the axis is $J_1(j_{1,1} r/R)$, with time constant $R^2/(\nu\, j_{1,1}^2) = 68.1$ seconds. The run reads 0.2186 at $t = 70$, which is about right.

So Newton's sentence needs one amendment. Viscosity by itself is not what kills a vortex. A vortex dies when there is viscosity and a boundary through which angular momentum can leave. In an unbounded domain viscosity only spreads vorticity around. That was still enough to refute Descartes: finite universe or infinite, a planet riding a vortex cannot reproduce Kepler's exponent.

Before you write your next vortex test case#

Three practical consequences follow.

Do not use a free vortex for accuracy verification. Outside the core the viscous term is identically zero, so any error there is convection-scheme error, not viscous-discretization error. To verify the viscous term, resolve a Lamb–Oseen core on the grid and measure the growth of $r_c(t)$.

Do use it to measure numerical damping. Turned around, the same property makes an excellent diagnostic. Seed a free vortex, and any velocity drop outside the core is your code's numerical dissipation rather than physics. Upwind-family schemes are caught immediately.

Domain size changes the answer. In problems where a vortex has to survive — wingtip vortices, rotor wakes — an outer boundary placed too close behaves like a wall. Compare $R^2/(\nu\, j_{1,1}^2)$ against the physical time of interest before you fix the domain. If the run uses a turbulence model, $\nu$ becomes $\nu_t$ and that time constant shortens by orders of magnitude.