乱流は三本の物差しを同時に使う — Kolmogorovカスケードと-5/3則

当てはめ窓を 10/L と 0.1/η の外まで開くと、−5/3 であるべき傾きが −1.4 まで寝ます。

1941年、モスクワで書かれた一冊のノート#

戦争が始まる二か月前、Kolmogorovは三本の短い論文を学術誌に送りました。合わせて十六ページに満たない分量です。風洞も、測定データも、図版も一枚もありません。あったのは次元解析と二つの統計的仮定だけでした。

そのノートの一行が八十年後のCFDを縛り続けています。十分高いReynolds数では、エネルギーがスケール間にどう分布するかは一本の長さだけで決まる — Kolmogorovのミクロスケール $\eta$ です。 この記事では、その一行がどう生まれ、なぜ格子解像度とDNS(Direct Numerical Simulation、直接数値計算)の費用を決めるのかを追います。最後に合成データに対して-5/3の傾きを当てはめます。

time speed のスライダーを上げると、砕ける間隔だけが縮み、大渦から小渦への順番はそのままです。

Each eddy lives a turnover time τ ∝ ℓ^(2/3), then splits into smaller ones. The smallest eddies (gray) dissipate as heat.

time speed: 1.0×

大渦(シアン)は長く生きます — $\tau \sim L^{2/3}$ 。中間渦(オレンジ)はもっと早く砕け、小渦(ピンク)は灰色(粘性散逸)へほぼ瞬時に渡されます。各渦が生きている間に下流へ流すエネルギー流束は、すべてのスケールで等しい — それが慣性小領域の定義です。

一つの流れに三本の物差しを詰め込む

乱流の一かたまりを広げてみると、三つの長さが見えてきます。

積分スケール $L$ — 最大渦の大きさです。普通は流れの幾何(管径、煙突の幅)で決まります。エネルギーを蓄える場所です。
Taylorミクロスケール $\lambda$ — 二点速度自己相関関数が原点で持つ曲率から定義されます。大スケールと小スケールの仲介役です。
Kolmogorovミクロスケール $\eta$ — 粘性が運動エネルギーを熱に擦り変える最小の物差しです。分子平均自由行程よりは遥かに大きいですが、 $L$ に比べれば桁違いに小さいです。

三本の物差しが同じ流れの中に同居しています。Reynolds数が大きくなると、 $L$ と $\eta$ の隔たりが広がります。その隔たりが慣性小領域(inertial subrange)です。

次元解析が橋を架ける

Kolmogorovの仮定は二行で済みます。

十分小さなスケールでは、統計量は粘性 $\nu$ とエネルギー散逸率 $\varepsilon$ のみに依存する。
慣性小領域では粘性も落ち、統計量は $\varepsilon$ だけに依存する。

仮定1から最小スケールが組み立てられます。 $\nu$ の単位は $\mathrm{m^2/s}$ 、 $\varepsilon$ は $\mathrm{m^2/s^3}$ 。長さを取り出すには:

\eta = \left(\frac{\nu^3}{\varepsilon}\right)^{1/4}, \quad u_\eta = (\nu \varepsilon)^{1/4}, \quad \tau_\eta = \left(\frac{\nu}{\varepsilon}\right)^{1/2}

$\eta$ が長さ、 $u_\eta$ が速度、 $\tau_\eta$ が時間。単位だけで系が閉じます。

仮定2はもっと深く切り込みます。慣性小領域のエネルギースペクトル $E(k)$ は $\varepsilon$ と波数 $k$ にしか依存しません。 $E$ の単位は $\mathrm{m^3/s^2}$ 。次元を揃えると:

E(k) = C\, \varepsilon^{2/3}\, k^{-5/3}

$C$ はKolmogorov定数、実験的に $C \approx 1.5$ 。-5/3則は仮定と単位だけから出てきます。測定点を一つも使っていません。

スケール分離 — Reが広げる隔たり#

大渦の見積もり $\varepsilon \sim u'^3 / L$ (エネルギー注入が散逸と等しい定常状態)を使うと:

