「蓋を動かすだけで、箱の中に宇宙が広がる」

そんな言葉を聞いたら、ちょっと興味が湧きませんか？実はこれ、流体力学の世界では**「蓋駆動キャビティ流れ」**と呼ばれる、CFD（数値流体力学）の“Hello World”のような問題なんです。

「え、流体ってそんなに複雑なの？」
そう、あなたが思っているよりずっと奥が深い。しかも、ブラウザ上でリアルタイムにその挙動を観察できるツールがあるんです。

▶ 2D Navier-Stokes 流れシミュレーター

今回は、このシミュレーターを題材に、**「渦度-ストリーム関数法」**の本質から、実際の実装コード、そして数値実験までを一気に駆け抜けます。

ざっくり本質——「渦」で考えると、なぜ楽なのか？

流体の運動を支配するのは、Navier-Stokes方程式という、人類がまだ完全には解けていない偏微分方程式です。でも、2次元に限れば、ある「裏技」が使えます。

それが渦度-ストリーム関数法。

「速度を直接追いかけるのではなく、『渦（渦度）』と『流れの道筋（ストリーム関数）』を計算する」

これの何が嬉しいかというと：

圧力を計算しなくて済む（計算量が激減）
質量保存則（連続の式）を自動的に満たす（数値的に安定）
渦の生成・消滅・移動が直感的に理解できる

つまり、**「流れの“模様”」**を描くことに特化した手法なんです。

数式で理解する——3つの式が全てを支配する

この手法のコアは、たった3つの式に集約されます。

1. 渦度輸送方程式（渦の生成・拡散・移流）

\frac{\partial \omega}{\partial t} + u\frac{\partial \omega}{\partial x} + v\frac{\partial \omega}{\partial y} = \nu \nabla^2 \omega

日本語訳：
「ある場所の渦度 $\omega$ は、時間とともに変化する。その変化は、

第1項〜第2項（移流項）：流体が流れに乗って渦を運ぶ効果
右辺（拡散項）：粘性によって渦が広がっていく効果
の合計で決まる」

ここで $\nu$ は動粘性係数。レイノルズ数 $Re$ は、この移流と拡散のバランスを決めるパラメータです。

2. ポアソン方程式（渦度からストリーム関数を求める）

\nabla^2 \psi = -\omega

日本語訳：
「渦度 $\omega$ がわかれば、この式を解くことでストリーム関数 $\psi$ が求まる。ストリーム関数の等高線が、そのまま『流線』になる」

3. 速度の定義（ストリーム関数から速度を復元）

u = \frac{\partial \psi}{\partial y}, \quad v = -\frac{\partial \psi}{\partial x}

日本語訳：
「ストリーム関数 $\psi$ を $y$ 方向に微分すると $x$ 方向速度 $u$ が、$x$ 方向に微分すると $y$ 方向速度 $v$ が得られる」

つまり、この3式を順番に解けば、圧力を一切計算せずに流れ場が完全に記述できる——これが渦度-ストリーム関数法の神髄です。

コードで実装する——JavaScriptで動く2D Navier-Stokesソルバー

ここからが本番。上記の理論を、実際に動くコードに落とし込みます。

以下のコードは、シミュレーターのコアロジックを30行に凝縮したものです。差分法（中心差分＋前方差分）で離散化し、SOR法でポアソン方程式を解いています。

// 2D Navier-Stokesソルバー（渦度-ストリーム関数法）の核心部分
// グリッドサイズ N x N, 時間刻み dt, 動粘性係数 nu

// 1. 渦度輸送方程式（移流項＋拡散項）の更新
function updateVorticity(w, psi, N, dt, nu, dx, dy) {
    const wNew = new Array(N*N).fill(0);
    for (let j = 1; j < N-1; j++) {
        for (let i = 1; i < N-1; i++) {
            const idx = i + j * N;
            // 速度をストリーム関数から計算
            const u =  (psi[i + (j+1)*N] - psi[i + (j-1)*N]) / (2*dy);
            const v = -(psi[(i+1) + j*N] - psi[(i-1) + j*N]) / (2*dx);
            
            // 移流項（風上差分で安定化）
            const adv_x = u > 0 ? 
                u * (w[idx] - w[(i-1)+j*N]) / dx :
                u * (w[(i+1)+j*N] - w[idx]) / dx;
            const adv_y = v > 0 ?
                v * (w[idx] - w[i+(j-1)*N]) / dy :
                v * (w[i+(j+1)*N] - w[idx]) / dy;
            
            // 拡散項（中心差分）
            const diff = nu * (
                (w[(i+1)+j*N] - 2*w[idx] + w[(i-1)+j*N]) / (dx*dx) +
                (w[i+(j+1)*N] - 2*w[idx] + w[i+(j-1)*N]) / (dy*dy)
            );
            
            // 時間発展（オイラー陽解法）
            wNew[idx] = w[idx] + dt * (-adv_x - adv_y + diff);
        }
    }
    return wNew;
}

