把 v 项和 v² 项相加，两端都对，中间却偏了 42% — Prony 两项阻力律与 Re=60

同样的水，同样的玻璃珠，不同的阻力律

往水里丢两颗玻璃珠。一颗直径 0.05 mm，另一颗 10 mm。小珠受到的阻力与速度成正比。大珠受到的阻力与速度的平方成正比。水相同，材料相同，只有大小不同。

18 世纪末的法国学者分别发现了这两条规律。随后他们把两者相加，写成 $F = a v^2 + b v$。本文用数字检验这个和在哪里成立、在哪里失效。读完后你会知道两项交接主导权时的 Reynolds 数（惯性力与黏性力之比），以及简单相加偏差最大的区间。

下面分成五个问题来走。

Q1. 1775 到 1800 年间，阻力到底是 v 还是 v²#

1752 年，达朗贝尔遇到了一个悖论：理想流体中物体受到的阻力为零。理论走不通，实验开了路。1775 年，法国科学院组建了流体阻力实验组，具体工作由 Bossut 负责。

Bossut 在水槽中拖动船模。精密测量表明，阻力与速度的平方成正比。与此同时，Du Buat 测量了管道内的流动。他看到的摩擦阻力与速度直接成正比。两组实验给出了不同的答案。

Prony 用更精密的实验把两者统一起来：阻力可以写成 $a v^2 + b v$。Coulomb 把它整理为两类阻力。与速度平方成正比的是形状阻力（form drag，由物体形状造成的压力差），与速度成正比的是摩擦阻力（friction drag）。这一划分后来成为 Navier 和 Stokes 研究黏性流动的基础。

谁对呢？两人都对。船模又大又快，管内流动又慢又窄。两组实验看到的是同一物理的两端。区分这两端的数字就是下一个问题。

Q2. 两项从何而来#

把球所受的阻力写成两项：

$\mu$ 是黏性系数，$\rho$ 是流体密度，$d$ 是球的直径，$C_\infty$ 是高速时的阻力系数（球约为 0.4）。

第一项是 Stokes 阻力，是完全没有惯性的慢速流动（蠕动流）的精确解。第二项是 Newton 阻力，假设撞上前表面的流体动量全部转化为压力即可得到。

需要注意一点。Stokes 阻力 $3\pi\mu d v$ 并不是纯粹的表面摩擦，其中 2/3 是摩擦，1/3 是压差阻力。因此比起 Coulomb 的"摩擦对形状"，"黏性主导对惯性主导"是更准确的叫法。

在下面的模拟中亲手操作一下。放下玻璃珠，速度上升时两种阻力随之增长。两者之和等于重量的那一点就是终端速度（terminal velocity，加速停止时的速度）。

把直径滑块从 0.02 mm 拖到 20 mm。小珠时蓝色直线（b·v）与重量线相交，大珠时橙色抛物线（a·v²）与之相交。把流体换成 glycerol，连 10 mm 的珠子也会移到 Du Buat 一侧。

Q3. 主导权何时交接#

求两项相等时的速度即可。

$Re$ 是以球直径为特征长度的 Reynolds 数。取 $C_\infty = 0.4$，得 $Re = 60$。

这个结果让 Bossut 和 Du Buat 达成了和解。决定阻力律的既不是流体种类，也不是速度本身，而是 $Re$。管内慢速流动的 $Re$ 很小，船模的 $Re$ 达到几十万。前面模拟里把水换成 glycerol，归根结底也是在移动 $Re$。

"连接两个极限"的同类问题在别处也出现过。在颗粒传热的 Knudsen 过渡区中，连续介质与自由分子两个极限是用调和平均连接的。在 Wood 声速中，压缩率按体积分数相加。极限怎么混合，会改变结果。Prony 选择的是简单相加。

Q4. 那么，直接相加就够了吗#

来验证一下。在水中从静止释放六颗玻璃珠（$\rho_p = 2500$ kg/m³），用两项阻力积分运动方程。

