CFDメッシュ独立性確認（Grid Independence Study）

理論と物理

概要

先生、CFDのメッシュ独立性確認って、構造FEMでやる格子収束確認と基本的に同じことですか？

原理は同じだよ。どちらも「メッシュをもっと細かくしても解がほとんど変わらない状態」を確認する作業だ。Richardson外挿でメッシュサイズ$h \to 0$の解を推定して、GCI（Grid Convergence Index）で格子誤差を定量化するのも共通。

じゃあ構造でやったことがあれば、CFDでもそのまま使えますか？

注意点がいくつか違うんだ。CFDでは3つの追加条件がある。第一に、y+が乱流モデルと整合しているか。これが合っていないと、メッシュをどんなに細かくしても正しい解に収束しない。第二に、CFDはナビエ-ストークス方程式の非線形性が強いから、観測される収束次数$p$が理論値（2次精度スキームなら$p=2$）より低くなることが多い。第三に、反復収束の残差を十分に下げないと、離散化誤差と反復誤差が混ざってGCIの意味がなくなる。

構造よりもハマりポイントが多いんですね…。体系的な手順ってあるんですか？

ある。Celikら（2008）の論文 "Procedure for Estimation and Reporting of Uncertainty due to Discretization in CFD Applications"（Journal of Fluids Engineering）が、CFDでのGCI適用の事実上の標準手順だ。ASME Journal投稿のレビュー基準にもなっている。この5段階の手順を押さえれば、実務でも論文でも使えるよ。

Richardson外挿法

Richardson外挿って、ざっくりどういう考え方なんですか？

離散解$f(h)$をTaylor展開すると、真の解$f_{\text{exact}}$との関係はこうなる：

$h$は代表セルサイズ、$p$は離散化スキームの精度次数、$C$は未知定数だ。例えば2次精度の中心差分なら理論上$p=2$だから、メッシュを半分にすると誤差は約$1/4$になる。この関係式から、メッシュサイズの異なる2つ以上の解を使って$f_{\text{exact}}$を推定するのがRichardson外挿の基本だよ。

でも$p$と$C$の2つが未知なら、方程式が2つ要りますよね。だからメッシュが3水準必要になるわけですか？

その通り。3水準のメッシュ（細$h_1$、中$h_2$、粗$h_3$、ただし$h_1 < h_2 < h_3$）で解$f_1, f_2, f_3$を得て、連立方程式を解くことで観測収束次数$p$を推定できるんだ。これが2水準だと$p$を仮定するしかなくなるから、「メッシュを細かくしたら変化が小さかったのでOK」という主観的な判断になってしまう。

GCIの定義と意味

GCIの式を教えてください。構造FEMで見たものと同じですか？

式そのものは同じだよ。細メッシュと中メッシュ間のGCIはこうなる：

ここで各変数の意味はこうだ：

• $e_a^{21} = (f_1 - f_2)/f_1$：細メッシュ解$f_1$と中メッシュ解$f_2$の相対差
• $r_{21} = h_2 / h_1$：格子細分化比（$r > 1$）
• $p$：観測収束次数（後述の式で算出）
• $F_s$：安全係数。3水準の外挿を使う場合は$F_s = 1.25$、2水準で$p$を仮定する場合は$F_s = 3.0$

$F_s = 1.25$と$3.0$でかなり違いますね。やっぱり3水準やるべきということですか？

そう。$F_s = 3$は「$p$を仮定したことの不確かさ」に対するペナルティだ。たった1水準余分にやるだけで安全係数が$1/2.4$になるんだから、3水準でやらない理由はないよ。例えば自動車の外部空力で抗力係数$C_d$のGCIが1.5%以下なら、風洞試験との比較に耐える精度と言えるだろう。

Celikの5段階手順

Celikの手順って具体的にどんなステップですか？

Celikら（2008, JFE 130(7)）の5段階手順を説明しよう。

Step 1. 代表セルサイズ$h$の定義

3Dの場合：

$N$は総セル数、$\Delta V_i$は各セルの体積だ。2Dなら面積の平方根を使う。

