Visual Basic 6.0 というと、業務アプリのフォームにボタンを並べる言語、という印象があるかもしれません。そこで私は、あえてこの古い言語で流体シミュレーションを書いてみました。対称翼型 NACA 0012 を大きな迎角（デフォルトで 25°）で流れの中に置き、翼のまわりに発生する剥離と渦を、渦度の大きさで色付けして表示するプログラムです。CFD パッケージのスクリーンショットでよく見る、あの Vorticity Magnitude の図をイメージしています。
描画も当時のやり方そのままです。OpenGL も GPU も使わず、ピクセルの配列を StretchDIBits で PictureBox に転送するだけ。それでも、前縁での剥離、せん断層の巻き上がり、翼の後ろに伸びる渦列まで、ちゃんと現れてくれます。1998 年当時の私の PC なら泣いていたと思います。
ソースコードは公開しています: https://github.com/makarov-mm/vb6-turbulence
何を解いているのか
解いているのは 2 次元の非圧縮 Navier–Stokes 方程式です。手法は J. Stam の "Stable Fluids"（1999 年）として知られる方法で、リアルタイム流体シミュレーションの定番です。1 ステップは次の 4 段階からなります。
半ラグランジュ移流
各格子点について、速度場に沿って軌跡を時間の逆向きにたどり、たどり着いた場所の値を双線形補間で拾ってきます。この方法の利点は無条件安定であることです。時間刻みを大きくしてもシミュレーションが発散しません。VB6 のように計算が速くない環境では、これは大きな助けになります。
粘性拡散
5 点ラプラシアンによる陽的な拡散です。ここで動粘性係数が効いてきます。UI のスライダーで粘性を下げるほど、流れは「乱流らしく」なっていきます。
渦度閉じ込め（Vorticity Confinement）
半ラグランジュ移流は安定な代わりに数値拡散が大きく、放っておくと渦がどんどんなまって、後流がのっぺりした絵になってしまいます。そこで Steinhoff / Fedkiw の渦度閉じ込めを入れています。速度場から渦度を計算し、渦の中心に向かう力を少しだけ加え戻すことで、数値拡散で失われた渦のエネルギーを補います。係数を 0 にすれば無効化でき、0.1 以上にすると後流が激しく「沸き立つ」ようになります。
圧力投影
非圧縮性を保つための段階です。ポアソン方程式 lap(p) = div(V) をガウス・ザイデル法（デフォルトで 24 反復）で解き、得られた圧力の勾配を速度から差し引きます。これで速度場はほぼ発散なしの状態に戻ります。反復回数を増やすほど非圧縮性は良くなりますが、そのぶん遅くなります。
翼と境界条件
翼はビットマップではなく、NACA 00xx 系列の翼厚分布の式（厚み 12%）から解析的に生成した固体セルのマスクとして表現しています。迎角のスライダーを動かすと、マスクはその場で作り直されます。固体セルでは流体の侵入を禁止する条件を課しています。
計算領域の境界は次のとおりです。左端は流入で、わずかなランダムな揺らぎを加えています。この揺らぎがせん断層の不安定性に「火をつけ」、渦の巻き上がりを引き起こします。右端は自由流出、上下の壁は自由すべり条件です。
RANS / URANS / DDES との関係
CFD パッケージの比較図では、同じ翼まわりの流れでも手法によって絵がまったく違って見えます。
- RANS は定常的な平均化で、乱流はすべてモデルの中に「隠れて」います。見えるのは滑らかな後流だけです
- URANS は非定常計算で、剥離やカルマン渦列といった大きな渦は見えますが、小さなスケールは依然としてモデル化されています
- DDES / LES では大きな渦を格子で直接解像し、モデル化するのは格子以下のスケールだけです。最も情報量の多い絵になります。
このソルバーは非定常で、明示的な乱流モデルを持ちません。大きな渦は格子で直接解像され、格子以下のスケールの役割はスキームの数値散逸が担います（いわゆる ILES 的なアプローチです）。そのため得られる絵は URANS / DDES のパネルに最も近く、前縁剥離、せん断層の渦への巻き上がり、翼後方の非定常な渦列が観察できます。もちろん、これは 220×110 格子・低レイノルズ数・2 次元の「おもちゃのモデル」であって、工学的な CFD ではありません。
VB6 ならではの話
ソースはプレーンな ASCII・CRLF 改行なので、どのロケールの Windows でも文字化けせずに開けます。VB6 の IDE で Turbulence.vbp を開いて F5 を押せば動きますが、IDE 上では p-code 実行になるため計算はかなり遅めです。本気で動かすなら File → Make Turbulence.exe でネイティブ EXE にビルドしてください。プロジェクト設定の Compile タブで最適化を有効にし、配列境界チェックとオーバーフローチェックを切ると、ソルバーは目に見えて速くなります。
操作できるパラメータ
- Start / Pause / Reset — 実行、停止、流れ場のクリア
- 迎角（0–40°）— 翼のマスクをその場で再構築
- 粘性 — 下げるほど乱流的に
- 流速 — 一様流の速さ
- トレーサー粒子 — 速度場に流される緑色の受動粒子
コード側では、frmMain.frm 冒頭の定数で格子サイズ（GRID_NX, GRID_NY）、時間スケール、粒子数を、Form_Load で圧力ソルバーの反復回数と渦度閉じ込めの強さを調整できます。
制約について
正直に書いておくと、手元に実物の VB6 IDE がなかったため、このコードを実際の VB6 でコンパイルして確認することはできていません。コードは厳密に VB6 の文法で書かれており、数値スキームは同一の実装の Python 移植版で一対一に検証しています（安定性、剥離の発生、迎角 25° での渦列を確認済みです）。また、繰り返しになりますが、これは学習とデモのための 2D トイモデルであって、検証済みの CFD コードではありません。
まとめ
Navier–Stokes 方程式を解くのに、最新の言語も GPU も必須ではありません。半ラグランジュ移流、圧力投影、渦度閉じ込め。この 3 つの部品を組み合わせれば、四半世紀前の言語でも、翼のまわりで渦が生まれて流れていく様子を目の前で見ることができます。ソースコードは MIT ライセンスで公開していますので、パラメータをいじって遊んでみてください。