リッド駆動キャビティ流れ — CFDベンチマーク検証ガイド

ざっくり言うと、CFDコードの検証で最も広く使われるベンチマーク問題だ。正方形の箱（キャビティ）の上壁だけが一定速度で水平にスライドすると、内部に再循環流れが生まれる。この流れ場を計算して、Ghia, Ghia & Shin (1982) が高精度差分法で求めた参照解と比較することで、ソルバーが正しく実装されているかチェックする。

いい質問だ。3つの理由がある。第一に、形状が極めて単純（正方形、壁だけ）なのでメッシュ生成の曖昧さがない。第二に、非自明な流れ場が得られる——再循環渦や角部の二次渦など、対流と粘性の競合を含む。管内流のような平行流とは違って、離散化スキームの質が結果に如実に表れる。第三に、Ghiaの論文がRe=100から10000まで7条件の速度プロファイルを数値テーブルとして公開してくれたおかげで、定量的比較が容易にできるんだ。

2次元・非圧縮・定常のNavier-Stokes方程式だ。連続の式と運動量方程式をセットで解く：

境界条件は以下の通りだ：

Ghiaの原論文がまさにそれだ。流れ関数 $\psi$ と渦度 $\omega$ を使って、圧力項を消去する：

Ghiaは $129 \times 129$ の等間隔格子で多重格子法を用いて計算した。この手法は圧力方程式が不要だから2D問題に向いているが、3D拡張が難しい。現在の商用CFDソルバーは速度-圧力形式（primitive variable法）が主流だ。

Re数が流れ場の性質を決定的に支配する。段階的に見ていこう：

2Dの場合、Re=10000まではGhiaの定常解が得られている。ただし収束は非常に遅い。Erturkら(2005)の研究では、Re=10000での定常解に $512 \times 512$ 格子で数十万ステップの反復が必要だった。Re=20000あたりからHopf分岐が起こり、周期的な流れに遷移するとされている。3Dの場合は遥かに低いRe数（およそRe=1000付近）で非定常性が表れることもある。

Ghia et al.は、キャビティ中心を通る 垂直線上の水平速度 $u$ と 水平線上の垂直速度 $v$ を17点ずつ表で提示している。Re=1000での垂直中心線上の $u$ 速度プロファイルが最もよく使われるデータだ：

この表を見ると、$y \approx 0.17$ 付近で $u/U \approx -0.38$ の最小値を取る。これは一次渦の下側で逆流が最も強い位置だ。自分のCFD結果をこのプロファイルと重ねてプロットし、全ポイントで相対誤差1%以内を達成できれば、ソルバーの実装は信頼に値する。

もちろん。Re=1000での一次渦中心の流れ関数値は $\psi_{\min} = -0.1189$ だ。これはスカラー値1つで検証できるので便利だが、速度プロファイル全体の比較のほうがより厳密な検証になる。

定常解を求めるなら SIMPLE（Semi-Implicit Method for Pressure-Linked Equations）が定番だ。各反復で圧力補正方程式を1回解き、速度場を更新する。亜緩和係数を適切に設定すれば安定に収束する。Re=1000程度なら速度の亜緩和0.7、圧力の亜緩和0.3くらいが出発点になる。

PISO（Pressure-Implicit with Splitting of Operators）は非定常解析向けだ。各時間ステップで圧力補正を2回行うことで、時間精度を確保する。Re=10000以上で非定常解が必要な場合に使う。キャビティ流れでRe=1000の定常解を求めるだけなら、SIMPLEのほうが計算効率が良い。

キャビティ流れではスキーム選択の影響が顕著に出る。具体例を見せよう。Re=1000で同じメッシュ（64×64）を使った場合：

1次風上法は数値拡散（人工粘性）が大きく、32×32では7%以上の誤差が出る。実務ではRe=1000の場合、最低でも2次精度スキーム＋64×64以上を使うべきだ。QUICKスキーム（3次精度上流内挿）なら64×64でも十分な精度が得られる。

検証目的には不適切だ。ただし、収束性が非常に良いので、初期解の推定や、2次精度スキームが発散する場合のデバッグ用には使える。「まず1次風上で回して流れ場の概形を確認し、その結果を初期値として2次精度に切り替える」という段階的アプローチは実務でもよく使われる。

Re数による。Re=100〜400なら等間隔で十分だが、Re=1000以上では壁面近傍に急な速度勾配が生じるから、壁面に向かって幾何級数的に格子を密にする不等間隔メッシュが有効だ。典型的には、膨張比（expansion ratio）1.05〜1.1で壁面の最初のセル厚を全体の1/200程度にする。

少なくとも3水準のメッシュ（例えば32×32, 64×64, 128×128）で同じ問題を解き、注目する物理量（$\psi_{\min}$、中心線速度など）の変化を確認する。Richardsonの外挿法で収束次数を評価し、GCI（Grid Convergence Index）を算出するのが標準的だ。