Step 2. 3水準のメッシュで解を求める

細分化比$r = h_{\text{coarse}}/h_{\text{fine}}$は理想的には$r \geq 1.3$。$r = 2$（セル数8倍）がベストだが、計算コストが厳しければ$r = 1.3$〜$1.5$（セル数2〜3倍）でもOK。3つの解$f_1, f_2, f_3$を得る。

Step 3. 観測収束次数$p$を求める

ここで$\varepsilon_{21} = f_2 - f_1$、$\varepsilon_{32} = f_3 - f_2$。$r_{21}$と$r_{32}$が等しい（等比細分化）なら$q=0$となり、簡略版の式が使える：

Step 4. Richardson外挿値$f_{\text{ext}}^{21}$を求める

Step 5. 誤差推定量を報告する

• 相対近似誤差：$e_a^{21} = |(f_1 - f_2)/f_1|$
• 相対外挿誤差：$e_{\text{ext}}^{21} = |(f_{\text{ext}}^{21} - f_1)/f_{\text{ext}}^{21}|$
• $\text{GCI}_{\text{fine}}^{21} = \frac{1.25 \, e_a^{21}}{r_{21}^{p} - 1}$

$p$の式が暗黙的（$p$が両辺にある）ですけど、どうやって解くんですか？

反復で解く。$p$の初期値に簡略版の値を入れて、$q(p)$を更新して$p$を再計算、というのを数回繰り返せば収束する。Pythonなら scipy.optimize.fsolve で一発だ。等比細分化（$r_{21} = r_{32}$）なら$q=0$だから反復不要、そのまま簡略版が使えるよ。

y+と乱流モデルの整合性

y+がメッシュ独立性確認とどう関係するのか、ピンと来ていないんですが…。

いい質問だ。y+は壁面第1層セルの無次元距離で、こう定義される：

$y$は壁面からの第1層セル中心までの距離、$u_\tau$は摩擦速度、$\nu$は動粘性係数だ。ポイントは、乱流モデルごとに要求されるy+の範囲が決まっていること：

• 壁関数使用（$k$-$\varepsilon$標準など）：$y^+ = 30\text{〜}300$
• 壁面解像（$k$-$\omega$ SST、SA、LES）：$y^+ < 1$（理想は$y^+ \approx 1$）
• ブレンド壁面処理（SSTのautomatic wall treatment）：$y^+ < 5$が望ましい

もし壁関数モードなのにy+が5くらいだったら、メッシュ独立性確認で何が起こるんですか？

壁関数は対数則（$y^+ > 30$の領域）を前提にしているから、$y^+ = 5$は緩衝層（buffer layer）の真ん中で、対数則も粘性底層の線形則も成り立たない「物理的に不整合」なゾーンだ。そこでは壁面せん断応力が大きく間違う。メッシュを細かくするほどy+が下がって壁関数の適用範囲から外れ、解が発散的な挙動を示したり、非単調な収束になったりする。結果としてGCIの前提（単調収束）が崩れるんだ。

つまりメッシュ独立性確認の「前提条件」として、まずy+が乱流モデルの想定範囲にあることを確認しないといけないんですね。

その通り。現場で多いのは「壁関数で計算したけど一部の壁でy+が3くらいになっていて、それに気づかずメッシュスタディをやって『収束しない！』と悩む」パターンだね。体温計をちゃんと挟まずに測って「おかしい」と言っているようなものだ。

数値解法と実装

代表セルサイズの定義

先生、代表セルサイズ$h$って、具体的にはどう計算するんですか？ 非構造格子だとセルのサイズがバラバラですよね。

Celikの手順ではグローバル平均を使う。3Dなら全セルの平均体積の立方根：

実務的には「総セル数$N$が分かれば、計算領域の全体積$V_{\text{total}}$から$h = (V_{\text{total}}/N)^{1/3}$で近似的に求められる」ので、ソルバーの出力にセル数さえあればすぐ計算できるよ。

