クーラン数 シミュレーター 戻る
CFD 基礎シミュレーター

クーラン数 シミュレーター — 陽的時間積分の安定性

クーラン数 $\mathrm{CFL} = u\,\Delta t/\Delta x$ をインタラクティブに計算します。移流速度 $u$、時間ステップ $\Delta t$、空間刻み $\Delta x$、最大許容 $\mathrm{CFL}_\mathrm{max}$ を変えて、安定・境界・不安定の判定、最大許容 $\Delta t$、最大許容 $u$、信号伝播の模式図を同時に可視化します。

パラメータ設定
移流速度 u
m/s
時間ステップ Δt
ms
空間刻み Δx
mm
最大許容 CFL_max

CFL 条件:$\mathrm{CFL} = u\,\Delta t/\Delta x \le \mathrm{CFL}_\mathrm{max}$。代表値は Forward Euler / Lax-Wendroff で $\mathrm{CFL}_\mathrm{max} = 1.0$、RK4 で約 2.06。

計算結果
CFL 数
安定性
最大許容 Δt
最大許容 u
1D グリッドと信号伝播

青枠=空間セル(幅 $\Delta x$)/黄矢印=1 ステップで信号が進む距離 $u\,\Delta t$/緑=CFL ≤ 1 の安定域/赤=CFL > 1 でセルを跨ぐ不安定域。矢印が 1 セルを超えると陽解法が破綻します。

安定性マップ(CFL vs 安定指標)

横軸=CFL [0, 3]/縦軸=安定指標 $\mathrm{CFL}_\mathrm{max} - \mathrm{CFL}$/緑帯=安定域(CFL < 0.9·CFL_max)/黄帯=境界域/赤帯=不安定域(CFL > CFL_max)/黄●=現在の CFL。

理論・主要公式

クーラン数(CFL 数)の定義:

$$\mathrm{CFL} = \frac{u\,\Delta t}{\Delta x}$$

陽的スキームの安定条件と最大許容 $\Delta t$・$u$:

$$\mathrm{CFL} \le \mathrm{CFL}_\mathrm{max},\quad \Delta t_\mathrm{max} = \frac{\mathrm{CFL}_\mathrm{max}\,\Delta x}{u},\quad u_\mathrm{max} = \frac{\mathrm{CFL}_\mathrm{max}\,\Delta x}{\Delta t}$$

安定指標(負なら不安定):

$$\mathrm{margin} = \mathrm{CFL}_\mathrm{max} - \mathrm{CFL}$$

$u$ は移流速度 [m/s]、$\Delta t$ は時間ステップ [s]、$\Delta x$ は空間刻み [m]、$\mathrm{CFL}_\mathrm{max}$ はスキーム固有の安定限界(Forward Euler・Lax-Wendroff で 1.0、RK4 で約 2.06)。陽的時間積分では $\mathrm{CFL} \le \mathrm{CFL}_\mathrm{max}$ が安定の必要条件です。

クーラン数とは

🙋
CFD の本で「CFL ≤ 1 にしないと発散する」って書いてあるんですけど、CFL って結局何ですか? なぜ 1 で線が引かれるんですか?
🎓
クーラン数 $\mathrm{CFL} = u\,\Delta t/\Delta x$ は「1 ステップで信号が進む距離 $u\,\Delta t$ が、セル幅 $\Delta x$ の何倍にあたるか」という比だよ。本ツールの既定値($u=10$ m/s、$\Delta t=1.00$ ms、$\Delta x=10.0$ mm)だと、$u\,\Delta t = 10 \cdot 10^{-3} = 10$ mm でちょうど 1 セル分。CFL = 10/10 = 1.000 になる。これが「境界(marginal)」で、これを超えると陽解法では破綻するんだ。
🙋
$\Delta t$ を 2 ms に上げると、CFL が 2.000 になって赤い「不安定」に変わりました。実際の解析だとどうなるんですか?
🎓
数値解が指数関数的に発散して、波の振幅が桁違いに大きくなる。「振動が止まらない」「Inf や NaN が出る」「結果が物理的にあり得ない値になる」といった現象が起きる。原因は、Forward Euler の増幅因子が $|1 - \mathrm{CFL}| > 1$ になり、毎ステップで誤差が増幅されてしまうから。von Neumann 安定性解析という線形理論で厳密に証明されていて、1928 年に Courant・Friedrichs・Lewy の 3 人が論文で示したのが起源だ。
🙋
$\Delta x$ を 20 mm に倍にすると CFL が 0.500 になって緑の「安定」になりました。逆に細かいメッシュにすると安定限界が厳しくなるんですね…
🎓
そうそう、これが「CFL の呪縛」って言われる課題で、解像度を上げたい($\Delta x$ を小さくしたい)と必ず $\Delta t$ も小さくしないといけない。例えば $\Delta x$ を 1/2 にしたら $\Delta t$ も 1/2、計算量は時間方向で 2 倍、空間方向で 2 倍、合計 4 倍に増える。3D だと $2^4 = 16$ 倍。だから DNS(直接数値計算)は莫大な計算資源が必要で、LES や RANS でモデル化するメリットが大きい。
🙋
CFL_max を 2.06 に上げると、CFL = 1.000 のまま緑(安定)になりました。これが RK4 ですか?
🎓
そう、RK4 は 4 段階の中間ステップで誤差を相殺する高次スキームで、実効的な安定限界が約 2.06 まで広がる。圧縮性流れの DG-FEM や Spectral Element でも RK4 系がよく使われる。ただし計算コストは Forward Euler の 4 倍、Lax-Wendroff の 2 倍。安定領域の広さと 1 ステップあたりのコストのバランスでスキームを選ぶのが実務だ。MUSCL + RK4 とか TVD-RK3 とか、組み合わせの工夫も多い。

