連続の方程式は、流体力学で最も基本的な支配方程式の一つです。「質量は生成も消滅もしない」という質量保存則を数学的に表現したものであり、ナビエ-ストークス方程式と並んで流体解析の基盤を成します。
本記事では、検査体積を用いた積分形の導出から微分形への変換、非圧縮性条件の意味、そして断面積変化と流速の関係をPythonで可視化します。
本記事の内容
- 検査体積と質量保存
- 連続の方程式の積分形
- 微分形への導出
- 非圧縮性流体の条件
- 断面積変化と流速の関係
- Pythonでの可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
検査体積
系と検査体積
流体力学では、保存則を適用するために2つの考え方を用います。
系（system）: 同じ流体粒子の集合を追跡する考え方（ラグランジュ的）。
検査体積（control volume, CV）: 空間に固定された領域を設定し、その領域を出入りする流体を考える方法（オイラー的）。
検査体積は空間に固定されているため、境界面（検査面, control surface, CS）を通じて流体が出入りします。
連続の方程式の導出
積分形
空間に固定された検査体積 $V$ を考えます。検査面 $S$ の外向き法線ベクトルを $\bm{n}$ とします。
検査体積内の全質量は、
$$ M = \int_V \rho \, dV $$
質量保存則は、「検査体積内の質量の時間変化率 = 検査面を通じた質量流入率」と表現されます。
$$ \frac{\partial}{\partial t} \int_V \rho \, dV + \oint_S \rho \bm{v} \cdot \bm{n} \, dS = 0 $$
ここで $\rho \bm{v} \cdot \bm{n} \, dS$ は、検査面の微小面 $dS$ を通じて単位時間に流出する質量です。$\bm{v} \cdot \bm{n} > 0$ のとき流出、$\bm{v} \cdot \bm{n} < 0$ のとき流入を表します。
微分形への変換
ガウスの発散定理を用いて面積分を体積分に変換します。
$$ \oint_S \rho \bm{v} \cdot \bm{n} \, dS = \int_V \nabla \cdot (\rho \bm{v}) \, dV $$
これを積分形に代入すると、
$$ \int_V \left[ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bm{v}) \right] dV = 0 $$
この等式が任意の検査体積 $V$ に対して成り立つためには、被積分関数がゼロでなければなりません。
$$ \boxed{\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bm{v}) = 0} $$
これが 連続の方程式の微分形 です。
成分表示
直交座標系 $(x, y, z)$ で速度を $\bm{v} = (u, v, w)$ とすると、
$$ \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} + \frac{\partial (\rho w)}{\partial z} = 0 $$
展開形
$\nabla \cdot (\rho \bm{v}) = \rho \nabla \cdot \bm{v} + \bm{v} \cdot \nabla \rho$ を用いて展開すると、
$$ \frac{\partial \rho}{\partial t} + \bm{v} \cdot \nabla \rho + \rho \nabla \cdot \bm{v} = 0 $$
左辺の最初の2項は物質微分 $D\rho/Dt$ に他ならないため、
$$ \frac{D\rho}{Dt} + \rho \nabla \cdot \bm{v} = 0 $$
これは「流体粒子に乗って観測した密度の変化率は、速度の発散に比例する」ことを表しています。
非圧縮性流体の条件
非圧縮性の定義
非圧縮性流体 とは、流体粒子を追跡したときに密度が変化しない流体です。
$$ \frac{D\rho}{Dt} = 0 $$
この条件を連続の方程式に代入すると、
$$ \boxed{\nabla \cdot \bm{v} = 0} $$
すなわち、非圧縮性流体では 速度場の発散がゼロ です。
成分表示では、
$$ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0 $$
この条件は「流体要素の体積が変化しない」ことを意味します。水などの液体の流れでは多くの場合この近似が適用可能です。
1次元定常流の連続の方程式
断面積変化と流速
1次元の定常流（$\partial / \partial t = 0$）において、流管の断面積 $A(x)$ が変化する場合を考えます。
断面 1（面積 $A_1$, 速度 $v_1$, 密度 $\rho_1$）と断面 2（面積 $A_2$, 速度 $v_2$, 密度 $\rho_2$）の間で質量保存を適用すると、
$$ \rho_1 A_1 v_1 = \rho_2 A_2 v_2 $$
非圧縮性（$\rho_1 = \rho_2$）の場合、
$$ \boxed{A_1 v_1 = A_2 v_2} $$
これは 体積流量の保存 を表しています。断面積が半分になると流速は2倍になります。
体積流量と質量流量
体積流量 $Q$ と質量流量 $\dot{m}$ は次のように定義されます。
$$ Q = A v \quad [\mathrm{m^3/s}] $$
$$ \dot{m} = \rho A v \quad [\mathrm{kg/s}] $$
非圧縮性定常流では $Q = \text{const.}$、一般の定常流では $\dot{m} = \text{const.}$ です。
Pythonでの可視化
import numpy as np
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# --- (1) 断面積変化と流速の関係 ---
x = np.linspace(0, 10, 500)
# 絞り管: 断面積が変化する
A = 1.0 - 0.6 * np.exp(-((x - 5)**2) / 2) # 最小で A=0.4
Q = 1.0 # 体積流量 [m^3/s]
v = Q / A # 流速
ax1_top = axes[0]
ax1_bot = ax1_top.twinx()
ax1_top.fill_between(x, 0, A, alpha=0.3, color='blue', label='Area $A(x)$')
ax1_top.plot(x, A, 'b-', lw=2)
ax1_bot.plot(x, v, 'r-', lw=2, label='Velocity $v(x)$')
ax1_top.set_xlabel('Position $x$ [m]', fontsize=12)
ax1_top.set_ylabel('Cross-sectional area $A$ [m$^2$]', fontsize=12, color='b')
ax1_bot.set_ylabel('Flow velocity $v$ [m/s]', fontsize=12, color='r')
ax1_top.set_title('Continuity: $Av$ = const.', fontsize=13)
lines1, labels1 = ax1_top.get_legend_handles_labels()
lines2, labels2 = ax1_bot.get_legend_handles_labels()
ax1_top.legend(lines1 + lines2, labels1 + labels2, fontsize=10)
# --- (2) 2次元非圧縮性流れの発散チェック ---
xg = np.linspace(-3, 3, 30)
yg = np.linspace(-3, 3, 30)
Xg, Yg = np.meshgrid(xg, yg)
# 非圧縮性流れの例: 渦流れ
r2 = Xg**2 + Yg**2 + 1e-10
u_vortex = -Yg / r2
v_vortex = Xg / r2
# 発散を計算（数値微分）
dx = xg[1] - xg[0]
dy = yg[1] - yg[0]
du_dx = np.gradient(u_vortex, dx, axis=1)
dv_dy = np.gradient(v_vortex, dy, axis=0)
div_v = du_dx + dv_dy
speed = np.sqrt(u_vortex**2 + v_vortex**2)
u_norm = u_vortex / (speed + 1e-10)
v_norm = v_vortex / (speed + 1e-10)
axes[1].quiver(Xg, Yg, u_norm, v_norm, speed, cmap='viridis', alpha=0.8)
axes[1].set_xlabel('$x$', fontsize=12)
axes[1].set_ylabel('$y$', fontsize=12)
axes[1].set_title('Vortex flow ($\\nabla \\cdot \\mathbf{v} = 0$)', fontsize=13)
axes[1].set_aspect('equal')
axes[1].grid(True, alpha=0.3)
# --- (3) 流管の概念図 ---
# 2次元ノズル形状の流線
xn = np.linspace(0, 10, 400)
yn = np.linspace(-2, 2, 400)
Xn, Yn = np.meshgrid(xn, yn)
# ノズル流れ（収束）
A_func = 2.0 - 1.2 * np.exp(-((Xn - 5)**2) / 3)
u_nozzle = 1.0 / A_func
v_nozzle = Yn * 1.2 * (Xn - 5) / 3 * np.exp(-((Xn - 5)**2) / 3) / A_func**2
axes[2].streamplot(Xn, Yn, u_nozzle, v_nozzle, density=2.0, color='blue', linewidth=0.8)
# ノズル壁を描画
wall_top = 1.0 - 0.6 * np.exp(-((xn - 5)**2) / 3)
wall_bot = -wall_top
axes[2].fill_between(xn, wall_top, 2.5, color='gray', alpha=0.3)
axes[2].fill_between(xn, -2.5, wall_bot, color='gray', alpha=0.3)
axes[2].plot(xn, wall_top, 'k-', lw=2)
axes[2].plot(xn, wall_bot, 'k-', lw=2)
axes[2].set_xlabel('$x$', fontsize=12)
axes[2].set_ylabel('$y$', fontsize=12)
axes[2].set_title('Stream tube in a nozzle', fontsize=13)
axes[2].set_ylim(-2, 2)
axes[2].set_aspect('equal')
plt.tight_layout()
plt.show()
左のグラフでは、絞り管における断面積と流速の関係を示しています。連続の方程式 $Av = \text{const.}$ により、断面積が最小の部分で流速が最大になることが確認できます。
中央のグラフでは、渦流れの速度場を示しています。この流れは $\nabla \cdot \bm{v} = 0$ を満たす非圧縮性流れです。
右のグラフでは、ノズル内の流管を流線で可視化しています。ノズルが絞られる部分で流線が密になり、流速が増加する様子がわかります。
まとめ
本記事では、連続の方程式（質量保存則）を導出しました。
- 積分形: $\frac{\partial}{\partial t}\int_V \rho \, dV + \oint_S \rho \bm{v} \cdot \bm{n} \, dS = 0$
- 微分形: $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bm{v}) = 0$
- 非圧縮性条件: $\nabla \cdot \bm{v} = 0$ — 速度場の発散がゼロ
- 1次元定常流: $A_1 v_1 = A_2 v_2$ — 断面積が小さいと流速が大きい
次のステップとして、以下の記事も参考にしてください。