直径读错1%，流量就差4% — 半个网格的壁面与Poiseuille的常数

黏性格式再完美，只要把壁面放进半个网格，流量就会以四倍的幅度偏掉。

验证算例差了4%，格式却没有问题#

圆管层流是少数几个有解析解的验证算例。把新写的黏性项离散放上去，流量低了4%。先怀疑格式是合理的顺序。可是把网格加密一倍，误差并没有减半，而是原封不动地留着。不收敛的误差不是离散误差，是几何误差。

问题出在壁面所站的位置。半径往里错1%，流量就掉4%。这篇文章追一下这个倍数4从哪里来，以及同样的指数4在1838年被Poiseuille在玻璃管实验里反过来用了什么。

流量挂在直径的四次方上

充分发展的圆管流动里，剩下的只有一条轴向动量平衡。

\frac{\mu}{r}\frac{d}{dr}\left(r\,\frac{du}{dr}\right) = \frac{dp}{dx}

$\mu$ 是黏性系数， $r$ 是从中心轴量起的径向坐标， $u$ 是轴向速度。在中心取 $du/dr = 0$ 、在壁面取 $u = 0$ ，就得到一条抛物线。

u(r) = \frac{\Delta p}{4\mu L}\left(R^{2} - r^{2}\right)

$\Delta p$ 是长度 $L$ 区间上的压降， $R$ 是管半径。把它在截面上积分就是流量。

Q = \int_{0}^{R} u(r)\,2\pi r\,dr = \frac{\pi\,\Delta p\,D^{4}}{128\,\mu L}

$D = 2R$ 是管直径。积分里多乘了一次 $r\,dr$ ， $R^{2}$ 于是变成 $R^{4}$ 。四次方的来历只有这一处。

这个指数一旦作用在误差上，就变成倍数。

\frac{\delta Q}{Q} = 4\,\frac{\delta D}{D} + \mathcal{O}\!\left(\left(\frac{\delta D}{D}\right)^{2}\right)

$\delta D$ 是直径误差， $\delta Q$ 是由此产生的流量误差。壁面位置误差在其他验证量里通常只是一阶出现。只有流量会四倍奉还。

在下面的模拟里亲手操作一下。

dD/D-2.5%dp/dx1.0

Drag dD/D to 1 % and watch the red bar: the wall moved by 1 %, the flow rate moved by 4.06 %. The tracer counters at the right are an independent measurement — they know nothing about the formula, and they still settle on 0.904 (currently 1.000 after 0 particles). Switch to wall pinned at cell centre and lower N: a mesh of 10 cells loses 19 % of the flow with a perfectly correct viscous scheme.

用 dD/D 滑块只移动下方管子的壁面，看右侧的红色条和计数器。计数器是一次独立的测量，它不知道公式，只是数着穿过出口的粒子，却仍然收敛到 $(1+e)^{4}$ 。切到 wall pinned at cell centre 并把N降到10，在黏性格式完全正确的情况下，流量少了19%。

1822年Navier留下的空位，1838年Poiseuille填上的指数#

Navier是桥梁工程师。流体阻力在实测里明明白白，欧拉方程里却写不出来——他被这一点抓住了。1822年，在弹性力学的数学基础尚未理清之前，他以速度拉普拉斯项的形式导出黏性项，接到欧拉方程上。在达朗贝尔之后停滞了半个多世纪的流体力学，就是从这里重新动起来的。

同一时期的École Polytechnique，Cauchy是教授，Coriolis在授课。旋转坐标系的能量方程出自水轮实验，也是在这所学校的同一片空气里。1816年学校一度关闭，刚入学的Poiseuille改了方向，去了医学院。

他成了医生，去测血流。血压如何随血管直径变化，是他的问题。为了模仿最细的血管，他自己拉玻璃管，最细的一根内径0.015 mm，是头发丝的五分之一。1838年他发表的结果是这个形状。

Q = K''\,\frac{\Delta p\,D^{4}}{L}, \qquad K'' = \frac{\pi}{128\,\mu}

右边的 $K''$ 是后来才填上的。Poiseuille本人并没有黏度这个概念，他只是把 $K''$ 写成一个常数。实验确定下来的不是常数，而是指数。

用Python量一下"半个网格的壁面"值多少#

把上面的动量平衡用径向有限体积解了一遍。壁面无滑移条件加在壁面上的情形，和干脆钉在最后一个网格中心的情形，并排跑。

import math
 
MU, GRAD_P, RADIUS = 1.0e-3, 100.0, 1.0e-3   # Pa*s, Pa/m, m
 
 
def poiseuille_q(radius, mu=MU, grad_p=GRAD_P):
    """解析流量  Q = pi*G*R^4/(8 mu)"""
    return math.pi * grad_p * radius ** 4 / (8.0 * mu)
 
 