よくある質問

多次元では各方向の CFL の和(または対角項)が制約されます。2D 移流では $\mathrm{CFL}_x + \mathrm{CFL}_y = u\,\Delta t/\Delta x + v\,\Delta t/\Delta y \le 1$、3D では $\mathrm{CFL}_x + \mathrm{CFL}_y + \mathrm{CFL}_z \le 1$ が陽解法 1 次精度の典型的条件です。等方メッシュで $u = v = w$ の場合、3D の $\Delta t$ は 1D の 1/3 まで小さくしないといけません。次元分割(Strang Splitting)を使えば各方向独立に評価できますが、それでも一番厳しい方向で制約されます。$\mathrm{CFL}_\mathrm{max}$ は高次スキーム(Lax-Wendroff、MacCormack、RK4)でやや緩和されますが、本質的に多次元では 1 を超えるのが難しいです。
対流拡散方程式 $\partial \phi/\partial t + u\,\partial \phi/\partial x = \nu\,\partial^2 \phi/\partial x^2$ では 2 つの安定条件が同時に必要で、$\mathrm{CFL} = u\,\Delta t/\Delta x \le 1$ と拡散数 $d = \nu\,\Delta t/\Delta x^2 \le 0.5$(陽的 1 次精度)の厳しい方が支配します。さらに対流+拡散の組み合わせ条件として $\mathrm{CFL}^2 \le 2d$ と書ける場合もあり、移流支配(高 Re)では CFL、拡散支配(低 Re)では d で制約されます。実用 CFD では拡散項を陰的に、対流項を陽的に解く半陰解法(IMEX)も多用され、対流の CFL だけ気にすればよくなります。
圧縮性 Euler/Navier-Stokes では、特性速度として $|u| + a$($a$ は局所音速)を使い、$\mathrm{CFL} = (|u| + a)\,\Delta t/\Delta x \le 1$ になります。亜音速流れ($M = u/a \ll 1$)では音速 $a$ が支配的で、低 Mach 流れほど $\Delta t$ が小さくなるのが圧縮性 CFD の悩みです。例えば空気中の $u = 10$ m/s、$a = 340$ m/s では特性速度は 350 m/s、$\Delta x = 10$ mm なら $\Delta t \le 28.6$ μs と非常に小さい値になります。低 Mach 流れには Preconditioning(Turkel、Choi-Merkle)や圧縮性/非圧縮性ハイブリッドソルバが用いられます。
毎ステップ全セルで最大 CFL を計算し、目標 CFL_target(典型は 0.5〜0.8)に対して次ステップの $\Delta t_\mathrm{new} = \mathrm{CFL}_\mathrm{target}/\mathrm{CFL}_\mathrm{current} \cdot \Delta t_\mathrm{current}$ を更新するのが基本です。OpenFOAM の adjustTimeStep 機能、ANSYS Fluent の CFL-based time step、SU2 のシステムなど主要 CFD ソルバは同じ思想で動いています。注意点として、$\Delta t$ を毎ステップ大きく変えると 2 次精度の時間積分(BDF2 など)の精度が落ちるので、増減率を 1 ステップで 1.2 倍までに制限することが多いです。また、衝撃波が発生したり境界条件が急変したりすると CFL が急増するので、上限・下限のクリップも実装します。

実世界での応用

