渦度方程式 — 数値解法と実装

2次元非圧縮流れでは渦度-流れ関数法（vorticity-stream function method）が古典的な手法だ。渦度方程式と流れ関数のPoisson方程式を連立させる。

速度は $u = \frac{\partial \psi}{\partial y}$、$v = -\frac{\partial \psi}{\partial x}$ で復元する。圧力が未知数から消えるのがこの手法のメリットだよ。

2次元では未知数が $\omega$ と $\psi$ の2つだけで、圧力-速度連成の問題（SIMPLE法などの反復が不要）を回避できる。ただし3次元への拡張は渦度が3成分ベクトルになるため複雑で、実用的にはあまり使われない。

有限差分法（FDM）での典型的な離散化を見てみよう。均一格子間隔 $h$ の中心差分を用いると、Poisson方程式は

移流項は中心差分だと格子Reynolds数 $Re_h = |u|h/\nu > 2$ でウィグルが発生する。これを避けるために風上差分やQUICK法を使うことが多い。

FVMではセル界面での渦度フラックスを評価する。移流フラックスには2次精度風上スキーム（Second Order Upwind）やTVDスキーム、拡散フラックスには中心差分を使う。OpenFOAMのscalarTransportFoamソルバーが渦度の移流拡散をFVMで解く実装の参考になるよ。

代表的な時間積分スキームをまとめよう。

実用的にはCrank-Nicolsonで拡散項を陰的に、移流項をAdams-Bashforthで陽的に扱う半陰的手法がよく使われる。

もちろん。移流のCFL条件は $\Delta t < h/|u_{max}|$ だ。拡散の安定性条件は $\Delta t < h^2/(4\nu)$ で、高Reynolds数では移流のCFLが支配的、低Reynolds数では拡散の制約が厳しくなる。

壁面での渦度は直接与えられないから、流れ関数から間接的に導出する。代表的な手法はThomの公式だ。壁面（$j=0$）で滑りなし条件 $u=v=0$ と $\psi=0$ を使うと

これは2次精度で、反復法の中でPoisson方程式と交互に更新する必要がある。

Woods の公式は $\omega_{i,0} = -\frac{3\psi_{i,1}}{h^2} - \frac{\omega_{i,1}}{2}$ で3次精度になる。またBriley の公式は壁面近傍の渦度分布をTaylor展開して4次精度を達成する。精度が上がるほど収束が遅くなる傾向があるので、実用上はThomの公式で十分な場合が多い。

その通りで、ここがボトルネックになる。代表的な解法は以下だ。

• SOR法: 実装が簡単。緩和係数 $\omega_{SOR} \approx 2 - \pi h$ が最適
• ADI法: 交互方向陰解法。各方向の三重対角行列をThomas法で解く
• マルチグリッド法: $O(N)$ の計算量で最も高速。粗いグリッドで低周波誤差を効率的に除去
• FFT法: 矩形領域で一様格子の場合に $O(N \log N)$ で直接解ける

格子数が $100 \times 100$ 程度ならSOR法で十分だ。$1000 \times 1000$ 以上ならマルチグリッド法が圧倒的に速い。OpenFOAMのfvMatrixではGAMGソルバーがマルチグリッド法に対応している。