局所的な評価量（壁面摩擦など）を見ているのに、グローバル平均で$h$を定義していいんですか？

鋭い指摘だ。厳密には、評価量が依存する領域のメッシュを均一に細分化している場合にのみグローバル$h$が適切になる。例えば壁面近傍だけ細かくして遠方は粗いままだと、細分化比$r$が実効的に変わってしまう。ベストプラクティスは全領域を相似的に細分化する（ジオメトリの比率を保ったまま全セル数を均一に増やす）ことだ。FluentならMesh > Controls > Sizingの全パラメータを一括でスケーリングする。

観測収束次数の計算

実際のCFD問題で$p$はどのくらいの値になるんですか？ 2次精度スキームなら$p=2$？

理想的にはそうだけど、実務では$p = 0.5$〜$3.0$くらいの範囲に散らばることが多い。非線形問題では$p < 2$になるケースが多いんだ。目安としては：

• $p \approx 1.5$〜$2.5$：正常。2次精度の漸近域にいると判断できる
• $p < 0.5$：メッシュがまだ粗すぎて漸近域に入っていない、またはバグがある
• $p > 3.0$：偶然の相殺が起きている可能性。もう1水準追加して確認すべき
• $p < 0$（$\varepsilon_{32}/\varepsilon_{21} < 0$）：非単調収束。標準手順が適用できない

$p$が1.8くらいだったとして、それでGCIは計算していいんですか？

全く問題ない。Celik手順の利点は、$p$を理論値に固定せず、観測値を使うことで、非線形問題でもGCIが信頼できる誤差推定になる点だ。$p=2$を仮定して計算するより遥かに正直な推定になるよ。

外挿値とGCI算出

具体的な数字で見てみたいです。例えば管内流の圧力損失で、3水準のメッシュの結果が$\Delta P_1 = 245.3$Pa、$\Delta P_2 = 248.1$Pa、$\Delta P_3 = 258.7$Paだとしたら？

等比細分化$r = 2.0$の場合を計算してみよう。

• $\varepsilon_{21} = f_2 - f_1 = 248.1 - 245.3 = 2.8$
• $\varepsilon_{32} = f_3 - f_2 = 258.7 - 248.1 = 10.6$
• $p = \ln(10.6/2.8)/\ln(2.0) = \ln(3.786)/\ln(2.0) = 1.331/0.693 = 1.92$
• $f_{\text{ext}}^{21} = (2.0^{1.92} \times 245.3 - 248.1)/(2.0^{1.92} - 1) = (3.784 \times 245.3 - 248.1)/2.784 = 244.3$ Pa
• $e_a^{21} = |245.3 - 248.1|/245.3 = 1.14\%$
• $\text{GCI}_{\text{fine}}^{21} = 1.25 \times 0.0114 / (2.0^{1.92} - 1) = 0.01425/2.784 = 0.51\%$

GCI 0.51%！ かなり良さそうですね。これで「メッシュ独立」と言えますか？

グローバルな圧力損失としてはGCI < 1%で十分だ。ただし複数の評価量で同時に確認することが大事だよ。圧力損失が収束していても壁面摩擦係数や局所的な剥離点位置がまだ動いているかもしれない。最低でも「グローバル量（力係数や流量）」と「局所量（壁面y+、特定断面の速度分布）」の両方でGCIを報告すべきだ。

収束判定の落とし穴

メッシュ独立性確認で見落としがちなポイントってありますか？

よくある落とし穴を3つ挙げよう：

• 反復残差の収束不足：粗メッシュでは残差$10^{-4}$で収束しても、細メッシュでは同じレベルに達するまで反復数が増える。全メッシュで同じ残差レベル（少なくとも$10^{-5}$）まで収束させないと、反復誤差が離散化誤差を汚染する。
• 非定常効果の無視：「定常解析でやったけど残差が下がりきらない」場合、物理的に非定常な流れ（カルマン渦など）の可能性がある。この場合は非定常解析に切り替えて、時間平均値でGCIを計算する必要がある。
• 境界層メッシュの不整合：全体を細分化しても、インフレーション層の成長率や層数を変えないとy+が変わってしまう。壁面近傍は第1層高さ$\Delta y_1$を指定して成長率で制御するため、全体の細分化比$r$と壁面近傍の細分化比が乖離することがある。

