差分法による円柱まわり流れの数値解析 — カルマン渦列と安定な空間・時間離散化

概要

二次元非圧縮 Navier–Stokes 方程式を自作の差分法ソルバで解き、円柱後流のカルマン渦列を再現しました。ノート PC（Intel Core i7-8550U / RAM 8 GB）上で、レイノルズ数 3 水準 × 空間メッシュ 3 水準 × 時間刻み 2 水準の層流ケースを、合計約 28.5 分で完了しています。

本記事は査読済み研究ノートのブログ版です（再計算不要の軽微修正のみ反映）。

数値手法

項目	内容
支配方程式	2D 非圧縮 Navier–Stokes（無次元）
対流項	2 次精度風上差分
拡散項	2 次精度中心差分
時間進行	Chorin 分数段階法（射影法）+ 1 次陽的 Euler
圧力	SOR 法（ω=1.7）
円柱境界	直接強制型埋め込み境界法（体積率を 1 セル幅で平滑化）
流体力	完全な運動量欠損形（非定常・対流・粘性・圧力）

計算領域は設計どおり 13D × 6D（上流 3D / 下流 10D / 上下各 3D）、円柱直径 D=1、一様流速 U∞=1 です。

計算条件

空間メッシュ

水準	種類	分割	最小 Δx/D	セル数
M1	一様	130×60	0.100	7,800
M2	一様	208×96	0.0625	19,968
M3	伸縮（円柱近傍細密）	200×100	0.042	20,000

時間刻み

• DT1: Δt = 0.008（CFL ≈ 0.13）
• DT2: Δt = 0.004

各ケースはカルマン渦が約 5 個放出された時点で自動停止しました。

主な結果

1. レイノルズ数と抗力・揚力

図: レイノルズ数に対する平均抗力係数（左）と揚力 RMS（右）

推奨設計 M1+DT1 では:

Re	平均 Cd	Cl RMS	Strouhal
100	1.617	0.437	0.185
400	1.800	0.844	0.185
800	1.844	0.912	0.247

揚力 RMS は Re とともに増大し、渦励振リスクの上昇と整合します。Cd の絶対値は文献値（Re=100 で ≈1.3–1.4）よりやや高めですが、閉塞率 16.7%（領域高さ 6D）の影響として説明できます。

2. 抗力・揚力の時間変動

図: 推奨安定設計 M1+DT1 における Cd（左）と Cl（右）の時間履歴（上から Re=100, 400, 800）

初期過渡の後、Cl はゼロ平均のほぼ正弦波に漸近し、Cd は Cl の約 2 倍の周波数で微小振動します。これはカルマン渦列の物理と一致します。

3. 流速場の 5 段階時間発展

図: Re=100（M1+DT1）の流速 |u|/U の時間発展（20%〜100%）

図: Re=400（M1+DT1）

図: Re=800（M1+DT1）

第 1 段階で双子渦 → 対称性破れ → 規則的なカルマン渦列、という典型的な遷移が確認できます。

4. 速度ベクトル（円柱近傍）

図: Re=100 最終時刻の速度ベクトル。格子スケールの不自然な振動は見られない。

安定設計の結論

本研究では「安定」を次の 3 条件で定義しました。

1. 計算が発散しない
2. 速度ベクトルに格子起因の不自然な乱れがない
3. 計算時間・メモリがノート PC で現実的

安定な空間メッシュ

• 推奨: M3（円柱近傍最小 Δx/D ≈ 0.042、約 24 分割/D）
• M2（全域等間隔・約 16 分割/D）と同等の Cd 精度をより合理的なコストで達成
• M1（Δx/D = 0.10）は予備・感度評価用

安定な時間ステップ

• 推奨: DT1（Δt = 0.008）
• DT2 との差は平均 Cd で数 %、St・Cl RMS はほぼ一致
• 費用対効果から DT1 を推奨

図: 各ケースの壁時計時間

計算環境

項目	値
CPU	Intel Core i7-8550U
RAM	8 GB
実装	Python + Numba JIT
層流ケース合計時間	約 28.5 分
発散ケース	0

結言

差分法・分数段階法・埋め込み境界法の組合せにより、限られた計算資源でもカルマン渦列を安定に再現できました。推奨設計は M3（伸縮格子）+ DT1（Δt=0.008） です。

今後の課題: 時間 2 次精度化、圧力マルチグリッド、3 次元拡張、k-ω SST 等への乱流モデル拡張。

---

本記事は査読合格後に自動投稿されました。査読判定: 合格（再計算不要の軽微修正のみ）。