ナビエ・ストークス方程式（Navier-Stokes equations, NS方程式）は、粘性流体の運動を支配する基本方程式です。連続の方程式と組み合わせることで、ほぼ全ての流体現象を記述できます。
本記事の内容
- NS方程式の導出（微小流体要素への力の釣り合い）
- 各項の物理的意味
- 非圧縮性NS方程式
- 無次元化とレイノルズ数
- Pythonでの可視化
前提知識
導出
ニュートンの第2法則
微小流体要素に対して $\bm{F} = m\bm{a}$ を適用します。流体粒子の加速度は物質微分で表されます：
$$ \frac{D\bm{v}}{Dt} = \frac{\partial \bm{v}}{\partial t} + (\bm{v} \cdot \nabla)\bm{v} $$
作用する力
- 圧力: $-\nabla p$（圧力勾配に逆らう方向）
- 粘性力: $\mu \nabla^2 \bm{v}$（ニュートン流体、非圧縮性の場合）
- 外力: $\rho \bm{g}$（重力など）
NS方程式
$$ \boxed{\rho\left(\frac{\partial \bm{v}}{\partial t} + (\bm{v} \cdot \nabla)\bm{v}\right) = -\nabla p + \mu \nabla^2 \bm{v} + \rho \bm{g}} $$
左辺は慣性項（流体粒子の加速度 × 密度）、右辺は圧力勾配力 + 粘性力 + 体積力です。
成分表示（直交座標、非圧縮性）
$x$ 成分：
$$ \rho\left(\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z}\right) = -\frac{\partial p}{\partial x} + \mu\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}\right) + \rho g_x $$
$y$, $z$ 成分も同様の構造です。
無次元化とレイノルズ数
特性長さ $L$、特性速度 $U$ で無次元化すると：
$$ \frac{\partial \bm{v}^*}{\partial t^*} + (\bm{v}^* \cdot \nabla^*)\bm{v}^* = -\nabla^* p^* + \frac{1}{Re}\nabla^{*2}\bm{v}^* $$
ここでレイノルズ数 $Re = \rho UL/\mu$ が唯一の無次元パラメータです。$Re \gg 1$ で慣性力支配（乱流）、$Re \ll 1$ で粘性力支配（ストークス流れ）となります。
特殊ケース
| 条件 | 結果 |
|---|---|
| 非粘性（$\mu = 0$） | オイラー方程式 |
| 定常・非粘性・非回転 | ベルヌーイの定理 |
| $Re \ll 1$ | ストークス方程式 $\nabla p = \mu\nabla^2\bm{v}$ |
Pythonでの可視化
import numpy as np
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# (1) ストークス流れ（球まわりの低Re流れ）
theta = np.linspace(0, 2*np.pi, 100)
r_vals = np.linspace(1.0, 4.0, 20)
R_sphere = 1.0
ax = axes[0]
for r in r_vals:
# ストークスの流れ関数から流線
x_line = r * np.cos(theta)
y_line = r * np.sin(theta)
# 速度場（ストークス解）
vr = (1 - 3*R_sphere/(2*r) + R_sphere**3/(2*r**3)) * np.cos(theta)
vt = -(1 - 3*R_sphere/(4*r) - R_sphere**3/(4*r**3)) * np.sin(theta)
xg = np.linspace(-4, 4, 40)
yg = np.linspace(-4, 4, 40)
X, Y = np.meshgrid(xg, yg)
R = np.sqrt(X**2 + Y**2)
R = np.where(R < R_sphere, np.nan, R)
THETA = np.arctan2(Y, X)
Vr = (1 - 3*R_sphere/(2*R) + R_sphere**3/(2*R**3)) * np.cos(THETA)
Vt = -(1 - 3*R_sphere/(4*R) - R_sphere**3/(4*R**3)) * np.sin(THETA)
U = Vr * np.cos(THETA) - Vt * np.sin(THETA)
V = Vr * np.sin(THETA) + Vt * np.cos(THETA)
ax.streamplot(X, Y, U, V, density=1.5, color='blue', linewidth=0.8)
circle = plt.Circle((0, 0), R_sphere, color='gray', alpha=0.5)
ax.add_patch(circle)
ax.set_xlim(-4, 4); ax.set_ylim(-4, 4); ax.set_aspect('equal')
ax.set_title('Stokes flow around sphere ($Re \\ll 1$)', fontsize=13)
ax.set_xlabel('$x$'); ax.set_ylabel('$y$'); ax.grid(True, alpha=0.3)
# (2) NS方程式の各項のスケール比較
Re_range = np.logspace(-2, 6, 200)
inertia = np.ones_like(Re_range) # 慣性項（基準=1）
viscous = 1.0 / Re_range # 粘性項
pressure = np.ones_like(Re_range) # 圧力項
axes[1].loglog(Re_range, inertia, 'b-', lw=2, label='Inertia')
axes[1].loglog(Re_range, viscous, 'r-', lw=2, label='Viscous ~ 1/Re')
axes[1].axvline(1, color='gray', ls='--', alpha=0.5)
axes[1].set_xlabel('Reynolds number $Re$', fontsize=12)
axes[1].set_ylabel('Relative magnitude', fontsize=12)
axes[1].set_title('Balance of NS terms', fontsize=13)
axes[1].legend(fontsize=11); axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
まとめ
- NS方程式: $\rho D\bm{v}/Dt = -\nabla p + \mu\nabla^2\bm{v} + \rho\bm{g}$
- 物質微分 $D/Dt$ が流体粒子の加速度を表す
- レイノルズ数 $Re$ が慣性力と粘性力の比を決定
- $Re \ll 1$ でストークス流れ、$Re \gg 1$ で慣性支配
次のステップとして、以下の記事も参考にしてください。