$m$ 是球的质量，$V$ 是体积，$g$ 是重力加速度。忽略附加质量（added mass，球带动的周围流体的惯性）。

把终端速度与参考值比较。参考值取 Schiller–Naumann 关联式 $C_D = \frac{24}{Re}(1 + 0.15 Re^{0.687})$，已知它在 $Re < 1000$ 时与实验相差在几个百分点以内。

import math
 
# 玻璃珠在水中沉降：Du Buat（线性）对 Bossut（平方）阻力
RHO_P, RHO, MU, G = 2500.0, 1000.0, 1.0e-3, 9.81
C_INF = 0.4          # 两项律中的平方（形状）阻力系数
 
def weight(d):
    return (RHO_P - RHO) * math.pi * d**3 / 6 * G
 
def prony_drag(d, v):
    lin = 3 * math.pi * MU * d * v                 # Du Buat / Stokes:  b v
    quad = C_INF * math.pi / 8 * RHO * d**2 * v**2  # Bossut / Newton:  a v^2
    return lin, quad
 
def vt_prony(d):
    a = C_INF * math.pi / 8 * RHO * d**2
    b = 3 * math.pi * MU * d
    return (-b + math.sqrt(b * b + 4 * a * weight(d))) / (2 * a)
 
def cd_schiller(re):
    return 24 / re * (1 + 0.15 * re**0.687) if re < 1000 else 0.44
 
def vt_reference(d):
    lo, hi = 1e-9, 10.0
    for _ in range(200):
        v = 0.5 * (lo + hi)
        re = RHO * v * d / MU
        f = cd_schiller(re) * 0.5 * RHO * v * v * math.pi * d * d / 4
        lo, hi = (v, hi) if f < weight(d) else (lo, v)
    return 0.5 * (lo + hi)
 
def settle(d, steps_per_tau=400):
    m = RHO_P * math.pi * d**3 / 6
    vt = vt_prony(d)
    tau = vt / ((1 - RHO / RHO_P) * G)       # 加速的时间尺度
    dt = tau / steps_per_tau
    v, t = 0.0, 0.0
    while v < 0.99 * vt:
        lin, quad = prony_drag(d, v)
        v += dt * (weight(d) - lin - quad) / m
        t += dt
    return t, tau, vt
 
print(" d[mm]   Re_t    vt[mm/s]  ref[mm/s]  err%   lin%   t99/tau  Stokes/vt  Newton/vt")
for dmm in [0.05, 0.1, 0.3, 1.0, 3.0, 10.0]:
    d = dmm * 1e-3
    t99, tau, vt = settle(d)
    lin, quad = prony_drag(d, vt)
    ref = vt_reference(d)
    v_st = weight(d) / (3 * math.pi * MU * d)
    v_nw = math.sqrt(weight(d) / (C_INF * math.pi / 8 * RHO * d * d))
    re = RHO * vt * d / MU
    print(f"{dmm:6.2f} {re:7.2f} {vt*1e3:10.2f} {ref*1e3:10.2f} {100*(vt-ref)/ref:6.1f} "
          f"{100*lin/(lin+quad):6.1f} {t99/tau:8.2f} {v_st/vt:10.2f} {v_nw/vt:10.2f}")
re_x = 24 / C_INF  # 3 pi mu d v = C_inf pi/8 rho d^2 v^2
print(f"crossover Re where lin = quad: 24/C_inf = {re_x:.1f}")

 d[mm]   Re_t    vt[mm/s]  ref[mm/s]  err%   lin%   t99/tau  Stokes/vt  Newton/vt
  0.05    0.10       2.04       1.98    2.9   99.8     4.59       1.00      24.27
  0.10    0.81       8.07       7.29   10.6   98.7     4.55       1.01       8.68
  0.30   17.16      57.21      40.21   42.3   77.8     3.92       1.29       2.12
  1.00  193.49     193.49     145.95   32.6   23.7     2.93       4.22       1.14
  3.00 1121.20     373.73     365.75    2.2    5.1     2.70      19.69       1.03
 10.00 6973.63     697.36     667.76    4.4    0.9     2.65     117.23       1.00