// 2. ポアソン方程式 ∇²ψ = -ω をSOR法で解く
function solveStreamFunction(omega, psi, N, dx, dy, maxIter=1000, tol=1e-6) {
    const omegaFactor = 1.8; // SOR加速係数
    for (let iter = 0; iter < maxIter; iter++) {
        let maxDiff = 0;
        for (let j = 1; j < N-1; j++) {
            for (let i = 1; i < N-1; i++) {
                const idx = i + j * N;
                // 5点差分の離散化
                const psiOld = psi[idx];
                psi[idx] = (1 - omegaFactor) * psiOld + omegaFactor * (
                    (psi[(i+1)+j*N] + psi[(i-1)+j*N]) / (dx*dx) +
                    (psi[i+(j+1)*N] + psi[i+(j-1)*N]) / (dy*dy) +
                    omega[idx]
                ) / (2/(dx*dx) + 2/(dy*dy));
                maxDiff = Math.max(maxDiff, Math.abs(psi[idx] - psiOld));
            }
        }
        if (maxDiff < tol) break; // 収束判定
    }
    return psi;
}

// 3. 境界条件の設定（蓋駆動キャビティ）
function applyBC(w, psi, N, dx, dy, U_lid) {
    // 壁面の渦度（トンプソン条件）
    for (let i = 0; i < N; i++) {
        // 上壁（蓋）：速度U_lidで移動
        w[i + (N-1)*N] = -2 * (psi[i + (N-1)*N] - psi[i + (N-2)*N]) / (dy*dy) 
                         + 2 * U_lid / dy;
        // 下壁：固定壁（速度0）
        w[i + 0*N] = -2 * (psi[i + 1*N] - psi[i + 0*N]) / (dy*dy);
    }
    for (let j = 0; j < N; j++) {
        // 左壁：固定壁
        w[0 + j*N] = -2 * (psi[1 + j*N] - psi[0 + j*N]) / (dx*dx);
        // 右壁：固定壁
        w[(N-1) + j*N] = -2 * (psi[(N-2) + j*N] - psi[(N-1) + j*N]) / (dx*dx);
    }
    // ストリーム関数の境界条件（壁面で0）
    for (let i = 0; i < N; i++) {
        psi[i + 0*N] = 0;      // 下壁
        psi[i + (N-1)*N] = 0;  // 上壁
    }
    for (let j = 0; j < N; j++) {
        psi[0 + j*N] = 0;      // 左壁
        psi[(N-1) + j*N] = 0;  // 右壁
    }
}

このコードのポイント：

渦度の移流項には風上差分を使用（数値振動を抑制）
ポアソン方程式はSOR法で高速収束（加速係数1.8が経験的に最適）
壁面の渦度はトンプソン条件で計算（速度境界条件から渦度を算出）

数値例で確かめる——Re=100の流れ場を実際に計算

理論とコードだけではピンと来ない。実際の数値で確認しましょう。

条件設定：

グリッドサイズ：N = 64（65×65点）
キャビティサイズ：1.0 × 1.0（無次元化）
蓋の速度：U_lid = 1.0
レイノルズ数：Re = 100（動粘性係数 ν = 0.01）
時間刻み：dt = 0.001

定常状態に達した後の流れ場の特徴：

中心渦の位置：キャビティのほぼ中央（x=0.5, y=0.6付近）
最大ストリーム関数値：ψ_max ≈ 0.1（Ghiaらのベンチマークデータと一致）
速度分布：蓋付近でU_lid=1.0から急激に減速

具体的な数値（中心断面のx方向速度）：

