ポアズイユ流れ：解析解の導出と有限差分法による検証

平板が2枚、間隔は $2h$ です。圧力は $x$ 方向にだけ下がり、完全発達なら $v=0$ で $u$ は $y$ だけの関数です。円管も同じ名前（Poiseuille flow）で呼ばれますが、下の式は2次元チャンネルです。解がこの放物線から外れたら、粘性項か no-slip が壊れています。

−dp/dx と μ のスライダーを動かすと、画面下の $u_{\max}$ だけが変わり、曲線の形はそのままです。

描かれている曲線は $u/u_{\max} = 1 - (y/h)^2$ です。Re のスライダーは隅の数字だけを動かします。$h=1$ のとき、ピークの表示は −dp/dx の値を $2\mu$ で割ったものです。

問題の設定

形状および境界条件

2枚の無限平行平板間の完全発達した層流を考えます。

• チャンネル半幅：$h$（平板は $y = \pm h$ に位置）
• 流れ方向：$x$（無限に長いと仮定）
• 境界条件：両壁面で no-slip $\Rightarrow$ $u(\pm h) = 0$
• 駆動源：一定の圧力勾配 $\dfrac{dp}{dx} = -G < 0$（$G > 0$ は定数）

支配方程式

完全発達した流れなので $\partial/\partial x = 0, v = 0$ です。$x$ 方向のナビエ・ストークス方程式は以下のようになります。

$\rho \left( u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} \right) = -\frac{dp}{dx} + \mu \frac{\partial^2 u}{\partial y^2}$

完全発達の条件（$\partial u/\partial x = 0, v = 0$）を適用すると：

$\mu \frac{d^2 u}{dy^2} = \frac{dp}{dx} = -G$

解析解の導出

ステップ1：常微分方程式（ODE）への単純化#

$\frac{d^2 u}{dy^2} = -\frac{G}{\mu} \equiv -C \quad (C > 0)$

ステップ2：1回目の積分#

$\frac{du}{dy} = -C y + A_1$

ステップ3：2回目の積分#

$u(y) = -\frac{C}{2} y^2 + A_1 y + A_2$

ステップ4：境界条件の適用#

$u(+h) = 0$：

$u(-h) = 0$：

(1) $-$ (2) より：$0 = 2 A_1 h \Rightarrow A_1 = 0$（対称性）

(1) $+$ (2) より：$0 = -C h^2 + 2A_2 \Rightarrow A_2 = \dfrac{C h^2}{2}$

最終的な解析解

$\boxed{u(y) = \frac{G}{2\mu}\bigl(h^2 - y^2\bigr)}$

チャンネル中心（$y=0$）での最大速度：

$u_{\max} = \frac{G h^2}{2\mu}$

物理的解釈

速度プロファイル

解析解は**放物線（parabola）**分布となります。壁面の摩擦が流れの断面全体に均一に伝わり、中心で最大、壁面で0になります。

無次元速度 $u/u_{\max} = 1 - (y/h)^2$ は、チャンネルの半幅 $h$ や粘度 $\mu$ に関わらず一定となる普遍的なプロファイルです。

壁面せん断応力

$\tau_w = \mu \left.\frac{du}{dy}\right|_{y=h} = \mu \cdot \frac{G}{\mu}(-h) = -G h$

（大きさ：$\tau_w = G h$） 圧力勾配が大きいほど、またチャンネルが広いほど、壁面摩擦は大きくなります。

体積流量

$Q = \int_{-h}^{h} u \, dy = \frac{G}{2\mu} \int_{-h}^{h} (h^2 - y^2) \, dy = \frac{G}{2\mu} \cdot \frac{4h^3}{3} = \frac{2 G h^3}{3\mu}$

これは ハーゲン・ポアズイユの法則 の2D版と呼ばれます。流量は $h^3$ に比例するため、チャンネル幅が2倍になると流量は8倍に増加します。

コードの既定値 $h=G=\mu=1$ を同じ式へ入れると、$u_{\max}=1/2$、$Q=2/3$ です。高さ $2h$ で割ると平均速度は $1/3$ です。ピークと平均の比は $3/2$ なので、平均は中心の $2/3$ です。$1-(y/h)^2$ を $-h$ から $h$ まで積分し、$2h\,u_{\max}$ で割ると同じ $2/3$ が出ます。比のなかで $h$ と $\mu$ は消えます。