実践ガイド

5段階実践フロー

実際の業務では、どういう手順でやればいいですか？

実務で使える5ステップをまとめよう：

1. 評価量の決定：何を見て「メッシュ独立」と判断するかを最初に決める（抗力係数、出口温度、圧力損失など）。少なくともグローバル量1つ＋局所量1つ。
2. ベースメッシュの作成：まず中程度のメッシュで計算を回し、y+の分布を確認する。乱流モデルの要件を満たしていなければ先にインフレーション層を修正する。
3. 3水準メッシュの作成：ベースメッシュを基準に、$r \approx 1.3$〜$2.0$の範囲で粗・細メッシュを作る。全領域を相似的に細分化すること。
4. 3水準すべてを十分に収束させて解く：残差を$10^{-5}$以下に。評価量のモニターが十分にフラットになっていることも確認。
5. Celik手順でGCIを算出・報告：$p$、$f_{\text{ext}}$、$e_a$、$\text{GCI}$をテーブルにまとめる。

評価量の選び方

評価量って、抗力係数みたいなグローバル量だけじゃダメなんですか？

ダメではないけど、不十分だ。例えば円柱まわりの流れで抗力係数$C_d$が収束していても、後流の渦構造や剥離点位置がまだ変動しているかもしれない。推奨される評価量の組み合わせを示そう：

メッシュ細分化比の設計

$r = 1.3$と$r = 2.0$だと、計算コストがかなり違いますよね？

3Dだとセル数は$r^3$倍になるから：

• $r = 1.3$：セル数は$1.3^3 = 2.2$倍。3水準で粗の$2.2^2 = 4.8$倍
• $r = 1.5$：セル数は$1.5^3 = 3.4$倍。3水準で粗の$3.4^2 = 11.4$倍
• $r = 2.0$：セル数は$2.0^3 = 8$倍。3水準で粗の$8^2 = 64$倍

Celikは$r \geq 1.3$を推奨している。$r < 1.3$だとメッシュ間の差が小さすぎて丸め誤差に埋もれるリスクがある。$r = 1.5$が「信頼性とコストのバランス」として一番人気だね。

Fluent・STAR-CCM+・OpenFOAMでの実装

各ソフトでメッシュ独立性確認をやるときのコツってありますか？

ソフトごとの実践ポイントをまとめよう。

Ansys Fluent：Mesh > Adapt > Uniformで均一細分化が可能。ただしこれは六面体メッシュを各方向2分割するため$r=2$固定。非構造メッシュではMeshing側でサイズパラメータを一括スケーリングする方が確実。Fluent 2024R2以降はAdjoint Solverベースのメッシュ適応機能もあるが、GCI計算には手動の均一細分化が必要。

Simcenter STAR-CCM+：Parts > Mesh Operationsで「Base Size」を変えるだけで全メッシュが相似的に再生成される。これがSTAR-CCM+の最大の利点で、3水準のメッシュ作成が非常に楽。さらにJavaマクロで「ベースサイズを変えて再メッシュ → 解く → 結果取得」を自動化できる。GCIまで一気通貫で計算するマクロを書いている人もいるよ。

OpenFOAM：blockMeshを使っているなら、blockMeshDictのセル数指定を変えるだけで細分化比を制御できる。snappyHexMeshの場合はmaxLocalCells、maxGlobalCells、refinementLevelを調整する。GCI計算はPythonスクリプトで行うのが一般的で、postProcessingディレクトリから評価量を読み取って$p$とGCIを自動計算できる。

品質保証チェックリスト

報告書やレビューで「ちゃんとメッシュ独立性確認やった？」と聞かれたときに備えて、チェックリストが欲しいです。

