泊肃叶流动：解析解推导与有限差分验证

两块平板，间距 $2h$。压力只沿 $x$ 下降。充分发展时 $v=0$，$u$ 只剩 $y$。圆管也叫同一个名字（Poiseuille flow），下面的式子却是二维通道。解一旦离开这条抛物线，粘性项或无滑移就错了。

拨动 −dp/dx 或 μ 滑块，画面下方的 $u_{\max}$ 会变，曲线形状不动。

画布上的曲线是 $u/u_{\max} = 1 - (y/h)^2$。Re 滑块只改角落里的数字。$h=1$ 时，峰值读数是 −dp/dx 的值除以 $2\mu$。

问题设置

几何形状与边界条件

考虑两个无限大平行平板之间的完全发展层流。

• 通道半宽：$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$ 方向的 Navier–Stokes 方程为：

$\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$

解析解推导

第一步：简化为常微分方程 (ODE)#

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

第二步：一次积分

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

第三步：二次积分

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

第四步：应用边界条件

$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）**分布。壁面摩擦均匀地传递到整个流动截面，导致中心流速最大，壁面流速为零。

无量纲速度 $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}$

这被称为 Hagen-Poiseuille 定律 的 2D 版本。由于流量与 $h^3$ 成正比，通道宽度增加一倍，流量将增加八倍。

把代码默认值 $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$ 加倍，幂次就分开了。$u_{\max}$ 按 $h^2$ 走，变成 4 倍；$Q$ 按 $h^3$ 走，变成 8 倍；壁面应力大小 $Gh$ 只变成 2 倍。若把应力也当成跟着流量的倍数走，壁面值会被读成四倍大。转到圆管，把这条抛物线的动能除以平均速度，系数会变。层流圆管里该系数锁在 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$ 进行离散化。在均匀网格 $\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）线性方程组：

$\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}")

收敛结果

“

注：由于泊肃叶流动是二次多项式分布，因此使用二阶中心差分格式可以获得接近**机器精度（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$ 的二次式三阶导数为 0，二阶中心差分的截断误差就是 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;         // 壁面无滑移
    }
    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 以上即可确保二阶精度。
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 不是这些比值算出来的。在二次多项式上，中心差分的截断误差是 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