$h$ だけを2倍にすると次数が分かれます。$u_{\max}$ は $h^2$ なので4倍、$Q$ は $h^3$ なので8倍、壁面応力の大きさ $Gh$ は2倍です。流量が8倍になったから応力も同じ倍率だとすると、壁の値を4倍大きく読みます。断面を円管に移し、この放物線の運動エネルギーを平均速度で割ると係数が変わります。層流の円管でその係数が2.0000に固定される計算は 速度水頭の前に付く係数 にあります。

Pythonによる解析解の可視化#

import numpy as np
import matplotlib.pyplot as plt
 
# パラメータ
h   = 1.0    # チャンネル半幅 [m]
G   = 1.0    # 圧力勾配 [Pa/m]
mu  = 1.0    # 粘性係数 [Pa·s]
 
# 格子
y = np.linspace(-h, h, 400)
u_exact = G / (2 * mu) * (h**2 - y**2)
u_max   = G * h**2 / (2 * mu)
 
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
 
# 速度プロファイル (横方向)
axes[0].plot(u_exact / u_max, y / h, "k-", lw=2)
axes[0].fill_betweenx(y / h, u_exact / u_max, alpha=0.15)
axes[0].axvline(0, color="gray", lw=0.8, ls="--")
axes[0].set_xlabel(r"$u / u_{\max}$")
axes[0].set_ylabel(r"$y / h$")
axes[0].set_title("ポアズイユ速度プロファイル")
axes[0].set_xlim(-0.05, 1.1)
 
# せん断応力分布
tau = mu * np.gradient(u_exact, y)
axes[1].plot(tau / G / h, y / h, "r-", lw=2, label=r"$\tau / \tau_w$")
axes[1].axvline(0, color="gray", lw=0.8, ls="--")
axes[1].set_xlabel(r"$\tau / \tau_w$")
axes[1].set_ylabel(r"$y / h$")
axes[1].set_title("せん断応力分布 (線形)")
axes[1].legend()
 
plt.tight_layout()
plt.savefig("poiseuille_profile.png", dpi=150)
plt.show()

有限差分法による数値検証

$d^2u/dy^2 = -G/\mu$ を2次中心差分で離散化します。一様格子 $\Delta y = 2h/(N-1)$：

$\frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta y)^2} = -\frac{G}{\mu}$

境界条件：$u_0 = u_{N-1} = 0$

これを三対角行列（tridiagonal matrix）の連立方程式として書くと：

$\begin{pmatrix} -2 & 1 \\ 1 & -2 & 1 \\ & \ddots & \ddots & \ddots \\ & & 1 & -2 \end{pmatrix} \mathbf{u} = -\frac{G (\Delta y)^2}{\mu} \mathbf{1}$

Pythonでの実装#

import numpy as np
 
def solve_poiseuille_fd(N, h=1.0, G=1.0, mu=1.0):
    """
    N : 内部格子点数 (境界を除く)
    returns (y, u_num, u_exact, L2_error)
    """
    dy = 2 * h / (N + 1)
    y_inner = np.linspace(-h + dy, h - dy, N)
 
    # 三対角係数行列
    diag  = -2 * np.ones(N)
    upper =  np.ones(N - 1)
    lower =  np.ones(N - 1)
    A = (np.diag(diag) + np.diag(upper, k=1) + np.diag(lower, k=-1))
 
    rhs = -G * dy**2 / mu * np.ones(N)
 
    u_num = np.linalg.solve(A, rhs)
 
    # 境界を含む全格子
    y_full  = np.concatenate([[-h], y_inner, [h]])
    u_full  = np.concatenate([[0.0], u_num, [0.0]])
    u_exact = G / (2 * mu) * (h**2 - y_full**2)
 
    L2 = np.sqrt(np.mean((u_full - u_exact)**2))
    return y_full, u_full, u_exact, L2
 
# 格子収束テスト
Ns = [4, 8, 16, 32, 64, 128, 256]
errors = []
for N in Ns:
    _, _, _, err = solve_poiseuille_fd(N)
    errors.append(err)
 
# 収束次数
for i in range(1, len(Ns)):
    ratio = errors[i-1] / errors[i]
    order = np.log2(ratio)
    print(f"N={Ns[i]:4d}  L2={errors[i]:.3e}  収束次数={order:.2f}")

収束結果

“

備考: ポアズイユ流れは2次多項式であるため、2次中心差分を用いると**マシン精度（machine precision）**レベルの正確さが得られます。実際の誤差は $\Delta y^2$ に比例するのではなく、浮動小数点演算の丸め誤差レベルで飽和します。