\frac{\eta}{L} = \mathrm{Re}^{-3/4}, \quad \frac{\lambda}{L} \sim \mathrm{Re}^{-1/2}, \quad \frac{\tau_\eta}{T} \sim \mathrm{Re}^{-1/2}

ここで $\mathrm{Re} = u'L/\nu$ です。Reynolds数が10倍になると、 $\eta/L$ は $10^{-3/4} \approx 0.18$ 倍に縮みます。すべてのスケールを格子で解こうとすると、セル数は $\mathrm{Re}^{9/4}$ で爆発します。実験室の風洞 $\mathrm{Re} = 10^4$ では $L/\eta \approx 1{,}000$ 、3D DNSは $\sim 10^9$ セル。航空機の翼 $\mathrm{Re} = 10^7$ なら $10^{15}$ セル — 現在のスーパーコンピュータでも届きません。

当てはめ窓の十進幅は、境界からそのまま出ます。 $k$ が $10/L$ から $0.1/\eta$ までなら、幅は $\log_{10}(L/\eta)-2$ decade です。 $\mathrm{Re}=10^{4}$ で $L/\eta\approx 1{,}000$ なら、幅は 1 decade です。十倍ごとに $\eta/L$ は $0.18$ 倍です。その縮小を二回かけた $\mathrm{Re}=10^{6}$ では $\eta/L$ が約 $0.032$ 倍になり、 $L/\eta$ は約 $3\times 10^{4}$ まで開きます。窓幅は約 2.5 decade です。スライダーの $\mathrm{Re}=10^{6}$ で直線が 4 decade 近く見えるのは、窓の外の峰と粘性の尾まで目盛りに入るからです。

コードの定数で $u'=1$ 、 $L=1$ なら $\varepsilon=u'^{3}/L=1$ です。 $f_L$ と $f_\eta$ を掛ける前の慣性区間は $E(k)=1.5\,k^{-5/3}$ です。コメントの当てはめ傾き −1.661 と期待値 −1.667 の差は 0.006 です。雑音 10% を載せても、窓を守れば三桁が残ります。

この隔たりが管内で最初に開く過程は、水筋が乱れて境界層が層を成す臨界へ続きます。公称 Reynolds が大きくなるほど、 $L$ と $\eta$ のあいだに入れる格子は $\mathrm{Re}^{9/4}$ で増え、 $10^{4}$ の三次元はすでに $10^{9}$ セルです。

エネルギーが流れる道

大渦が一つ回転しています。粘性は弱すぎてそのエネルギーを直接削れません。代わりに大渦は小渦に砕け、小渦はさらに小さな渦に砕けます。砕ける時間は渦の回転時間(turnover time) $\tau_\ell \sim \ell/u_\ell$ です。この破砕が、エネルギー散逸率 $\varepsilon$ という一つの数値ですべてのスケールを束ねます。

-5/3傾きと両端の罠#

$E(k)$ を両対数で描くと、中央に傾き-5/3の直線帯が現れます。しかし両端では別の物理が支配します。

$k \sim 1/L$ 付近 — エネルギー含有域。幾何が決める最大渦が住む場所。 $E(k)$ がピークを打ちます。
$k \sim 1/\eta$ 付近 — 粘性散逸域。 $E(k)$ が指数的に落ちます。粘性がエネルギーを飲み干します。

Reynolds数のスライダーを上げると、-5/3の直線帯が長くなります。

Model energy spectrum E(k) on log-log axes (Pope 2000 form).

Reynolds number: Re = 10^5.0 ≈ 1.0e+5 show -5/3 reference line

L/η = 5.62e+3 · η = 1.78e-4 · λ (Taylor) = 1.22e-2

$\mathrm{Re} = 10^3$ では-5/3帯はほとんど見えません。 $\mathrm{Re} = 10^6$ になると、ほぼ四桁にわたる直線が開きます。実験的には、一桁でも綺麗な-5/3傾きが取れれば、その流れは十分に発達した乱流とみなされます。

Python — 合成スペクトルに-5/3を当てはめる#

実験で-5/3傾きをどう抽出するか。両対数軸で直線当てはめ、それだけです。合成データで試してみましょう。

import numpy as np
 
C_K = 1.5
c_L, p0 = 6.78, 2
beta_e, c_eta = 5.2, 0.40
 
