阻力系数削掉 9%，液滴却慢了 35% — 变形与蒸发修正该加在哪里

蒸发修正压低了 $C_D$，但减速取决于 $C_D/d$。直径一旦开始缩小，省下的那点系数就又还了回去。

1983 年，燃烧液滴面前的阻力系数对不上#

Renksizbulut 和 Yuen 测量了高温气流中蒸发液滴的阻力。实测阻力始终低于球形关联式给出的预测值。而当时液滴仍然保持着球形。真正的原因是，从表面喷出的蒸气把边界层往外顶开了。

同一时期，人们也报告了方向相反的偏差。置于高速气流中的液滴会被压扁，一旦压扁，阻力就比球形更大。O'Rourke 和 Amsden 在 1987 年用 TAB 模型总结的正是这一侧。

所以喷雾程序里的液滴阻力系数通常是三层的：一个球形关联式，再乘上一个变形修正和一个蒸发修正。本文把这三层逐层剥开，在相同条件下测量每一层各自把轨迹推动了多少。先说结论：把系数削得最狠的那个模型，最先停下来。

三个关联式在同一个位置改出不同方向

起点是球形液滴的阻力系数，这里采用 Schiller–Naumann 一系的关联式。

$C_{D,\mathrm{sph}} = \frac{24}{Re}\left(1 + \tfrac{1}{6} Re^{2/3}\right), \qquad Re < 1000$

$Re = \rho_g u_{\mathrm{rel}} d / \mu_g$ 是用液滴直径和相对速度构造的雷诺数。前面的 $24/Re$ 是 Stokes 阻力，括号里则是惯性修正。

挂在它上面的两个修正，性质并不相同。

修正	相乘因子	方向	依据
无 (球形)	$1$	—	刚性球，表面无物质传递
TAB 变形	$1 + 2.632\,y$	阻力增大	迎风面积变大
蒸发 (R–Y)	$(1+B_M)^{-0.2}$	阻力减小	蒸气把边界层推开

$y$ 是 TAB 模型的无量纲变形量。 $y = 0$ 表示球形， $y = 1$ 判定为破碎。 $B_M$ 是 Spalding 物质传递数(把表面与远场的蒸气质量分数之差无量纲化后的量)，也就是蒸发源项与 $d^2$ 定律里讨论过的那个数。

变形究竟是怎么长起来的，下面直接动手调一调。

relative speed25 m/sdiameter40 umgas density5.0

y 0.000y_ss 0.000C_D x1.00We_r 0.00Re 0

Raise the speed and watch the drop flatten into a disc while the amber bar climbs past the grey sphere value. The oscillation overshoots the dashed steady line by almost 2x, so the peak drag is roughly twice what an equilibrium estimate would give. Turn liquid damping off and the ring never settles. Push far enough and y hits 1 — TAB stops and calls it breakup.

把相对速度和直径调大，液滴会被压成圆盘，系数条也会越过灰色的球形基准值。要观察的是：虚线画出的稳态值 $y_{ss}$ ，实际的 $y$ 究竟冲过头多少。

变形加在系数上，而不是迎风面积上

TAB 把液滴看成一个阻尼谐振子。气体动压推它，表面张力把它拉回，液体粘性让它衰减。

\frac{d^2 y}{dt^2} = \frac{C_F}{C_b}\frac{\rho_g u_{\mathrm{rel}}^2}{\rho_l r^2} - \frac{C_k \sigma}{\rho_l r^3} y - \frac{C_d \mu_l}{\rho_l r^2}\frac{dy}{dt}

$r$ 是液滴半径， $\sigma$ 是表面张力， $\mu_l$ 是液体粘性。常数取 $C_F = 1/3$ 、 $C_k = 8$ 、 $C_d = 5$ 、 $C_b = 1/2$ 。这个方程一路推到破碎的过程，在Weber 数与 TAB 破碎里讲过。

把所有导数项置零，就得到稳态变形量。

$y_{ss} = \frac{C_F}{C_b C_k}We_r = \frac{We_r}{12}, \qquad We_r = \frac{\rho_g u_{\mathrm{rel}}^2 r}{\sigma}$

这里有一个工程实践中经常出错的地方。把 $y_{ss}$ 直接塞进阻力修正的代码很常见。可这个振子几乎没有阻尼。正庚烷 40 μm 液滴的阻尼比只有 $10^{-2}$ 量级。从静止出发的 $y$ 会冲过 $y_{ss}$ ，一直涨到接近两倍。

也就是说，用稳态值算出的阻力只有真实第一次振荡的一半。出问题的不是系数本身，而是取用这个系数的时刻。

蒸发把边界层顶开，从而削掉系数

Renksizbulut–Yuen 关联式把蒸发处理成一个乘性修正。