def solve_pipe_fv(ncell, wall_at_face=True, radius=RADIUS, mu=MU, grad_p=GRAD_P):
    """轴对称充分发展速度分布，环形网格ncell个的有限体积"""
    dr = radius / ncell
    rf = [i * dr for i in range(ncell + 1)]          # 面的半径
    lo, dg, up, rhs = [0.0] * ncell, [0.0] * ncell, [0.0] * ncell, [0.0] * ncell
    for i in range(ncell):
        rhs[i] = -grad_p * (rf[i + 1] ** 2 - rf[i] ** 2) / 2.0
        if i > 0:
            w = mu * rf[i] / dr
            lo[i], dg[i] = w, dg[i] - w
        if i < ncell - 1:
            e = mu * rf[i + 1] / dr
            up[i], dg[i] = e, dg[i] - e
        else:
            if wall_at_face:                          # 无滑移条件加在壁面上
                dg[i] -= mu * rf[i + 1] / (dr / 2.0)
            else:                                     # 钉在网格中心的无滑移条件
                lo[i], dg[i], up[i], rhs[i] = 0.0, 1.0, 0.0, 0.0
    for i in range(1, ncell):                         # Thomas算法
        m = lo[i] / dg[i - 1]
        dg[i] -= m * up[i - 1]
        rhs[i] -= m * rhs[i - 1]
    u = [0.0] * ncell
    u[-1] = rhs[-1] / dg[-1]
    for i in range(ncell - 2, -1, -1):
        u[i] = (rhs[i] - up[i] * u[i + 1]) / dg[i]
    q = sum(u[i] * math.pi * (rf[i + 1] ** 2 - rf[i] ** 2) for i in range(ncell))
    return q
 
 
def fit_slope(xs, ys):
    """log y 对 log x 的最小二乘斜率"""
    lx = [math.log(x) for x in xs]
    ly = [math.log(y) for y in ys]
    n = len(lx)
    mx, my = sum(lx) / n, sum(ly) / n
    num = sum((lx[i] - mx) * (ly[i] - my) for i in range(n))
    den = sum((lx[i] - mx) ** 2 for i in range(n))
    return num / den
 
 
_seed = 20260901
 
 
def unit_normal():
    """架在LCG上的Box-Muller，保证在哪里跑都是同一组数"""
    global _seed
    out = []
    for _ in range(2):
        _seed = (1103515245 * _seed + 12345) % (2 ** 31)
        out.append((_seed + 0.5) / 2 ** 31)
    return math.sqrt(-2.0 * math.log(out[0])) * math.cos(2 * math.pi * out[1])
 
 
def measured_exponent(dmin, dmax, ntube, noise, repeat=200):
    """Poiseuille的实验: 流量取真值，直径带着相对误差读出"""
    slopes = []
    for _ in range(repeat):
        ds, qs = [], []
        for k in range(ntube):
            f = k / (ntube - 1)
            d_true = dmin * (dmax / dmin) ** f
            qs.append(poiseuille_q(d_true / 2.0))
            ds.append(d_true * (1.0 + noise * unit_normal()))
        slopes.append(fit_slope(ds, qs))
    mean = sum(slopes) / len(slopes)
    sd = math.sqrt(sum((s - mean) ** 2 for s in slopes) / len(slopes))
    return mean, sd
 
 
print("[1] wall half a cell off  (R = 1.000 mm)")
print(" N   Q_face/Q_exact   Q_centre/Q_exact   (1-1/2N)^4")
qex = poiseuille_q(RADIUS)
for n in (10, 20, 40, 80):
    a = solve_pipe_fv(n, True) / qex
    b = solve_pipe_fv(n, False) / qex
    print(f"{n:3d}      {a:.4f}            {b:.4f}          {(1-1/(2*n))**4:.4f}")
 
print()
print("[2] 1% error on D vs the resulting error on Q")
for e in (0.005, 0.01, 0.02):
    print(f"  dD/D = {e*100:4.1f}%  ->  dQ/Q = {((1+e)**4-1)*100:5.2f}%")
 
print()
print("[3] exponent fitted from 12 tubes, 200 repeats")
print(" range of D            noise on D    exponent (mean +- sd)")
for (lo_d, hi_d, tag) in ((0.10e-3, 0.30e-3, "0.10 - 0.30 mm"),
                          (0.015e-3, 0.60e-3, "0.015- 0.60 mm")):
    for nz in (0.0, 0.01, 0.03):
        m, s = measured_exponent(lo_d, hi_d, 12, nz)
        print(f" {tag}        {nz*100:4.1f}%       {m:.3f} +- {s:.3f}")

[1] wall half a cell off  (R = 1.000 mm)
 N   Q_face/Q_exact   Q_centre/Q_exact   (1-1/2N)^4
 10      1.0100            0.8100          0.8145
 20      1.0025            0.9025          0.9037
 40      1.0006            0.9506          0.9509
 80      1.0002            0.9752          0.9752
 
