乱流粘性モデル — Boussinesq、混合長、k-ε、k-ω

第5章で閉じなかったレイノルズ応力を「巨大な粘性」としてモデル化する工学的トリック。混合長から産業標準の2方程式モデルまでの系譜を整理する。

はじめに

第5章でRANS方程式を導出して我々は袋小路にぶつかった — レイノルズ応力 −ui′uj′‾-\overline{u'_i u'_j} という新しい未知数が六つ現れたのに、方程式は依然として四つ(連続方程式1 + 運動量方程式3)のままだった。これがRANSの閉包問題(closure problem) である。本章では、一世紀にわたって工学者が用いてきた最も単純な解法、「粘性のように扱おう」をたどる。本章を終えると、OpenFOAMのケースファイルで kEpsilon、kOmegaSST といった名前を見たときに、それがどのような仮定を置き、どのような限界をもつかを一文で答えられるようになる。

本論 1 — Boussinesq仮説: 乱流を粘性のように扱う

ブシネスク(Boussinesq、ブシネスク)は1877年、レイノルズ応力が平均ひずみ速度に比例するという単純な仮定を提案した。分子粘性で粘性応力がひずみ速度に比例する(τ=μ ∂u/∂y\tau = \mu \, \partial u / \partial y)のと同様に、乱流でも同じ形の関係が成り立つと仮定するわけだ:

−ui′uj′‾=νt(∂uˉi∂xj+∂uˉj∂xi)−23kδij-\overline{u'_i u'_j} = \nu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) - \frac{2}{3} k \delta_{ij}

ここで νt\nu_t (nu-t、渦粘性または乱流動粘性)は新たに導入された係数、k=12ui′ui′‾k = \tfrac{1}{2}\overline{u'_i u'_i} は単位質量あたりの乱流運動エネルギー(turbulent kinetic energy)、δij\delta_{ij} はクロネッカーのデルタ(Kronecker delta) — i=ji=jのとき1、それ以外0となる記号である。最後の項 23kδij\tfrac{2}{3} k \delta_{ij} はトレース(trace)補正で、初読では「i=ji=jのとき式が破綻しないように付け足された項」程度に理解すれば十分だ。

ここで決定的な違いに注意したい — 分子粘性 ν\nu は流体の物性(空気なら約 1.5×10−5 m2/s1.5 \times 10^{-5} \text{ m}^2/\text{s} で固定)であるのに対し、渦粘性 νt\nu_t は流れの性質である。同じ空気でも静止状態では νt=0\nu_t = 0、強い乱流では νt\nu_t が ν\nu の数百〜数千倍になりうる。つまり νt\nu_t は未知数であり、これをどう決めるかが以降すべてのモデルの出発点となる。

この仮定の利点は明確だ — レイノルズ応力の六つの未知数がスカラー一つ νt\nu_t に圧縮される。閉包問題は「νt\nu_t をどう与えるか」という一つの問題に集約される。限界も明確だ — 応力とひずみ速度が比例しない流れ(強い曲率を伴う流れ、回転、剥離など)では仮定そのものが破綻する。それでも産業CFDの8割以上が今もBoussinesqに依拠する。

本論 2 — Prandtlの混合長モデル

プラントル(Prandtl、プラントル)が1925年に提案した混合長モデル(mixing length model) は最も単純な νt\nu_t の決め方である:

νt=ℓm2∣∂uˉ∂y∣\nu_t = \ell_m^2 \left| \frac{\partial \bar{u}}{\partial y} \right|

ここで ℓm\ell_m (ell-m、混合長)は「流体塊が他の塊と混合するまでに平均的に移動する距離」という直観的意味をもつ。気体分子が衝突までに平均的に進む平均自由行程(mean free path)に着想を得た概念だ。