ここで $f_1, f_2, f_3$ はそれぞれ密・中・粗メッシュの結果、$r$ はメッシュ密度比（通常2）、$F_s$ は安全係数（3水準なら1.25）だ。2次精度スキームを使っていれば、収束次数 $p$ はおよそ2になるはずだ。$p$ が2から大きくかけ離れていたら、何かおかしい。

実はOpenFOAMにはチュートリアルケースとして cavity が同梱されている。$FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity にあるから、まずはそこから始めるといい。ただしデフォルトは20×20メッシュでRe=10と非常に低いので、検証用途にはカスタマイズが必要だ。

Re=1000でやるための設定ポイント：

divSchemes に Gauss linearUpwind grad(U) を使えば2次精度上流内挿になる。Gauss linear（中心差分）でもいいが、粗いメッシュでは振動することがある。Gauss upwind（1次風上）は検証用途にはおすすめしない。laplacianSchemes は Gauss linear corrected で2次精度の拡散項になる。

Fluentでの手順はこうだ：

Fluentなら Solution > Results > Plots > XY Plot で、垂直中心線（$x=0.5$）上の $u$ 速度、水平中心線（$y=0.5$）上の $v$ 速度をプロットする。「Write to File」で数値データをCSV出力すれば、ExcelやPythonでGhiaのデータと重ね合わせプロットできる。

以下の5ステップで検証するのが標準的な手順だ：

複数のRe数で検証できるとさらに良い。Re=100で正しくても、Re=1000で対流スキームの問題が顕在化することがある。最低でも Re=100とRe=1000の2条件 で検証することを勧める。

まともなソルバーなら、適切な設定で Ghia値との誤差1%以内 に収まるはずだ。差が出るのはむしろ「使い方」の問題で、スキーム選択・メッシュ品質・収束判定の設定で結果が変わる。ソルバー間クロスチェックの実例を見せよう：

有限要素法（FEM）ベースのCOMSOLやFEniCSでは、Taylor-Hood要素（P2速度＋P1圧力）を使えば比較的粗いメッシュでも高精度が得られる。これは速度場の多項式次数が高いからだ。有限体積法（FVM）のソルバーは同等の精度を得るのにより細かいメッシュが必要になるが、計算コストあたりの効率は高い。

目的に応じて使い分けるべきだ：

3Dキャビティ（立方体、リッドが上面全体でスライド）では、スパン方向の端壁の影響が加わる。2Dでは見えなかったTaylor-Görtler風の縦渦が端壁近傍に発生し、流れ場が大幅に複雑になる。Ku, Hirsh & Taylor (1987) やKoseff & Street (1984) の風洞実験データが3D検証の参照として使われる。

128×128の2Dで約16,000セルに対し、128×128×128の3Dでは約200万セルになる。メモリも計算時間も桁違いに増える。3D検証は最初のステップとしてはヘビーすぎるので、まず2Dで検証を通してから3Dに進むのがセオリーだ。

2Dの場合、Re=10000〜20000付近でHopf分岐が起こり、定常解が不安定になる。周期的な渦の振動が発生し、さらにRe数を上げるとカオス的な挙動に遷移する。3Dでは遥かに低いRe数（Re=1000〜3000程度）で3次元不安定性が現れる。

高Re数のキャビティ流れを解く手法としては：

LBMはNavier-Stokes方程式を直接解くのではなく、粒子の衝突・移流を格子上でシミュレーションする手法だ。キャビティ流れはLBMのベンチマークとしても定番で、D2Q9モデル（2次元9速度モデル）でGhia値を良好に再現できることが知られている。

LBMの利点は並列化の容易さだ。各格子点の計算が局所的なので、GPUとの相性が非常に良い。近年はNVIDIAが推進するGPUベースのLBMソルバー（例えばWaLBerla、Palabos）でキャビティ流れの高速計算が報告されている。ただし高Re数での安定性確保にはMRT（多緩和時間）モデルなどの工夫が必要だ。

よくあるケースだ。原因を順番にチェックしよう：

64×64で2〜3%の誤差は実はメッシュ不足として正常な範囲だ。心配なら128×128に上げてみて、誤差が1%未満になれば問題ない。GCIを計算してメッシュ依存性を定量的に評価するのが確実だ。

高Re数では対流が強く、収束が困難になる。以下の対策を試してみてくれ：

角部の速度不連続は数学的な特異点で、正常な挙動だ。上壁の角では速度が $(U,0)$ から $(0,0)$ に瞬間的に変化する。厳密解では圧力と渦度が無限大に発散する。メッシュを細かくすると、その近傍の圧力値はどんどん大きくなっていく——これは「メッシュ収束しない」ように見えるが、特異点の近傍ではそれが正しい。

実用上の対処法：

その通りだ。キャビティ流れの検証は、角部の局所的な振る舞いではなく、全体的な流れ場の再現性で判定する。中心線速度プロファイル、一次渦の位置と強度、二次渦の存在——これらが参照データと一致していれば、ソルバーの実装は信頼に値する。