メッシュ細密化比とGCI収束確認の手法
理論と物理
概要 — メッシュ収束確認とは
先生、メッシュ収束確認って何段階のメッシュで比較すればいいんですか? 「3つ試せ」って言われたんですが、なぜ2つじゃダメなんでしょう。
最低3段階だ。メッシュを2つ比較しただけだと「結果が近いから大丈夫」としか言えない。3つあれば収束の速度(次数 p)を推定でき、メッシュサイズ→0の極限解を外挿で予測できる。これがRichardson外挿であり、その誤差バンドを定量化するのがGCI(Grid Convergence Index)だ。
なるほど、2点じゃ直線しか引けないけど、3点あれば曲線のカーブも分かるって感じですね。
そういうこと。ASME V&V 20(計算流体力学・固体力学の検証・妥当性確認ガイドライン)でも、この3メッシュ法が標準手順として規定されている。論文やプロジェクトの報告書で「メッシュ依存性はありません」と書くには、GCIの数値を示すのが国際的な共通言語になっているんだ。
メッシュ収束確認(Mesh Convergence Study)は、離散化誤差を定量的に評価し、解の信頼性を保証するV&Vプロセスの根幹をなす。その核心は以下の3要素に集約される。
- 細密化比 r — 粗メッシュと細メッシュの代表サイズの比。$r \geq 1.3$(推奨 $r = \sqrt{2} \approx 1.41$)
- Richardson外挿 — 3つ以上のメッシュ結果からメッシュサイズ→0の極限解を推定する
- GCI — 離散化不確かさを95%信頼区間として数値化する
メッシュ細密化比 r の定義
細密化比って、メッシュをどれくらい細かくするかの「倍率」ですよね? $r = 2$ にすれば要素を半分のサイズにできるから簡単そうですけど。
構造格子なら $r = 2$ は確かに簡単だけど、要素数が3D で 8倍になる。計算コストが爆発するんだ。だからASME V&V 20では $r = \sqrt{2} \approx 1.41$ を推奨している。要素数は約2.8倍で済む。逆に $r < 1.3$ だと、メッシュ間の差が小さすぎて丸め誤差に埋もれてしまい、観測次数 $p$ の推定が不安定になる。
メッシュ細密化比(grid refinement ratio)の定義は次の通りである。
ここで $h$ は代表メッシュサイズ(representative grid size)である。常に $r > 1$ が成り立つよう、分子を粗メッシュ、分母を細メッシュとする。
r の推奨範囲と根拠
| $r$ の値 | 3D要素数の増加率 | 判定 |
|---|---|---|
| $r = 2.0$ | $\times 8$ | 理想的だが計算コスト大 |
| $r = \sqrt{2} \approx 1.41$ | $\times 2.8$ | ASME V&V 20推奨 |
| $r = 1.3$ | $\times 2.2$ | 最低ライン |
| $r < 1.3$ | — | 非推奨($p$推定が不安定) |
非構造メッシュの代表サイズ
構造格子なら要素サイズが均一だから $h$ は明確ですけど、テトラメッシュみたいな非構造メッシュだとどうやって $h$ を定義するんですか?
一番標準的なのは、全要素の体積の平均から等価長さを逆算する方法だ。$d$ は空間次元数で、2Dなら面積、3Dなら体積を使う。
$N$ は要素数、$V_i$ は各要素の体積(2Dの場合は面積)、$d$ は空間次元数(2 or 3)。この定義により、非構造メッシュでもスカラー値として代表サイズを得ることができる。
じゃあ非構造メッシュの場合、細密化比 $r$ は要素数の比から計算できるんですか?
そうだ。等方的に細密化した場合、要素数比と $r$ の関係はこうなる。
例えば3Dモデルで $N_{\text{fine}} = 800{,}000$、$N_{\text{coarse}} = 200{,}000$ なら、$r = (800{,}000 / 200{,}000)^{1/3} = 4^{1/3} \approx 1.587$ となり、十分な細密化比が確保されている。
Richardson外挿
Richardson外挿って、要は「メッシュを無限に細かくしたら答えはいくつになるか」を予測する方法ですよね? 具体的にはどう計算するんですか?
離散化された解 $f_h$ は、真の解 $f_{\text{exact}}$ に対して次のTaylor展開で表せる。
ここで $p$ は離散化スキームの理論次数、$g_p$ は未知係数である。高次項を無視すると、3つのメッシュ($h_1 < h_2 < h_3$, $f_1, f_2, f_3$)から外挿値が得られる。
ここで $r_{21} = h_2 / h_1$ であり、$\hat{p}$ は3つの解から推定した観測次数(後述)である。
つまり最も細かいメッシュの結果 $f_1$ をさらに補正して、真の解に近づけるわけですね。
その通り。ただし高次項を切り捨てているので、外挿値は「推定」であって真値ではない。だからこそ不確かさをGCIで定量化するんだ。
観測次数 p の算出
等比の細密化比($r_{21} = r_{32} = r$)の場合、観測次数 $\hat{p}$ は次式で直接計算できる。
$f_1$: 細メッシュの解、$f_2$: 中メッシュの解、$f_3$: 粗メッシュの解。$\hat{p}$ が理論次数 $p_{\text{formal}}$(例:2次要素なら2)に近ければ、解は漸近範囲(asymptotic range)に入っていると判断できる。
もし $r_{21} \neq r_{32}$(等比じゃない)場合はどうなりますか? 実務だとメッシャーの都合でぴったり等比にならないことが多いと思うんですが。
いい質問だね。その場合は陽的に解けないので、次の非線形方程式を反復法(固定点反復やニュートン法)で解く。
ここで $\varepsilon_{32} = f_3 - f_2$、$\varepsilon_{21} = f_2 - f_1$、$s = \text{sign}(\varepsilon_{32}/\varepsilon_{21})$ である。通常10回以内の反復で収束する。
GCI(Grid Convergence Index)
GCIの式を見ると安全係数 $F_s$ がありますよね。これはどういう意味ですか?
$F_s$ は信頼度を上げるための「ふかし係数」だ。Roacheの原論文では2メッシュ法で $F_s = 3$(保守的)、3メッシュ法で $F_s = 1.25$ を推奨している。3メッシュ法なら観測次数 $\hat{p}$ で外挿精度の整合性を確認済みだから、安全係数を小さくできるわけだ。
細メッシュ側のGCI($\text{GCI}_{21}^{\text{fine}}$)は次式で定義される。
ここで $\varepsilon$ は相対誤差である。
| 方法 | 安全係数 $F_s$ | 意味 |
|---|---|---|
| 2メッシュ法 | 3.0 | $\hat{p}$ 未知のため保守的 |
| 3メッシュ法 | 1.25 | $\hat{p}$ 推定済みで95%信頼区間相当 |
具体例で計算してみたいです。例えば3つのメッシュで最大応力がそれぞれ 248.3 MPa、251.6 MPa、260.1 MPa だったら?
$f_1 = 248.3$(細)、$f_2 = 251.6$(中)、$f_3 = 260.1$(粗)、$r = \sqrt{2}$ とする。
$\varepsilon_{21} = 251.6 - 248.3 = 3.3$、$\varepsilon_{32} = 260.1 - 251.6 = 8.5$
$\hat{p} = \ln(8.5/3.3) / \ln(\sqrt{2}) = \ln(2.576)/0.3466 \approx 2.73$
$\varepsilon = (251.6 - 248.3)/248.3 = 0.01329$
$\text{GCI}_{21}^{\text{fine}} = 1.25 \times 0.01329 / (1.414^{2.73} - 1) = 0.01661 / 1.634 \approx 0.0102 = 1.02\%$
つまり細メッシュの結果 248.3 MPa に対して $\pm 1.0\%$ の離散化不確かさがある、と報告できる。
漸近範囲チェック
漸近範囲に入っているかどうかは、以下の指標で確認する。
この値が1に近い(目安: $0.95 \leq \text{AR} \leq 1.05$)なら、解は漸近範囲にあり、GCIとRichardson外挿の信頼性が高い。1から大きく外れる場合は、メッシュがまだ十分に細かくないか、モデルに特異点が含まれている可能性がある。
漸近範囲に入っていない場合はどうすれば?
選択肢は3つだ。(1) さらに細かいメッシュを追加して再計算、(2) 評価点を特異点から離す(例:応力集中点ではなく断面平均応力を使う)、(3) より高次の要素を使って漸近範囲に早く到達させる。実務では(2)が最もコストパフォーマンスが良いことが多い。
数値解法と実装
3メッシュ法の手順
ASME V&V 20に基づく3メッシュ法の標準手順を以下に示す。
- 3水準のメッシュを生成: 細密化比 $r \geq 1.3$(推奨 $\sqrt{2}$)で粗・中・細の3段階。メッシュトポロジーは共通に保つ
- 同一条件で解析実行: 境界条件・材料・ソルバー設定を完全に統一する
- 評価量を抽出: 着目する物理量(最大応力、平均温度、圧力損失など)の値 $f_1, f_2, f_3$ を記録
- 代表メッシュサイズ $h_1, h_2, h_3$ を計算: 構造格子ならセルサイズ、非構造なら前述の体積平均法
- 観測次数 $\hat{p}$ を計算
- Richardson外挿値 $f_{\text{ext}}$ を計算
- $\text{GCI}_{21}^{\text{fine}}$ と $\text{GCI}_{32}^{\text{coarse}}$ を計算
- 漸近範囲チェック AR を計算: 1に近ければOK
手順は分かりました。でも構造格子と非構造格子で「同じトポロジーを保つ」って、非構造メッシュだと難しくないですか?
確かに非構造メッシュでは厳密に同一トポロジーにはできない。だから「グローバルメッシュサイズを一律に $1/r$ 倍にする」というアプローチが現実的だ。ローカルリファインメントは3水準とも同じ領域・同じ比率で適用し、メッシュ生成の「レシピ」を統一する。
非均一細密化比への拡張
$r_{21} \neq r_{32}$ の場合、観測次数 $\hat{p}$ は前述の非線形方程式を反復法で解く必要がある。実務的には以下の固定点反復で十分な精度が得られる。
初期値 $\hat{p}_0$ は等比の場合の公式で計算し、$|\hat{p}_{k+1} - \hat{p}_k| < 10^{-6}$ まで反復する。
複数評価量の扱い
最大応力ではGCIが1%だけど、たわみでは5%みたいなケースもあり得ますよね? どっちの結果を報告すべきですか?
全部報告するのが正解だ。GCIは評価量ごとに独立に計算するもので、「モデル全体でGCI 1%」とは言えない。局所量(応力集中点の応力)と全体量(合計反力、最大たわみ)では収束性が全然違う。局所量ほど収束が遅いのが普通だ。設計判断に直結する量を最低3つは評価して、それぞれのGCIを報告するのがベストプラクティスだよ。
実践ガイド
メッシュ収束確認のワークフロー
| ステップ | 作業内容 | 判定基準 |
|---|---|---|
| 1. メッシュ生成 | 粗・中・細の3水準($r \geq 1.3$) | $r$ の確認、メッシュ品質指標(アスペクト比 < 5等) |
| 2. 解析実行 | 同一境界条件・ソルバー設定で3ケース | 全ケースが正常収束 |
| 3. $\hat{p}$ 計算 | 観測次数の算出 | $0 < \hat{p} \leq p_{\text{formal}} + 1$ |
| 4. GCI計算 | $\text{GCI}_{21}^{\text{fine}}$ の算出 | 設計マージンとの整合性 |
| 5. AR チェック | 漸近範囲の確認 | $0.95 \leq \text{AR} \leq 1.05$ |
| 6. 報告 | GCI値とメッシュ情報を記録 | 第三者が再現可能な詳細度 |
よくある失敗と対策
先輩から「メッシュ収束確認やったけどGCIが50%になった」って相談されたんですが、何が原因でしょう?
GCIが大きすぎる場合の原因トップ3はこれだ。
- 特異点での評価 — 応力集中点(ノッチ先端、荷重点)では理論解が無限大に発散するため、メッシュを細かくしても永久に収束しない。対策:特異点から離れた位置で評価するか、断面平均応力を使う
- 細密化比が小さすぎる — $r < 1.1$ だと丸め誤差と離散化誤差の区別がつかない。対策:$r \geq 1.3$ を確保
- 解析が未収束 — 非線形解析でソルバーの収束判定が甘いと、各メッシュの解自体が不正確になる。対策:残差を十分小さく設定する
特異点の話、実務ではめちゃくちゃ重要ですね。角部の応力は無限大だから、どれだけメッシュを細かくしても値が上がり続ける…。
そう。CAEの初心者がよく陥るのが「最大応力で収束確認」をしてしまうケースだ。特にフィレットのないシャープエッジのモデルだと、最大応力はメッシュを細かくするたびに単調増加して収束しない。評価量の選び方がメッシュ収束確認の成否を左右するんだ。
報告書への記載事項
メッシュ収束確認の結果を報告する際は、以下の情報を必ず含める。
- 各メッシュの要素数・代表サイズ $h$
- 細密化比 $r_{21}$, $r_{32}$
- 評価量(物理量名、位置、方向)
- 各メッシュの評価量の値 $f_1, f_2, f_3$
- 観測次数 $\hat{p}$
- Richardson外挿値 $f_{\text{ext}}$
- $\text{GCI}_{21}^{\text{fine}}$(パーセント表示)
- 漸近範囲チェック AR
- メッシュ品質指標(最小アスペクト比、最大スキューネスなど)
ソフトウェア別の対応
主要CAEツールでのメッシュ収束機能
| ツール | 自動メッシュ収束 | GCI計算 | 備考 |
|---|---|---|---|
| Ansys Mechanical | Convergence Tool(適応メッシュ細密化) | 手動 or ACT Extension | Workbench内でメッシュサイズのパラメトリック掃引が可能 |
| Abaqus | Python scripting + Adaptive Remeshing | 手動(Pythonスクリプト) | abaqus.rpy でメッシュサイズをパラメータ化して自動化 |
| COMSOL Multiphysics | Parametric Sweep + Adaptive Mesh | 手動 | メッシュサイズをパラメータとしたスイープが容易 |
| OpenFOAM | snappyHexMesh のレベル指定 | 手動 | pyFoam等でGCI計算を自動化するスクリプトが公開されている |
| Ansys Fluent | Adaptive Mesh Refinement (AMR) | 手動 or Journal | 解適応型メッシュ細密化と組み合わせ可能 |
| Simcenter STAR-CCM+ | Mesh Pipeline + Design Manager | 手動 | Design Managerでメッシュサイズの自動掃引が可能 |
どのソフトも「GCI計算は手動」なんですね。自動でやってくれるツールはないんですか?
GCIの計算自体は式に値を代入するだけだから、Excelやスクリプトで十分対応できる。むしろ重要なのは「メッシュの生成レシピを統一して3ケースを効率よく回す仕組み」の方で、これはパラメトリックスタディ機能があるツールが有利だ。Ansys WorkbenchのDesign Points機能やCOMSOLのParametric Sweepを使えば、メッシュサイズをパラメータにして一括実行できるよ。
先端技術
適応メッシュとGCIの融合
解適応型メッシュ細密化(Adaptive Mesh Refinement: AMR)は、誤差推定量に基づいて局所的にメッシュを細密化する手法であり、GCIとは本質的に異なるアプローチである。しかし近年、両者を組み合わせた手法が提案されている。
AMRだと局所的にメッシュが異なるので、「代表サイズ $h$ の定義」が難しそうですね。
その通りで、AMRを使う場合はGCIの適用にいくつか制約がある。一つのアプローチは「AMRの最大レベルを固定して、ベースメッシュのサイズだけを $r$ 倍に変えて3水準作る」方法だ。これならGCIの枠組みをそのまま使える。もう一つは、AMR+事後誤差推定量(ZZ推定量やSPR法)を使って、GCIを使わずに離散化誤差を直接推定する方法だ。
最小二乗GCI(多メッシュ法)
4つ以上のメッシュが利用可能な場合、$\hat{p}$ と $g_p$ を最小二乗法で推定することで、外挿の安定性が大幅に向上する。特に観測次数が安定しない(漸近範囲に入っていない)場合に有効であり、Eqaら (2014) の研究で体系化されている。
$M \geq 4$ のメッシュ水準で上記を解くことで、3メッシュ法よりもロバストな $\hat{p}$ と $f_{\text{ext}}$ が得られる。計算コストは増えるが、高信頼性が要求される認証解析(certification analysis)で活用されている。
トラブルシューティング
観測次数が発散する
$\hat{p}$ が 8 とか 15 とか、理論次数をはるかに超える値になったんですが、これって何が起きてるんですか?
$\hat{p}$ が異常に大きい場合、$\varepsilon_{21}$ と $\varepsilon_{32}$ の比が1に非常に近いことを意味する。つまり中メッシュと細メッシュの差がほとんどないのに、粗メッシュとの差だけが大きい状態だ。原因は丸め誤差領域に入っているか、中と細のメッシュが実質的に同程度の解像度であることが多い。対策:$r$ を大きくするか、有効桁数を増やす(倍精度→4倍精度)。
p が負値になる
$\hat{p} < 0$ は $\varepsilon_{32}$ と $\varepsilon_{21}$ の符号が異なること(振動的収束)を示す。この場合、単調収束の仮定が崩れているためGCIは適用できない。
振動的収束が起きる典型例は、CFDで対流スキームの次数が不十分な場合や、応力の評価点がメッシュ境界に近い場合だ。評価点を変えるか、上位の離散化スキームに切り替えることで解決することが多い。
単調収束しない
細→中→粗の順に結果が単調に変化しないんですが、これだとRichardson外挿は使えませんよね?
使えない。非単調収束の場合は以下を確認しよう。
- メッシュの品質: 細メッシュにした結果、高アスペクト比やスライバー要素が増えていないか
- 境界条件の一貫性: メッシュ変更時に荷重面や拘束面の定義がずれていないか
- 物理モデルの整合性: 接触やタービュレンスモデルがメッシュサイズに依存していないか
- 丸め誤差の支配: 解の差が $10^{-12}$ オーダーなら倍精度の限界に達している
これらを全てクリアしても非単調なら、4メッシュ以上を使った最小二乗法アプローチに切り替えるのが安全だ。
関連トピック
なった
詳しく
報告