低い $N$ でも放物線は式に張り付きます。粘性項が死んでいるかはここで分かります。次数の列はその判定に使えません。

離散化の文は $\Delta y = 2h/(N-1)$ と書いてあります。表は別です。$N=4,8,16,32,64$ の間隔は 0.4000、0.2222、0.1176、0.0606、0.0308 です。$h=1$ なら $2/(N+1)$ と一致します。$2/5=0.4$、$2/9 \approx 0.2222$、$2/17 \approx 0.1176$、$2/33 \approx 0.0606$、$2/65 \approx 0.0308$ です。コードは両壁を未知数から外すので、内部 $N$ 点の間隔は $2h/(N+1)$ です。$N-1$ は境界点まで未知数に数えた個数です。

$L_2$ の五つは $2.78\times 10^{-16}$、$1.67\times 10^{-16}$、$2.22\times 10^{-16}$、$2.78\times 10^{-16}$、$2.22\times 10^{-16}$ と上下します。隣り合う比の二進対数はおよそ 0.7、次の二つは負、最後はおよそ 0.3 です。表の ~2.0 とは違います。$y$ の2次式は3階導関数が0なので、2次中心差分の打ち切り誤差は0です。この列では次数を測れません。

解析結果の比較プロット

import matplotlib.pyplot as plt
 
fig, ax = plt.subplots(figsize=(6, 5))
 
colors = plt.cm.viridis(np.linspace(0.2, 0.9, 4))
Ns_plot = [4, 8, 16, 64]
 
for N, c in zip(Ns_plot, colors):
    y, u_num, _, _ = solve_poiseuille_fd(N)[:3]
    _, _, u_exact, _ = solve_poiseuille_fd(N)
    ax.plot(u_num, y, "o--", color=c, ms=4, label=f"FD N={N}")
 
# 解析解
y_fine = np.linspace(-1, 1, 400)
ax.plot(0.5 * (1 - y_fine**2), y_fine, "k-", lw=2, label="解析解")
 
ax.set_xlabel(r"$u$ (無次元)")
ax.set_ylabel(r"$y/h$")
ax.set_title("FD vs 解析解の比較")
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig("poiseuille_fd_vs_exact.png", dpi=150)
plt.show()

OpenFOAM設定のヒント#

ポアズイユ流れを OpenFOAM の icoFoam で設定する際の主要パラメータ：

0/U — 速度境界条件#

boundaryField
{
    walls
    {
        type    noSlip;         // 壁面 no-slip
    }
    inlet
    {
        type    fixedValue;
        value   uniform (1 0 0); // または parabolicVelocity
    }
    outlet
    {
        type    zeroGradient;
    }
}

0/p — 圧力境界条件#

boundaryField
{
    inlet   { type fixedValue; value uniform 1; }
    outlet  { type fixedValue; value uniform 0; }
    walls   { type zeroGradient; }
}

constant/transportProperties#

nu              nu [ 0 2 -1 0 0 0 0 ] 0.01;  // 動粘性係数の調整

発達長さと中心速度

1. 完全発達条件：チャンネル長が入口の発達長さ（$L_{dev} \approx 0.05 \, Re \cdot h$）よりも十分に長いこと
2. 格子独立性：$y$ 方向の格子数が16以上であれば2次精度を確保
3. 検証指標：中心流速 $u_c = G h^2 / (2\mu)$ と壁面せん断応力 $\tau_w = G h$ を比較

L2が1e-16なら次数欄を見るか#

$N=4$ の時点で $L_2$ はすでに $2.78\times 10^{-16}$ です。格子を64まで増やしても、その欄は $10^{-16}$ 台で上下するだけです。~2.0 はその比から出ていません。2次多項式の上では中心差分の打ち切り誤差が0なので、この表の次数で合格を押す判断は外れます。$\mu \nabla^2 \mathbf{u}$ が生きているかは、中心速度と壁面応力で見ます。

照合する数字は $Gh^2/(2\mu)$ と $Gh$ です。次数が実際に分かれるデータは、不連続を含む解の側にあります。同じコードが0.79次と2.00次を同時に出した場合 は、次数欄がスキーム次数とは限らないことを別に測っています。残差が丸めに張り付いているここでは、格子を足しても次数は出ません。

“

参考文献

• Batchelor, G. K. An Introduction to Fluid Dynamics. Cambridge, 1967.
• Ferziger, J. H., Perić, M. Computational Methods for Fluid Dynamics. Springer, 2002.
• OpenFOAM Documentation — icoFoam tutorial: cavity