直観は次のようなものだ。ある位置 yy で流体塊が距離 ℓm\ell_m だけ上方に移動したとすれば、その塊は元の位置の速度を保って周囲と衝突する。二つの位置の平均速度差は約 ℓm⋅∂uˉ/∂y\ell_m \cdot \partial \bar{u}/\partial y であり、これがまさに速度変動 u′u' の典型的大きさになる。この推論を応力の定義に代入すると上の式が得られる。

肝心なのは ℓm\ell_m が位置に依存する点だ。壁面近くでは流体塊が壁の向こうに行けないので ℓm→0\ell_m \to 0、したがって νt→0\nu_t \to 0 となる。自由流れでは ℓm\ell_m は流れの幅に比例する。壁面近くの最も簡単な仮定はカルマン(Kármán、カルマン)定数 κ≈0.41\kappa \approx 0.41 を用いて:

ℓm(y)=κ⋅min⁡(y,2h−y)\ell_m(y) = \kappa \cdot \min(y, 2h - y)

両側の壁で距離0、チャネル中央で最大となる。混合長モデルは単純で未知パラメータがなく、壁近傍の流れに対する適合性が良いため、今でも境界層コードや教科書例題で生き残っている。限界は — 流れごとに ℓm\ell_m を人間が予め決める必要があることだ。剥離、再付着、衝突噴流のように幾何が複雑な流れでは、どの ℓm\ell_m が正しいか事前にはわからない。

本論 3 — k-ε 2方程式モデル

混合長の限界を解くには、νt\nu_t を人間が与えるのではなく方程式に解かせる必要がある。k-εモデル(Launder & Spalding, 1974)は二つの追加輸送方程式を導入する:

  • kk — 乱流運動エネルギー、単位質量あたりの変動運動エネルギー 12ui′ui′‾\tfrac{1}{2} \overline{u'_i u'_i}。単位は m²/s²。乱流の「強さ」 を表す。
  • ε\varepsilon (epsilon、イプシロン) — 乱流散逸率(dissipation rate)。単位時間・単位質量あたり粘性によって熱に散逸する運動エネルギーで、単位は m²/s³。乱流が小さな渦に砕けて消える速さ を表す。

両者で次元解析を行うと、粘性と同じ次元 m²/s をもつ組合せは k2/εk^2/\varepsilon ただ一つで、これに補正定数 Cμ≈0.09C_\mu \approx 0.09 を掛けた形が定義式となる:

νt=Cμk2ε\nu_t = C_\mu \frac{k^2}{\varepsilon}

kk と ε\varepsilon はそれぞれ独自の輸送方程式(生成 = 平均流が乱流を生み出す量、散逸 = ε\varepsilon そのもの、拡散 = νt\nu_t で広がる量)を解くことで格子点ごとに値を得る。本書ではこれらの方程式は導出せず、「この二量を解けば νt\nu_t が自動的に決まる」という事実だけ押さえて先に進む。

k-εは自由流れ(噴流、後流、大気境界層など)では安定で精度が高い。一方で壁近傍の粘性下層(viscous sublayer) では仮定が破綻するため、補正用の「壁関数(wall function)」が必要になる。産業的には1990年代まで最も広く用いられ、今も外部空気力学、換気、燃焼などで標準である。

本論 4 — k-ω モデルとSST

ウィルコックス(Wilcox、ウィルコックス)のk-ωモデル(1988)は二つ目の方程式を ε\varepsilon ではなく比散逸率 ω\omega (omega、オメガ)について解く。ω≈ε/(Cμk)\omega \approx \varepsilon / (C_\mu k) の関係で定義され、単位は 1/s すなわち周波数の次元をもつ。渦粘性は:

νt=kω\nu_t = \frac{k}{\omega}

この単純な変更が与える最大の利点は壁近傍の挙動だ。ω\omega 方程式は壁面で自然に発散する境界条件をもち、壁関数なしに粘性下層まで積分できる。一方の限界は — 自由流れの ω\omega 境界値に結果が敏感な点である。入口条件をわずかに変えただけで後流位置が大きく変わる事例が報告されてきた。

