ラージエディシミュレーション(LES) — 大きな渦は解き、小さな渦はモデル化する

DNSとRANSの中間策 — 空間フィルタで大きな渦だけを直接解き、格子より小さな渦はSGSモデルに任せるLESの数学と直観。

はじめに

第10章のDNSは全ての渦スケールを解くので正確だが、格子点数が Re9/4\mathrm{Re}^{9/4} で爆発し、産業用の高レイノルズ数流れにはほぼ使えなかった。第5章のRANSは全ての変動を一度に平均してしまうので速いが、渦そのものの時間発展を失う。本章を読み終えると、読者は LES(Large Eddy Simulation、ラージエディシミュレーション) がこの二つの両極端の間でどう折り合いをつけるか — 「大きな渦は解き、小さな渦はモデル化する」という一文の意味を — 数式と一行のコードで説明できるようになる。

本論 1 — 空間フィルタ:信号から小スケールを切り落とす操作

LESの出発点は 時間平均ではなく空間フィルタ である。あるフィルタ関数 GG (フィルタカーネル)とフィルタ幅 Δ\Delta (大文字デルタ、格子幅と同程度のスケール)を決めると、フィルタ済み速度 u~i\tilde{u}_i (「u チルダ」と読む)は次のように定義される:

u~i(x⃗,t)=∫G(x⃗−x⃗′) ui(x⃗′,t) dx⃗′\tilde{u}_i(\vec{x}, t) = \int G(\vec{x} - \vec{x}') \, u_i(\vec{x}', t) \, d\vec{x}'

言葉で言えば「ある点の速度を、その周囲 Δ\Delta の範囲で重み付き平均した値」である。Δ\Delta より大きい渦はほぼそのまま通過し、Δ\Delta より小さい渦は平均にならされて消える。実務でよく使われる GG には三種類あって — ボックスフィルタ(box filter、一定区間内の一様平均)、ガウシアンフィルタ(Gaussian filter、ガウス釣鐘形の重み)、シャープスペクトラルフィルタ(sharp-spectral filter、フーリエ空間で特定波数以上を切るフィルタ) — のいずれかを選ぶ。

第5章のレイノルズ分解 ui=uˉi+ui′u_i = \bar{u}_i + u'_i と形は似て見えるが本質は違う — レイノルズ分解は 時間平均 なので平均された場は時間に対して変わらないか非常にゆっくりとしか変わらないのに対し、LESのフィルタ済み場 u~i\tilde{u}_i は依然として 時間方向に揺らぐ 大きな渦の動力学をそのまま保つ。

本論 2 — フィルタ済みナビエ–ストークスとSGS応力

非圧縮ナビエ–ストークス方程式にフィルタを適用してみよう。時間微分・圧力勾配・粘性項は、フィルタと微分が交換できるおかげで、すべてフィルタ記号を内側に入れた形でそのまま残る。第5章のRANS導出と同様、問題は 対流項 uj∂ui/∂xju_j \partial u_i / \partial x_j である。

対流項にフィルタをかけると二つの速度の積が現れるが、積のフィルタは一般にフィルタの積に等しくない: uiuj~≠u~iu~j\widetilde{u_i u_j} \ne \tilde{u}_i \tilde{u}_j。この差を SGS応力(Subgrid-Scale stress、格子下スケール応力) と定義する:

τijsgs=uiuj~−u~iu~j\tau^{sgs}_{ij} = \widetilde{u_i u_j} - \tilde{u}_i \tilde{u}_j

このテンソルは格子より小さなスケールの渦が運ぶ運動量交換をまとめた量である。言い換えれば「自分が直接解いている大きな渦場 u~i\tilde{u}_i だけでは表せない、Δ\Delta 以下の領域のあらゆる未解決効果」が一つのテンソルに圧縮されている。したがってLESの閉じ込め問題は この τijsgs\tau^{sgs}_{ij} を u~i\tilde{u}_i のみの関数としてどう書くか に帰着する — 形は第5章の閉じ込め問題と同じだが、モデル化の対象が「全スペクトルの平均効果(RANS)」から「格子以下の小さな渦の効果(LES)」に絞られている。

本論 3 — スマゴリンスキーモデルとその後継

