杀死涡的不是粘性 — 牛顿的1/r与Lamb–Oseen涡核

验证粘性求解器时常用的初始条件是自由涡。给定速度场 $v_\theta = \Gamma/2\pi r$，然后推进时间。既然有粘性，涡应该慢慢变弱。可是半径0.5处的速度，在初始扩散时间的50倍之后，小数点后四位仍然纹丝不动。代码没有错。本文讨论这是为什么，以及真正杀死涡的东西是什么。答案的一半在1687年就已经印出来了。

1687年，《原理》第二卷瞄准的是什么#

《自然哲学的数学原理》共三卷。第一卷是力的定律，第三卷是万有引力。物理学家几乎不引用的第二卷是流体力学。它的靶子很明确：笛卡尔的涡旋（vortex）宇宙论。

笛卡尔把行星绕太阳的运动解释为被流体涡旋携带着旋转。宇宙充满细微物质，这些物质形成巨大的旋转流动。牛顿用流体力学反驳了这幅图景。流体运动必然伴随阻力（resistance），所以没有外力时涡旋终将消亡。

要让这个论证成立，必须把阻力定量化。于是在第二卷里，牛顿设定流体的阻力与剪切率成线性比例。

$\tau$ 是剪应力，$du/dy$ 是速度梯度，比例系数 $\mu$ 就是粘度。这是流体粘度第一次被数学地定义。后来人们发现不遵守这一线性关系的流体反而更多，于是把遵守的那一类称为牛顿流体。

从一个圆柱推出的1/r#

牛顿的推导今天读来仍然现代。一根无限长的圆柱在粘性流体中以恒定角速度旋转。定常状态下，穿过半径为 $r$ 的圆柱面向外传递的扭矩处处相等。否则角动量会在两个这样的面之间的流层里堆积。

柱坐标下纯旋转流动的剪应力是 $\tau_{r\theta} = \mu\, r\, \frac{d}{dr}(v_\theta/r)$。它不是简单的 $\mu\, dv_\theta/dr$。刚体旋转（$v_\theta \propto r$）不应有剪切，而这个形式自动满足该条件。单位长度的扭矩等于应力乘以力臂 $r$ 再乘以周长 $2\pi r$。

解这个常微分方程会得到两项。

第一项是刚体旋转，第二项是自由涡。如果外边界在无穷远且那里流体静止，则 $A=0$。剩下的就是 $v_\theta \propto 1/r$。这就是牛顿在第二卷得到的结果，也正是今天作为 Taylor–Couette 流动定常解讲授的那个公式。

开普勒要求的指数是1/2#

反驳到这里就完成了。如果笛卡尔的涡旋携带行星，那么这个涡旋的速度分布必须重现开普勒第三定律。公转周期满足 $T \propto r^{3/2}$，所以速度满足 $v = 2\pi r / T \propto r^{-1/2}$。

牛顿的流体涡旋给出 $r^{-1}$。开普勒要求 $r^{-1/2}$。两个指数不同。粘性流体中的定常涡旋无法产生行星轨道。

在下面的旋钮里亲手移动这个指数。

n = −1 时辐条保持笔直，这就是刚体旋转。推到 n = 1，辐条被卷得很紧，可右下角的 "viscous force" 条却降到了绿色。在开普勒的 n = 0.5 处两条都是红的。观察点就在这里 — 要维持那条剖面，必须有人持续施力。

有应力，却没有力

开头的疑问在这里解开。把 $v_\theta = C\,r^{-n}$ 代入旋转分量的粘性项。

最后那个 $-v_\theta/r^2$ 项来自曲率。直角坐标的拉普拉斯算子里没有对应项。整理括号之后剩下 $n^2-1$。

这个系数在 $n = \pm 1$ 处恰好为0。$n = -1$ 是刚体旋转，$n = +1$ 就是自由涡。也就是说，自由涡上一点粘性力也没有。

原因不是剪应力为0。应力本身是 $\tau_{r\theta} = -\mu\Gamma/\pi r^{2}$，在涡核附近相当大。发生的事情是：一个流体微元内侧面受到的扭矩与外侧面受到的扭矩恰好抵消。力不是应力，而是应力的散度。自由涡正是那条散度为0的特殊剖面。

同样的话也可以用涡量来说。$\omega = \frac{1}{r}\frac{d(r v_\theta)}{dr}$，而 $v_\theta = \Gamma/2\pi r$ 时 $r v_\theta$ 是常数，所以在 $r>0$ 的任何地方 $\omega = 0$。不可压缩流动的粘性力可以写成 $-\nu\,\nabla\times\boldsymbol{\omega}$。没有涡量，就没有粘性力。

