天気予報が数日先の風や雨を予測できるのは、なぜでしょうか。航空機の翼が揚力を発生するメカニズムを、実機を作らずに検証できるのはなぜでしょうか。その答えの中心にあるのが数値流体力学（Computational Fluid Dynamics, CFD）です。
流体の運動はナビエ・ストークス方程式という非線形偏微分方程式に支配されています。この方程式は、ごく限られた特殊な条件を除いて解析解（紙と鉛筆で得られる厳密解）が存在しません。しかし現実の工学・科学では、複雑な形状の周りの流れ、乱流、熱と物質の輸送など、解析解が得られない問題こそが重要です。CFDは、この支配方程式をコンピュータで近似的に解くことで、流体現象を定量的に予測する手法です。
CFDの応用範囲は驚くほど広いです。
- 航空宇宙: 翼型周りの空力解析、ロケットエンジン燃焼室内の流れ、再突入カプセルの空力加熱予測
- 自動車: 車体の空気抵抗低減、エンジン冷却系の最適化、車室内の空調設計
- 気象・海洋: 数値天気予報、台風の進路予測、海流シミュレーション
- 医療: 血管内の血流解析、人工心臓弁の設計、呼吸器内の気流シミュレーション
- エネルギー: 風力タービンの翼設計、原子炉内の冷却材流動解析、建築物の換気設計
本記事では、CFDの基礎を「支配方程式 → 離散化 → 求解」の3ステップで体系的に解説します。空間離散化では差分法と有限体積法を、時間積分では陽的・陰的スキームを取り上げ、安定性を保証するCFL条件についても説明します。最後に、CFDの古典的なベンチマーク問題である2次元リッド駆動キャビティ流れをPythonでスクラッチ実装し、流れの物理を自分の目で確認します。
本記事の内容
- CFDの基本的な考え方（支配方程式 → 離散化 → 求解）
- 空間離散化: 差分法（前進差分・後退差分・中心差分）とその精度
- 時間積分: 陽的スキームと陰的スキームの特徴
- 安定性条件（CFL条件）の直感的理解と数学的導出
- 有限体積法の基本概念と保存則
- 格子生成の基礎（構造格子と非構造格子）
- Pythonによるスクラッチ実装: 1次元移流方程式、2次元拡散方程式、リッド駆動キャビティ流れ
- 結果の可視化と物理的解釈
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ナビエ・ストークス方程式 — 流体の運動を支配する基本方程式
- 連続の式 — 質量保存則の微分形
- レイノルズ数 — 慣性力と粘性力の比を表す無次元数
- 流体の性質 — 密度・粘性係数・圧縮性の基礎
- 偏微分方程式 — 偏微分の定義と計算
- テイラー展開 — 関数の近似と打ち切り誤差の理解に必要
CFDの基本的な考え方
流体のシミュレーションとは何か
日常的な直感として、「水が障害物にぶつかると二手に分かれ、下流で渦を巻く」という経験を持っている方は多いでしょう。しかし「どのくらいの速度で」「どのような形の渦が」「いつ発生するのか」を正確に答えるには、流体の運動方程式を解く必要があります。
流体の振る舞いは、ナビエ・ストークス方程式と連続の式によって記述されます。非圧縮性のニュートン流体（水や低速の空気がこれに該当します）の場合、支配方程式は次のようになります。
連続の式（質量保存）:
$$ \begin{equation} \nabla \cdot \bm{u} = 0 \end{equation} $$
ナビエ・ストークス方程式（運動量保存）:
$$ \begin{equation} \frac{\partial \bm{u}}{\partial t} + (\bm{u} \cdot \nabla)\bm{u} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \bm{u} \end{equation} $$
ここで $\bm{u}$ は速度ベクトル、$p$ は圧力、$\rho$ は流体の密度、$\nu$ は動粘性係数です。左辺の第2項 $(\bm{u} \cdot \nabla)\bm{u}$ は移流項と呼ばれ、流体が自分自身を運ぶ非線形効果を表しています。右辺第1項は圧力勾配による力、第2項は粘性による拡散力です。
この方程式は非線形偏微分方程式であり、一般には解析解を求めることができません。そこで、空間と時間を有限個の離散点に分割し、連続的な微分方程式を代数方程式の集合に変換して数値的に解く――これがCFDの基本戦略です。
CFDの3ステップ
CFDの手順は、どのような手法を用いる場合でも、大きく以下の3つのステップに分けられます。
ステップ1: 支配方程式の選定
対象とする流れの特性に応じて、適切な支配方程式を選びます。非圧縮性流れならナビエ・ストークス方程式、圧縮性流れならオイラー方程式またはその粘性拡張、乱流を扱うならRANS（Reynolds-Averaged Navier-Stokes）やLES（Large Eddy Simulation）の方程式を用います。
ステップ2: 離散化
連続的な偏微分方程式を、有限個の代数方程式に変換します。空間の離散化には差分法、有限体積法、有限要素法（FEM）の3つが代表的です。時間の離散化には陽的スキームや陰的スキームを用います。
ステップ3: 求解と後処理
離散化によって得られた連立方程式を反復法や直接法で解きます。得られた数値解を流線、等値面、ベクトル場などで可視化し、物理的な妥当性を検証します。
この3ステップの中で、CFDの精度と効率を決定づけるのがステップ2の「離散化」です。以下では、まず空間の離散化手法として最も基本的な差分法から解説していきます。
差分法の基礎
テイラー展開と差分近似
差分法の核心は、微分を「近傍の関数値の差の商」で近似することにあります。この近似がどの程度正確なのかを理解するために、テイラー展開が不可欠です。
1次元の関数 $f(x)$ を点 $x_i$ の周りでテイラー展開すると、次のようになります。
$$ \begin{equation} f(x_i + \Delta x) = f(x_i) + \Delta x \frac{df}{dx}\bigg|_{x_i} + \frac{(\Delta x)^2}{2!} \frac{d^2 f}{dx^2}\bigg|_{x_i} + \frac{(\Delta x)^3}{3!} \frac{d^3 f}{dx^3}\bigg|_{x_i} + \cdots \end{equation} $$
同様に、$x_i – \Delta x$ の方向に展開すると次のようになります。
$$ \begin{equation} f(x_i – \Delta x) = f(x_i) – \Delta x \frac{df}{dx}\bigg|_{x_i} + \frac{(\Delta x)^2}{2!} \frac{d^2 f}{dx^2}\bigg|_{x_i} – \frac{(\Delta x)^3}{3!} \frac{d^3 f}{dx^3}\bigg|_{x_i} + \cdots \end{equation} $$
これらの式を操作することで、さまざまな差分近似を導くことができます。
前進差分
式(3)を $df/dx$ について解くと、次のようになります。
$$ \begin{equation} \frac{df}{dx}\bigg|_{x_i} = \frac{f(x_i + \Delta x) – f(x_i)}{\Delta x} – \frac{\Delta x}{2} \frac{d^2 f}{dx^2}\bigg|_{x_i} – \cdots \end{equation} $$
2次以上の項を打ち切ると、前進差分が得られます。
$$ \begin{equation} \frac{df}{dx}\bigg|_{x_i} \approx \frac{f_{i+1} – f_i}{\Delta x} \end{equation} $$
ここで $f_i = f(x_i)$ という簡略記法を使いました。打ち切り誤差は $O(\Delta x)$ であり、これを1次精度と呼びます。格子間隔 $\Delta x$ を半分にすると誤差も半分になるという性質です。
後退差分
式(4)を同様に変形すると、後退差分が得られます。
$$ \begin{equation} \frac{df}{dx}\bigg|_{x_i} \approx \frac{f_i – f_{i-1}}{\Delta x} \end{equation} $$
これも1次精度です。前進差分が「今の点と右隣」の情報を使うのに対し、後退差分は「今の点と左隣」の情報を使います。
中心差分
式(3)から式(4)を引き算してみましょう。
$$ \begin{equation} f(x_i + \Delta x) – f(x_i – \Delta x) = 2 \Delta x \frac{df}{dx}\bigg|_{x_i} + \frac{2(\Delta x)^3}{3!} \frac{d^3 f}{dx^3}\bigg|_{x_i} + \cdots \end{equation} $$
$df/dx$ について解くと、中心差分が得られます。
$$ \begin{equation} \frac{df}{dx}\bigg|_{x_i} \approx \frac{f_{i+1} – f_{i-1}}{2\Delta x} \end{equation} $$
打ち切り誤差は $O((\Delta x)^2)$ であり、2次精度です。格子間隔を半分にすると誤差は4分の1になります。前進・後退差分に比べて、同じ格子間隔でより正確な近似が得られるため、多くのCFDコードで中心差分が基本的な近似として使われています。
2階微分の差分近似
拡散項 $\nabla^2 f$ に現れる2階微分も差分で近似する必要があります。式(3)と式(4)を足し算すると、1階微分の項が打ち消し合い、次の式が得られます。
$$ \begin{equation} f_{i+1} + f_{i-1} = 2f_i + (\Delta x)^2 \frac{d^2 f}{dx^2}\bigg|_{x_i} + O((\Delta x)^4) \end{equation} $$
$d^2f/dx^2$ について解くと、2階微分の中心差分近似が得られます。
$$ \begin{equation} \frac{d^2 f}{dx^2}\bigg|_{x_i} \approx \frac{f_{i+1} – 2f_i + f_{i-1}}{(\Delta x)^2} \end{equation} $$
これも2次精度です。この式は「右隣と左隣の平均値が中心の値からどれだけずれているか」を測っており、物理的には拡散（均一化しようとする作用）を表現していることが直感的にも理解できます。
ここまでで、1階微分と2階微分を離散的に近似する方法を学びました。次に、時間方向の離散化について見ていきましょう。時間方向の離散化は、計算の安定性と精度に直結する極めて重要な選択です。
時間積分スキーム
陽的スキーム（Explicit scheme）
時間発展する偏微分方程式を数値的に解くとき、時間方向の離散化が必要です。最も直感的な方法は、「現在の時刻の情報だけを使って、次の時刻の値を直接計算する」というアプローチです。これを陽的スキームと呼びます。
1次元の移流方程式を例に考えましょう。
$$ \begin{equation} \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 \end{equation} $$
この方程式は、波形が速度 $c$ で右に移動することを表しています。時間微分を前進差分、空間微分を後退差分で近似すると、次のような離散式が得られます。
$$ \begin{equation} \frac{u_i^{n+1} – u_i^n}{\Delta t} + c \frac{u_i^n – u_{i-1}^n}{\Delta x} = 0 \end{equation} $$
ここで上付き添字 $n$ は時間ステップ番号です。$u_i^{n+1}$ について解くと次のようになります。
$$ \begin{equation} u_i^{n+1} = u_i^n – c \frac{\Delta t}{\Delta x}(u_i^n – u_{i-1}^n) \end{equation} $$
この式の右辺には時刻 $n$ の値しか含まれていません。つまり、時刻 $n$ の情報がすべてわかっていれば、連立方程式を解くことなく、各点の次の時刻の値を独立に計算できます。これが陽的スキームの最大の利点です。実装が簡単で、並列計算にも適しています。
陰的スキーム（Implicit scheme）
一方、「次の時刻の情報も使って離散化する」アプローチが陰的スキームです。同じ移流方程式に対して、空間微分を時刻 $n+1$ で評価すると次のようになります。
$$ \begin{equation} \frac{u_i^{n+1} – u_i^n}{\Delta t} + c \frac{u_i^{n+1} – u_{i-1}^{n+1}}{\Delta x} = 0 \end{equation} $$
この式には時刻 $n+1$ の値が複数の項に現れています。したがって、全ての空間点について連立方程式を立てて同時に解く必要があります。計算の手間は陽的スキームより大きくなりますが、後述するように安定性の面で大きな優位性があります。
クランク・ニコルソン法
陽的と陰的の「いいとこ取り」をする手法がクランク・ニコルソン法です。空間微分を時刻 $n$ と $n+1$ の平均で評価します。
$$ \begin{equation} \frac{u_i^{n+1} – u_i^n}{\Delta t} = -c \cdot \frac{1}{2}\left(\frac{u_i^n – u_{i-1}^n}{\Delta x} + \frac{u_i^{n+1} – u_{i-1}^{n+1}}{\Delta x}\right) \end{equation} $$
この方法は時間方向に2次精度を達成し、かつ無条件安定（後述のCFL条件に縛られない）という優れた特性を持ちます。拡散方程式の数値解法では特に広く使われています。
では、陽的スキームが「不安定になる」とは具体的にどういうことなのでしょうか。ここで、CFDの最も重要な概念の一つであるCFL条件について見ていきましょう。
CFL条件と安定性
不安定性の直感
数値計算が「不安定になる」とは、本来滑らかであるべき解が時間ステップを進めるにつれて振動し、最終的に発散してしまう現象です。なぜこのようなことが起こるのでしょうか。
直感的に考えてみましょう。速度 $c$ で移動する波を追いかけるとき、1時間ステップ $\Delta t$ の間に波は距離 $c \Delta t$ だけ進みます。一方、格子間隔は $\Delta x$ です。もし $c \Delta t > \Delta x$ であれば、波は1時間ステップの間に格子点を1つ以上飛び越えてしまいます。飛び越えた先の情報は差分式に含まれていないので、波の動きを正しく追跡できず、数値的な振動や発散が起こるのです。
CFL条件の数学的な意味
Courant-Friedrichs-Lewy（CFL）条件は、陽的スキームが安定であるための必要条件を与えます。1次元の移流方程式に対しては、次の条件が成り立つ必要があります。
$$ \begin{equation} \text{CFL} = \frac{c \Delta t}{\Delta x} \leq 1 \end{equation} $$
この無次元数 $c\Delta t / \Delta x$ をCFL数（クーラン数）と呼びます。CFL数が1を超えると計算は発散します。
拡散方程式 $\partial u / \partial t = \alpha \partial^2 u / \partial x^2$ に対する陽的スキームでは、安定条件はさらに厳しくなります。
$$ \begin{equation} \frac{\alpha \Delta t}{(\Delta x)^2} \leq \frac{1}{2} \end{equation} $$
格子間隔 $\Delta x$ の2乗に比例して時間刻み $\Delta t$ を小さくしなければならないため、高解像度の計算では膨大な時間ステップ数が必要になります。これが陰的スキームが好まれる理由の一つです。陰的スキームは（理論上は）無条件安定であり、CFL条件に制限されずに大きな時間刻みを取ることができます。ただし、精度の観点からは時間刻みを無限に大きくしてよいわけではないことに注意が必要です。
フォン・ノイマン安定性解析
CFL条件を厳密に導出する方法として、フォン・ノイマンの安定性解析があります。この方法では、数値解の誤差を波数 $k$ のフーリエモードで展開し、各モードの増幅率を調べます。
離散解を $u_i^n = G^n e^{jk i \Delta x}$ の形に書くと（$G$ は増幅因子、$j$ は虚数単位）、陽的前進差分法（式(15)）の増幅因子は次のようになります。
$$ \begin{equation} G = 1 – \sigma(1 – e^{-jk\Delta x}) \end{equation} $$
ここで $\sigma = c\Delta t / \Delta x$ はCFL数です。安定性の条件は全ての波数 $k$ に対して $|G| \leq 1$ が成り立つことです。$|G|^2$ を計算すると次のようになります。
$$ \begin{equation} |G|^2 = 1 – 2\sigma(1 – \sigma)(1 – \cos(k\Delta x)) \end{equation} $$
全ての $k$ に対して $|G|^2 \leq 1$ とするには、$\sigma(1 – \sigma) \geq 0$、すなわち $0 \leq \sigma \leq 1$ が必要です。これがCFL条件の数学的な導出です。
ここまでで、差分法による空間離散化と時間積分の基礎を理解しました。しかし、差分法には「保存則を離散レベルで厳密に満たす保証がない」という弱点があります。この問題を克服する手法として、実用的なCFDコードで最も広く使われている有限体積法を次に見ていきましょう。
有限体積法の基本概念
差分法との違い
差分法は「微分を差分で近似する」手法でした。これに対して有限体積法（Finite Volume Method, FVM）は、「保存則の積分形を直接離散化する」手法です。
流体力学の支配方程式は本質的に保存則です。質量は創られも消えもせず、運動量の変化は力に等しく、エネルギーの変化は仕事と熱の出入りに等しい。有限体積法は、この保存則を離散レベルでも厳密に満たすように設計されています。これが有限体積法がCFDの主流となった最大の理由です。
保存則の積分形
1次元のスカラー保存則を例に考えましょう。
$$ \begin{equation} \frac{\partial u}{\partial t} + \frac{\partial F(u)}{\partial x} = 0 \end{equation} $$
ここで $u$ は保存量（例えば密度）、$F(u)$ はフラックス（保存量の流れ）です。この方程式を空間の区間 $[x_{i-1/2}, x_{i+1/2}]$ で積分します。
$$ \begin{equation} \frac{d}{dt}\int_{x_{i-1/2}}^{x_{i+1/2}} u \, dx + F(u(x_{i+1/2})) – F(u(x_{i-1/2})) = 0 \end{equation} $$
セル平均値を $\bar{u}_i = \frac{1}{\Delta x}\int_{x_{i-1/2}}^{x_{i+1/2}} u \, dx$ と定義すると、次のように書けます。
$$ \begin{equation} \frac{d\bar{u}_i}{dt} = -\frac{F_{i+1/2} – F_{i-1/2}}{\Delta x} \end{equation} $$
この式は物理的に明快な意味を持っています。「セル内の保存量の時間変化率は、セルの右面から出ていくフラックスと左面から入ってくるフラックスの差に等しい」ということです。あるセルの右面から出ていくフラックス $F_{i+1/2}$ は、隣のセルの左面から入ってくるフラックスと完全に同じです。したがって、隣接するセル間でフラックスがキャンセルされ、離散レベルでも保存則が厳密に満たされます。
セル界面のフラックス評価
有限体積法の核心は、セル界面でのフラックス $F_{i+1/2}$ をどのように評価するかにあります。最も単純な方法は、セル平均値からの補間（線形再構成）を使う方法です。
1次精度（定数再構成）: セル内の値が一定だと仮定し、上流側のセルの値を使います（上流差分法に相当）。数値拡散が大きいですが安定です。
2次精度（線形再構成）: セル内の値を線形分布で近似し、界面での値を補間で求めます。精度は上がりますが、不連続面の近くで数値振動が生じる可能性があります。
実用的なCFDコードでは、上流側の情報を多く取り入れつつ精度を維持する「上流型スキーム」や、TVD（Total Variation Diminishing）制限関数を使って振動を抑えるスキームが広く使われています。
有限体積法は有限要素法（FEM）と並んでCFDの二大手法です。FEMは変分原理に基づき形状関数を使って近似するのに対し、FVMは保存則の積分形を直接離散化します。流体力学ではFVMの方が主流ですが、構造解析ではFEMが主流です。それぞれの手法には得意・不得意があり、問題に応じて使い分けることが重要です。
次に、離散化を行うための「舞台」となる計算格子について見ていきましょう。格子の質はCFDの結果に直結するため、格子生成の基礎を理解しておくことは重要です。
格子生成の基礎
計算格子とは
CFDで空間を離散化するとき、計算領域を有限個のセル（要素）に分割します。この分割のパターンを計算格子（メッシュ）と呼びます。格子は「物理的な空間をコンピュータが扱える数値データに変換する橋渡し」です。
格子の品質は計算結果の精度に直結します。格子が粗すぎれば流れの詳細を捉えられず、格子が歪んでいれば数値誤差が大きくなります。一方、格子を細かくしすぎれば計算コストが膨大になります。適切な格子を生成することは、CFDにおける最も重要な技術の一つです。
構造格子と非構造格子
計算格子は大きく2種類に分けられます。
構造格子（Structured grid）: 格子点が規則的に並んでいるもので、(i, j, k) のインデックスで全ての格子点にアクセスできます。直交格子（デカルト格子）や曲線座標系に沿った格子がこれに当たります。格子点の配列が規則的なため、メモリ効率が良く、高次精度のスキームを適用しやすいという利点があります。反面、複雑な形状への適合が難しいという欠点があります。
非構造格子（Unstructured grid）: 三角形（2次元）や四面体（3次元）などの要素を自由に配置する格子です。格子点の接続関係をテーブルで管理する必要があるため、メモリ使用量は増えますが、複雑な形状にも柔軟に対応できます。実用的なCFDソフト（OpenFOAMなど）の多くは非構造格子を採用しています。
格子の品質指標
良い格子の条件として、以下の点が挙げられます。
- 直交性: 隣接セルを結ぶ線とセル界面が直交に近いほど、フラックスの評価精度が高い
- 滑らかさ: 隣接するセルのサイズ比が急激に変化しないこと。一般に、拡大率は1.2倍以下が推奨される
- アスペクト比: セルの縦横比が極端にならないこと。境界層の近くでは壁面垂直方向に薄いセルを使うことが多いですが、それでもアスペクト比が大きすぎると精度が落ちます
- 解像度: 流れの変化が激しい領域（壁面近傍、衝撃波、剪断層など）ではセルを細かくする
本記事のPython実装では、最もシンプルな等間隔の直交格子（構造格子）を使います。これは実装が簡単で、差分法や有限体積法の概念を学ぶのに最適です。
ここまでで、CFDの理論的な基盤は一通り押さえました。いよいよ、これらの概念をPythonで実際に動かしてみましょう。まず最も単純な1次元移流方程式から始め、段階的に複雑な問題へ進みます。
Python実装1: 1次元移流方程式の差分法
問題設定
最初のPython実装として、1次元の移流方程式を差分法で解きます。移流方程式は波動現象の最も基本的なモデルであり、差分法の性質を理解するのに最適です。
$$ \begin{equation} \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0, \quad x \in [0, 2], \quad c = 1 \end{equation} $$
初期条件として正弦波を与え、この波が速度 $c = 1$ で右に移動する様子をシミュレーションします。厳密解は $u(x, t) = u_0(x – ct)$ であり、初期波形がそのまま右に平行移動するだけです。数値解がこの厳密解をどの程度再現できるかを確認します。
以下のコードでは、前進差分・後退差分（1次上流差分法）と、CFL条件を満たす場合・満たさない場合の比較を行います。
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ設定 ---
nx = 81 # 格子点数
dx = 2.0 / (nx - 1) # 格子間隔
c = 1.0 # 移流速度
nt = 50 # 時間ステップ数
# CFL条件を満たす場合 (CFL = 0.5) と満たさない場合 (CFL = 1.5)
cfl_values = [0.5, 1.5]
labels = ["CFL = 0.5 (安定)", "CFL = 1.5 (不安定)"]
x = np.linspace(0, 2, nx)
# 初期条件: 正弦波パルス
u0 = np.sin(2 * np.pi * x) * (x <= 1.0)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
for idx, cfl in enumerate(cfl_values):
dt = cfl * dx / c # CFL数からdtを逆算
u = u0.copy()
for n in range(nt):
un = u.copy()
# 1次上流差分法 (c > 0 なので後退差分)
for i in range(1, nx):
u[i] = un[i] - c * dt / dx * (un[i] - un[i - 1])
u[0] = un[0] - c * dt / dx * (un[0] - un[-1]) # 周期境界
# 厳密解
u_exact = np.sin(2 * np.pi * (x - c * nt * dt)) * ((x - c * nt * dt) % 2 <= 1.0)
axes[idx].plot(x, u0, "k--", label="初期条件", linewidth=1)
axes[idx].plot(x, u_exact, "b-", label="厳密解", linewidth=2)
if np.max(np.abs(u)) < 10: # 発散していなければプロット
axes[idx].plot(x, u, "r-o", label="数値解", markersize=3, linewidth=1)
else:
axes[idx].set_title(f"{labels[idx]}\n(数値解は発散)", fontsize=13)
axes[idx].set_xlabel("x", fontsize=12)
axes[idx].set_ylabel("u", fontsize=12)
axes[idx].set_title(labels[idx], fontsize=13)
axes[idx].legend(fontsize=10)
axes[idx].set_ylim(-2, 2)
axes[idx].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("advection_cfl_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
上のコードの実行結果から、2つの重要な事実が確認できます。
- CFL = 0.5 の場合（左図）: 数値解は厳密解と概ね一致していますが、波形が少し平坦化しています。これは1次上流差分法に固有の数値拡散（数値粘性）と呼ばれる現象で、1次精度スキームの宿命です。打ち切り誤差の主要項が拡散型であるため、波の振幅が時間とともに減衰していきます
- CFL = 1.5 の場合（右図）: 数値解は発散し、もはや物理的に意味のある結果ではなくなります。CFL条件を満たさないと計算が破綻するという事実を、このシンプルな例で直接確認できます
この結果は、先に述べたフォン・ノイマンの安定性解析と完全に整合しています。CFL数が1を超えると増幅因子 $|G| > 1$ となるモードが存在し、誤差が指数的に成長するのです。
1次元の移流方程式で差分法の基本的な振る舞いを確認しました。次は、拡散現象の数値解法に進みましょう。拡散方程式はナビエ・ストークス方程式の粘性項に直接対応するため、CFDにとって不可欠な要素です。
Python実装2: 2次元拡散方程式
問題設定
2つ目の実装として、2次元の拡散方程式（熱伝導方程式）を解きます。
$$ \begin{equation} \frac{\partial u}{\partial t} = \alpha \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \end{equation} $$
物理的には、ある領域内の温度分布（または濃度分布）が時間とともに均一化していく過程を表しています。初期条件として中央にホットスポット（高温領域）を配置し、それが拡散していく様子をシミュレーションします。
2次元の拡散方程式を陽的オイラー法（FTCS: Forward Time, Central Space スキーム）で離散化すると、次のようになります。
$$ \begin{equation} u_{i,j}^{n+1} = u_{i,j}^n + \alpha \Delta t \left(\frac{u_{i+1,j}^n – 2u_{i,j}^n + u_{i-1,j}^n}{(\Delta x)^2} + \frac{u_{i,j+1}^n – 2u_{i,j}^n + u_{i,j-1}^n}{(\Delta y)^2}\right) \end{equation} $$
2次元の場合の安定条件は次のようになります。
$$ \begin{equation} \alpha \Delta t \left(\frac{1}{(\Delta x)^2} + \frac{1}{(\Delta y)^2}\right) \leq \frac{1}{2} \end{equation} $$
$\Delta x = \Delta y$ の等方格子の場合、$\alpha \Delta t / (\Delta x)^2 \leq 1/4$ となり、1次元の場合よりもさらに厳しい制約です。
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ設定 ---
nx, ny = 51, 51 # 格子点数
Lx, Ly = 2.0, 2.0 # 領域サイズ
dx = Lx / (nx - 1)
dy = Ly / (ny - 1)
alpha = 0.01 # 熱拡散率
dt = 0.2 * dx**2 / (2 * alpha) # 安定条件を十分に満たすdt
nt = 500 # 時間ステップ数
x = np.linspace(0, Lx, nx)
y = np.linspace(0, Ly, ny)
X, Y = np.meshgrid(x, y)
# 初期条件: 中央にガウス分布のホットスポット
u = np.exp(-((X - 1.0)**2 + (Y - 1.0)**2) / 0.05)
# 結果を保存するリスト (初期, 中間, 最終)
snapshots = [u.copy()]
snap_times = [0]
# --- 時間発展ループ ---
for n in range(nt):
un = u.copy()
# 内部点に対してFTCSスキームを適用
u[1:-1, 1:-1] = (un[1:-1, 1:-1]
+ alpha * dt / dx**2 * (un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[:-2, 1:-1])
+ alpha * dt / dy**2 * (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, :-2]))
# 境界条件: ディリクレ (u = 0)
u[0, :] = 0; u[-1, :] = 0
u[:, 0] = 0; u[:, -1] = 0
if n == nt // 4:
snapshots.append(u.copy())
snap_times.append(n * dt)
elif n == nt - 1:
snapshots.append(u.copy())
snap_times.append(n * dt)
# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
titles = [f"t = {snap_times[0]:.3f}", f"t = {snap_times[1]:.3f}", f"t = {snap_times[2]:.3f}"]
for ax, snap, title in zip(axes, snapshots, titles):
im = ax.contourf(X, Y, snap, levels=30, cmap="hot")
ax.set_title(title, fontsize=13)
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_aspect("equal")
plt.colorbar(im, ax=ax, shrink=0.8)
plt.suptitle("2D Diffusion Equation (FTCS scheme)", fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig("diffusion_2d.png", dpi=150, bbox_inches="tight")
plt.show()
上のグラフから、拡散方程式の基本的な振る舞いが確認できます。
- 初期（$t = 0$）: 中央に鋭いガウス分布のホットスポットが存在します。温度勾配はホットスポットの縁で最大です
- 中間時刻: ホットスポットは周囲に広がり、ピーク値が低下しています。温度の高い領域から低い領域へ熱が流れ、分布が均一化に向かう様子がわかります。これは拡散方程式の本質的な性質です
- 最終時刻: さらに拡散が進み、温度分布はほぼ均一に近づいています。境界条件として $u = 0$ を課しているため、最終的には全領域で $u \to 0$ に収束します
等温線（等高線）が同心円状に広がっていることから、この拡散過程が等方的（方向によらず同じ速さで拡散する）であることも視覚的に確認できます。これは拡散係数 $\alpha$ が空間的に一様で、格子が等方的（$\Delta x = \Delta y$）であることに対応しています。
1次元の移流と2次元の拡散を個別に扱ってきました。いよいよ、これらを統合した本格的な流体シミュレーション――リッド駆動キャビティ流れに挑戦しましょう。
Python実装3: 2次元リッド駆動キャビティ流れ
問題設定
リッド駆動キャビティ流れ（Lid-driven cavity flow）は、CFDの最も有名なベンチマーク問題の一つです。正方形の箱（キャビティ）の上面（リッド）だけが一定速度で水平に動き、残りの3面は静止壁（滑りなし条件）です。この単純な設定にもかかわらず、レイノルズ数が上がるにつれて、主循環渦、コーナー渦、さらには非定常渦構造といった豊かな流体現象が現れます。
Ghia et al. (1982) による数値解がベンチマークとして広く使われており、自分のコードの検証に利用できます。
非圧縮性ナビエ・ストークス方程式の離散化
2次元非圧縮性流れの支配方程式を成分で書くと次のようになります。
$x$ 方向運動量方程式:
$$ \begin{equation} \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x} + \nu\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \end{equation} $$
$y$ 方向運動量方程式:
$$ \begin{equation} \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y} + \nu\left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) \end{equation} $$
連続の式:
$$ \begin{equation} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \end{equation} $$
ここで $u$, $v$ はそれぞれ $x$, $y$ 方向の速度成分、$p$ は圧力、$\nu$ は動粘性係数です。
圧力ポアソン方程式法
非圧縮性流れの数値解法で最も重要な課題は、「速度と圧力をどのように連成させるか」です。連続の式は速度場に対する拘束条件（発散がゼロ）であり、圧力はこの拘束条件を満たすように決まります。
本実装では、圧力ポアソン方程式法（Chorin のプロジェクション法の変形）を用います。手順は以下の通りです。
ステップ1: 運動量方程式の右辺（移流項 + 拡散項）を計算し、仮の速度場を求める
ステップ2: 連続の式を満たすように、圧力に関するポアソン方程式を解く
$$ \begin{equation} \nabla^2 p = \frac{\rho}{\Delta t}\left(\frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y}\right) \end{equation} $$
ここで $u^*$, $v^*$ はステップ1で得られた仮速度です。
ステップ3: 得られた圧力場を使って速度を補正し、発散ゼロの速度場を得る
$$ \begin{equation} u^{n+1} = u^* – \frac{\Delta t}{\rho}\frac{\partial p}{\partial x}, \quad v^{n+1} = v^* – \frac{\Delta t}{\rho}\frac{\partial p}{\partial y} \end{equation} $$
このプロセスを各時間ステップで繰り返すことで、非圧縮性流れの時間発展を追跡します。
境界条件
キャビティ流れの境界条件は次の通りです。
- 上面（リッド）: $u = U_{\text{lid}}$, $v = 0$（リッドが一定速度で右に動く）
- 下面・左面・右面: $u = 0$, $v = 0$（滑りなし条件、静止壁）
- 圧力: 壁面で $\partial p / \partial n = 0$（ノイマン条件）。ただし圧力の絶対値は不定なので、1点を $p = 0$ に固定します
実装コード
以下のコードは、上記のアルゴリズムをPythonで実装したものです。等間隔の構造格子上で差分法を適用し、圧力ポアソン方程式はヤコビ反復法で解きます。
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ設定 ---
nx, ny = 41, 41 # 格子点数
Lx, Ly = 1.0, 1.0 # キャビティサイズ [m]
dx = Lx / (nx - 1)
dy = Ly / (ny - 1)
rho = 1.0 # 密度 [kg/m^3]
nu = 0.01 # 動粘性係数 [m^2/s]
U_lid = 1.0 # リッド速度 [m/s]
dt = 0.001 # 時間刻み [s]
nt = 5000 # 時間ステップ数
nit = 50 # 圧力ポアソン方程式の反復回数
# レイノルズ数
Re = U_lid * Lx / nu
print(f"Reynolds number: Re = {Re:.0f}")
# 格子座標
x = np.linspace(0, Lx, nx)
y = np.linspace(0, Ly, ny)
X, Y = np.meshgrid(x, y)
# 変数の初期化
u = np.zeros((ny, nx)) # x方向速度
v = np.zeros((ny, nx)) # y方向速度
p = np.zeros((ny, nx)) # 圧力
def pressure_poisson(p, rhs, dx, dy, nit):
"""ヤコビ反復法で圧力ポアソン方程式を解く"""
for _ in range(nit):
pn = p.copy()
p[1:-1, 1:-1] = (
((pn[1:-1, 2:] + pn[1:-1, :-2]) * dy**2
+ (pn[2:, 1:-1] + pn[:-2, 1:-1]) * dx**2
- rhs[1:-1, 1:-1] * dx**2 * dy**2)
/ (2 * (dx**2 + dy**2))
)
# 圧力の境界条件 (ノイマン条件)
p[:, -1] = p[:, -2] # 右壁
p[:, 0] = p[:, 1] # 左壁
p[0, :] = p[1, :] # 下壁
p[-1, :] = 0 # 上面 (基準圧力)
return p
def cavity_flow(nt, u, v, p, dx, dy, dt, rho, nu, nit):
"""リッド駆動キャビティ流れのメインループ"""
for n in range(nt):
un = u.copy()
vn = v.copy()
# --- 移流項 + 拡散項の右辺を計算 ---
# 移流項 (中心差分)
dudx = (un[1:-1, 2:] - un[1:-1, :-2]) / (2 * dx)
dudy = (un[2:, 1:-1] - un[:-2, 1:-1]) / (2 * dy)
dvdx = (vn[1:-1, 2:] - vn[1:-1, :-2]) / (2 * dx)
dvdy = (vn[2:, 1:-1] - vn[:-2, 1:-1]) / (2 * dy)
# 拡散項 (中心差分)
d2udx2 = (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, :-2]) / dx**2
d2udy2 = (un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[:-2, 1:-1]) / dy**2
d2vdx2 = (vn[1:-1, 2:] - 2*vn[1:-1, 1:-1] + vn[1:-1, :-2]) / dx**2
d2vdy2 = (vn[2:, 1:-1] - 2*vn[1:-1, 1:-1] + vn[:-2, 1:-1]) / dy**2
# 仮速度場
u[1:-1, 1:-1] = (un[1:-1, 1:-1]
- dt * (un[1:-1, 1:-1] * dudx + vn[1:-1, 1:-1] * dudy)
+ dt * nu * (d2udx2 + d2udy2))
v[1:-1, 1:-1] = (vn[1:-1, 1:-1]
- dt * (un[1:-1, 1:-1] * dvdx + vn[1:-1, 1:-1] * dvdy)
+ dt * nu * (d2vdx2 + d2vdy2))
# --- 境界条件 (速度) ---
u[0, :] = 0 # 下壁
u[-1, :] = U_lid # 上面 (リッド)
u[:, 0] = 0 # 左壁
u[:, -1] = 0 # 右壁
v[0, :] = 0 # 下壁
v[-1, :] = 0 # 上面
v[:, 0] = 0 # 左壁
v[:, -1] = 0 # 右壁
# --- 圧力ポアソン方程式の右辺 ---
rhs = np.zeros_like(p)
rhs[1:-1, 1:-1] = (rho / dt
* ((u[1:-1, 2:] - u[1:-1, :-2]) / (2 * dx)
+ (v[2:, 1:-1] - v[:-2, 1:-1]) / (2 * dy)))
# --- 圧力を解く ---
p = pressure_poisson(p, rhs, dx, dy, nit)
# --- 速度の補正 ---
u[1:-1, 1:-1] -= dt / rho * (p[1:-1, 2:] - p[1:-1, :-2]) / (2 * dx)
v[1:-1, 1:-1] -= dt / rho * (p[2:, 1:-1] - p[:-2, 1:-1]) / (2 * dy)
# 速度境界条件の再適用
u[0, :] = 0; u[-1, :] = U_lid; u[:, 0] = 0; u[:, -1] = 0
v[0, :] = 0; v[-1, :] = 0; v[:, 0] = 0; v[:, -1] = 0
return u, v, p
# --- シミュレーション実行 ---
u, v, p = cavity_flow(nt, u, v, p, dx, dy, dt, rho, nu, nit)
print("シミュレーション完了")
このコードの構造を解説します。メインループ cavity_flow は各時間ステップで以下の処理を行っています。
- 移流項と拡散項の計算: 速度場の空間微分を中心差分で計算し、仮速度場を求めます。移流項 $(\bm{u} \cdot \nabla)\bm{u}$ は非線形項であり、速度自身が速度を輸送するという物理を表しています
- 速度境界条件の適用: 上面にリッド速度 $U_{\text{lid}}$ を、他の壁面にゼロ速度を設定します
- 圧力ポアソン方程式の求解: 仮速度場の発散から圧力の右辺を構成し、ヤコビ反復法で圧力場を求めます。反復回数
nitは精度と計算速度のトレードオフです - 速度の補正: 圧力勾配を用いて速度を修正し、連続の式（発散ゼロ条件）を近似的に満たす速度場を得ます
レイノルズ数 $\text{Re} = U_{\text{lid}} L / \nu = 100$ の設定では、流れは定常状態に収束します。$\text{Re}$ を上げると非定常な流れが現れますが、格子解像度と時間刻みも適切に調整する必要があります。
次に、得られた結果を可視化して物理的な解釈を行いましょう。
Python実装4: 結果の可視化（流線・速度場・圧力場）
流線と速度場の可視化
CFDの結果を可視化する際に最も重要なのは、流線と速度場です。流線は各点で速度ベクトルに接する曲線であり、流れのパターンを直感的に理解するために不可欠です。
import numpy as np
import matplotlib.pyplot as plt
# 前のコードブロックで得た u, v, p, X, Y を使用
# (再実行する場合は前のコードブロックを先に実行してください)
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
# --- (1) 流線図 ---
ax = axes[0, 0]
speed = np.sqrt(u**2 + v**2)
strm = ax.streamplot(X, Y, u, v, color=speed, cmap="coolwarm",
density=2, linewidth=1.5, arrowsize=1.2)
ax.set_title("Streamlines", fontsize=14)
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_aspect("equal")
plt.colorbar(strm.lines, ax=ax, label="Speed [m/s]")
# --- (2) 速度ベクトル場 ---
ax = axes[0, 1]
skip = 2 # 表示を間引く
ax.quiver(X[::skip, ::skip], Y[::skip, ::skip],
u[::skip, ::skip], v[::skip, ::skip],
speed[::skip, ::skip], cmap="coolwarm", scale=20)
ax.set_title("Velocity vectors", fontsize=14)
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_aspect("equal")
# --- (3) 圧力場 ---
ax = axes[1, 0]
cp = ax.contourf(X, Y, p, levels=30, cmap="RdBu_r")
ax.set_title("Pressure field", fontsize=14)
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_aspect("equal")
plt.colorbar(cp, ax=ax, label="Pressure [Pa]")
# --- (4) 中心線上の速度プロファイル ---
ax = axes[1, 1]
mid_x = nx // 2
mid_y = ny // 2
ax.plot(u[:, mid_x], y, "b-o", markersize=3, label="u at x = 0.5")
ax.plot(x, v[mid_y, :], "r-s", markersize=3, label="v at y = 0.5")
ax.set_xlabel("Velocity [m/s]", fontsize=12)
ax.set_ylabel("Position", fontsize=12)
ax.set_title("Centerline velocity profiles", fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.suptitle(f"Lid-Driven Cavity Flow (Re = {Re:.0f}, {nx}x{ny} grid)",
fontsize=15, y=1.01)
plt.tight_layout()
plt.savefig("cavity_flow_results.png", dpi=150, bbox_inches="tight")
plt.show()
上の4つのグラフから、リッド駆動キャビティ流れの物理的特徴を読み取りましょう。
-
流線図（左上）: キャビティ内に大きな時計回りの主循環渦が形成されていることがわかります。この渦はリッドの運動によって駆動されています。渦の中心はキャビティの幾何学的中心よりもやや上方に位置しており、これは Re = 100 における典型的な特徴です。また、右下と左下のコーナーにごく小さな二次渦が存在する兆候が見られます。レイノルズ数を上げると、これらのコーナー渦がより明確になります
-
速度ベクトル場（右上）: ベクトルの大きさと方向から、リッド（上面）近傍で速度が最大であること、壁面近傍で速度がゼロに減衰すること（滑りなし条件）がはっきりと確認できます。渦の中心付近では速度がほぼゼロになっています
-
圧力場（左下）: 圧力は上面の右上コーナーで高く、左上コーナーで低くなっています。これは直感的に理解できます。リッドが右に動くことで、流体は右壁にぶつかって圧縮され（高圧）、左壁から離れて膨張します（低圧）。この圧力差が渦循環を維持する駆動力の一部となっています
-
中心線速度プロファイル（右下）: 垂直中心線（$x = 0.5$）上の $u$ 速度は、上端で $U_{\text{lid}} = 1.0$、下端でゼロであり、中間に負の領域があります。これは主循環渦の下部で流れが左向き（リッドと逆方向）に戻ることを意味しています。水平中心線（$y = 0.5$）上の $v$ 速度は、渦の左右で符号が変わり、渦の循環構造を反映しています
渦度場の可視化
流れの渦構造をさらに詳しく見るために、渦度を計算して可視化します。2次元の渦度は速度場の回転として定義されます。
$$ \begin{equation} \omega = \frac{\partial v}{\partial x} – \frac{\partial u}{\partial y} \end{equation} $$
import numpy as np
import matplotlib.pyplot as plt
# 渦度の計算 (中心差分)
omega = np.zeros_like(u)
omega[1:-1, 1:-1] = ((v[1:-1, 2:] - v[1:-1, :-2]) / (2 * dx)
- (u[2:, 1:-1] - u[:-2, 1:-1]) / (2 * dy))
fig, ax = plt.subplots(figsize=(8, 7))
levels = np.linspace(-5, 5, 41)
cp = ax.contourf(X, Y, omega, levels=levels, cmap="RdBu_r", extend="both")
ax.contour(X, Y, omega, levels=levels, colors="k", linewidths=0.3, alpha=0.5)
ax.set_title(f"Vorticity field (Re = {Re:.0f})", fontsize=14)
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_aspect("equal")
plt.colorbar(cp, ax=ax, label=r"$\omega$ [1/s]")
plt.tight_layout()
plt.savefig("cavity_vorticity.png", dpi=150, bbox_inches="tight")
plt.show()
渦度場のグラフからは、以下の物理的特徴が読み取れます。
- リッド直下に強い負の渦度帯: リッドが右に動くことで生じる剪断層が、強い渦度を生成しています。渦度の絶対値はリッド近傍で最大であり、これは速度勾配 $\partial u / \partial y$ がリッド直下で最も大きいことに対応しています
- 壁面近傍の渦度: 左壁と右壁に沿って符号の異なる渦度が分布しており、壁面境界層の存在を示しています
- 渦中心付近のゼロ渦度: 主循環渦の中心では渦度がゼロに近くなります。これは剛体回転に近い領域では渦度が一定（この場合は弱い）であることを反映しています
収束性の確認
最後に、格子解像度を変えてシミュレーションを行い、解が格子に依存しないこと（格子収束性）を確認することはCFDにおいて極めて重要です。以下のコードでは、異なる格子点数での中心線速度プロファイルを比較します。
import numpy as np
import matplotlib.pyplot as plt
def run_cavity(nx_run, ny_run, nt_run, nit_run=50):
"""指定した格子点数でキャビティ流れを計算する"""
dx_r = 1.0 / (nx_run - 1)
dy_r = 1.0 / (ny_run - 1)
dt_r = 0.001
u_r = np.zeros((ny_run, nx_run))
v_r = np.zeros((ny_run, nx_run))
p_r = np.zeros((ny_run, nx_run))
for n in range(nt_run):
un = u_r.copy()
vn = v_r.copy()
# 移流 + 拡散
dudx = (un[1:-1, 2:] - un[1:-1, :-2]) / (2*dx_r)
dudy = (un[2:, 1:-1] - un[:-2, 1:-1]) / (2*dy_r)
dvdx = (vn[1:-1, 2:] - vn[1:-1, :-2]) / (2*dx_r)
dvdy = (vn[2:, 1:-1] - vn[:-2, 1:-1]) / (2*dy_r)
d2udx2 = (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, :-2]) / dx_r**2
d2udy2 = (un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[:-2, 1:-1]) / dy_r**2
d2vdx2 = (vn[1:-1, 2:] - 2*vn[1:-1, 1:-1] + vn[1:-1, :-2]) / dx_r**2
d2vdy2 = (vn[2:, 1:-1] - 2*vn[1:-1, 1:-1] + vn[:-2, 1:-1]) / dy_r**2
u_r[1:-1, 1:-1] = (un[1:-1, 1:-1]
- dt_r*(un[1:-1, 1:-1]*dudx + vn[1:-1, 1:-1]*dudy)
+ dt_r*nu*(d2udx2 + d2udy2))
v_r[1:-1, 1:-1] = (vn[1:-1, 1:-1]
- dt_r*(un[1:-1, 1:-1]*dvdx + vn[1:-1, 1:-1]*dvdy)
+ dt_r*nu*(d2vdx2 + d2vdy2))
u_r[0, :] = 0; u_r[-1, :] = U_lid; u_r[:, 0] = 0; u_r[:, -1] = 0
v_r[0, :] = 0; v_r[-1, :] = 0; v_r[:, 0] = 0; v_r[:, -1] = 0
rhs = np.zeros_like(p_r)
rhs[1:-1, 1:-1] = (rho/dt_r
* ((u_r[1:-1, 2:] - u_r[1:-1, :-2])/(2*dx_r)
+ (v_r[2:, 1:-1] - v_r[:-2, 1:-1])/(2*dy_r)))
for _ in range(nit_run):
pn = p_r.copy()
p_r[1:-1, 1:-1] = (
((pn[1:-1, 2:]+pn[1:-1, :-2])*dy_r**2
+(pn[2:, 1:-1]+pn[:-2, 1:-1])*dx_r**2
-rhs[1:-1, 1:-1]*dx_r**2*dy_r**2)
/ (2*(dx_r**2+dy_r**2)))
p_r[:, -1] = p_r[:, -2]
p_r[:, 0] = p_r[:, 1]
p_r[0, :] = p_r[1, :]
p_r[-1, :] = 0
u_r[1:-1, 1:-1] -= dt_r/rho*(p_r[1:-1, 2:]-p_r[1:-1, :-2])/(2*dx_r)
v_r[1:-1, 1:-1] -= dt_r/rho*(p_r[2:, 1:-1]-p_r[:-2, 1:-1])/(2*dy_r)
u_r[0, :] = 0; u_r[-1, :] = U_lid; u_r[:, 0] = 0; u_r[:, -1] = 0
v_r[0, :] = 0; v_r[-1, :] = 0; v_r[:, 0] = 0; v_r[:, -1] = 0
return u_r, v_r, p_r
# 異なる格子解像度で計算
grid_sizes = [21, 41, 61]
fig, ax = plt.subplots(figsize=(8, 6))
for ng in grid_sizes:
print(f"Computing with {ng}x{ng} grid...")
u_g, v_g, p_g = run_cavity(ng, ng, nt_run=5000)
y_g = np.linspace(0, 1, ng)
mid = ng // 2
ax.plot(u_g[:, mid], y_g, "-o", markersize=3,
label=f"{ng}x{ng} grid")
ax.set_xlabel("u velocity at x = 0.5", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_title(f"Grid convergence study (Re = {Re:.0f})", fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("cavity_grid_convergence.png", dpi=150, bbox_inches="tight")
plt.show()
格子収束性のグラフから、重要な知見が得られます。
- 21×21格子: 速度プロファイルは大まかな形状を捉えていますが、壁面近傍やリッド直下の急峻な速度勾配を十分に解像できていません。格子が粗すぎると、数値拡散が大きく物理的な詳細が失われます
- 41×41格子: 主要な流れの特徴が概ね再現されており、定性的にも定量的にも合理的な結果が得られています
- 61×61格子: 41×41の結果との差が小さくなっており、格子収束に近づいていることがわかります。理想的には、さらに細かい格子（81×81, 129×129）でも計算を行い、解がほぼ変化しなくなる（格子独立な解に到達する）ことを確認すべきです
格子収束性の検証は、CFDにおいて必須の手順です。格子に依存する結果は信頼できません。実際のCFD研究では、少なくとも3種類の格子解像度で計算を行い、リチャードソン外挿などの手法で格子独立解を推定します。
差分法の高次精度スキーム
なぜ高次精度が必要か
ここまでの実装で使ってきたのは2次精度の中心差分と1次精度の上流差分でした。実用的なCFDでは、計算コストと精度のバランスから、4次以上の高次精度スキームが使われることがあります。
高次精度スキームの考え方は、テイラー展開をより多くの項まで考慮し、より多くの近傍格子点の値を使って差分を構成するというものです。例えば、4次精度の中心差分は次のようになります。
$$ \begin{equation} \frac{df}{dx}\bigg|_i \approx \frac{-f_{i+2} + 8f_{i+1} – 8f_{i-1} + f_{i-2}}{12\Delta x} \end{equation} $$
打ち切り誤差は $O((\Delta x)^4)$ です。2次精度では格子間隔を半分にすると誤差が $1/4$ になりますが、4次精度では $1/16$ になります。つまり、同じ精度を達成するのに、格子点数を大幅に減らせる可能性があります。
移流項の風上スキーム
移流項 $u \partial f / \partial x$ の離散化では、中心差分は安定性の問題を抱えることがあります（チェッカーボード不安定性）。そのため、流れの方向を考慮した風上スキーム（upwind scheme）が広く使われています。
1次精度の風上差分は先ほどの1次元移流方程式で使いましたが、数値拡散が大きいという問題がありました。この問題を軽減するために、2次精度以上の風上型スキームが開発されています。代表的なものには、QUICK（Quadratic Upstream Interpolation for Convective Kinematics）スキームや、MUSCL（Monotonic Upstream-centered Scheme for Conservation Laws）スキームがあります。
これらの高次精度スキームを使う際には、不連続面や急峻な勾配の近くで数値振動が生じるギブス現象に注意が必要です。TVD（Total Variation Diminishing）やWENO（Weighted Essentially Non-Oscillatory）といった手法は、高次精度を維持しつつ振動を抑える工夫を施しています。
高次精度スキームと安定化手法の選択は、対象とする流れの特性に大きく依存します。ここまでの議論で、CFDの基本的なツールキットは一通り揃いました。最後に、実用的なCFDソフトウェアとの接点について触れておきましょう。
実用的なCFDソフトウェアとの接点
オープンソースCFDツール
本記事ではPythonでスクラッチ実装を行いましたが、実用的なCFD解析では専用のソフトウェアを使うのが一般的です。代表的なオープンソースCFDソフトウェアとしては以下があります。
- OpenFOAM: 最も広く使われているオープンソースCFDツールキット。有限体積法を基盤とし、非圧縮性/圧縮性流れ、乱流モデル、多相流、燃焼など幅広い物理をカバーしています
- SU2: スタンフォード大学で開発された、航空宇宙分野向けのCFDソルバー。形状最適化機能が充実しています
- FEniCS: 有限要素法ベースの汎用PDEソルバーフレームワーク。Pythonインターフェースを持ち、弱形式からの自動コード生成が特徴です
スクラッチ実装の意義
商用・オープンソースの優れたツールが存在する中で、なぜスクラッチ実装が重要なのでしょうか。
第一に、ブラックボックスの中身を理解するためです。本記事で実装した圧力ポアソン方程式法のアルゴリズムは、OpenFOAMなどの中核部分に相当します。自分で書くことで、「なぜ圧力を反復で解く必要があるのか」「CFL条件を満たさないと何が起こるのか」を体感として理解できます。
第二に、結果の妥当性を判断する力を養うためです。CFDソフトウェアは正しい設定で使えば正確な結果を出しますが、設定を間違えれば物理的に無意味な結果を出力します。基礎を理解していなければ、その結果が正しいのかどうかを判断できません。
第三に、新しいアルゴリズムを開発するための基盤となります。研究の最前線では、既存のソフトウェアにない手法を開発する必要がしばしばあります。その際に、離散化手法の基本原理を深く理解していることが不可欠です。
まとめ
本記事では、数値流体力学（CFD）の基礎を体系的に解説しました。
- CFDの3ステップ: 支配方程式の選定 → 離散化 → 求解と後処理。この枠組みはどのようなCFD手法にも共通する普遍的な構造です
- 差分法: テイラー展開を使って微分を差分で近似する手法です。前進差分・後退差分は1次精度、中心差分は2次精度であることを導出しました
- 時間積分: 陽的スキームは実装が簡単ですがCFL条件による制約を受けます。陰的スキームは無条件安定ですが連立方程式の求解が必要です。クランク・ニコルソン法は両者の利点を兼ね備えています
- CFL条件: 陽的スキームの安定性を保証する条件であり、「波が1時間ステップで格子点を飛び越えてはならない」という直感的な意味を持ちます。フォン・ノイマンの安定性解析から数学的に導出しました
- 有限体積法: 保存則の積分形を直接離散化する手法であり、離散レベルでの保存性が厳密に保たれます。実用的なCFDコードの大多数がこの手法を採用しています
- リッド駆動キャビティ流れ: CFDの古典的なベンチマーク問題をPythonでスクラッチ実装し、主循環渦、圧力分布、渦度場などの物理的特徴を可視化・考察しました
- 格子収束性: 格子解像度を変えて解の収束を確認することがCFDでは必須であることを示しました
CFDは流体力学の理論と計算科学の融合であり、本記事で扱った内容はその入口にすぎません。今後の学習として、以下のトピックに進むことをお勧めします。
次のステップとして、以下の記事も参考にしてください。
- ナビエ・ストークス方程式 — CFDの支配方程式の詳しい導出と物理的意味
- 境界層理論 — 壁面近傍の流れの薄い領域の理論。CFDでの格子設計に直結します
- 渦度と循環 — 渦の定量的な記述。キャビティ流れの渦構造の理解を深めます
- FEM入門 — 有限要素法による離散化手法。構造解析やマルチフィジックスでのCFDとの連成に重要です
- 流線・流跡線 — 流れの可視化手法の分類と物理的意味