$C_D = \frac{24}{Re}\left(1 + 0.2\,Re^{0.63}\right)(1+B_M)^{-0.2}$

关键在于指数只有 $-0.2$ ，很小。 $B_M = 0.6$ 时因子是 $0.91$ ， $B_M = 2$ 时是 $0.80$ 。即便是剧烈燃烧的液滴，系数也不过削掉 20% 出头。

原理与膜理论(film theory)相同。从表面流出的蒸气通量让边界层变厚，速度梯度变缓，壁面剪切减小，阻力随之下降。传热的 Nusselt 数上也挂着同样形式的修正。

用 Python 从同一个喷嘴打出四个液滴#

向静止气体中打入一个液滴，追踪 1 ms。控制方程是一维的。

$\frac{du_d}{dt} = -\frac{3}{4}\frac{\rho_g}{\rho_l}\frac{C_D}{d}\,u_d|u_d|$

四个算例只有阻力系数不同。A 是球形，B 是球形乘上 TAB 变形，C 只加蒸发修正，D 在蒸发修正之上再加 $d^2$ 定律的收缩。

import math
 
RHO_G, MU_G = 4.98, 3.3e-5      # 700 K, 10 bar 空气
RHO_L, SIG, MU_L = 620.0, 0.0145, 3.0e-4   # 正庚烷
D0, U0 = 40e-6, 25.0            # 破碎后的液滴、喷射速度
K_EVAP = 1.0e-6                 # d^2 定律常数 [m^2/s]
CF, CK, CD_T, CB = 1.0 / 3.0, 8.0, 5.0, 0.5   # TAB 常数
 
def cd_sphere(re):
    if re < 1e-8:
        return 0.0
    return 24.0 / re * (1.0 + re ** (2.0 / 3.0) / 6.0) if re < 1000.0 else 0.424
 
def cd_distorted(re, y):
    return cd_sphere(re) * (1.0 + 2.632 * max(0.0, min(1.0, y)))
 
def cd_blowing(re, bm):
    if re < 1e-8:
        return 0.0
    return 24.0 / re * (1.0 + 0.2 * re ** 0.63) * (1.0 + bm) ** -0.2
 
def tab_oscillate(y, ydot, u, d, dt):
    r = 0.5 * d
    acc = (CF / CB) * RHO_G * u * u / (RHO_L * r * r) \
        - (CK * SIG / (RHO_L * r ** 3)) * y \
        - (CD_T * MU_L / (RHO_L * r * r)) * ydot
    ydot += acc * dt
    y += ydot * dt
    return y, ydot
 
def run_droplet(law, evaporating, bm=0.6, t_end=1.0e-3, dt=5e-8):
    u, x, d, t = U0, 0.0, D0, 0.0
    y, ydot, ymax = 0.0, 0.0, 0.0
    cd_sum, cd0, n = 0.0, None, 0
    while t < t_end and u > 1e-6:
        re = RHO_G * u * d / MU_G
        if law == 'sphere':
            cd = cd_sphere(re)
        elif law == 'tab':
            y, ydot = tab_oscillate(y, ydot, u, d, dt)
            ymax = max(ymax, y)
            cd = cd_distorted(re, y)
        else:
            cd = cd_blowing(re, bm)
        if cd0 is None:
            cd0 = cd
        cd_sum += cd
        n += 1
        u += -0.75 * RHO_G * cd * u * u / (RHO_L * d) * dt
        x += u * dt
        if evaporating:
            d = math.sqrt(max(1e-18, d * d - K_EVAP * dt))
        t += dt
    return dict(x_mm=x * 1e3, u=u, d_um=d * 1e6, cd0=cd0,
                cd_mean=cd_sum / n, ymax=ymax)
 
CASES = [
    ('A sphere            ', 'sphere', False),
    ('B sphere x TAB      ', 'tab', False),
    ('C blowing (B_M=0.6) ', 'blow', False),
    ('D blowing + d2-law  ', 'blow', True),
]
 
print('Re0 = %.0f   We_r = %.2f   t = 1.0 ms' %
      (RHO_G * U0 * D0 / MU_G, RHO_G * U0 * U0 * 0.5 * D0 / SIG))
print('case                  Cd(t=0)  <Cd>    u(1ms)   x(1ms)   d(1ms)')
base = None
for name, law, ev in CASES:
    r = run_droplet(law, ev)
    if base is None:
        base = r
    print('%s  %6.3f  %6.3f  %6.2f   %6.3f   %5.1f  (%+.1f%% x)' %
          (name, r['cd0'], r['cd_mean'], r['u'], r['x_mm'], r['d_um'],
           100.0 * (r['x_mm'] / base['x_mm'] - 1.0)))
 