[2] 1% error on D vs the resulting error on Q
  dD/D =  0.5%  ->  dQ/Q =  2.02%
  dD/D =  1.0%  ->  dQ/Q =  4.06%
  dD/D =  2.0%  ->  dQ/Q =  8.24%
 
[3] exponent fitted from 12 tubes, 200 repeats
 range of D            noise on D    exponent (mean +- sd)
 0.10 - 0.30 mm         0.0%       4.000 +- 0.000
 0.10 - 0.30 mm         1.0%       3.996 +- 0.033
 0.10 - 0.30 mm         3.0%       3.993 +- 0.098
 0.015- 0.60 mm         0.0%       4.000 +- 0.000
 0.015- 0.60 mm         1.0%       4.000 +- 0.009
 0.015- 0.60 mm         3.0%       3.995 +- 0.028

第一张表的第二列，每加密一倍网格，误差就降到四分之一。这是二阶收敛。第三列不是这样。误差只减半，并且与第四列的 $(1-1/2N)^{4}$ 吻合到小数点后三位。

这个吻合就是诊断。把无滑移条件钉在网格中心的那一刻，计算解的是半径 $R - \Delta r/2$ 的管子。不是离散误差，是另一根管子。而那半个网格在流量上变成四倍。N = 20时半径小2.5%，流量小9.75%。

把壁面放在节点之间的一半处，并不是有限体积独有的问题。格子Boltzmann的bounce-back同样把壁面立在两个节点的中间。无论哪一种，都要先确定"代码认为壁面在哪里"。

指数4能在噪声里活下来吗#

Poiseuille的处境和我们相反。他不是知道 $D$ 去预测 $Q$ ，而是要从测到的 $Q$ 里读出指数。而直径恰恰是最难测的量。想一想把内径0.015 mm的玻璃管读到1%精度是什么工作。

第三张表就是那个实验。只在直径上加相对误差，把 $\log Q$ – $\log D$ 的斜率重新拟合200次。管径范围是0.10–0.30 mm时，1%的噪声让斜率抖动 $\pm 0.033$ 。把范围拓宽到0.015–0.60 mm，同样的噪声、同样的管数，抖动降到 $\pm 0.009$ 。

原因用一条回归估计就能看见。

\sigma_{n} \approx \frac{n\,\sigma_{\delta D/D}}{\sqrt{N}\;\sigma_{\ln D}}

$n = 4$ 是真指数， $N$ 是管数， $\sigma_{\ln D}$ 是这些管子在 $\log D$ 轴上散开的宽度。窄范围给出 $4 \times 0.01 / (\sqrt{12} \times 0.331) = 0.035$ ，宽范围给出 $0.010$ 。表里的0.033和0.009正是它们。

指数的精度不是靠把尺子做精来换的，是靠把 $\log D$ 轴上的杠杆拉长来换的。Poiseuille一路做到头发丝五分之一的粗细，理由就在这里。

noise on D3.0%Dmax/Dminx3tubes12

Leave the bench running with Dmax/Dmin at 3: the histogram spreads out and the fitted exponent wanders (mean 4.000, sd 0.000 over 0 fits). Now push the range to 40 without touching the noise — the same clumsy ruler, the same number of tubes, and the histogram collapses onto 4. The exponent is bought with the lever arm in log D, not with a better ruler.

先把 Dmax/Dmin 放在3，看着直方图铺开，然后噪声不动，只把范围推到40。尺子没换，管数没变，分布却收拢到4。

在Stokes把 $K''$ 换成 $\pi/128$ 之后#

Poiseuille的 $K''$ 变成 $\pi/128\mu$ ，不是一次单位整理。这是Navier方程里的 $\mu$ 走进实验留下的那个空常数位子的事件。Stokes从Navier方程里数学地导出这条定律，由此确认：在管子里测到的量，和喂给方程的物性，是同一个东西。

从那以后这个式子的方向反了过来。它不再是测指数的实验，而成了测 $\mu$ 的黏度计。CGS的黏度单位poise来自他的名字。方程经过Poiseuille的实验和Stokes的补强，确立为Navier–Stokes方程，并在2000年成为克雷研究所七大难题之一。

同样的结构在我们的代码里反复出现。壁面函数常数、有效直径、接触角——实验给出形式，理论填上系数。当形式像 $D^{4}$ 这样陡峭时，几何会比系数更早主导误差。

圆管验证里我最先看的东西

流量对不上时先去翻格式，是在扔时间。顺序是这样的。

先把网格加密一倍，看误差是否降到四分之一。不降，就不是离散的问题。接着把误差比与 $(1-1/2N)^{4}$ 这类几何因子对照。对上了，就是壁面位置。最后用流量误差除以直径误差。商接近4，说明一个半径就解释了全部。

Poiseuille无法知道直径，所以读指数；我们知道指数，所以能反推直径。同一个 $D^{4}$ ，从两头各用一次。