总压下降1,080 Pa的地方，耗散却是零 — 当伯努利常数横跨流线时

总压横跨流线变化不是损失，而是涡度。只有总压沿流线下降时，才真的有损失发生。

1738年，父亲把自己书的出版年份提前了六年#

Daniel Bernoulli 在1738年出版了《Hydrodynamica》。扉页上写着“Johann 之子”。他想跟父亲和解。父亲 Johann 却把几乎相同的内容另外写成《Hydraulica》出版，还向出版商施压，把发行年份印成1732年。这样看起来就比儿子早了六年。

一个半世纪之后，Horace Lamb 把这件事理清楚了：对欧拉方程积分，得到的就是伯努利定理。父子二人争夺的两本书，其实是同一个方程的两副面孔。不过这个“积分”附带着一个条件，而在 CFD 后处理中让人误读总压云图的原因，恰恰就是这个条件。

本文用一个兰金涡来量这个条件。先把结论写在前面：在旋转的涡核内部，总压下降了 1,080 Pa，粘性耗散却恰好为零。而在总压完全平坦的外区，耗散反而不为零。两者刚好相反。

在下面的模拟里，把探针直接移到涡核内外试试。

core a0.50 mOMEGA60 1/sprobe r0.25 m

u = 0.00 m/sp = 0.0 Pap0 = 0.0 Pad(p0)/dr = 0.0rho u omega = 0.0phi = 0.0000 W/m^3

Drag probe r across the dashed core edge. Inside the core the two right-hand numbers stay equal and nonzero, so p0 climbs with r while the red element rotates without changing shape. Outside, both go to zero — p0 is perfectly flat — yet the purple element keeps shearing, which is where a viscous fluid would actually lose energy.

拖动 probe r，看右侧图中绿色的 $p_0$ 曲线从哪里开始变平。同时看画布左侧的方块什么时候保持形状，什么时候被压垮。这两个观察互相矛盾的位置，就是本文要讲的东西。

对欧拉方程积分就得到伯努利 — 沿哪个方向?

密度恒定、无体积力的定常无粘流动，其运动方程是这样的。