tab = run_droplet('tab', False)
print('TAB peak distortion y_max = %.3f  ->  Cd multiplier %.2fx' %
      (tab['ymax'], 1.0 + 2.632 * tab['ymax']))

Re0 = 151   We_r = 4.29   t = 1.0 ms
case                  Cd(t=0)  <Cd>    u(1ms)   x(1ms)   d(1ms)
A sphere               0.910   1.708    3.36    9.261    40.0  (+0.0% x)
B sphere x TAB         0.910   2.233    2.66    7.490    40.0  (-19.1% x)
C blowing (B_M=0.6)    0.828   1.537    3.68    9.732    40.0  (+5.1% x)
D blowing + d2-law     0.828   2.096    2.18    8.868    24.5  (-4.3% x)
TAB peak distortion y_max = 0.632  ->  Cd multiplier 2.66x

变形修正把贯穿距离缩短了 19.1%，因为第一次振荡就把系数抬到了 2.66 倍。只加蒸发修正则相反，距离增加 5.1%。到这里为止，结果与各关联式的符号完全一致。

系数最低的那个为什么最先停下

问题出在 D 上。它用的是和 C 完全相同的阻力系数公式，初始系数同样是 0.828。可是 1 ms 之后的速度只有 2.18 m/s，是四个里最低的，比 A 的 3.36 m/s 慢了 35%。

再看一眼减速方程，答案就出来了。右端有 $C_D/d$ 。在阻力系数被削掉 9% 的同时，直径从 40 μm 缩到 24.5 μm，减少了 39%。分母这一侧动得大得多。

从物理上讲是这样：阻力正比于面积( $d^2$ )，惯性正比于质量( $d^3$ )，两者之比是 $1/d$ 。液滴越干瘪，同样的阻力就造成越大的减速。这与响应时间 $\tau_d \simeq \rho_l d^2 / 18\mu_g$ 按 $d^2$ 缩小是同一件事。

$\frac{a}{m} \propto \frac{d^2}{d^3} = \frac{1}{d}$

也就是说，蒸发同时具有压低阻力系数和放大减速这两种效果。前者被锁在 $-0.2$ 的指数里，后者直接挂在直径上。只要算得稍久一点，后者就会赢。

下面把四个液滴从同一个喷嘴同时打出去。

injection u025 m/sd040 umB_M0.60K1.00 mm2/s

A 0.00 mm (--)B 0.00 mm (-100.0%)C 0.00 mm (-100.0%)D 0.00 mm (-100.0%)

Push B_M up: lane C's drag coefficient drops and it pulls ahead of the grey sphere. Now push K up and lane D uses that same reduced coefficient yet falls behind — the shrinking diameter multiplies the deceleration faster than blowing removes it. Setting K to 0 collapses D onto C, which is the cleanest way to see which of the two effects is doing the work.

调高 $B_M$ ，C 泳道会超过灰色的球形基准。在这个状态下再调高 $K$ ，只有使用同一套系数的 D 泳道会掉队。把 $K$ 调回 0，D 又会与 C 精确重合，于是两种效果究竟哪一个在起作用，一次就分离开了。

什么时候该打开哪个关联式

情况	变形修正	蒸发修正	理由
$We_r < 0.5$ 、低温	关	关	$y_{ss} < 0.04$ ，修正不到 4%
柴油喷雾近场	开	开	相对速度最大，变形占主导
喷雾远场	关	开	减速后 $We$ 下降，只剩蒸发
固体颗粒	关	关	既无变形也无物质传递
近超临界	关	谨慎	$\sigma \to 0$ 使 TAB 失去意义， $B_M$ 发散

最后一行在实际工作中经常出问题。表面张力趋于 0 时，TAB 的回复项消失， $y$ 会发散。在那个区域应当关掉变形修正，转向扩散界面类的模型。

如果 $We$ 已经足够低、阻力却仍然对不上，那要查的就不是修正项，而是别的地方。比如壁面附近出现阻力危机与边界层分离那样的情形， $Re$ 本身已经跑出了关联式的适用范围。

换掉一个关联式，整套验证算例都会跟着动

在这次计算里，贯穿距离在 -19.1% 到 +5.1% 之间来回摆动。网格没动，时间步没动，湍流模型也没动。全部改动只是把乘在单个液滴上的两个无量纲因子开开关关。

所以在喷雾验证与实验对不上时，先去加密网格并不是正确的顺序。要先确认的是：阻力系数用的是哪个公式，其中的 $y$ 和 $B_M$ 取了什么值， $y$ 是按稳态代入的还是真的解了振荡方程。如果直径随时间在缩小，那就别只看 $C_D$ ，把 $C_D/d$ 输出来看。只盯着系数，液滴为什么先停下这件事到最后也看不出来。