y=0.9（蓋近く）：u ≈ 0.95（ほぼ蓋と同じ速度）
y=0.5（中央）：u ≈ -0.2（逆流！）
y=0.1（底面近く）：u ≈ 0.05（ほぼ静止）

「中央で逆流が発生する」——これが蓋駆動キャビティ流れの特徴です。蓋が流体を右に引きずることで、中央付近では循環流が生まれ、左向きの流れが発生します。

シミュレーターで遊ぶ——パラメータを変えて現象を体感

さあ、いよいよ実際のシミュレーターで実験です。

▶ 2D Navier-Stokes 流れシミュレーター

実験1：Re=100 → 層流の教科書通りの渦

設定：Re=100, N=64, Δt=0.001
観察結果：中央に1つの大きな渦が鎮座。渦度コンターは同心円状にきれいに広がる。
なぜ？：粘性が強いため、渦の拡散が移流に勝る。流れは完全に層流で安定。

実験2：Re=1000 → 二次渦の出現

設定：Re=1000, N=128（解像度を上げる！）
観察結果：
- 中心渦が左上に移動（x≈0.35, y≈0.75）
- 左下の角に小さな二次渦が出現（反時計回り）
- 右下の角にも小さな二次渦が出現（時計回り）
なぜ？：慣性力が強くなると、主流渦が角部に二次的な循環を誘起する。これがMoffattの渦と呼ばれる現象。

実験3：Re=5000 → 三次渦の発生と非定常性

設定：Re=5000, N=256, Δt=0.0005（CFL条件を厳守）
観察結果：
- 中心渦はさらに左上に移動
- 左下の二次渦の中に、さらに小さな三次渦が発生！
- 流れ場が時間的に振動し始める（非定常性の兆候）
なぜ？：高Reでは粘性による安定化が追いつかず、渦の生成・消滅が活発になる。ただし、これは2次元の限界——実際の3次元では乱流に遷移する。

現場でハマるポイント——数値流体力学の落とし穴

1. 「メッシュ依存性」を軽視すると大事故

「N=32でRe=1000を計算したら、二次渦が見えなかった」

これはよくある話。グリッドが粗すぎると、数値拡散が物理的な粘性を上回り、本来あるべき渦構造が消えてしまいます。

解決策：

必ずメッシュ依存性調査を行う（N=32, 64, 128, 256と変えて結果が変わらないことを確認）
経験則：1つの渦を解像するのに最低10グリッド必要

2. 「高Re＝現実的」は大間違い

「Re=10000で計算した結果を、そのまま実機設計に使おう」

これは危険。このシミュレーターは2次元・層流を仮定しています。現実の高Re流れは3次元乱流で、渦構造が根本的に異なります。

正しい使い方：

Re≤1000：定量的に信頼できる
Re=1000〜5000：定性的な傾向把握に使う
Re>5000：**「乱流の入り口」**としての参考値

3. 「圧力が見えない」ことの意味

渦度-ストリーム関数法は圧力を計算しません。そのため：

翼の揚力や配管の圧力損失は計算できない
圧力が必要な場合はプリミティブ変数法（SIMPLE法など）を使う

覚えておくべきこと：
この手法は「流れのパターン」を見るためのもの。圧力が必要な場面では、別のソルバーが必要です。

まとめ——「見えない流れ」を「見える化」する力

今回の記事でお伝えしたかったこと、それは：

渦度-ストリーム関数法は、2次元流体の“最適解”
圧力を計算せず、流れの構造だけを効率的に捉える——教育用から研究用まで幅広く使われる理由です。
レイノルズ数は「流れの性格」を決める
Re=100の穏やかな層流から、Re=5000の複雑な渦構造まで、たった1つのパラメータで世界が変わります。
数値計算の“常識”を身につける
メッシュ依存性、CFL条件、数値拡散——これらの概念を知っているだけで、シミュレーション結果の信頼性が大きく変わります。

そして何より、**「実際に動かして体感する」**ことの大切さ。

▶ 2D Navier-Stokes 流れシミュレーター
ブラウザで即動作、登録不要。Reスライダーを動かすだけで、流体の世界が目の前に広がります。

「見えない流れ」を「見える化」する——それが数値流体力学の力であり、このシミュレーターがあなたに提供する体験です。

NovaSolverでは700以上の工学シミュレーターを無料公開中 👉 一覧はこちら