• 3水準以上のメッシュで解を得たか
• 全メッシュで同じ離散化スキーム・同じ乱流モデルを使ったか
• 全メッシュで反復残差が$10^{-5}$以下に達しているか
• 全メッシュでy+が乱流モデルの要求範囲に入っているか
• 細分化比$r \geq 1.3$を満たしているか
• 観測収束次数$p$を計算し、理論値との比較を報告したか
• GCIをグローバル量と局所量の両方で報告したか
• 非単調収束の場合は対処を明記したか
• $\text{GCI}_{\text{coarse}}/(\text{GCI}_{\text{fine}} \times r^p) \approx 1$（漸近域チェック）を確認したか

ソフトウェア比較

自動メッシュ独立性確認機能の比較

メッシュ独立性確認を自動でやってくれるソフトってあるんですか？

自動GCI計算まで一体化したツールはまだ少ないけど、各ソフトのメッシュ細分化機能と自動化の容易さを比較しよう：

スクリプト化の容易さ

GCI計算まで自動化するなら、どのソフトが一番楽ですか？

個人的にはOpenFOAMが一番楽だと思う。理由は3つ：

1. テキストベースの辞書ファイルだからsedやPythonで一括置換できる
2. postProcessingに評価量が自動出力されるから結果の取得が簡単
3. シェルスクリプトでメッシュ生成→解析→結果取得→GCI計算まで一気通貫で書ける

STAR-CCM+はJavaマクロでGUIの全操作が自動化でき、特にBase Size変更→再メッシュ→再計算が1つのマクロで完結するのが強い。FluentはJournalファイルでTUIコマンドを自動実行できるけど、メッシュ再生成は別ツール（ANSYS Meshing）なのでワークフローが少し煩雑になりがちだ。

import numpy as np
f1, f2, f3 = 245.3, 248.1, 258.7  # 細→粗
r = 2.0
eps21, eps32 = f2 - f1, f3 - f2
p = np.log(abs(eps32/eps21)) / np.log(r)
f_ext = (r**p * f1 - f2) / (r**p - 1)
ea21 = abs((f1 - f2) / f1)
gci = 1.25 * ea21 / (r**p - 1)
print(f"p={p:.2f}, f_ext={f_ext:.1f}, GCI={gci*100:.2f}%")

先端技術

適応格子細分化（AMR）との関係

AMR（Adaptive Mesh Refinement）を使えば、メッシュ独立性確認は不要になりますか？

ならない。AMRは解のグラデーションが大きい領域を自動的に細分化するけど、「解がメッシュに依存していないか」を保証するものではないんだ。AMRは計算効率の向上手段であって、精度の保証手段ではない。AMRを使った場合でも、最終的な結果に対してGCIベースの確認を行うべきだよ。

でもAMRだとメッシュが非均一に細分化されますよね。Celik手順の前提と矛盾しませんか？

いいポイントだ。AMRの局所細分化とCelik手順の均一細分化は確かに矛盾する。実務的な解決策としては、(1) AMRのrefinement level上限を変えて3水準を作る、(2) 代表セルサイズ$h$の計算で重み付き平均を使う、(3) 評価量が局在する領域のメッシュ密度でGCIを計算する、といった工夫がある。どれも完全ではないけどね。

LES/DNS特有の課題

LES（Large Eddy Simulation）やDNS（Direct Numerical Simulation）でもGCIは使えますか？

ここがCFDのメッシュ独立性確認で最も難しい部分だ。LESではフィルタ幅（≒メッシュサイズ）が物理モデルの一部になっている。つまりメッシュを変えると「解いている方程式自体が変わる」んだ。だからRANSとは根本的に事情が異なる。

LESでのメッシュ感度評価のアプローチとしては：

• Pope criterion：解像されている乱流運動エネルギーの割合が80%以上であることを確認する（$M = k_{\text{resolved}}/(k_{\text{resolved}} + k_{\text{sgs}}) > 0.8$）
• エネルギースペクトル：$-5/3$乗則の慣性小領域がメッシュのカットオフ周波数より十分手前で観測されるか
• 2点相関：空間相関長がセルサイズより十分大きいか

