渦を殺すのは粘性ではない — ニュートンの1/rとLamb

自由渦の粘性項は三つの項が打ち消し合います。コアの外で速度が落ちるなら、それは数値散逸です。

粘性ソルバの検証によく使われる初期条件が自由渦です。速度場 $v_\theta = \Gamma/2\pi r$ を与えて時間を進めます。粘性があるので渦はゆっくり弱まるはずです。ところが半径0.5の速度は、初期拡散時間の50倍が過ぎても小数第4位まで変わりません。コードが間違っているわけではありません。

指数ダイヤルの n を 1 に押すと、右下の「viscous force」バーが緑に落ちます。自由渦の側です。

exponent n in v ~ r^-n1.00

lap difference 0.00 · |tau| 0.00 · |force| 0.00

Slide to n = -1: the spoke stays straight and both bars go green — rigid rotation. Slide to n = 1: the spoke winds up hard, the stress bar is red, but the force bar is green. That is the free vortex. Now stop at n = 0.5, the only exponent that reproduces Kepler’s 3/2 slope on the plot — both bars are red there, so no viscous fluid can hold that profile without something driving it.

ケプラーの n = 0.5 では両方のバーが赤です。そのプロファイルを保つには、誰かが力を入れ続ける必要があります。

1687年の円柱解が、すでに半分を言っています。時計を握っているのはコアと壁です。

1687年、『プリンキピア』第2巻が狙ったもの#

『プリンキピア』は全3巻です。第1巻は力の法則、第3巻は万有引力です。物理学者がほとんど引用しない第2巻が流体力学です。標的は明確でした。デカルトの渦動説（vortex 宇宙論）です。

デカルトは、惑星が太陽のまわりを回るのを、流体の渦に乗って運ばれる現象として説明しました。宇宙は微細な物質で満たされ、その物質が巨大な回転流をなしているという描像です。ニュートンはこの描像を流体力学で反駁しました。流体の運動には必ず抵抗（resistance）が伴うので、外力がなければ渦はいずれ消滅するという論法です。

論法を成立させるには抵抗を定量化する必要がありました。そこで第2巻でニュートンは、流体の抵抗がせん断率に線形に比例すると置きます。

\tau = \mu \frac{du}{dy}

$\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$ を掛けた値です。

T(r) = 2\pi r^{2}\,\tau_{r\theta} = 2\pi \mu\, r^{3} \frac{d}{dr}\!\left(\frac{v_\theta}{r}\right) = \text{const}

この常微分方程式を解くと二つの項が出ます。

v_\theta(r) = A\,r + \frac{B}{r}

前の項は剛体回転、後ろの項は自由渦です。外側境界が無限遠にあり、そこで流体が静止していれば $A=0$ です。残るのは $v_\theta \propto 1/r$ です。これがニュートンが第2巻で得た結果であり、今日 Taylor–Couette 流れの定常解として習うまさにその式です。強制回転と自由渦が一つの円の中に並ぶ Rankine 渦も、この二つの項の断片です。

ケプラーが要求する指数は1/2#

ここで反駁が完成します。デカルトの渦が惑星を運ぶなら、その渦の速度分布はケプラーの第3法則を再現しなければなりません。公転周期が $T \propto r^{3/2}$ なので、速度は $v = 2\pi r / T \propto r^{-1/2}$ です。

ニュートンの流体渦は $r^{-1}$ を与えます。ケプラーは $r^{-1/2}$ を要求します。指数が違います。粘性流体の定常渦は惑星軌道を作れません。冒頭のダイヤルでは n = −1 のスポークがまっすぐで、n = 1 は巻きついても粘性力のバーだけが緑です。

応力はあるのに力がない

ここで冒頭の疑問が解けます。 $v_\theta = C\,r^{-n}$ を回転成分の粘性項に代入してみましょう。

\nu\left(\frac{d^{2}v_\theta}{dr^{2}} + \frac{1}{r}\frac{dv_\theta}{dr} - \frac{v_\theta}{r^{2}}\right) = \nu\,C\,(n^{2}-1)\,r^{-n-2}

最後の $-v_\theta/r^2$ 項は曲率のために付きます。直交座標のラプラシアンには対応する項がありません。括弧を整理すると $n^2-1$ が残ります。

この係数は $n = \pm 1$ でちょうど0になります。 $n = -1$ は剛体回転です。 $n = +1$ が自由渦です。つまり自由渦には粘性力がまったく働きません。

ケプラーの指数 $n=1/2$ では $n^2-1=1/4-1=-3/4$ です。0ではありません。係数が消えるのは $n=1$ だけです。後で出力する r=0.10 の三項 318.81、−159.40、−159.30 は、この0のせいで和が $+9.98\times 10^{-2}$ に留まります。318.81 に対する相対の大きさは約 $3\times 10^{-4}$ で、2次差分の残りです。同じ場所のせん断応力 −0.0318 は0ではありません。応力はあり、発散だけがありません。

せん断応力が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 渦が出ます。

v_\theta(r,t) = \frac{\Gamma}{2\pi r}\left(1 - e^{-r^{2}/4\nu t}\right)