crossover Re where lin = quad: 24/C_inf = 60.0

在浏览器里直接运行上面的代码。两颗相同的珠子并排下沉，一颗受 Prony 两项阻力，另一颗受 Schiller–Naumann 参考阻力。

默认 d = 0.30 mm 时得到 err 42.3%、lin 77.8%、t99/τ 3.92，与输出块一致。右侧曲线是整个直径范围上的误差。把 C∞ 调到 0.6，大珠一侧的误差变为负值；降低 ρp，同一直径下 Re_t 变小，峰值向右移动。

输出里能读出三件事。第一，两端吻合得很好：0.05 mm 时 2.9%，3 mm 时 2.2%。第二，只用一项会在另一端崩溃。只用 Stokes，10 mm 珠快了 117 倍；只用 Newton，0.05 mm 珠快了 24 倍。第三，问题在中间：0.3 mm（$Re_t = 17$）偏 42%，1 mm（$Re_t = 193$）偏 33%。

加速时间也随阻力律变化。只有线性阻力时，速度按指数趋近，$t_{99}/\tau = \ln 100 = 4.61$。只有平方阻力时，速度按 $\tanh$ 趋近，$t_{99}/\tau = \operatorname{artanh} 0.99 = 2.65$。输出中的 4.59 和 2.65 就是这两个极限，中间直径落在两者之间。

Q5. 为什么中间会偏#

把两项律改写成阻力系数，原因就显现了。

两式在 $Re \to 0$ 时都趋于 $24/Re$，差别在括号里的第二项。

代入 $Re = 17$。Prony 给出 $1.41 + 0.40 = 1.81$，Schiller–Naumann 给出 $1.41 \times 2.05 = 2.90$。阻力比参考值小 38%。重量相同，Prony 珠必须落得更快才能平衡，结果就是终端速度偏大 42%。

缺少的阻力从哪里来？$Re$ 超过 1 后，前后对称的流动被打破。球后形成慢速尾流，在 $Re \approx 20$ 附近表面流动分离，出现回流区。这部分尾流压差阻力由黏性和惯性共同产生，既不干净地属于 $v$，也不属于 $v^2$。Prony 的加法里没有这个混合项的位置。

所以当两项大小相近时，也就是 $Re \approx 60$ 附近，相加的偏差最大。历史上两组实验各自正确，原因也一样：它们都在这个区间之外。

两项律今天仍然活着的地方

尽管存在失效区间，$a v^2 + b v$ 的形式并没有消失。多孔介质的 Ergun 方程和 Darcy–Forchheimer 定律正是这个形式，管道损失的 Prony 公式也同根同源。在多孔介质中，即使孔隙尺度的 $Re$ 跨越很宽的范围，两项之和也相当准确，因为系数是按实验拟合的。

Euler–Lagrange 颗粒程序则不同。单个颗粒的 $Re$ 往往正好落在这个问题区间（11000）：喷雾液滴、流化床颗粒、沉降颗粒都是如此。在这里使用 $24/Re + 0.4$，终端速度会偏 3040%，沉降时间和停留时间也偏同样多。

为程序选择阻力模型时，可以按以下顺序检查。

1. 先提取颗粒 $Re$ 的分布范围。小于 1 用 Stokes，大于 1000 用常数 $C_D$ 即可。
2. 若大量质量落在 1~1000 区间，使用 Schiller–Naumann 这类带中间项的关联式。
3. 跑一个单颗粒沉降测试，与上表中 0.3 mm 和 1 mm 两行比较。若出现 30% 误差，说明代码里藏着两项相加。

Bossut 的船和 Du Buat 的管子，是在不同 $Re$ 下看同一片水。两者之间的空白，单靠加法填不满。