def pope_spectrum(k, eps, L, eta):
    """Pope のモデルスペクトル(等方等質乱流)。"""
    kL = k * L
    kEta = k * eta
    f_L = (kL / np.sqrt(kL**2 + c_L)) ** (5/3 + p0)
    f_eta = np.exp(-beta_e * ((kEta**4 + c_eta**4)**0.25 - c_eta))
    return C_K * eps**(2/3) * k**(-5/3) * f_L * f_eta
 
# 設定
L_int = 1.0      # 積分長さ
u_rms = 1.0      # rms 速度
Re = 1e5
nu = u_rms * L_int / Re
eps = u_rms**3 / L_int
eta = (nu**3 / eps) ** 0.25
 
# 波数格子上の擬似測定(乗法的雑音付き)
rng = np.random.default_rng(0)
k = np.logspace(-1, np.log10(1.0 / eta) - 0.5, 80)
E_true = pope_spectrum(k, eps, L_int, eta)
E_meas = E_true * np.exp(0.10 * rng.standard_normal(k.size))   # ~10% 対数正規雑音
 
# 慣性域だけに -5/3 を当てはめる  (10/L < k < 0.1/eta)
mask = (k > 10 / L_int) & (k < 0.1 / eta)
slope, _ = np.polyfit(np.log(k[mask]), np.log(E_meas[mask]), 1)
print(f"L/eta = {L_int / eta:.1f}")
print(f"fitted slope = {slope:+.3f}  (expected -1.667)")
# L/eta = 17783.3
# fitted slope = -1.661  (expected -1.667)

80点の雑音込みデータでも当てはめ傾きは-1.66、理論値-5/3と三桁一致します。雑音を30%まで上げても傾きは-1.59~-1.74の範囲に留まります — 慣性域の頑健さを示しています。

ここで一つ質問。なぜマスクに $k > 10/L$ と $k < 0.1/\eta$ を入れたのか。両端の $f_L$ 、 $f_\eta$ が直線を曲げてしまうからです。実測でも当てはめ窓を狭く取らないと、傾きが-1.4まで落ちます。「私のデータは-5/3を示さない」という報告の半分は、窓の取り方が原因です。

格子が届かない場所で生き延びる

DNSの費用が見えたので、迂回路が読みやすくなります。

RANS(Reynolds-Averaged Navier–Stokes) — 時間平均だけを解き、変動はモデルで閉じます。 $\eta$ を解かないので、格子コストは $\mathrm{Re}^0$ で済みます。産業の主力ですが、異方性レイノルズ応力には弱いです。
LES(Large Eddy Simulation) — 大渦のみ解き、格子以下の小渦をSGS(サブグリッドスケール)モデルで閉じます。壁モデルの有無で格子コストは $\mathrm{Re}^{0.5}$ ~ $\mathrm{Re}^1$ 。DNSとRANSの中間です。
DNS — $\eta$ まで解き、 $\tau_\eta$ まで時間積分します。コストは $\mathrm{Re}^{3}$ (3D + 時間)。学術的な道具です。

三手法ともに同じ-5/3線の上に立っています。RANSは慣性領域全体をモデルに置き換え、LESはその一部だけを解き、DNSはすべてを解く。どこで切るかがモデリングの判断です。

当てはめ窓を両端まで開くと傾きが −1.4 まで寝る#

$\eta/L = \mathrm{Re}^{-3/4}$ が窓の幅を決めます。マスクを $10/L < k < 0.1/\eta$ の中に置かないと、 $f_L$ と $f_\eta$ が両端を曲げ、流束 $\varepsilon$ だけが残すべき傾きは −1.4 まで下がります。コードのコメントにある −1.661 は、その窓を守ったときの数字です。

$\mathrm{Re} = 10^4$ でも三次元はすでに重く、 $\mathrm{Re} = 10^7$ では格子が $10^{15}$ まで行きます。その下で切ってサブグリッドに渡すとき、係数が負に漏れる場所は Germano 恒等式で動的係数が負になった箇所に書いてあります。カスケードの流束を信じることと、切ったあとの係数をそのまま使うことは、別の判断です。