流体力学の基礎 — ナビエ・ストークス方程式から乱流モデルまで完全解説
流体力学は、液体や気体の流れを支配する物理法則の体系です。CAEにおける流体解析（CFD: Computational Fluid Dynamics）は、自動車のエアロダイナミクス、航空機翼周りの流れ、熱交換器設計、建築物の風荷重評価、電子機器の冷却設計など、現代のエンジニアリングのあらゆる場面で活用されています。本稿では、流体力学の数学的基礎から、実務で不可欠な乱流モデルの選び方、CFDソルバーの内部アルゴリズムまでを体系的に解説します。
1. 流体の基本性質
連続体仮説と流体粒子
先生、流体って「水分子が集まったもの」ですよね。CFDって分子1個1個を計算してるわけじゃないですよね？
そうそう、そこ大事なポイント。流体力学は連続体仮説を使う。1 cm³の空気の中には約2.7×10¹⁹個もの分子があるわけだけど、工学的なスケールで見ると「連続的に分布する物質」として扱っても十分精度が出る。CFDのメッシュセル1つが1mm³でも、その中には分子が何兆個もいるからね。
ただ、マイクロ流体デバイス（チャネル幅が数μm以下）や真空に近い極低圧では連続体仮説が破綻する。そのときはクヌーセン数 $Kn = \lambda/L$（$\lambda$：平均自由行程, $L$：代表長さ）でチェックする。$Kn > 0.01$ くらいから注意が必要だ。
普通の自動車や飛行機の流れなら連続体でOKってことですね。「流体粒子」っていうのは何ですか？分子と違うんですよね？
流体粒子は「無限小だけど十分多くの分子を含む微小体積要素」のことだ。分子より大きく、流れの構造より小さいスケールを想像してほしい。この流体粒子に対して、位置・速度・圧力・温度などの状態量が定義できる。それが連続体力学の出発点になる。
粘性の物理的意味
粘性は「流体の層間で運動量を伝達する性質」です。ニュートン流体では、せん断応力 $\tau$ はせん断速度 $du/dy$ に比例します（ニュートンの粘性法則）：
ここで $\mu$ は動力粘性係数（粘度）[Pa·s]、$u$ は流速、$y$ は壁面からの距離です。また動粘性係数 $\nu = \mu/\rho$ [m²/s]（$\rho$：密度）もCFDでよく使われます。温度依存性が強く、液体では温度上昇で粘度は下がり、気体では上がります。
| 流体 | 動力粘性係数 μ [Pa·s]（20℃） | 動粘性係数 ν [m²/s] | 備考 |
|---|---|---|---|
| 空気 | 1.81×10⁻⁵ | 1.51×10⁻⁵ | 温度で大きく変化 |
| 水 | 1.00×10⁻³ | 1.00×10⁻⁶ | 基準値として使いやすい |
| 軽油（20℃） | 3.0×10⁻³ | 3.4×10⁻⁶ | 温度依存性が強い |
| グリセリン（20℃） | 1.49 | 1.19×10⁻³ | 水の1490倍の粘性 |
| エンジンオイル SAE30（40℃） | ~0.1 | ~1.1×10⁻⁴ | 高温で急激に低下 |
| 水銀（20℃） | 1.53×10⁻³ | 1.13×10⁻⁷ | 密度が大きいため ν が小さい |
圧縮性 vs 非圧縮性
空気って圧縮できますよね。でもCFDで「非圧縮性」として解く、っていうのを見るんですが、それっていつ使えるんですか？
実務の目安はマッハ数 $Ma < 0.3$。流速が音速（空気中約340 m/s）の30%未満なら、密度変化が1%未満なので非圧縮として扱える。時速100 kmの車（約28 m/s、$Ma \approx 0.08$）は余裕で非圧縮。逆にジェット戦闘機の巡航（$Ma \approx 0.85$）や超音速ロケットは圧縮性を無視したら全然ダメだ。
非圧縮性にすると何がうれしいかというと、エネルギー方程式を独立して解けるし、密度が定数になるから計算が格段に楽になる。SIMPLE法みたいな圧力ベースの解法が使えるのもそのおかげだ。
マッハ数は $Ma = V/a$（$V$：流速, $a = \sqrt{\gamma RT}$：音速）で定義されます。$\gamma$ は比熱比（空気：1.4）、$R$ は気体定数（空気：287 J/(kg·K)）、$T$ は絶対温度です。
| マッハ数域 | 名称 | 密度変化 | 推奨アプローチ | 代表例 |
|---|---|---|---|---|
| $Ma < 0.3$ | 非圧縮領域 | <2% | 非圧縮性ソルバー（SIMPLE等） | 自動車・建屋・低速ファン |
| $0.3 < Ma < 0.8$ | 亜音速圧縮性 | 2〜20% | 圧縮性ソルバー or 密度補正 | 旅客機巡航・高速ファン |
| $0.8 < Ma < 1.2$ | 遷音速 | 大きい・衝撃波の可能性 | 圧縮性ソルバー必須 | 翼型トランソニック設計 |
| $Ma > 1.2$ | 超音速 | 衝撃波・膨張波 | 密度ベースソルバー（Roe/AUSM） | ロケット・超音速機 |
ニュートン流体 vs 非ニュートン流体
$\tau = \mu \, du/dy$ が成立するのがニュートン流体（水・空気・多くの低分子有機溶媒）です。血液・ポリマー溶液・グリースなど多くの実用流体は非ニュートン挙動を示します。
| モデル | 構成則 | 代表例 | CFDでの扱い |
|---|---|---|---|
| ニュートン | $\tau = \mu \dot{\gamma}$ | 水・空気・エタノール | 標準NS方程式そのまま |
| べき乗則（Power Law） | $\tau = K \dot{\gamma}^n$（$n<1$：擬塑性，$n>1$：膨張性） | ケチャップ・ポリマー溶融体 | 有効粘度 $\mu = K\dot{\gamma}^{n-1}$ として代入 |
| Bingham塑性体 | $\tau = \tau_0 + \mu_p \dot{\gamma}$（$\tau > \tau_0$ のとき流動） | グリース・コンクリート・泥流 | 降伏応力 $\tau_0$ を考慮したモデル |
| Herschel-Bulkley | $\tau = \tau_0 + K\dot{\gamma}^n$ | 血液・ヘアジェル | Power Law + 降伏応力の複合 |
| Carreau-Yasuda | $\mu = \mu_\infty + (\mu_0-\mu_\infty)[1+(\lambda\dot{\gamma})^a]^{(n-1)/a}$ | 血液・歯磨き粉 | Ansys Fluent等で直接入力可 |
2. 連続の方程式（質量保存）
積分形と微分形の導出
任意の検査体積 $V$ に対して、質量保存則を適用します。検査体積内の質量変化率 = 流入質量フラックス − 流出質量フラックス：
ガウスの発散定理（$\oint_S \mathbf{F}\cdot\hat{n}\,dS = \int_V \nabla\cdot\mathbf{F}\,dV$）を適用すると微分形が得られます：
これを展開すると：$\frac{\partial \rho}{\partial t} + \rho\nabla\cdot\mathbf{v} + \mathbf{v}\cdot\nabla\rho = 0$、すなわち物質微分を使って：
非圧縮性流体（$\rho = \text{const}$、すなわち $D\rho/Dt = 0$）では：
$\nabla \cdot \mathbf{v} = 0$ っていう条件、よく出てきますよね。これ、具体的に何を言ってるんですか？
ざっくり言うと「どこにも流体が溜まらないし、どこからも湧き出ない」という条件だ。例えば断面積が変わる管を想像して。細くなったら流速が速くなる、太くなったら遅くなる。その関係を式で表したのが $A_1 v_1 = A_2 v_2$（連続の式の積分形）で、それの微分版が $\nabla \cdot \mathbf{v} = 0$ だ。
CFDでは、この条件が満たされていないと「質量が消えたり生えたり」してエラーになる。収束判定に使う残差の中で continuity residual が大きいときは、大抵メッシュか境界条件の設定ミスだよ。
管の断面積変化と流速の関係（実例）
直径 $D_1 = 50$ mm の管が $D_2 = 25$ mm に縮小する場合、入口流速 $v_1 = 2$ m/s のとき出口流速は：
断面積が1/4になると流速は4倍になります。ベルヌーイの定理（$p + \frac{1}{2}\rho v^2 + \rho g z = \text{const}$）と組み合わせると、その圧力低下量 $\Delta p = \frac{1}{2}\rho(v_2^2 - v_1^2) = \frac{1}{2} \times 1000 \times (64-4) = 30\,000$ Pa ≈ 0.3 atm も求められます。
3. ナビエ・ストークス方程式（運動量保存）
物質微分の概念
NS方程式に $D/Dt$ っていう記号が出てくるんですが、普通の偏微分 $\partial/\partial t$ と何が違うんですか？
これ、流体力学で最初に躓くポイントだよね。$\partial/\partial t$ は「ある固定した場所で時間変化を見る」オイラー描像。$D/Dt$ は「流体粒子に乗って動きながら時間変化を見る」ラグランジュ描像。
例えば川の水温を考えると：$\partial T/\partial t$ は「橋の上から同じ場所の水温変化を測る」で、$DT/Dt$ は「ゴムボートに乗って流されながら水温変化を感じる」イメージだ。この2つの関係が $\frac{D}{Dt} = \frac{\partial}{\partial t} + (\mathbf{v} \cdot \nabla)$ で、$(\mathbf{v} \cdot \nabla)$ の項が「移流」を表す。
じゃあ CFD のコードでは $D/Dt$ をどうやって離散化するんですか？
FVMでは $\partial/\partial t$ の部分を時間差分（後退差分・クランク・ニコルソン等）で離散化して、移流項 $(\mathbf{v} \cdot \nabla)\phi$ は対流フラックスとしてセル境界面でのフラックスを評価する。この移流フラックスの評価方法が「対流スキーム」で、1次風上・QUICK・リミター付き高次スキームなどがある。後でまた詳しく説明するよ。
NS方程式の導出と物理的意味
ニュートンの第2法則（$F = ma$）を流体粒子に適用します。流体粒子に作用する力は、(1) 圧力勾配力、(2) 粘性力（せん断応力の発散）、(3) 体積力（重力など）です。ニュートン流体の応力テンソル $\sigma_{ij}$ を使って：
これを運動量保存則に代入し、非圧縮性（$\nabla\cdot\mathbf{v}=0$）と定粘性係数を仮定すると、ナビエ・ストークス方程式（NS方程式）が得られます：
展開した形（x成分）：
各項の物理的意味を整理します：
| 項 | 数式 | 物理的意味 | CFDで注意する点 |
|---|---|---|---|
| 非定常項 | $\rho \partial\mathbf{v}/\partial t$ | 流体粒子の運動量の時間変化 | 定常解析ではゼロとして省略 |
| 慣性項（移流項） | $\rho(\mathbf{v} \cdot \nabla)\mathbf{v}$ | 流れによる運動量の輸送 | 非線形！数値不安定の主因 |
| 圧力勾配項 | $-\nabla p$ | 圧力差による駆動力 | 非圧縮では圧力方程式を別途解く |
| 粘性項 | $\mu \nabla^2 \mathbf{v}$ | 粘性による運動量拡散（滑らかにする） | 高 Re では小さい→数値拡散に注意 |
| 体積力 | $\rho \mathbf{g}$ | 重力・浮力（熱対流で重要） | Boussinesq近似か完全圧縮性で扱う |
慣性項が「非線形」って強調されてましたけど、それって何が問題なんですか？
非線形ってつまり $\mathbf{v}$ と $\mathbf{v}$ のかけ算項があるってこと。線形方程式なら重ね合わせが使えて解きやすいんだけど、非線形だとそうはいかない。この慣性項が乱流・渦・カルマン渦列みたいな複雑な現象を生み出す根本原因だ。CFDで反復計算（アンダーリラクセーション）が必要なのも、この非線形移流項の扱いのせいだよ。
逆に Re が非常に小さい（$Re \ll 1$）ストークス流れでは慣性項を無視できて「クリーピング流れ」になる。これは線形方程式だから厳密解が得られる。例えばマイクロ流路の生体流体解析や潤滑膜はこの近似が使えるよ。
無次元化とレイノルズ数の導出
代表長さ $L$、代表速度 $U$、動圧 $\rho U^2$ で無次元化（アスタリスク量は無次元）すると：
ここでレイノルズ数は唯一の無次元パラメータとして現れます：
無次元化の重要な示唆：同じ $Re$ 数の流れは幾何学的に相似なら流れパターンが相似になります。これが風洞実験の縮小模型試験が成立する根拠です。例えば実物の1/10模型で実物の10倍の流速を使えば同じ $Re$ 数が得られ、圧力係数などの無次元量は一致します。
4. エネルギー方程式
圧縮性流れや熱流体連成では、エネルギー方程式を連続式・運動量式と連立して解きます。全エンタルピー $H = h + V^2/2$（$h$：比エンタルピー, $V = |\mathbf{v}|$）を使った保存形：
非圧縮性・定物性・粘性散逸を無視した場合は、温度のみの移流拡散方程式になります：
ここで $c_p$ は定圧比熱 [J/(kg·K)]、$k$ は熱伝導率 [W/(m·K)] です。この方程式は温度場が速度場に依存する一方向的な連成（速度→温度のみ）ですが、自然対流やBoussinesq近似では密度が温度に依存するため双方向連成になります。
エネルギー方程式って、CFDでは毎回解くんですか？それとも「解かなくていい」ケースがある？
実務では「等温解析」として温度方程式をオフにするケースが多い。例えば自動車ボディのエアロダイナミクスだと「とりあえず風圧・ドラッグ係数だけ知りたい」なら非圧縮・等温でOKだ。でも電子機器の冷却（CPU排熱）、エンジン燃焼室、超音速流れは必ずエネルギー方程式が要る。
注意点は、Ansys FluentやOpenFOAMでは「Energy: OFF」を忘れると温度が全然無関係な値になってるのに気づかないことがある。密度が温度依存（理想気体 $\rho = p/(RT)$）のときは特に要注意だよ。自動車のエンジンルーム熱流動みたいな「高温ガス＋強制対流」を解くときは必ずエネルギー方程式をONにしておこう。
粘性散逸項の重要性
高速流れや高粘性流体では粘性散逸（$\Phi = \boldsymbol{\tau}:\nabla\mathbf{v}$）が重要になります。
ここで $e_{ij}$ はひずみ速度テンソルの成分です。高粘性ポリマーの押出成形や、マッハ数が高い超音速流れでは粘性散逸による温度上昇（空力加熱）を無視できません。大気圏再突入カプセルが高温になるのはこの現象です。
5. レイノルズ数と流れのスケール
層流・遷移・乱流の判定
| 流れの状態 | 円管内（水力直径 $D_h$） | 平板上（距離 $x$） | 特徴 |
|---|---|---|---|
| 層流 | $Re_D < 2300$ | $Re_x < 5 \times 10^5$ | 流線が整然、圧力損失∝流速の1乗 |
| 遷移流 | $2300 \sim 4000$ | $5\times10^5 \sim 10^6$ | 不安定、解析的な予測が困難 |
| 乱流 | $Re_D > 4000$ | $Re_x > 10^6$ | 乱れた渦運動、圧力損失∝流速の1.75〜2乗 |
実際の機器って、だいたいどのくらいの Re 数なんですか？感覚がつかめなくて。
実例を挙げると：
- 毛細血管（$D \approx 10\,\mu\text{m}$）: $Re \approx 0.001$ — 完全層流、慣性は無視できる Stokes 流れ
- シャワーヘッドのノズル: $Re \approx 500\sim2000$ — 層流〜遷移域
- 家庭の水道管（φ13 mm, 2 L/min の水）: $Re \approx 3000$ — 遷移〜乱流域
- 自動車ボディ（100 km/h, 代表長さ 4 m）: $Re \approx 1.1\times10^7$ — 完全乱流
- 大型旅客機の翼（巡航速度, 翼弦長 7 m）: $Re \approx 5\times10^8$ — 極めて高い Re 乱流
- パイプラインポンプの羽根車: $Re \approx 10^6\sim10^7$ — 乱流
実務で使うCFDのほとんどは乱流域だから、乱流モデルの選択が超重要になる。
境界層の発達
平板上の境界層厚さ $\delta(x)$ は流れ方向距離 $x$ とともに成長します：
壁面摩擦係数 $C_f$（局所）：
乱流境界層は層流より厚く成長し、壁面せん断応力も大きいため摩擦抵抗が増大します。一方、乱流境界層の方が逆圧力勾配（adverse pressure gradient）に対して剥離しにくいという特性があります。ゴルフボールのディンプル（小さな窪み）は意図的に乱流境界層を形成して剥離点を後退させ、抵抗を減らす工夫です。
6. 乱流の基礎と乱流モデル
レイノルズ平均（RANS）と乱流変動
乱流では速度が時間的に変動します。変数を時間平均値と変動成分に分解します（レイノルズ分解）：
これをNS方程式に代入して時間平均を取るとRANS方程式が得られます：
$-\rho\overline{u'_i u'_j}$ がレイノルズ応力テンソルです。対称テンソルなので独立成分は6つ（3D）。これが未知量として新たに現れるため、方程式を閉じるための「乱流モデル」が必要になります（閉合問題）。
「閉合問題」って何ですか？未知数が増えちゃって困る、ってことですか？
まさにその通り。RANS方程式にはレイノルズ応力という新しい未知量（3D空間で6成分）が増えてしまう。でも方程式の数は増えていない。これを「方程式が閉じていない」という。乱流モデルの役割は、このレイノルズ応力を既知の平均量で近似すること。
一番シンプルな近似が渦粘性仮説（Boussinesq近似）で、$-\rho\overline{u'_i u'_j} \approx \mu_t\left(\partial\overline{u}_i/\partial x_j + \partial\overline{u}_j/\partial x_i\right) - \frac{2}{3}\rho k \delta_{ij}$ と書いて、乱流粘性係数 $\mu_t$ を追加の輸送方程式で求める。$k = \frac{1}{2}\overline{u'_i u'_i}$ は乱流運動エネルギーだ。
乱流エネルギーカスケードとコルモゴロフスケール
コルモゴロフの理論によれば、乱流は大きなスケール（エネルギー注入スケール $L$）から小さなスケール（コルモゴロフスケール $\eta$）へとエネルギーが連鎖的に伝達されます（エネルギーカスケード）。
ここで $\varepsilon$ は乱流エネルギー散逸率 [m²/s³] です。高 Re 流れではスケール比 $L/\eta \sim Re^{3/4}$ となり、全スケールを解像するDNSは計算量 $\propto Re^{9/4}$ が必要です。
| 手法 | 解像スケール | 計算コスト | 対象 Re 数 |
|---|---|---|---|
| DNS（直接数値シミュレーション） | コルモゴロフスケール $\eta$ まで | $\propto Re^{9/4}$（非常に高） | $Re < 10^4$（基礎研究） |
| LES（大渦シミュレーション） | 慣性小領域より大きいスケール | $\propto Re^{13/7}$（高） | $Re \sim 10^5\sim10^6$（詳細解析） |
| RANS（レイノルズ平均） | 全スケール平均化 | 低〜中 | 任意（実用的） |
| DES（分離渦シミュレーション） | 壁近傍RANS + 自由域LES | 中〜高 | $Re \sim 10^6\sim10^7$ |
主要乱流モデルの比較
| モデル | 輸送方程式 | 得意な流れ | 苦手な流れ | 実務用途 |
|---|---|---|---|---|
| Spalart-Allmaras（SA） | 1本（$\tilde{\nu}$） | 翼周り添付境界層 | 剥離・逆圧力勾配 | 航空機外部空力 |
| Standard k-ε | 2本（$k$, $\varepsilon$） | 完全発達管内・噴流 | 壁面近傍・剥離 | 建屋換気・煙突プルーム |
| RNG k-ε | 2本（RNG導出） | 渦・曲率がある流れ | 強い剥離 | 混合タンク・室内換気 |
| Realizable k-ε | 2本（可実現条件付き） | 噴流・剪断流・チャンネル流 | 強い回転 | ターボ機械初期検討 |
| Standard k-ω | 2本（$k$, $\omega$） | 壁面近傍・逆圧力勾配 | 自由流依存性（外部流で問題） | 境界層詳細解析 |
| k-ω SST (Menter) | 2本（ブレンド関数付き） | 壁面近傍 + 自由流の両立、剥離流れ | 強い浮力効果（自然対流） | 最汎用。自動車・翼・ターボ機械 |
| RSM（レイノルズ応力） | 7本（全応力成分） | 強い異方性・スワール流れ | 収束が遅い・複雑 | 旋回燃焼器・サイクロン |
| LES（Smagorinsky等） | 非定常+SGSモデル | 大スケール渦の詳細、音響 | 壁近傍高 Re（コスト過大） | 騒音・燃焼詳細解析 |
| DES/DDES | SA or SST ベース | 大規模剥離（車体後流等） | 収束判断が難しい | 自動車後流・ブラフボディ |
k-ω SST を「最汎用」って言いましたけど、なんでこれが人気なんですか？k-ε じゃダメなんですか？
k-ε の弱点は壁面近傍だ。正確に解くには壁関数（wall function）に頼る必要があって、y+が 30〜300 の範囲に入っていないと誤差が大きい。逆に k-ω は壁面近傍は得意だけど、自由流（壁から遠い領域）で初期乱流条件に過度に依存してしまう問題がある。
SST（Shear Stress Transport）は「壁面近傍では k-ω、自由流では k-ε ライク」というブレンディング関数 $F_1$ で両者の長所を組み合わせ、さらに乱流せん断応力の過大評価を補正するリミターを付けたモデルだ。1994年にMenterが提案した後、航空・自動車・ターボ機械でデファクトスタンダードになった。「何を使えばいいかわからない」なら SST を選ぶのが最初の一手だよ。
y⁺（yプラス）の意味と壁面メッシュ設計
$y^+$ は壁面からの距離を「粘性サブ層スケール」で無次元化したものです：
壁面近傍の速度プロファイル（対数則）：
| 乱流モデル / アプローチ | 推奨 y⁺ | 壁面処理 | 備考 |
|---|---|---|---|
| k-ε（標準壁関数） | 30〜300 | 壁関数による補間 | 粗いメッシュでOK。コスト低 |
| k-ε（Enhanced Wall Treatment） | y⁺ ≈ 1 でも可 | 自動切り替え | Ansys Fluent の推奨設定 |
| k-ω SST（低 Re 解法） | y⁺ ≈ 1 | 粘性サブ層まで解像 | 精度重視。計算コスト高 |
| k-ω SST（壁関数） | y⁺ = 30〜300 | 壁関数 | 大規模解析での妥協策 |
| LES / DES | y⁺ ≈ 1〜5 | 壁面解像 | 非常に細かいメッシュが必要 |
y⁺ が推奨範囲から外れると、壁面せん断応力・熱伝達係数が大幅に誤差を持ちます。解析前に y⁺計算ツール を使って第一層メッシュ高さを必ず確認しましょう。
7. 圧縮性流れと衝撃波
等エントロピー流れとよどみ点条件
摩擦・熱交換のない等エントロピー流れでは、完全気体の状態量はマッハ数のみで決まります：
添字 0 はよどみ点（全量）条件、添字なしは静条件です。$\gamma = 1.4$（空気）、$Ma = 1$ でスロート条件 $T^*/T_0 = 2/(\gamma+1) = 0.833$、$p^*/p_0 = 0.528$ が得られます。
デ・ラバルノズルと面積比
ノズル内で超音速流れを生成するには、流路が一度収縮し（スロート）から拡大する形状（デ・ラバルノズル、収縮拡大ノズル）が必要です。スロート面積 $A^*$ との面積比 $A/A^*$ とマッハ数の関係：
ロケットのノズルってこの形ですよね。でも、面積を広げると超音速になるってのが直感的に理解できないんですが…亜音速では「細くなると速くなる」なのに、なぜ逆になるんですか？
これは「音速を超えると連続の式が逆に働く」という現象だ。亜音速では $d\rho/\rho \ll dA/A$ なので、管が細くなると流速が上がる（体積保存 ≒ 断面積保存）。でも超音速では $d\rho/\rho$ が $dA/A$ より大きくなる——密度の減少が速いから、面積が広がっても連続性を保つために流速がさらに上がらないといけない。
スロート（のど部）で $Ma = 1$ になると「チョーキング（閉塞）」が起きて、上流の圧力を上げてもスロートを通過する質量流量は増えない。ロケットエンジンのベル形ノズルはまさにこの原理を利用して、燃焼ガスを超音速まで加速して推力を最大化してるんだよ。
垂直衝撃波の関係式（ランキン・ユゴニオ）
垂直衝撃波前後（添字 1：上流, 2：下流）の関係式：
衝撃波は数平均自由行程程度の薄い不連続面で、圧力・密度・温度が急上昇し、エントロピーが増加します（不可逆過程）。CFDでの衝撃波捕捉には数値拡散制御が重要で、密度ベースソルバー（Roe法・AUSM法・JST法など）と TVD スキームが必要です。
| $Ma_1$（上流） | $Ma_2$（下流） | $p_2/p_1$ | $T_2/T_1$ | $\rho_2/\rho_1$ |
|---|---|---|---|---|
| 1.0 | 1.000 | 1.000 | 1.000 | 1.000 |
| 1.5 | 0.701 | 2.458 | 1.320 | 1.862 |
| 2.0 | 0.577 | 4.500 | 1.688 | 2.667 |
| 3.0 | 0.475 | 10.33 | 2.679 | 3.857 |
| 5.0 | 0.415 | 29.00 | 5.800 | 5.000 |
8. 数値流体力学（CFD）の概要
有限体積法（FVM）による離散化
CFDで最もよく使われるのは有限体積法（FVM: Finite Volume Method）です。計算領域をセル（コントロールボリューム）に分割し、各セルに対して保存則の積分形を適用します。セル中心 $P$ とその隣接セル $N$ の間の境界面でフラックスを評価します。
$\phi$ は保存変数（速度成分・温度・$k$・$\varepsilon$ 等）、$\dot{m}_f$ は境界面での質量フラックス、$\Gamma_\phi$ は拡散係数、$A_f$ は面積、$S_\phi$ はソース項です。FVMは保存則をセル単位で厳密に満たす（保存性が高い）のが特長です。
圧力速度連成：SIMPLE法
非圧縮性のNS方程式って、圧力の方程式がないですよね？CFDではどうやって圧力を求めてるんですか？
鋭い！非圧縮では圧力は「連続の式 $\nabla \cdot \mathbf{v} = 0$ を満たすための制約」として決まる。これが圧力速度連成問題だ。SIMPLE（Semi-Implicit Method for Pressure-Linked Equations）法のやり方はこう：
- 仮の圧力場 $p^*$ を使って運動量方程式を解き、仮の速度場 $\mathbf{v}^*$ を得る
- $\mathbf{v}^*$ は連続の式を満たさないので、圧力補正量 $p'$ を解くポアソン方程式を立てて解く
- $p = p^* + \alpha_p p'$、$\mathbf{v} = \mathbf{v}^* + \mathbf{v}'$ でアンダーリラクセーション付きで更新
- 残差が十分小さくなるまで繰り返す
改良版のSIMPLEC（係数をより正確に計算）、非定常向けのPISO（Pressure-Implicit Splitting of Operators）もある。OpenFOAMでは PIMPLE（SIMPLE+PISO のハイブリッド）が非定常解析でよく使われるよ。
対流スキームの比較
| スキーム | 精度次数 | 安定性 | 数値拡散 | 振動 | 用途 |
|---|---|---|---|---|---|
| 1次風上差分（UDS/UD） | 1次 | ◎（無条件安定） | 大 | なし | 収束確認・初期解・ロバスト性重視 |
| 中心差分（CDS） | 2次 | △（高 Pe 数で不安定） | なし | あり | LES・DNS（数値拡散が最小限必要） |
| 線形風上（LUDS） | 2次 | ○ | 中 | 小 | 定常RANS標準的な選択 |
| QUICK（3次風上加重） | 3次 | ○ | 小 | わずかにあり | 精度重視のRANS |
| TVDスキーム（MUSCL + limiters） | 2〜3次 | ○ | 小 | 抑制 | 衝撃波・密度不連続を含む流れ |
| WENO（高次重み付き本質非振動） | 5次以上 | ○ | 非常に小 | なし | 超音速・爆発波・精密LES |
実務では最初に1次風上で収束させ、その後2次精度スキームに切り替えて精度を上げる「2段階解法」が定石です。いきなり高次スキームを使うと発散しやすくなります。
収束判定と残差の読み方
OpenFOAMで計算してると「残差」がいっぱい出てくるんですが、どの値になったら「収束した」って言えるんですか？
一般的な目安はこう：
- 速度成分 ($U_x, U_y, U_z$): $10^{-4}$ 未満が目標、$10^{-3}$ で工学的に許容されることも
- 連続残差 (continuity): $10^{-4}$ 未満。これが大きいと質量保存が崩れてる
- エネルギー（温度）: $10^{-7}$ 未満（数値が小さく収束しやすい）
- 乱流変数 ($k, \omega, \varepsilon$): $10^{-5}$ 未満
ただし残差だけ見てたら危ない。関心のある物理量（ドラッグ係数 $C_D$、圧力損失 $\Delta p$、熱流束 $q$）が変化しなくなることを必ず確認すること。残差は下がってるのに $C_D$ がフラフラしてる、というのは収束してない証拠だ。定量的なモニタリングを必ず設定しよう。
9. 実務でのCFDの注意点
境界条件の設定ミスと対策
CFD解析でもっとも多い失敗の一つが境界条件の誤設定です。計算が回っていても、物理的に意味のない境界条件を設定していると正しい結果は得られません。
| 境界 | よくあるミス | 正しい設定・確認ポイント |
|---|---|---|
| 入口（Inlet） | 乱流強度と乱流長さスケールをデフォルト（TI=5%）放置 | TI = 1〜5%（内部流路）または 0.1〜1%（外部流）。乱流長さスケールは水力直径の0.07倍程度 |
| 出口（Outlet） | 逆流が生じる場所に圧力出口を設定、逆流条件を無設定 | 逆流が予想される場合は backflow 乱流条件を設定。計算域を出口より先まで延長するのも有効 |
| 壁面（Wall） | 回転壁面の速度を0に設定してしまう（ターボ機械） | Moving Wall / MRF（Multiple Reference Frame）/ Sliding Mesh の適切な使用 |
| 対称面（Symmetry） | 非対称な流れ（カルマン渦）に対称条件を設定 | カルマン渦列・非対称剥離は symmetry 使用不可。半モデルでも非定常解析を |
| 熱境界 | 外壁を断熱（adiabatic）にしたまま熱解析。全温と静温の混同 | 測定値や設計仕様から熱流束または対流熱伝達係数と外気温を明確に設定 |
| 周期境界 | 周期長さが物理的な周期と合っていない | 圧力降下指定型（Fluent の Periodic BC）と流量指定型の使い分けに注意 |
メッシュ品質と数値拡散
「数値拡散」ってよく聞くんですが、これって何が問題なんですか？物理的な拡散と区別できるんですか？
数値拡散は「数値的な誤差が物理的な拡散のように見えてしまう現象」だ。1次風上スキームの截断誤差を展開すると、実は擬似的な拡散係数 $\nu_{\text{num}} \approx U\Delta x/2$ が加わっていることがわかる。これが大きいと、本来鋭いはずの温度フロントや渦が「ぼやけた」結果になる。特に流れ方向に対して斜めにメッシュが走っているとひどくなる。
見分け方はメッシュを細かくして結果が変わるかどうかだ。粗いメッシュと細かいメッシュで温度分布が大きく変わるなら数値拡散が支配している証拠。これをメッシュ収束確認（Grid Independence Study）という。実務では少なくとも2段階（粗・細）のメッシュで比較するのが原則だよ。
| 品質指標 | 良好な範囲 | 問題が出やすい値 | 影響 |
|---|---|---|---|
| アスペクト比（Aspect Ratio） | < 100（境界層以外） | > 1000 | 収束悪化・精度低下 |
| スキューネス（Skewness） | < 0.5 | > 0.85 | 離散化誤差増大・発散リスク |
| 直交性（Non-orthogonality） | < 40° (OpenFOAM定義) | > 70° | 圧力速度連成の誤差増大 |
| 体積比（Volume Change Ratio） | < 10 | > 100 | 急激な体積変化で局所誤差 |
物性値のデータソースと注意点
物性値の誤りは「計算が通る」が「答えが全然違う」という最悪のパターンを生みます。以下の注意点を守ってください。
- 温度依存性の考慮：空気の粘度は20℃〜100℃で約1.5倍変化します。温度範囲が広い解析では定物性仮定が不適切な場合があります
- 混合物の物性：冷却液（エチレングリコール水溶液）の物性は濃度・温度に大きく依存します。メーカーのデータシートを必ず確認してください
- 推奨データソース：
- NIST Webbook（https://webbook.nist.gov）: 気体・液体の熱物性を温度・圧力の関数で無料取得可能
- Engineering Toolbox: 実用的な範囲の物性値を素早く確認
- CRC Handbook of Chemistry and Physics: 実験値の信頼性が高い
先生、最後に「CFDを始める人が一番やりがちなミス」を3つ教えてください。
①検証なしでいきなり細かいメッシュを切る：まず粗いメッシュで動かして「物理的に妥当か」を確認し、収束したらメッシュを細かくするのが正解。最初から10倍のメッシュを切って、答えが合わなかったときの修正コストは大きい。
②デフォルト乱流強度のまま計算：入口条件の乱流強度（Turbulent Intensity）は結果に大きく影響するのに、デフォルト（5%）で放置する人が多い。実験や規格を参考に適切に設定すること。
③残差だけ見て収束と判断する：モニタリングポイントで $C_D$ や $\Delta p$ など関心量を監視すること。あと「計算が回った＝正解」ではない。必ずベンチマーク（Moody線図との比較、解析解との照合）で手法の妥当性を確かめてから本計算に進むこと。
インタラクティブCFDツール
ブラウザ上でリアルタイムに流体パラメータを計算・可視化
レイノルズ数・y⁺・境界層厚さ・圧力損失・ノズル流れなど、CFD前処理で必要な計算をワンクリックで。