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.

Push the exponent dial to n = 1 and the "viscous force" bar at lower right drops to green. That is the free vortex.

At Kepler's n = 0.5 both bars stay red. Holding that profile takes a continuous push.

The cylinder solution from 1687 already says half of this. The core and the wall hold the clock.

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 leading term is rigid rotation, the trailing term 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. A Rankine vortex that puts forced rotation and a free vortex in one circle is built from those same two terms.

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. On the dial above, n = −1 keeps the spoke straight, and n = 1 winds it while the viscous-force bar alone goes green.

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.

At the Kepler exponent $n=1/2$, $n^2-1=1/4-1=-3/4$. Not zero. The coefficient vanishes only at $n=1$. The three terms printed later at r=0.10, 318.81, −159.40 and −159.30, sum to only $+9.98\times 10^{-2}$ because of that zero. Relative to 318.81 that is about $3\times 10^{-4}$, the leftover of a second difference. The shear stress there, −0.0318, is not zero. Stress is present. Only its divergence is gone.

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}$.

Circulation is the reference. As $r \to \infty$ it is always exactly $\Gamma$. Viscosity spreads vorticity; it does not remove it. Without a wall that $\Gamma$ stays inside the domain.

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

Change ν and the wall switch, and the core radius moves together with the circulation bar at lower right.

Wall off, the peak falls while the circulation bar stays stuck at 1.000. Wall on, that bar drops, and the walled circulation 0.0342 at t=200 is where it ends.

The opening 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 time 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 columns with the peak 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.

Marching t from 0.02 to 1.00 is a factor of 50. The value 0.3183 at r=0.50 matches $\Gamma/(2\pi r)=1/(2\pi\cdot 0.5)$, and the change is +0.0%. Only r=0.02 moves, by −90.5%. The core radius grows from 0.0163 at t=0.05 to 0.0713 at t=1.00. Putting $\nu=10^{-3}$ and t=1 into $2.2418\sqrt{\nu t}$ gives 0.0709. The gap is 0.0004, smaller than the spacing dr=1/400=0.0025. v_peak*sqrt(t) stays between 1.6032 and 1.6057, inside 0.2%.

Put a wall around it and the same vortex dies#

The wall 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. Only 3.4% of the original circulation of 1 is left. Next to the unbounded 0.6367 at the same time, the wall has taken far more. 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 code's $j_{1,1}=3.8317$ at $\nu=10^{-3}$ and $R=1$ gives $3.8317^2\approx 14.68$, and $1/(10^{-3}\times 14.68)\approx 68.1$. 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.

A velocity drop outside the core is not viscosity#

Using the free vortex outside the core to check accuracy does not test the viscous discretization. The viscous term there is analytically zero. Error in that region is convection-scheme error. To see the viscous term, cover a Lamb–Oseen core and measure $r_c(t)=2.2418\sqrt{\nu t}$. At $\nu=10^{-3}$ and t=1 that is 0.0709, which matches the printed core radius 0.0713.

Turned around, a velocity drop outside the core is numerical dissipation, not the physics. Upwind schemes get caught there. In a run where the vortex has to survive, an outer boundary placed too close acts like a wall. Compare $R^2/(\nu\, j_{1,1}^2)=68.1$ seconds with the physical time before fixing the domain. With a turbulence model, $\nu_t$ replaces $\nu$ and that time constant gets much shorter. Same judgment as the Reynolds number, which swaps the roles of viscosity and inertia. Viscosity alone, with nowhere for angular momentum to leave, does not kill the vortex.