有限体積法の基礎
理論：保存則を「箱の収支」として書く
概要
CFDのソフトはほとんど有限体積法だと聞きました。構造解析で使う有限要素法と、何がそんなに違うんですか？
発想の出発点が違うんだ。有限要素法は微分方程式に重み関数を掛けて領域全体で積分する「弱形式」から出発する。有限体積法は、空間を小さな箱（検査体積、セル）に区切って、「箱の中身の増え方＝面から入ってくる量−出ていく量＋湧き出し」という収支をそのまま式にする。隣の箱とは同じ面を通る同じ量をやり取りするから、どこかで質量が勝手に消えたり増えたりしない。これが流体計算で好まれる最大の理由だよ。
積分形の輸送方程式
密度 $\rho$、速度 $\mathbf{u}$、拡散係数 $\Gamma$ の場で、任意の量 $\phi$（速度成分・温度・濃度など）の保存を検査体積 $V$ とその表面 $S$ について書くと次のようになる。
左から、時間変化・対流・拡散・生成項。セル $P$ について面 $f$ ごとの和に置き換えると、離散式は次の形になる（$N$ は面の向こうの隣接セル、$F_f$ は面を通る質量流量）。
式の形はシンプルですね。難しいのはどこなんですか？
未知数はセル中心の値なのに、式に出てくるのは面の値 $\phi_f$ と面の勾配なんだ。セル中心から面の値をどう補間するか——これが有限体積法の精度と安定性をほぼ決める。次の節からはその話が中心になる。
セル中心型と節点中心型
| 方式 | 未知数の位置 | 主な採用例 |
|---|---|---|
| セル中心型 | セルの中心。セル＝検査体積 | OpenFOAM, Ansys Fluent, Simcenter STAR-CCM+ |
| 節点中心型 | メッシュの節点。節点のまわりに双対セルを作る | Ansys CFX など |
熱と流れの計算を広めた一冊
有限体積法が工学の現場に広まったきっかけとしてよく挙げられるのが、S. V. Patankar の教科書『Numerical Heat Transfer and Fluid Flow』（1980年）だ。検査体積の収支から離散式を組み立てる手順と、圧力と速度を交互に修正する SIMPLE 法（Patankar と Spalding が1972年に発表）が、手で追える形で説明されていた。いま商用コードの内部で動いている考え方の多くは、この本の延長線上にある。
離散化：面の値をどう決めるか
拡散項：中心差分と非直交補正
拡散項の面勾配は、両側のセル中心の差 $(\phi_N-\phi_P)/|\mathbf{d}|$ で近似するのが基本。これはセル中心を結ぶベクトル $\mathbf{d}$ が面に垂直なときだけ正確で、斜めに交わる（非直交な）メッシュでは、セル中心の勾配を使った補正項を加える。補正項は陽的に扱うことが多く、非直交性が大きいと反復が不安定になりやすい。
対流項：スキームの選び方
| スキーム | 精度 | 性質 |
|---|---|---|
| 1次風上差分 | 1次 | 常に有界で安定。数値拡散が大きく、勾配がなまる |
| 中心差分 | 2次 | セルペクレ数が2を超えると振動する |
| 2次風上（線形風上） | 2次 | 上流側の勾配で外挿。実用の標準。わずかにオーバーシュートし得る |
| QUICK | 3次（一様格子） | 精度は高いが有界性は保証されない |
| TVD（リミッタ付き） | 2次（急変部で1次） | 急な変化の近くだけ風上に切り替えて振動を防ぐ |
セルペクレ数
対流と拡散のどちらが強いかは、セル幅 $\Delta x$ を代表長さとするペクレ数で判断できる。
中心差分では、隣のセルにかかる係数が $D - F/2$（$D=\Gamma/\Delta x$、$F=\rho u$）になる。$Pe_\Delta>2$ だとこれが負になって、「上流が上がると下流が下がる」という非物理的な関係が式に入る。これが振動の正体なんだ。
計算例：1次元の対流拡散
長さ $L=1$ m、$\rho=1$ kg/m³、$\Gamma=0.1$ kg/(m·s)、左端 $\phi=1$、右端 $\phi=0$ の定常対流拡散は厳密解をもつので、スキームの性質を確かめるのにちょうどいい。
流速 $u=2.5$ m/s（$Pe=25$）を5セルで解いた結果（セルペクレ数5）を示す。
| 位置 x [m] | 厳密解 | 中心差分 | 1次風上 |
|---|---|---|---|
| 0.1 | 1.0000 | 1.0356 | 0.9998 |
| 0.3 | 1.0000 | 0.8694 | 0.9987 |
| 0.5 | 1.0000 | 1.2573 | 0.9921 |
| 0.7 | 0.9994 | 0.3521 | 0.9524 |
| 0.9 | 0.9179 | 2.4644 | 0.7143 |
同じ条件でセル数を増やしたときの最大誤差は次のとおり。
| セル数 | セルペクレ数 | 中心差分 | 1次風上 |
|---|---|---|---|
| 5 | 5 | 1.547 | 0.204 |
| 10 | 2.5 | 0.537 | 0.158 |
| 20 | 1.25 | 0.160 | 0.120 |
| 40 | 0.625 | 0.044 | 0.079 |
| 80 | 0.31 | 0.012 | 0.048 |
| 160 | 0.16 | 0.003 | 0.026 |
5セルの中心差分、境界の値は0と1なのに2.46って……完全に壊れてますね。
そう、しかも収束計算としてはちゃんと解けているから、残差を見ても気づけない。一方で風上差分は振動しないけれど、セル数を2倍にしても誤差は約1.7倍しか減らない（1次精度）。中心差分はセルペクレ数が2を下回ると誤差が約4倍ずつ減る（2次精度）。40セル以上では中心差分のほうが正確になる。実務で2次の風上系スキームやTVDが標準なのは、この両方の欠点を避けたいからなんだ。
圧力と速度の連成
非圧縮流れでは圧力を直接決める方程式がない。そこで運動量式で仮の速度を求め、連続の式を満たすように圧力補正方程式を解いて速度と圧力を修正する、という反復を行う。定常計算では SIMPLE法（と改良版のSIMPLEC）、非定常計算では PISO法 がよく使われる。
- セル中心に速度と圧力を同居させる配置では、そのままだと圧力が市松模様に振動する。Rhie–Chow 補間で面の流束を作ることで防いでおり、主要な商用コードは内部でこれを行っている。
- SIMPLE 系では不足緩和が必要。圧力0.3・運動量0.7が広く使われる初期値で、発散するときは下げ、遅すぎるときは上げる。
実務：メッシュとソルバー設定
メッシュ品質の見方
| 指標 | 何が悪くなるか | 目安 |
|---|---|---|
| 非直交性 | 拡散項の精度と反復の安定性 | OpenFOAM の checkMesh は70°超を警告。補正の回数や制限を設定で調整する |
| 歪度（skewness） | 面の値の補間精度 | checkMesh は4超を警告 |
| 隣接セルの体積比 | 勾配の評価精度 | 急な寸法変化を避け、徐々に変える |
| アスペクト比 | 流れに沿わない方向では精度低下 | 境界層では大きくてよい（流れ方向に長いセル） |
しきい値はツールごとに異なるので、あくまで目安として使い、最終的には次の節の収束確認で判断する。
収束の判定
- 残差が3桁下がっただけで終わりにしない。圧力損失・揚力・出口温度など、知りたい量が反復で変わらなくなったことを確認する。
- 流入と流出の質量流量の差が、流量に対して十分小さい（例えば0.1%程度以下）ことを確認する。
- 最終結果は2次精度のスキームで出す。1次風上は立ち上げの安定化に使い、途中で切り替える。
主要ツールでの設定
| ツール | 対流スキームの指定 | 圧力・速度の連成 |
|---|---|---|
| OpenFOAM | system/fvSchemes の divSchemes に Gauss upwind、Gauss linearUpwind grad(U)、Gauss linear などを書く | ソルバーで選ぶ（simpleFoam＝SIMPLE、pimpleFoam＝PIMPLE） |
| Ansys Fluent | 空間離散化で First Order Upwind / Second Order Upwind / QUICK などを選ぶ | SIMPLE / SIMPLEC / PISO / Coupled |
| Simcenter STAR-CCM+ | 物理モデルで対流の離散化次数を選ぶ | 分離型（Segregated）と連成型（Coupled）のソルバー |
斜めの流れで現れる「偽の拡散」
1次風上差分の数値拡散は、流れがメッシュに対して斜め45°のときに最も大きくなることが知られている。拡散がまったくない流れで、斜めに走る鋭い境界面を計算しても、結果では境界がぼやけて広がってしまう。物理の拡散係数よりこの「偽の拡散」のほうが大きいと、結果は物理ではなくメッシュを表していることになる。混合や熱の広がりを評価する計算では特に注意が必要だ。
よくあるトラブルと検証
症状と原因
| 症状 | 考えられる原因 | 対策 |
|---|---|---|
| 境界値を超える値・縞模様 | セルペクレ数が大きいのに中心差分 | メッシュ細分化、または線形風上・TVDへ |
| 数十反復で発散 | 非直交性の高いセル、緩和係数が大きい | 該当セルを直す、緩和係数を下げる、1次で立ち上げる |
| 残差が一定値で止まる | 本当に非定常な流れ（はく離・渦放出）を定常で解いている | 監視量が振動しているなら非定常計算へ（CFDの時間積分法） |
| スキームを変えると結果が大きく変わる | メッシュが粗く、数値拡散が支配的 | メッシュ収束の確認を先に行う |
検証の進め方
- 上の1次元問題のように厳密解のある問題で、スキームの精度次数（メッシュ2倍で誤差が何分の1になるか）を確認する。
- 実際の形状では、少なくとも3段階のメッシュで注目量を比べ、Richardson 外挿や GCI（格子収束指数）で離散化誤差を見積もる。
- 実験値との比較は、離散化誤差を把握したうえで行う。メッシュ依存の結果を実験に合わせ込むと、別の条件で外れる。
自分でもペクレ数を変えて振動を再現してみたいです。
対流拡散方程式シミュレーターで流速と拡散係数を動かすと、中心差分が振動し始める境目を確かめられるよ。検査体積の組み立て方は1D有限体積法シミュレーター、できた三重対角の連立方程式の解き方はトーマス法（TDMA）シミュレーターで追える。スキームの詳しい比較は風上差分スキームの記事にまとめてある。
関連トピック
なった
詳しく
報告