直接数値シミュレーション (DNS) — なぜ研究専用なのか

モデルなしでナビエ–ストークス方程式を直接解くDNSのコストが、レイノルズ数の9/4乗で爆発する理由をナプキン計算で確かめる。

はじめに

第9章では、最大の渦 LL (エル、積分スケール)と、粘性によって運動エネルギーが熱に変換される最小の渦 η\eta (エータ、コルモゴロフスケール)が同時に存在することを確認した。本章では、その両端のスケールをまとめて格子に乗せようとした瞬間にコンピュータに何が起きるかを見積もる。読み終えたあと、「なぜ産業現場でDNSを使わないのか?」という問いに、ナプキン一枚の計算で答えられるようになる。

本論 1 — DNSの定義: モデルがない

DNS(Direct Numerical Simulation、直接数値シミュレーション)とは、ナビエ–ストークス方程式を、いかなる乱流モデルも用いずにそのまま解く方法である。第5章のRANSは平均を分離して −ui′uj′‾-\overline{u'_i u'_j} をモデルで閉じ、第11章のLESは大きな渦だけを解いて小さな渦をモデルで処理する。DNSはそのような閉包モデルを持たない。代わりに格子を η\eta まで細かく切って、エネルギーカスケードのあらゆるスケールを数値的に直接解像する。

そのためDNSは「近似のない数値実験」と呼ばれる。風洞実験と同等の資格で新しい乱流統計を生み出す研究ツールである。

本論 2 — 格子点数が Re9/4\mathrm{Re}^{9/4} で増える

押さえるべき式は一本だけだ。領域の一辺の長さを LL、最小セルサイズを η\eta とする。1次元方向に必要な格子点数は単純に

Lη\frac{L}{\eta}

となる。第9章で示したとおり、積分スケールのレイノルズ数 Re=UL/ν\mathrm{Re} = UL/\nu (ν\nu、ニュー、動粘性係数)について

Lη∼Re3/4\frac{L}{\eta} \sim \mathrm{Re}^{3/4}

が成り立つ。3次元では三方向すべてを解像する必要があるため、総格子点数 NN は

N∼(Lη)3∼Re9/4N \sim \left(\frac{L}{\eta}\right)^3 \sim \mathrm{Re}^{9/4}

となる。9/4 は 2.25 である。レイノルズ数が10倍になると格子は約178倍、100倍になると約31,623倍に膨らむ。

さらに時間積分の安定条件(CFL条件) Δt∝Δx/U∝η/U\Delta t \propto \Delta x / U \propto \eta / U が加わり、空間と時間を併せて見ると総コストは Re3\mathrm{Re}^3 に比例する。本文では空間格子だけを表で確かめる。

本論 3 — ナプキン計算: どのRe からスパコンが要るか

N∼Re9/4N \sim \mathrm{Re}^{9/4} という一行だけで、産業適用の可否はほぼ決まる。

Re\mathrm{Re}N∼Re9/4N \sim \mathrm{Re}^{9/4}現実的な解釈
10210^2∼104.5≈3×104\sim 10^{4.5} \approx 3 \times 10^4ノートPCで数分
10310^3∼106.75≈6×106\sim 10^{6.75} \approx 6 \times 10^6ワークステーションで一晩
10410^4∼109≈109\sim 10^{9} \approx 10^9国家スパコン1回分
10510^5∼1011.25≈2×1011\sim 10^{11.25} \approx 2 \times 10^{11}世界トップクラスでの数年プロジェクト
10610^6∼1013.5≈3×1013\sim 10^{13.5} \approx 3 \times 10^{13}現世代のハードウェアでは不可能

自動車の外部空力は Re∼107\mathrm{Re} \sim 10^7、旅客機の翼は Re∼108\mathrm{Re} \sim 10^8 の領域だ。表の最後の行から、さらに2〜3段下がってようやく産業問題に届く。だからDNSは産業コードの検証や新しいモデル用のデータ生成に使われるだけで、設計ループにはRANSとLESが入る。

Pythonで確かめる

先の表を自分で出力してみよう。numpy のみ、Re9/4\mathrm{Re}^{9/4} を計算するだけだ。

import numpy as np

# DNS の格子点数の見積もり: N ~ Re^(9/4)
# 9/4 = 2.25, 第9章の L/eta ~ Re^(3/4) を3次元に持ち上げた結果

Re_values = [1e2, 1e3, 1e4, 1e5, 1e6]
exponent = 9 / 4  # 2.25

print(f"{'Re':>10} | {'N ~ Re^(9/4)':>15}")
print("-" * 30)
for Re in Re_values:
    N = Re ** exponent
    print(f"{Re:>10.1e} | {N:>15.2e}")

# 参考: 時間積分の CFL 条件まで考慮すると総コストは Re^3 に比例する。
# 上の表は空間格子のみの見積もりで、実際の実行時間はさらに急に増える。

実行すると Re=104\mathrm{Re} = 10^4 で N≈3.16×109N \approx 3.16 \times 10^9 が出る。第1章で見た「Reという単一変数が流れの定性的状態を決める」という言葉は、コストの観点では単一変数が計算可能性の境界を決めるに置き換わる。

次章へ

DNSの格子コストを受け入れられないなら、選択肢は二つある。平均を取ってすべての変動をモデルで処理する(RANS、第5–6章)か、大きな渦だけ解いて小さな渦をモデルで処理する(LES)か。第11章: ラージエディシミュレーション (LES)では後者を定量的に扱う — DNSの Re9/4\mathrm{Re}^{9/4} がLESではどの指数まで下がるかが核心となる。