「蓋を動かすだけで、箱の中に宇宙が広がる」
そんな言葉を聞いたら、ちょっと興味が湧きませんか？実はこれ、流体力学の世界では**「蓋駆動キャビティ流れ」**と呼ばれる、CFD（数値流体力学）の“Hello World”のような問題なんです。
「え、流体ってそんなに複雑なの？」
そう、あなたが思っているよりずっと奥が深い。しかも、ブラウザ上でリアルタイムにその挙動を観察できるツールがあるんです。
今回は、このシミュレーターを題材に、**「渦度-ストリーム関数法」**の本質から、実際の実装コード、そして数値実験までを一気に駆け抜けます。
ざっくり本質——「渦」で考えると、なぜ楽なのか？
流体の運動を支配するのは、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（ほぼ静止）
「中央で逆流が発生する」——これが蓋駆動キャビティ流れの特徴です。蓋が流体を右に引きずることで、中央付近では循環流が生まれ、左向きの流れが発生します。
シミュレーターで遊ぶ——パラメータを変えて現象を体感
さあ、いよいよ実際のシミュレーターで実験です。
実験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以上の工学シミュレーターを無料公開中 👉 一覧はこちら