流れ関数 — 数値解法と実装

アルゴリズムの流れはこうだ。各タイムステップで以下を繰り返す。

1. 渦度方程式を時間積分して $\omega^{n+1}$ を求める

2. Poisson方程式 $\nabla^2 \psi^{n+1} = -\omega^{n+1}$ を解いて流れ関数を更新

3. 速度場を $u = \partial\psi/\partial y$、$v = -\partial\psi/\partial x$ で計算

4. 壁面渦度を更新してステップ1に戻る

均一格子間隔 $h$ のスタガード格子で、Poisson方程式の5点差分は

速度の中心差分は

中心差分だと格子Reynolds数 $Re_h = |u|h/\nu$ が2を超えると振動する。実用的な選択肢を比較しよう。

実際、QUICK（Quadratic Upstream Interpolation for Convective Kinematics）は多くの問題でバランスが良い。ただし単調性が保証されないので、急な不連続がある場合はTVD（Total Variation Diminishing）制限子を組み合わせるとよい。

そうだ。代表的な反復法の収束速度を比較しよう。$N \times N$ 格子の場合。

マルチグリッドのV-cycleは、粗い格子で低周波誤差を効率的に除去する。実装はやや複雑だが、$256 \times 256$ 以上の格子では他の手法と桁違いの速度差になるよ。PythonならPyAMGライブラリで簡単にAMG（Algebraic Multigrid）を使えるんだ。

代表的な境界条件の実装をまとめよう。

静止壁面: $\psi = 0$（Dirichlet）、$\omega_{wall}$ はThomの公式 $\omega_{i,0} = -2\psi_{i,1}/h^2$ で計算

移動壁面（速度 $U$ で移動）: $\psi = 0$、$\omega_{i,0} = -2(\psi_{i,1})/h^2 - 2U/h$

流入境界: 速度分布を指定し、$\psi$ を壁面から積分して求める。例えば一様流 $U_\infty$ なら $\psi(y) = U_\infty y$

流出境界: $\partial^2 \psi / \partial x^2 = 0$ や $\partial \omega / \partial x = 0$ のNeumann条件が一般的

自由面: $\psi = \text{const.}$、$\omega = 0$（せん断応力ゼロ）

よく変わるよ。流出境界が渦の存在する領域に近すぎると、非物理的な反射が起きて結果がおかしくなる。流出境界は関心領域から十分離すこと（最低でも特性長さの10倍以上）が鉄則だ。

Python + NumPyで渦度-流れ関数法のLid-driven cavity問題を解くなら、200行程度で書ける。核心部分は以下のような構造だ。

• SOR法でPoisson方程式を解く関数（約30行）
• 移流拡散方程式の時間積分関数（約20行）
• 壁面渦度の更新関数（約15行）
• メインの時間ループ（約20行）

Re=100のLid-driven cavityなら $64 \times 64$ 格子で十分な精度が得られる。計算時間はノートPCで数秒程度だよ。

その姿勢はとても良い。自分でコードを書くと、離散化スキームの選択が結果にどう影響するか体感できる。教科書で読むだけでは身につかない感覚が養えるんだ。