涡核铺开的速度决定了时钟

说没有涡量，准确地讲是除原点之外。环量 $\Gamma$ 必须存在于某处，在理想自由涡里它集中成原点上的一个δ函数。粘性真正干活的地方就在这里。

把环量 $\Gamma$ 集中在原点然后求解扩散方程，得到的就是 Lamb–Oseen 涡。

指数项造出涡核。在 $r \gg \sqrt{4\nu t}$ 处括号趋于1，回到自由涡。涡核半径按 $r_c = 2.2418\sqrt{\nu t}$ 增长，峰值速度按 $t^{-1/2}$ 下降。

关键量是环量。当 $r \to \infty$ 时环量永远精确等于 $\Gamma$。粘性只是把涡量摊开，并不消灭它。在下面的实验里拨动壁面开关。

壁面关闭时涡核鼓起、峰值下降，但右下角的环量条始终停在1.000。打开壁面，同一条就开始往下掉。要操作的是 $\nu$ 和壁面，要观察的是涡核半径以及两条环量历史之间的差别。

用Python数出的分半径衰减率#

把话变成数字。用前向欧拉法，在400个格心径向网格上推进轴对称旋转方程 $\partial_t v_\theta = \nu(\partial_{rr} v_\theta + \partial_r v_\theta / r - v_\theta / r^2)$。

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):
    """Lamb-Oseen涡的切向速度。t很小时收敛到裸的1/r涡。"""
    return GAMMA / (2 * np.pi * rr) * (1 - np.exp(-rr * rr / (4 * NU * t)))
 
def swirl_terms(v):
    """旋转拉普拉斯算子的三块。分开放是为了看清相消。"""
    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): 纯旋转中粘性项剩下的全部。"""
    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. 裸的1/r涡: 三个大项互相抵消
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. 真的推进一段时间，与解析解对照
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. 涡核按sqrt(nu t)铺开，峰值按1/sqrt(t)下降。环量不变
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. 在 r = R 立一堵壁，同一个涡就死了
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

第一块展示了相消。三项各自是300量级，和却是0.1。相对而言是 $3\times10^{-4}$，这是二阶差分的离散误差。解析上和精确为0。同一位置的剪应力则不为0。

第二块回答了开头的疑问。$r = 0.5$ 和 $r = 0.95$ 处的变化是0.0%。推进了初始时间的50倍也如此。只有 $r = 0.02$ 被削掉了90%。粘性只在涡核里干了活。

第三块最后两列印证了标度关系。v_peak*sqrt(t) 从1.603到1.606，固定在0.2%之内。涡核半径与 $2.2418\sqrt{\nu t} = 0.0709$ 相比是0.0713，也对得上。整个过程中环量一直是1.0000。

立一堵壁，同一个涡就死了

第四块是结论。在无界区域里，$t = 200$ 时 $r=0.9$ 处的环量仍有0.6367。这不是衰减，而是涡核已经鼓到了那个半径。到更大的半径去测，它仍然是1。

立了壁的那一侧是0.0342。95%以上消失了。壁做的事情是把角动量搬出系统。衰减时间由最慢的模态决定。满足壁面 $v_\theta = 0$ 与轴上正则性的模态是 $J_1(j_{1,1} r/R)$，时间常数为 $R^2/(\nu\, j_{1,1}^2) = 68.1$ 秒。计算在 $t = 70$ 时降到0.2186，大致吻合。

所以牛顿那句话需要修一处。杀死涡的不是粘性本身。当有粘性、而且有让角动量流出去的边界时，涡才会死。在无界区域里粘性只是把涡量摊开。但用来反驳笛卡尔已经足够了：无论宇宙有限还是无限，骑在涡旋上的行星都造不出开普勒的指数。

写下一个涡的测试算例之前

这个结果在工程上有三种用法。

不要用自由涡做精度验证。 涡核之外粘性项恰好为0，所以那片区域的误差是对流格式的误差，而不是粘性离散的误差。要验证粘性项，就得用网格充分覆盖 Lamb–Oseen 的涡核，并测量 $r_c(t)$ 的增长。

数值衰减要在这里量。 反过来，这个性质是极好的诊断。放入一个自由涡，如果涡核之外速度下降了，那不是物理，而是代码的数值耗散。迎风类格式在这里立刻现形。

计算域大小会改变答案。 在涡必须长时间存活的问题里 — 翼尖涡、旋翼尾流 — 外边界放得太近就会像壁一样起作用。定域之前，先把 $R^2/(\nu\, j_{1,1}^2)$ 与关心的物理时间比一比更稳妥。如果计算用了湍流模型，$\nu$ 的位置要换成 $\nu_t$，这个时间常数会短上几个数量级。