MMS: 非圧縮性ナビエ・ストークス方程式の検証
理論と物理
MMSとは何か
先生、MMS（製造解法）って「偽の解を作る」って聞いたんですけど、なんでわざわざ偽物を作るんですか？ それって意味あるんですか？
いい質問だ。MMS（Method of Manufactured Solutions）はコード検証の最強ツールなんだ。ポイントは「偽の解」じゃなくて「自分で選んだ解」ということ。
ナビエ・ストークス方程式みたいな複雑な方程式は、一般に解析解が存在しない。でもMMSなら、任意の滑らかな関数を「解だったことにする」ことで、ソルバーの実装が数学的に正しいかを厳密にチェックできる。
え、「解だったことにする」ってどういうことですか？ 方程式を満たさない関数を勝手に解って言っていいんですか？
そこがMMSのトリックだ。例えばナビエ・ストークス方程式の微分演算子を $\mathcal{L}$ とすると、本来は
$$ \mathcal{L}(u) = 0 $$を満たす $u$ を探す。でもMMSでは、自分で好きな $u_\text{mfg}$ を選んで、それを方程式に代入する。当然 $\mathcal{L}(u_\text{mfg}) \neq 0$ だから、その残差をソース項 $f$ として右辺に足す：
$$ \mathcal{L}(u_\text{mfg}) = f $$つまり「$f$ というソース項がある世界では、$u_\text{mfg}$ が厳密解」という問題を人工的に作るわけだ。ソルバーにこの問題を解かせて、$u_\text{mfg}$ と一致するかを確認する。
なるほど！ ソース項を逆算して辻褄を合わせるんですね。でも、それで本当にバグが見つかるんですか？
見つかる。なぜなら、メッシュを系統的に細かくしたときの収束次数（order of accuracy）を確認するからだ。2次精度のスキームで離散化しているなら、メッシュ幅 $h$ を半分にすると誤差は $1/4$ になるはず。これが崩れたら、コードのどこかにバグがある。
例えば自動車メーカーの外装CFDコードで、対流項の離散化に符号ミスがあったケースがある。通常のベンチマーク（Lid-driven cavity等）では「なんとなく合っている」結果が出ていたが、MMSで収束次数が1次に落ちていることが発覚して、バグが見つかった。
非圧縮性ナビエ・ストークス方程式
では、そのMMSを適用する対象の方程式を教えてください。非圧縮性ナビエ・ストークス方程式ってどんな形ですか？
非圧縮性NS方程式は、運動量保存式と連続の式（質量保存）の2つからなる。2次元で書くとこうだ：
運動量保存式（$x$ 方向）:
運動量保存式（$y$ 方向）:
連続の式:
各項の物理的意味
- 非定常項 $\partial u / \partial t$：流速の時間変化率。定常問題ではゼロ。
- 対流項 $(u \cdot \nabla)u$：流体自身が運動量を輸送する非線形効果。高レイノルズ数で支配的。例えば川の合流点で流れが加速する現象。
- 圧力勾配項 $-\nabla p / \rho$：圧力差が流体を押す力。ポンプの駆動力や翼表面の揚力の源。
- 粘性項 $\nu \nabla^2 u$：分子粘性による運動量の拡散。蜂蜜のようにゆっくり広がる流れ。
- ソース項 $f$：重力や外力。MMSではここに「辻褄合わせ」の項が入る。
仮定条件と適用限界
- 非圧縮（マッハ数 $\ll$ 0.3）：密度 $\rho$ が一定
- ニュートン流体：粘性応力が歪み速度に線形比例
- 等温：温度変化による物性変動を無視
- 連続体仮定：Knudsen数 $\ll$ 1（分子の平均自由行程がはるかに小さい）
製造解の構築
NS方程式にMMSを適用するとき、製造解はどう選ぶんですか？ 何でもいいんですか？
理論的には何でもいいが、実用上は守るべきルールがある：
- 連続の式を満たすこと：$\nabla \cdot \mathbf{u}_\text{mfg} = 0$ でないと、非圧縮性ソルバーに矛盾する入力を与えることになる
- 十分に滑らかであること：微分可能でないと、高次精度スキームの検証ができない
- 全ての項を「活性化」すること：例えば定数の製造解だと対流項も粘性項もゼロになってしまい、テストの意味がない
代表的なのがTaylor-Green渦だ。これは実はNS方程式の厳密解でもある（ソース項なしで成立する）が、MMSの説明にはうってつけだ。
Taylor-Green渦型の製造解（2D）:
この製造解が連続の式を満たすことは、どうやって確認するんですか？
偏微分を実際に計算すればいい：
$$ \frac{\partial u_\text{mfg}}{\partial x} = \pi\sin(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$ $$ \frac{\partial v_\text{mfg}}{\partial y} = -\pi\sin(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$足すとゼロ。完璧に連続の式を満たしている。この製造解は $u$ と $v$ の構造が「符号反転＋sin/cos入れ替え」になっているから、必ず非圧縮条件を満たすように設計されているんだ。
ソース項の導出
では、いよいよソース項の計算ですね。製造解をNS方程式に代入するんですよね？
その通り。各項を順に計算していく。$x$ 方向の運動量方程式について示そう。
まず非定常項：
$$ \frac{\partial u_\text{mfg}}{\partial t} = 2\pi^2\nu\cos(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$対流項（$u \cdot \partial u/\partial x + v \cdot \partial u/\partial y$）：
$$ u\frac{\partial u}{\partial x} = -\pi\cos(\pi x)\sin(\pi y)\cdot\pi\sin(\pi x)\sin(\pi y)\cdot e^{-4\pi^2\nu t} $$ $$ v\frac{\partial u}{\partial y} = \pi\sin(\pi x)\cos(\pi y)\cdot(-\pi)\cos(\pi x)\cos(\pi y)\cdot e^{-4\pi^2\nu t} $$粘性項：
$$ \nu\nabla^2 u = \nu\cdot 2\pi^2\cos(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$非定常項と粘性項が相殺し、Taylor-Green渦の場合はソース項 $f_x = 0$ になる。つまりこの製造解はNS方程式の厳密解でもある。
えっ、ソース項がゼロ？ それってMMSとして使えるんですか？
使える。$f=0$ でも「厳密解が分かっている」ことに変わりはないから、収束次数の検証は問題なくできる。ただし実務では、より一般的な製造解を使うことが多い。例えば：
$$ u_\text{mfg} = \sin(ax)\sin(by) + c $$ $$ v_\text{mfg} = \cos(ax)\cos(by) + d $$のように連続の式を満たさないものを選び、$a, b$ の係数を調整する方法もある。この場合はソース項が非ゼロになるから、ソース項の実装テストも兼ねられる。
ここでは教育的にTaylor-Green渦を使うが、実務ではもっと複雑な製造解（多項式＋三角関数の組み合わせ）を使って、全ての離散化項を活性化させるのが鉄則だ。
MMS手順のまとめ：
- 製造解 $u_\text{mfg}$ を選ぶ（連続の式を満たすように設計）
- NS方程式の微分演算子 $\mathcal{L}$ に代入してソース項 $f = \mathcal{L}(u_\text{mfg})$ を計算
- ソース項 $f$ をソルバーの右辺に追加し、$u_\text{mfg}$ を境界条件として設定
- 複数のメッシュ解像度で解を求め、$u_\text{mfg}$ との誤差を計算
- 収束次数 $p$ が離散化スキームの理論値と一致するか確認
収束次数の評価
収束次数って、具体的にどうやって数字を出すんですか？
メッシュを系統的に細分化して、各メッシュでの誤差 $e$ を計算する。2つのメッシュ間での観測収束次数は：
ここで $e_1, e_2$ はそれぞれメッシュ幅 $h_1, h_2$ での誤差ノルム。
具体例で見せてもらえますか？ 例えば2次精度スキームだとどうなりますか？
実際の数値例を見せよう。メッシュを $h$, $h/2$, $h/4$, $h/8$ と細分化した場合：
| メッシュ | 要素数 | $h$ | $L_2$ 誤差ノルム | 観測次数 $p$ | 判定 |
|---|---|---|---|---|---|
| Level 1 | $16 \times 16$ | $1/16$ | $3.87 \times 10^{-3}$ | -- | -- |
| Level 2 | $32 \times 32$ | $1/32$ | $9.72 \times 10^{-4}$ | 1.99 | PASS |
| Level 3 | $64 \times 64$ | $1/64$ | $2.43 \times 10^{-4}$ | 2.00 | PASS |
| Level 4 | $128 \times 128$ | $1/128$ | $6.08 \times 10^{-5}$ | 2.00 | PASS |
判定基準: 観測次数が理論値（2.0）の $\pm 0.1$ 以内なら PASS。系統的に低い場合はバグを疑う。
$L_2$ ノルムの計算式も書いておこう：
離散的には：
ここで $V_i$ はセル $i$ の体積（面積）。
数値解法と実装
有限体積法による離散化
CFDだと有限体積法（FVM）が主流ですよね。MMSのソース項はFVMだとどう扱うんですか？
FVMではセル平均値を扱うから、ソース項もセル体積で積分する。セル $i$ の離散化された運動量方程式は：
$$ \sum_f \phi_f u_f - \sum_f \nu (\nabla u)_f \cdot \mathbf{S}_f = -\frac{1}{\rho}\sum_f p_f \mathbf{S}_f + \int_{V_i} f \, dV $$ここで $\phi_f = \mathbf{u}_f \cdot \mathbf{S}_f$ は面フラックス、$\mathbf{S}_f$ は面法線ベクトル。ソース項の体積積分 $\int_{V_i} f \, dV$ は、セル中心の値 $f_i \cdot V_i$ で近似する（セル中心一点近似）。
一点近似で精度は大丈夫なんですか？
いい着眼点だ。セル中心一点近似は2次精度なので、2次精度スキームの検証なら問題ない。ただし3次以上の精度を検証したい場合は、Gauss求積でソース項を高精度に積分する必要がある。ここは見落としやすいポイントだから注意してほしい。
ソース項の実装
ソース項をソルバーにどうやって追加するんですか？
大きく2つの方法がある：
- 方法1: 明示的にソース項を追加 -- 運動量方程式の右辺にユーザー定義のソース項関数を追加。OpenFOAMなら
fvOptions、Fluentなら UDF のDEFINE_SOURCEマクロを使う。 - 方法2: 体積力として追加 -- 重力項と同じ形式でソース項を与える。実装が簡単だが、圧力項のソース項を別途処理する必要がある場合もある。
実務では方法1が一般的だ。ソース項の関数をハードコードするか、座標と時刻を引数にしたユーザー関数として実装する。
誤差ノルムの計算
誤差ノルムは $L_1$, $L_2$, $L_\infty$ の3種類を計算するのが望ましい：
3種類も計算する意味は何ですか？ $L_2$ だけじゃダメですか？
ダメとは言わないが、$L_\infty$（最大誤差）ノルムが崩れるケースがある。例えば境界近傍だけ精度が落ちているバグは、$L_2$ では平均化されて見えにくいが、$L_\infty$ なら一発で分かる。3ノルム全てで正しい収束次数が出れば、コードの信頼性はかなり高い。
SymPyによるソース項の自動導出
ソース項の手計算って、対流項の微分で間違えそうで怖いです...
実務ではSymPy（Pythonのシンボリック計算ライブラリ）で自動導出するのが鉄則だ。手計算は検算用にとどめる。特に非線形対流項 $(u \cdot \nabla)u$ の展開で符号ミスが頻発するから、シンボリック計算でクロスチェックすることを強く推奨する。参考コード：
import sympy as sp
x, y, t, nu = sp.symbols('x y t nu')
pi = sp.pi
# 製造解の定義
u = -sp.cos(pi*x) * sp.sin(pi*y) * sp.exp(-2*pi**2*nu*t)
v = sp.sin(pi*x) * sp.cos(pi*y) * sp.exp(-2*pi**2*nu*t)
p = -sp.Rational(1,4)*(sp.cos(2*pi*x)+sp.cos(2*pi*y))*sp.exp(-4*pi**2*nu*t)
# 連続の式の確認
div = sp.diff(u, x) + sp.diff(v, y)
print("div(u) =", sp.simplify(div)) # => 0
# x方向ソース項: f_x = du/dt + u*du/dx + v*du/dy + dp/dx - nu*(d2u/dx2 + d2u/dy2)
fx = (sp.diff(u,t) + u*sp.diff(u,x) + v*sp.diff(u,y)
+ sp.diff(p,x) - nu*(sp.diff(u,x,2) + sp.diff(u,y,2)))
print("f_x =", sp.simplify(fx)) # => 0 (Taylor-Green渦の場合)
実践ガイド
ステップバイステップ手順
実際にMMSをやるとき、最初から最後まで何をすればいいんですか？
以下の7ステップだ：
- 製造解を選定：連続の式を満たす滑らかな関数を選ぶ。全ての方程式項が活性化されることを確認
- ソース項を導出：SymPy等でシンボリックに計算。手計算とのクロスチェック推奨
- 計算領域と境界条件を設定：境界には製造解の値をDirichlet条件として与える
- 系統的メッシュ生成：少なくとも4水準（例: 16x16, 32x32, 64x64, 128x128）。メッシュ比は2倍ずつが望ましい
- 各メッシュで計算実行：同一のソルバー設定。残差を十分に収束させる
- 誤差ノルムを計算：$L_1$, $L_2$, $L_\infty$ の3種類。速度 $u$, $v$ と圧力 $p$ のそれぞれについて
- 収束次数を判定：$p = \log(e_1/e_2) / \log(h_1/h_2)$ が理論値と一致するか確認
境界条件は全部Dirichletでいいんですか？ 圧力の境界条件はどうするんですか？
速度はDirichlet（製造解の値を固定）でOK。圧力は少し注意が必要で、非圧縮性ソルバーでは圧力は1点を固定（reference pressure）するか、全境界でNeumann条件（圧力勾配を製造解から計算）にする。Dirichlet条件を全面に課すとoverconstrainedになることがあるから、ソルバーのマニュアルを確認してほしい。
OpenFOAMでの実装例
OpenFOAMだと具体的にどうやるんですか？
OpenFOAMでは icoFoam（非圧縮性層流ソルバー）に fvOptions でソース項を追加する方法が最も簡単だ。手順は：
blockMeshDictで $[0,1] \times [0,1]$ の正方形領域を生成0/Uに製造解の初期条件と境界条件をcodedFixedValueで設定constant/fvOptionsにcodedSourceでソース項を定義- 4水準のメッシュで実行し、
postProcess -func 'error'で誤差を収集
Taylor-Green渦の場合はソース項ゼロなので、fvOptions は不要。ただし汎用の製造解ではソース項の実装が必須になる。
Ansys Fluentでの実装
Fluentの場合は UDF（User Defined Function）を使う：
DEFINE_SOURCE(x_mom_source, ...)：$x$ 方向運動量のソース項DEFINE_PROFILE(u_bc, ...)：境界条件として製造解の速度を与えるDEFINE_INIT(mms_init, ...)：初期条件を製造解で設定
FluentではC言語でUDFを書く必要があるが、ソース項の数式はSymPyで導出した式をそのままCコードに変換できる（sp.ccode() 関数）。
品質保証チェックリスト
MMSを実行するとき、チェックすべきポイントをまとめてもらえますか？
以下を全項目クリアすれば、MMS検証は合格だ：
| # | チェック項目 | 合格基準 |
|---|---|---|
| 1 | 製造解が連続の式を満たすか | $\nabla \cdot \mathbf{u}_\text{mfg} = 0$（シンボリックに確認） |
| 2 | ソース項の導出が正しいか | SymPyと手計算で一致 |
| 3 | 境界条件が製造解と整合するか | 境界上で $u_\text{num} = u_\text{mfg}$ |
| 4 | 各メッシュで残差が十分に収束しているか | 残差 $< 10^{-10}$（離散化誤差を丸め誤差以下に） |
| 5 | $L_2$ ノルムの収束次数が理論値と一致するか | $|p_\text{obs} - p_\text{theory}| < 0.1$ |
| 6 | $L_\infty$ ノルムの収束次数も確認したか | 同上 |
| 7 | 速度 $u$, $v$ と圧力 $p$ の全てで確認したか | 全変数で合格 |
ソフトウェア比較
CFDソルバー別のMMS対応状況
ソルバーによってMMSのやりやすさって変わりますか？
かなり変わる。ソース項をどれだけ柔軟に追加できるかがカギだ。
| ソルバー | ソース項追加方法 | MMS向け機能 | 難易度 |
|---|---|---|---|
| OpenFOAM | fvOptions / codedSource | フルソースアクセス、カスタムソルバー作成可 | 低 |
| Ansys Fluent | UDF (DEFINE_SOURCE) | C言語UDF、Scheme scripting | 中 |
| STAR-CCM+ | Field Function / User Code | Java API、テーブル入力 | 中 |
| COMSOL | 方程式ベースモデリング | 任意PDEの右辺にソース項を直接記述 | 低 |
| FEniCS/Firedrake | 弱形式で直接定義 | Pythonでシンボリック定義、UFL記法 | 低 |
| SU2 | ソースコード修正 | MMS検証フレームワーク内蔵 | 中 |
OpenFOAMとCOMSOLが一番やりやすいんですね。商用ソルバーだとUDFの開発環境が必要になるのがハードル高そう...
その通り。ただ、最近のFluent 2024R2以降は Python UDF にも対応し始めたから、だいぶ楽になってきている。SU2はオープンソースで、MMS検証のリグレッションテストが公式にCI/CDに組み込まれている。コード検証を重視する文化が根付いているソルバーほど、MMS対応が充実しているね。
先端技術
3次元MMSへの拡張
ここまで2次元の話でしたが、実務のCFDは3次元ですよね。3Dへの拡張は難しいですか？
基本的な考え方は同じだが、連続の式 $\nabla \cdot \mathbf{u} = 0$ を3成分で満たす必要がある。3D Taylor-Green渦は：
$$ u = A\cos(ax)\sin(by)\sin(cz) $$ $$ v = B\sin(ax)\cos(by)\sin(cz) $$ $$ w = C\sin(ax)\sin(by)\cos(cz) $$ただし $Aa + Bb + Cc = 0$ を満たすように係数を選ぶ。ソース項の導出は3Dだと項数が爆発的に増えるから、SymPyの利用はほぼ必須になる。
非定常問題の時間精度検証
時間積分スキームの精度も同時に検証できる。手順は：
- 空間メッシュを十分に細かく固定（空間誤差が支配的にならないように）
- 時間刻み $\Delta t$ を系統的に変化させる（$\Delta t$, $\Delta t/2$, $\Delta t/4$...）
- 最終時刻での誤差を計算し、時間方向の収束次数を確認
Euler陰解法なら1次、Crank-Nicolsonなら2次、BDF2なら2次が期待される。空間と時間を独立に検証するのが重要で、同時に細分化すると、どちらの精度が崩れているか区別できなくなる。
乱流モデルへの適用
RANS乱流モデル（例えばSpalart-Allmaras）にもMMSは使えますか？
使える。むしろ乱流モデルこそMMSが威力を発揮する。なぜなら、乱流モデルの実装にはソース項（生成項・散逸項）が多数含まれ、バグが潜みやすいからだ。
Spalart-Allmarasモデルの場合、輸送方程式の乱流粘性 $\tilde{\nu}$ にも製造解を設定し、NS方程式と結合したソース項を導出する。Rumsey（NASA LaRC）による公開検証データがデファクトスタンダードになっている。
トラブルシューティング
収束次数が理論値より低い
先生、2次精度のはずなのに収束次数が1.5くらいしか出ないんですけど、これってバグですか？
バグの可能性もあるが、まずは以下を確認しよう：
- メッシュが十分にasymptotic rangeに入っているか？ 粗すぎるメッシュでは高次の打ち切り誤差項が支配的で、理論通りの収束次数が出ない。解像度を上げて再計算してみる
- 反復残差が十分に収束しているか？ 残差が $10^{-6}$ 程度だと、微細メッシュでは離散化誤差と反復誤差が同程度になり、収束次数が崩れる。$10^{-10}$ 以下まで落とす
- ソース項の実装で丸め誤差はないか？
floatではなくdoubleを使っているか確認 - 境界条件の離散化精度は？ 内部が2次でも、境界処理が1次だと全体の精度が落ちる
あ、残差を $10^{-6}$ で止めていました... それが原因かもしれません。
MMSでは「離散化誤差 >> 反復誤差」の状態を作らないと意味がない。実務CFDでは $10^{-6}$ で十分だけど、MMS検証では残差を極限まで落とす。これは初心者が最もハマるポイントだから覚えておいてほしい。
ソース項の不一致
ソース項を実装したのに、収束次数が0次（つまり精度なし）になりました。完全に壊れてます...
収束次数0次は「メッシュを細かくしても誤差が減らない」状態。これはほぼ確実にソース項の計算ミスだ。よくある原因：
- 符号ミス：対流項 $(u \cdot \nabla)u$ の展開で $+$ と $-$ を逆にしている
- 密度の扱い：$\rho = 1$ を仮定しているのに、ソース項では $\rho$ をかけている（またはその逆）
- 座標系の不一致：ソース項の $(x,y)$ とソルバー内部の座標系がずれている
- 時刻の不一致：非定常問題で、ソース項の評価時刻がタイムステップの始点なのか終点なのかが不整合
対処法は、まずソース項をゼロにして、製造解をNS方程式の厳密解（Taylor-Green渦）に限定してテストすること。これで正しい収束次数が出れば、ソース項の実装に問題があると切り分けられる。
境界条件の落とし穴
境界条件で気をつけるべきことは他にありますか？
いくつか典型的な落とし穴がある：
- 周期境界条件：Taylor-Green渦は周期境界と相性が良いが、製造解が本当に周期的でないと、境界で不連続が生じて収束次数が崩壊する
- 圧力基準点：非圧縮性ソルバーでは圧力は定数分の不定性がある。圧力のMMSでは、計算値と製造解の平均値を引いてから比較する必要がある
- 壁関数との干渉：乱流計算で壁関数を使っている場合、壁近傍のモデル化が MMS の厳密解と矛盾する。壁関数なし（resolve to wall）で検証すべき
MMSって奥が深いですね... でも、これをちゃんとやればCFDコードの信頼性が格段に上がることは分かりました！
その通り。MMSは「コードが正しいことの数学的証明」に最も近い手法だ。ASME V&V 20やAIAA Guide for V&Vでも、MMSはコード検証の推奨手法として明記されている。特に新しいソルバーを開発したり、既存ソルバーに新機能を追加したりしたときは、必ずMMSで検証してほしい。
関連トピック
なった
詳しく
報告