DNSの場合はKolmogorov長スケール$\eta$以下のメッシュを全領域で使えば理論上メッシュ独立だが、コストが$Re^{9/4}$に比例するため、高Re数では事実上不可能だ。

機械学習による格子感度予測

最近はAIで「このメッシュで足りるか」を予測する研究とかあるんですか？

ある。主な方向性は2つだ。1つ目は、粗メッシュの解と局所的なメッシュ特徴量（セルサイズ、アスペクト比、y+等）から、細メッシュでの解の変化量を予測するサロゲートモデル。これは「GCIを計算するための細メッシュ計算コストを削減する」アプローチだ。2つ目は、過去の格子収束スタディのデータベースから、類似ジオメトリ・類似流れ条件での必要メッシュ密度を推定する手法。ただし、どちらもまだ研究段階で、実務で「GCI計算の代替」として使えるレベルには達していない。当面はCelik手順が標準だと思っていい。

トラブルシューティング

収束次数$p$が異常に低い（$p < 0.5$）

計算したら$p = 0.3$とかになったんですけど、何がおかしいんでしょうか？

$p$が異常に低い場合の原因と対策：

• 漸近域に入っていない：粗メッシュが粗すぎる。ベースメッシュを細かくしてからやり直す。
• 反復残差が不十分：細メッシュの残差が$10^{-3}$程度で止まっていないか確認。残差を$10^{-6}$まで下げてみる。
• y+の不整合：メッシュ間でy+が壁関数の有効範囲から外れていないか。3水準すべてのy+ヒストグラムを確認。
• 離散化スキームの不一致：風上差分（1次精度）を使っていないか。全メッシュで2次精度以上の同じスキームを使用する。
• 特異点を含む評価量：角部や再付着点での圧力など、特異点を含む位置のデータを評価量にしている場合、$p$が低くなる。評価点を特異点から離す。

GCIが大きすぎる

GCIが15%とか出ちゃったんですけど、これは「メッシュ不足」ってことですか？

そう、現在の細メッシュでも離散化誤差が大きいことを意味する。対策は：

• メッシュをさらに細かくする：計算コストが許すなら最もストレートな方法。
• 高次精度スキームに変更：1次風上を2次に、2次をMUSCL/QUICK等に上げる。$p$が上がればGCIは劇的に下がる。
• 重要領域に集中的にメッシュを投入：壁面近傍やウェイク領域など、評価量に効くエリアを重点的に細分化。
• 評価量の変更を検討：局所的なピーク値（最大壁面せん断応力など）はGCIが大きくなりやすい。面積平均値に変更できないか検討。

非単調収束への対処

3水準の解が$f_1 = 10.2$、$f_2 = 10.5$、$f_3 = 10.3$みたいに振動しているんですが…。

$\varepsilon_{32}/\varepsilon_{21} < 0$だから非単調収束だね。$\ln$の中が負になるので$p$が計算できない。これはCelik手順の適用外ケースだ。対処法は：

1. 反復残差を確認：本当に収束しているか。定常解析で残差が振動していないか。
2. 物理的に非定常な流れでないか確認：カルマン渦やフラッタリングがある場合、定常解は存在しない。非定常解析に切り替えて時間平均値を使う。
3. 4水準に増やす：もう1つメッシュを追加して、3つの隣接ペアのうち単調収束するペアでGCIを計算する。
4. 最悪のケースで誤差を報告：$p$が求まらない場合は、$e_a^{21}$（相対近似誤差）をそのまま「離散化誤差の保守的推定」として報告する。AIAA-G-077はこのアプローチを許容している。

メッシュ独立性確認の全体像がよく分かりました。まずはy+を確認してから3水準で回して、Celik手順でGCIを報告する——この流れですね。

その通り。最後にもう1つだけ。メッシュ独立性確認はVerificationの一部であって、Validationの代わりにはならない。「方程式を正しく解けている」ことと「正しい方程式を解いている」ことは別問題だ。両方をやって初めて信頼できるCFD結果になるよ。