円管内の定常層流はハーゲン・ポアズイユ流れと呼ばれ、NS方程式の厳密解が得られる重要な基本問題です。速度は放物線分布をとり、流量は管径の4乗に比例します。
本記事の内容
- 放物線速度分布の導出
- 体積流量と平均流速
- ダルシー・ワイスバッハの式
- Pythonでの可視化
前提知識
速度分布の導出
仮定と支配方程式
半径 $R$ の円管内の非圧縮性定常流で、軸方向（$z$方向）のみに流速があると仮定します：$\bm{v} = (0, 0, u_z(r))$。
NS方程式の $z$ 成分は：
$$ 0 = -\frac{\partial p}{\partial z} + \mu\left(\frac{\partial^2 u_z}{\partial r^2} + \frac{1}{r}\frac{\partial u_z}{\partial r}\right) $$
圧力勾配は一定とします：$-\partial p / \partial z = \Delta p / L$。
$$ \frac{1}{r}\frac{d}{dr}\left(r\frac{du_z}{dr}\right) = -\frac{\Delta p}{\mu L} $$
積分
1回積分：
$$ r\frac{du_z}{dr} = -\frac{\Delta p}{2\mu L}r^2 + C_1 $$
$r = 0$ で速度勾配が有限（$C_1 = 0$）より：
$$ \frac{du_z}{dr} = -\frac{\Delta p}{2\mu L}r $$
2回積分：
$$ u_z = -\frac{\Delta p}{4\mu L}r^2 + C_2 $$
壁面 $r = R$ で $u_z = 0$（no-slip条件）：
$$ C_2 = \frac{\Delta p R^2}{4\mu L} $$
$$ \boxed{u_z(r) = \frac{\Delta p}{4\mu L}(R^2 – r^2)} $$
中心速度（$r = 0$）：$u_{\max} = \Delta p R^2 / (4\mu L)$
体積流量
$$ Q = \int_0^R u_z(r) \cdot 2\pi r \, dr = \frac{2\pi\Delta p}{4\mu L}\int_0^R (R^2 – r^2)r \, dr $$
$$ = \frac{\pi\Delta p}{2\mu L}\left[\frac{R^2 r^2}{2} – \frac{r^4}{4}\right]_0^R = \frac{\pi\Delta p R^4}{8\mu L} $$
$$ \boxed{Q = \frac{\pi \Delta p R^4}{8\mu L} = \frac{\pi \Delta p d^4}{128\mu L}} $$
流量は管径の4乗に比例します。平均流速は $\bar{u} = Q/(\pi R^2) = u_{\max}/2$ です。
ダルシー・ワイスバッハの式
管摩擦による圧力損失の一般式：
$$ \Delta p = f \cdot \frac{L}{d} \cdot \frac{\rho \bar{u}^2}{2} $$
層流（$Re < 2300$）での摩擦係数：
$$ \boxed{f = \frac{64}{Re_D}} $$
Pythonでの可視化
import numpy as np
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (1) 放物線速度分布
R = 0.025 # 管半径 [m]
r = np.linspace(-R, R, 200)
dp_L = 1000 # 圧力勾配 [Pa/m]
mu = 1e-3 # 水の粘性 [Pa·s]
u = dp_L / (4*mu) * (R**2 - r**2)
axes[0].plot(u, r*1e3, 'b-', lw=2.5)
axes[0].fill_betweenx(r*1e3, 0, u, alpha=0.2, color='blue')
axes[0].axhline(0, color='r', ls='--', alpha=0.5)
axes[0].set_xlabel('Velocity $u_z$ [m/s]', fontsize=12)
axes[0].set_ylabel('Radial position $r$ [mm]', fontsize=12)
axes[0].set_title('Hagen-Poiseuille velocity profile', fontsize=13)
axes[0].grid(True, alpha=0.3)
# (2) 流量 vs 管径（4乗則）
d_range = np.linspace(0.01, 0.10, 100)
dp = 10000; L = 10.0
Q_range = np.pi * dp * d_range**4 / (128 * mu * L)
axes[1].plot(d_range*100, Q_range*1e3, 'b-', lw=2.5)
axes[1].set_xlabel('Pipe diameter $d$ [cm]', fontsize=12)
axes[1].set_ylabel('Volume flow rate $Q$ [L/s]', fontsize=12)
axes[1].set_title('$Q \\propto d^4$ (Hagen-Poiseuille)', fontsize=13)
axes[1].grid(True, alpha=0.3)
# (3) 摩擦係数 f vs Re
Re = np.logspace(1.5, 4.5, 300)
f_lam = 64 / Re
f_turb = 0.316 * Re**(-0.25)
axes[2].loglog(Re[Re<2300], f_lam[Re<2300], 'b-', lw=2.5, label='Laminar: $f=64/Re$')
axes[2].loglog(Re[Re>4000], f_turb[Re>4000], 'r-', lw=2.5, label='Turbulent (Blasius)')
axes[2].axvline(2300, color='gray', ls='--', alpha=0.5)
axes[2].set_xlabel('$Re_D$', fontsize=12); axes[2].set_ylabel('$f$', fontsize=12)
axes[2].set_title('Darcy friction factor', fontsize=13)
axes[2].legend(); axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
まとめ
- 放物線速度分布: $u_z(r) = \frac{\Delta p}{4\mu L}(R^2 – r^2)$
- 体積流量: $Q = \pi\Delta p d^4 / (128\mu L)$（管径の4乗に比例）
- 層流摩擦係数: $f = 64/Re$
- 平均流速は最大流速の半分
次のステップとして、以下の記事も参考にしてください。