商用 CFD ソルバの時間ステップ設計:ANSYS Fluent、STAR-CCM+、OpenFOAM、SU2 など主要 CFD ソルバはすべて CFL 条件を内部で監視しています。陽的時間積分(コンプレッシブル流れ、爆発・衝撃解析)では CFL = 0.3〜0.7 を目標とし、陰的時間積分(インコンプレッシブル流れの定常解析)では CFL = 5〜100 を狙うのが標準的です。本ツールで $u = 10$ m/s、$\Delta x = 10$ mm のとき $\Delta t_\mathrm{max} = 1$ ms と計算できることから、実用解析では LES でも $\Delta t \approx 0.5$ ms 程度を CFL_target = 0.5 で設定する流れがイメージできます。

気象・海洋モデルのスーパーコンピューティング:WRF(米気象数値予報モデル)、MPAS、NICAM(日本の全球非静水力モデル)など気象モデルは典型的な陽解法 CFD で、CFL 条件が計算コストを決定します。水平 $\Delta x = 1$ km の予報モデルで風速 $u = 50$ m/s の場合、$\Delta t \le 20$ 秒に制限されます。地球シミュレータや「富岳」のような大型計算機を用いても 24 時間先までの予報計算に数時間かかるのは CFL 条件のためです。

地震波伝播・FDTD 電磁解析:FDTD(時間領域差分)法は電磁界解析・地震波伝播・音響解析で広く用いられ、CFL 条件は Courant 数として $c\,\Delta t \le \Delta x/\sqrt{N}$($N$ は次元、$c$ は伝播速度)の形で現れます。3D FDTD では $c\,\Delta t \le \Delta x/\sqrt{3} \approx 0.577\,\Delta x$ が必要で、光速 $c \approx 3 \times 10^8$ m/s、$\Delta x = 1$ mm のとき $\Delta t \le 1.93$ ps と極端に小さくなります。RF・アンテナ解析でメッシュを細かくできない大きな理由がここにあります。

衝撃波・爆発解析(ALE・SPH):ハイドロコード(LS-DYNA、AUTODYN、Abaqus/Explicit)の爆発・貫通解析では、特性速度に音速+粒子速度を採用するため、衝撃面付近で局所 CFL が急増しがちです。実務では CFL_target = 0.6〜0.7 を維持し、衝撃面追従の自動再分割(adaptive mesh refinement)と組み合わせて精度・コストのバランスを取ります。航空機の鳥衝突解析、自動車衝突、防爆ブロック設計などで標準的に使われる手法です。

よくある誤解と注意点

最も多い誤解は、「CFL ≤ 1 なら必ず精度が高い」と思い込むことです。CFL は安定の必要条件であって、精度の保証ではありません。CFL = 0.001 のような極端に小さい値だと、数値拡散が増大して鋭い波形が鈍ります。逆に CFL = 0.7〜1.0 の領域では Lax-Wendroff など 2 次精度スキームで波形保存性が最良になることが知られています。「安定 ≠ 精度」を意識し、メッシュ解像度($\Delta x$)と時間刻み($\Delta t$)を両方バランスよく設定する必要があります。

次に多いのが、「陰解法を使えば $\Delta t$ をいくら大きくしても良い」という思い込みです。Backward Euler や BDF2 は無条件安定ですが、$\Delta t$ が大きいと数値拡散が支配的になり、非定常流れの渦構造が消えてしまいます。LES で大規模渦を解像したい場合、$\Delta t$ は流れの最小時間スケール(コルモゴロフ時間スケール)以下に取らないと意味がありません。定常解析なら CFL = 100 でも問題ありませんが、過渡解析では CFL = 1〜5 程度に抑え、必要に応じて時間刻みの収束性検証(Δt を半分にして解が変わらないことを確認)を行うのが鉄則です。

最後に、「CFL は線形理論だから非線形問題では関係ない」という誤解です。Burgers 方程式や Euler 方程式のような非線形対流問題では、衝撃波形成に伴って局所的に特性速度が急増し、見かけの CFL がジャンプします。理論上の CFL = 1 を狙うと衝撃波付近で発散することが多く、実務では CFL = 0.3〜0.5 と保守的に設定するのが安全です。さらに圧縮性流れでは音速 $a$ も含めた特性 CFL = $(|u|+a)\,\Delta t/\Delta x$ を用い、本ツールの単純な定義 $u\,\Delta t/\Delta x$ とは区別する必要があります。線形 CFL は教育用ツールでのみ十分で、実 CFD では各ソルバの内部定義を必ず確認してください。