最も単純なSGS閉じ込めは、1963年にスマゴリンスキー(Smagorinsky)が提案した 渦粘性モデル である。格子以下の渦の効果が追加的な粘性のように働くと仮定し、その粘性の大きさを局所の解像歪み速度(resolved strain rate)の大きさに比例させる:

νsgs=(CsΔ)2∣S~∣,∣S~∣=2S~ijS~ij\nu_{sgs} = (C_s \Delta)^2 |\tilde{S}|, \qquad |\tilde{S}| = \sqrt{2 \tilde{S}_{ij} \tilde{S}_{ij}}

ここで S~ij=12(∂u~i/∂xj+∂u~j/∂xi)\tilde{S}_{ij} = \tfrac{1}{2}(\partial \tilde{u}_i / \partial x_j + \partial \tilde{u}_j / \partial x_i) はフィルタ済み速度場から計算される歪み速度テンソルで、CsC_s (スマゴリンスキー定数)は流れの種類に応じて 0.10.1–0.170.17 の値を取る — 等方性乱流では 0.170.17 付近、壁面境界層では 0.10.1 以下が普通である。CsC_s をユーザーがいちいち決めなければならないという欠点を補うために、1991年にジャーマノ(Germano)が提案した 動的スマゴリンスキー(dynamic Smagorinsky) モデルでは、解像場そのものから二種類のフィルタ幅を比較して CsC_s を時々刻々自動で決める。

この節約のおかげで、LESの格子点数はDNSの Re9/4\mathrm{Re}^{9/4} ではなく Re1.8\mathrm{Re}^{1.8} 程度 にしか増えない。Re=106\mathrm{Re} = 10^6 のときDNSは 1013.510^{13.5}、LESは 1010.810^{10.8} の格子点が必要となり、約500倍の差がつく — 同じスーパーコンピュータでもDNSは不可能だがLESは可能、という領域が存在することを意味する。産業用航空機の空力解析、自動車外部流れ、都市風環境シミュレーションがLESの主戦場である理由はここにある。

Pythonで確かめる

LESの最も単純なモデルケースであるボックスフィルタを1次元で直接適用してみよう。合成信号に短波長のノイズを乗せ、幅5のボックスフィルタで平均すると、フィルタ済み信号が大スケールのトレンドは保ちつつ短波長変動だけを抑えるのが一目で分かる。

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)

# 大きな渦(低周波サイン波) + 小さな渦(高周波ノイズ)
x = np.linspace(0, 10, 500)
large_scale = np.sin(x)                              # 大きな渦
small_scale = 0.4 * rng.standard_normal(x.size)      # 小さな渦
u = large_scale + small_scale

# 幅5のボックスフィルタを適用 — LESフィルタの1D版
kernel = np.ones(5) / 5
u_tilde = np.convolve(u, kernel, mode='same')

# 可視化:原信号とフィルタ済み信号を同じ軸に重ねて描く
plt.figure(figsize=(8, 4))
plt.plot(x, u, label='原信号 u (大きな渦 + 小さな渦)', alpha=0.5)
plt.plot(x, u_tilde, label='フィルタ済み ũ (LESが直接解く量)', linewidth=2)
plt.xlabel('x')
plt.ylabel('速度')
plt.title('ボックスフィルタ:大きな渦は残し、小さな渦は減衰させる')
plt.legend()
plt.tight_layout()
plt.show()

実行すると、原信号の曲線はノイズで上下に揺れるが、フィルタ済みの曲線はサイン波の大きな流れを滑らかに追跡する。LESがシミュレーションで直接解いているのはまさにこの滑らかな曲線で、二つの曲線の差に相当する情報 — つまり切り落とされた小さな渦 — がSGS応力 τijsgs\tau^{sgs}_{ij} に圧縮されてモデルへ渡される。

次章へ

DNSは正確だが高価で、RANSは速いが平均しか教えてくれず、LESはその中間で大きな渦の動力学を残す。産業現場ではこの三つのうちどれをいつ使うかを決めることが、最も重要なエンジニアリング判断である。第12章: 産業CFDにおけるモデル選択では、航空・自動車・建築・気象分野の代表事例を通じて「この流れにはどのモデルをどんな格子で使うか」を実際に選びながら、本書全体を総括する。