メンター(Menter、メンター)は1994年、両モデルの長所だけを取ったSST k-ω(Shear-Stress Transport k-ω)を提案した。中核となる発想は単純だ — 壁近傍ではk-ω、自由流れではk-εに自動で切り替えるブレンディング関数を導入する。SST k-ωは現在、航空宇宙・自動車・ターボ機械産業で事実上のデフォルトであり、第12章で産業CFDのモデル選定を扱う際に再登場する。

モデル強み弱み
混合長単純、未知パラメータ0ℓm\ell_m を人間が決める
k-ε自由流れで安定壁近傍で壁関数が必要
k-ω壁近傍で精度良好自由流れの境界値に敏感
SST k-ω両領域とも安定計算コストがやや増える

Pythonで確かめる

混合長モデルがチャネル断面でどのような νt\nu_t 分布を生むかをプロットする。チャネル半高さ h=1h=1、平均速度は放物線 uˉ(y)=1−((y−h)/h)2\bar{u}(y) = 1 - ((y-h)/h)^2 と仮定する(中央で最大1、両側の壁で0)。すると ∂uˉ/∂y=−2(y−h)/h2\partial \bar{u}/\partial y = -2(y-h)/h^2 となり、混合長は ℓm(y)=0.41⋅min⁡(y,2h−y)\ell_m(y) = 0.41 \cdot \min(y, 2h - y) である。結果の曲線は両側の壁で0、チャネル中央付近で最大となり、「壁近傍では乱流拡散が弱く中央で強い」という第1章の直観に一致する。

import numpy as np
import matplotlib.pyplot as plt

# チャネル形状: 半高さ h=1、全高さ 2h=2
h = 1.0
y = np.linspace(0.0, 2.0 * h, 401)

# 平均速度プロファイル(放物線): 中央で最大、壁で0
u_bar = 1.0 - ((y - h) / h) ** 2

# 速度勾配(解析的): du/dy = -2(y-h)/h^2
du_dy = -2.0 * (y - h) / (h ** 2)

# 混合長: l_m = κ · min(y, 2h - y)、κ = 0.41(カルマン定数)
kappa = 0.41
l_m = kappa * np.minimum(y, 2.0 * h - y)

# Prandtl混合長による渦粘性: ν_t = l_m^2 · |du/dy|
nu_t = (l_m ** 2) * np.abs(du_dy)

# 可視化
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(y, nu_t, label=r'$\nu_t = \ell_m^2\,|\partial \bar{u}/\partial y|$')
ax.set_xlabel('y (壁からの距離)')
ax.set_ylabel(r'$\nu_t$ (渦粘性)')
ax.set_title('混合長モデル: チャネル断面の渦粘性')
ax.axvline(h, color='gray', linestyle='--', alpha=0.5, label='チャネル中央')
ax.legend()
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig('nu_t_mixing_length.png', dpi=120)
print(f"最大 ν_t = {nu_t.max():.4f} (y = {y[nu_t.argmax()]:.2f})")

曲線を見ると νt\nu_t は両壁で0に落ち(ℓm→0\ell_m \to 0 のため)、中央付近でなだらかな山をなす。ちょうど y=hy=h では ∂uˉ/∂y=0\partial \bar{u}/\partial y = 0 となるため小さな谷ができ、山はその左右に二つに分かれる。この小さなディテールこそが「混合長モデルは対称軸で νt\nu_t が一点だけ0になる非物理的挙動を示す」という、よく知られた限界の視覚的証拠である。k-εのような2方程式モデルはこの欠点を自動的に補正する。

次章へ

ここまで流れ全体に同じモデルを一様に適用してきたが、実際の流れの最も興味深い部分は壁の近くに集中している。第7章: 境界層理論では粘性下層・バッファ層・対数則層という三つの構造を整理し、第6章のモデルたちが各層でどのように(あるいは、なぜ合わないのか)振る舞うかを見る。第7章を終えるとCFD結果を見たときに「壁近傍が合っているか」を自分で判定できるようになる。