指数項がコアを作ります。 $r \gg \sqrt{4\nu t}$ では括弧が1になり自由渦に戻ります。コア半径は $r_c = 2.2418\sqrt{\nu t}$ で成長し、最大速度は $t^{-1/2}$ で落ちます。

基準は循環です。 $r \to \infty$ で循環はつねに正確に $\Gamma$ です。粘性は渦度を広げるだけで消しません。壁がなければ、その $\Gamma$ は領域の中に残ります。

Pythonで数えた半径ごとの減衰率#

言葉にしたことを数字で確かめます。軸対称の回転方程式 $\partial_t v_\theta = \nu(\partial_{rr} v_\theta + \partial_r v_\theta / r - v_\theta / r^2)$ を、セル中心格子400個で前進オイラー積分します。

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

ν と壁のスイッチを変えると、コア半径と右下の循環バーが一緒に動きます。

nu x1e-32.0core at release0.080time rate1.0x

t 0.00 · r_core 0.000 · v_peak 0.000 · circulation 1.000

Leave the wall off and run it: the core swells, the peak drops, and the amber curve outside the core stays welded to the dashed 1/r line — the circulation bar never leaves 1.000. Now switch the wall on and watch the same bar fall. Raising nu speeds both up by the same factor, which is the point: viscosity sets the clock, the boundary decides whether there is anything to run down.

壁を切るとピークは下がっても循環バーは1.000に貼り付いたままです。壁を入れると同じバーが下がり、t=200 の壁側の循環 0.0342 がその終わりです。

最初のブロックが打ち消し合いを示しています。三項がそれぞれ300規模なのに、和は0.1です。相対的には $3\times10^{-4}$ で、これは2次差分の離散化誤差です。解析的にはちょうど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 を 0.02 から 1.00 まで進めると50倍です。r=0.50 の 0.3183 は $\Gamma/(2\pi r)=1/(2\pi\cdot 0.5)$ と一致し、変化は +0.0% です。動くのは r=0.02 の −90.5% だけです。コア半径は t=0.05 の 0.0163 から t=1.00 の 0.0713 まで育ちます。 $2.2418\sqrt{\nu t}$ に $\nu=10^{-3}$ 、t=1 を入れると 0.0709 です。差は 0.0004 で、格子間隔 dr=1/400=0.0025 より小さいです。v_peak*sqrt(t) は 1.6032 から 1.6057 まで、0.2% の内側です。

壁を立てると同じ渦が死ぬ

壁のブロックが結論です。無限領域では $t = 200$ でも $r=0.9$ の循環が0.6367残っています。これは減衰ではなく、コアがその半径まで膨らんだ結果です。より大きな半径で測れば依然として1です。

壁を立てたほうは0.0342です。最初の循環1の3.4%しか残っていません。同じ時刻の無限領域 0.6367 と比べると、壁が運び出した分のほうがずっと大きいです。壁がしていることは、角運動量を系の外へ運び出すことです。減衰時間は最も遅いモードが決めます。壁で $v_\theta = 0$ 、軸で正則という条件を満たすモードは $J_1(j_{1,1} r/R)$ で、時定数は $R^2/(\nu\, j_{1,1}^2) = 68.1$ 秒です。コードの $j_{1,1}=3.8317$ を $\nu=10^{-3}$ 、 $R=1$ に入れると $3.8317^2\approx 14.68$ で、 $1/(10^{-3}\times 14.68)\approx 68.1$ です。計算では $t = 70$ で0.2186まで落ちており、おおむね合っています。

ですからニュートンの一文はこう修正すべきです。渦を殺すのは粘性そのものではありません。粘性があり、角運動量が抜けていく境界があるときに渦は死にます。無限領域では粘性は渦度を広げるだけです。デカルトを反駁するにはそれで十分でした。宇宙が有限でも無限でも、渦に乗った惑星はケプラーの指数を作れません。

コアの外で速度が落ちるなら、それは粘性ではない

精度検証にコアの外の自由渦を使うと、粘性の離散化は検査されません。その領域の粘性項は解析的に0です。そこでの誤差は対流スキームの誤差です。粘性項を見るなら Lamb–Oseen のコアを格子で覆い、 $r_c(t)=2.2418\sqrt{\nu t}$ の成長を測ります。 $\nu=10^{-3}$ 、t=1 なら 0.0709 で、出力のコア半径 0.0713 と合います。

逆に、コアの外で速度が落ちるなら物理ではなく数値散逸です。風上差分はここで露見します。渦が長く残るべき計算で外側境界を近くに置くと、壁のように働きます。 $R^2/(\nu\, j_{1,1}^2)=68.1$ 秒を物理時間と比べてから領域を決めます。乱流粘性なら $\nu$ の位置に $\nu_t$ が入り、この時定数はずっと短くなります。粘性と慣性の役割を入れ替える Reynolds 数と同じ判断です。粘性だけで、角運動量の出口がなければ、渦は死にません。

渦を殺すのは粘性ではない — ニュートンの1/rとLamb–Oseenコア