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

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

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

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

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

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

$\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$ です。これがニュートンが第2巻で得た結果であり、今日 Taylor–Couette 流れの定常解として習うまさにその式です。

ケプラーが要求する指数は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 に押すとスポークは強く巻きつきますが、右下の「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で数えた半径ごとの減衰率#

言葉にしたことを数字で確かめます。軸対称の回転方程式 $\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

最初のブロックが打ち消し合いを示しています。三項がそれぞれ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 = 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$ が入り、この時定数は桁で短くなります。