(\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p

$\mathbf{u}$ 是速度， $p$ 是压力， $\rho$ 是密度。对左边使用矢量恒等式。

(\mathbf{u}\cdot\nabla)\mathbf{u} = \nabla\!\left(\frac{|\mathbf{u}|^2}{2}\right) - \mathbf{u}\times\boldsymbol{\omega}

$\boldsymbol{\omega} = \nabla\times\mathbf{u}$ 是涡度。把两式合并，再用总压 $p_0 = p + \tfrac{1}{2}\rho|\mathbf{u}|^2$ 归并，剩下的只有一行。

\nabla p_0 = \rho\,\mathbf{u}\times\boldsymbol{\omega}

这就是 Crocco 形式的欧拉方程(不可压、等熵条件)。只要右端不为零，总压就在空间中变化。但关键在于它沿哪个方向变化。

$\mathbf{u}\times\boldsymbol{\omega}$ 垂直于 $\mathbf{u}$ 。所以两边点乘 $\mathbf{u}$ ，右端就消失了。

\mathbf{u}\cdot\nabla p_0 = 0

沿着流线走， $p_0$ 永远是常数。有没有涡度都一样。这才是伯努利定理准确的适用范围。反过来，要让 $p_0$ 在所有位置都相同，就必须 $\mathbf{u}\times\boldsymbol{\omega} = 0$ ，这实质上意味着无旋流动。父子相争的两本书的差别不在这里，条件藏在 Lamb 后来补上的那条注里。

兰金涡：旋转的涡核与不旋转的外区

兰金涡就是这个条件在同一个流场里同时开启又关闭的例子。半径 $a$ 以内像刚体一样转，以外是自由涡。

u_\theta(r) = \begin{cases} \Omega r, & r < a \\ \Omega a^2 / r, & r \ge a \end{cases} \qquad \omega_z(r) = \begin{cases} 2\Omega, & r < a \\ 0, & r \ge a \end{cases}

$\Omega$ 是涡核的角速度， $u_\theta$ 是周向速度。压力由径向动量方程 $\mathrm{d}p/\mathrm{d}r = \rho u_\theta^2 / r$ 积分得到。在外区，压力下降的量与动压上升的量恰好抵消， $p_0$ 被钉在 $p_\infty$ 上。在涡核内部则抵消不了。

p_0(r) - p_\infty = \rho\Omega^2 (r^2 - a^2), \qquad r < a

$r = 0$ 处的亏损是 $\rho\Omega^2 a^2$ ，也就是涡核边缘速度的平方乘以 $\rho$ 。取 $\rho = 1.2$ 、 $\Omega = 60\,\mathrm{s^{-1}}$ 、 $a = 0.5\,\mathrm{m}$ ，则 $u_{\max} = 30$ m/s, 亏损为1,080 Pa。

用 Python 在同一半径上把两边都量一遍#

用中心差分求 $\nabla p_0$ ，再与 $\rho\,u_\theta\omega_z$ 对照。耗散函数 $\phi = \mu(2S_{r\theta})^2$ 也一并测量。

import math
 
RHO, MU = 1.2, 1.8e-5      # kg/m^3, Pa*s
OMEGA, A = 60.0, 0.5       # 1/s, m   -> u_max = 30 m/s
P_INF = 101325.0           # Pa
 
def rankine_velocity(r):
    return OMEGA * r if r < A else OMEGA * A * A / r
 
def rankine_spin(r):
    return 2.0 * OMEGA if r < A else 0.0
 
def rankine_pressure(r):
    if r >= A:
        return P_INF - 0.5 * RHO * (OMEGA * A * A / r) ** 2
    return P_INF - RHO * OMEGA**2 * A**2 + 0.5 * RHO * OMEGA**2 * r**2
 
def total_head(r):
    return rankine_pressure(r) + 0.5 * RHO * rankine_velocity(r) ** 2
 
def shear_dissipation(r):
    s = 0.0 if r < A else -OMEGA * A * A / r**2   # 2*S_rtheta
    return MU * s * s
 
def crocco_residual(r, h=1e-6):
    lhs = (total_head(r + h) - total_head(r - h)) / (2 * h)   # d(p0)/dr
    rhs = RHO * rankine_velocity(r) * rankine_spin(r)         # rho*(u x omega)_r
    return lhs, rhs
 
print("  r[m]   u[m/s]   p[Pa]      p0[Pa]     dp0/dr    rho*u*w    phi[W/m^3]")
for r in (0.10, 0.25, 0.40, 0.60, 1.00):
    lhs, rhs = crocco_residual(r)
    print("%6.2f %8.2f %10.1f %10.1f %9.1f %10.1f %11.4f"
          % (r, rankine_velocity(r), rankine_pressure(r), total_head(r),
             lhs, rhs, shear_dissipation(r)))
 
def bernoulli_gap(r1, r2):
    return total_head(r2) - total_head(r1)
 
print()
print("core  p0(0.00) - p0(%.2f) = %8.1f Pa" % (A, total_head(0.0) - total_head(A)))
print("outer p0(%.2f) - p0(1.00) = %8.1f Pa" % (A, total_head(A) - total_head(1.0)))
print("cross-streamline gap r=0.10 -> 0.40 : %8.1f Pa" % bernoulli_gap(0.10, 0.40))
print("along-streamline  gap r=0.40 -> 0.40 : %8.1f Pa" % bernoulli_gap(0.40, 0.40))
print("dissipation at r=0.25 (core)  : %.4f W/m^3" % shear_dissipation(0.25))
print("dissipation at r=0.60 (outer) : %.4f W/m^3" % shear_dissipation(0.60))

  r[m]   u[m/s]   p[Pa]      p0[Pa]     dp0/dr    rho*u*w    phi[W/m^3]
  0.10     6.00   100266.6   100288.2     864.0      864.0      0.0000
  0.25    15.00   100380.0   100515.0    2160.0     2160.0      0.0000
  0.40    24.00   100590.6   100936.2    3456.0     3456.0      0.0000
  0.60    25.00   100950.0   101325.0       0.0        0.0      0.0313
  1.00    15.00   101190.0   101325.0       0.0        0.0      0.0040
 
core  p0(0.00) - p0(0.50) =  -1080.0 Pa
outer p0(0.50) - p0(1.00) =      0.0 Pa
cross-streamline gap r=0.10 -> 0.40 :    648.0 Pa
along-streamline  gap r=0.40 -> 0.40 :      0.0 Pa
dissipation at r=0.25 (core)  : 0.0000 W/m^3
dissipation at r=0.60 (outer) : 0.0313 W/m^3

第四列和第五列在所有半径上都相同。Crocco 关系精确到小数位都对得上。而在涡核内部， $r = 0.10$ 与 $r = 0.40$ 之间的总压差是648 Pa。同一半径上两个不同角度之间则是0。沿流线用伯努利是对的，横跨流线用就会丢掉这648 Pa。

总压平坦的地方，耗散反而不为零

最后一列是反着的。总压下降1,080 Pa 的涡核内部耗散为零，总压完全平坦的外区耗散却不为零。

原因是耗散挂在变形上，而不是挂在旋转上。刚体旋转只是把流体微团转过去，并不把它压扁。应变率张量为零，所以 $\phi = 0$ 。自由涡正好相反：涡度为零，但 $u_\theta \propto 1/r$ , 内侧和外侧以不同的速度掠过。微团一直在被剪切。

上面模拟中的方块把这一点原样展示了出来。在涡核内(红色)正方形保持原样地旋转，在外区(紫色)则塌成平行四边形。旋转坐标系里什么守恒、什么不守恒，科里奥利力与转焓一文里讲过，这里的轴线也一样：转动和做功记在不同的账本上。

把五种流动放进同一张表

涡度和耗散是彼此独立的。四种组合全都真实存在。

流动	$\boldsymbol{\omega}$	沿流线的 $p_0$	横跨流线的 $p_0$	粘性耗散
均匀流	0	恒定	恒定	0
自由涡 ( $r > a$ )	0	恒定	恒定	> 0
兰金涡核 ( $r < a$ )	$2\Omega$	恒定	变化1,080 Pa	0
剪切入口，无粘	$\ne 0$	恒定	变化	0
粘性尾流	$\ne 0$	下降	变化	> 0

表中加粗的格子只有一个，这一点很重要。能称为损失的，只有 $p_0$ 沿流线减小的那一种。其余四行里 $p_0$ 在空间上的变化，全都是涡度的几何，并不是能量消失了。

当总压云图因为两种原因看起来一样时

工程实践中这个区分最容易失守的地方，是出口截面的总压云图。把大气边界层或充分发展的管流作为入口条件给进去，从第一层网格开始 $p_0$ 就沿 $y$ 方向变化。而损失此时仍是零。尾流造成的总压亏损画出来也是同样的云图。这一边才是真正的损失。

下面把两个通道放在同一色标下并排跑一遍。

inlet shear S0.15wake deficit D0.32speed1.0x

A spread 0 PaA loss 0 PaB spread 0 PaB loss 0 Pa

Press match outlet spread — the two outlet profiles now cover the same range of p0, so a contour plot of the outlet cannot tell them apart. The dots can: in A every parcel keeps the color it entered with, in B the parcels that pass the body change color on the way out. Only a color change along a streamline is a loss.

按下 match outlet spread，两个通道的出口 $p_0$ 范围就对齐了。也就是说，单看出口云图已经无法区分。改看粒子的颜色在流动过程中会不会变。上面通道的粒子颜色一路保持不变，只有下面通道的粒子在绕过物体时改变了颜色。

所以计算损失系数时，如果把基准值取成截面平均 $p_0$ ，那么剪切入口明明没有损失，却会算出非零的值。基准应当是那条流线自己的入口 $p_0$ 。即便用流量加权平均，也必须对入口和出口采用同一种平均方式，差值里剩下的才只有损失。

数值上还要再加一条。网格偏粗时，旋转区域里的 $p_0$ 会被人为地抹平。数值扩散一旦把涡度糊掉， $\nabla p_0 = \rho\,\mathbf{u}\times\boldsymbol{\omega}$ 的右端就变小，涡核的亏损会比实际更浅。在检查穿过涡核的网格时，总压亏损的深度可以当作涡度分辨率的一个指标。就像圆管验证中直径误差以四次方压在流量上那样，这里的误差同样是通过一个不显眼的量溜进来的。

达朗贝尔在1752年撞上的也是同一处#

Johann 至死没有接受牛顿的粘性理论，他的儿子 Daniel 和学生欧拉也一样。只用无粘理论求解绕物体的流动，阻力会算成零。这就是1752年达朗贝尔公布这个结果时学界陷入恐慌的原因。它与达朗贝尔同一时期在波动方程上撞到的问题同根同源：方程所允许的解，其范围到底该划到哪里。

用我们今天的说法重写一遍是这样。在无旋无粘流动中， $p_0$ 在整个区域是同一个常数，于是物体前后的压力对称，积分为零。要产生阻力，就必须在某处让 $p_0$ 沿流线下降。提供这个位置的，正是粘性以及它在壁面上生产的涡度。

所以打开总压云图时该问的，不是“下降了多少”。而是“是沿流线下降的，还是横跨流线下降的”。前者是损失，后者是涡度。