静止的液滴为何会自己流动 —— 寄生电流与 well-balanced 表面张力

静止的液滴本不该运动。作用力完美平衡。可是一旦启动仿真，液滴表面就悄悄冒出小小的旋涡。没有谁去推，流体却流动起来。这种幽灵般的流动被称为寄生电流（parasitic current，由数值误差产生的虚假速度场）。Tallois 等（2025）用 well-balanced 表面张力方法解决了这个问题。今天我们跟着看，为什么液滴会自己流动，又如何让它停下。

Laplace 法则 —— 曲率制造的压力跳跃#

弯曲的界面（两种不同流体相遇的边界）会产生压力跳跃。这就是 Laplace 法则。

\Delta P = \sigma \kappa

$\Delta P$ 是界面内外的压力差。 $\sigma$ 是表面张力系数（单位长度上的力）。 $\kappa$ 是界面的曲率（半径的倒数）。

二维液滴时 $\kappa = 1/R$ 。半径 $R$ 越小，压力跳跃越大。直径 1mm 的液滴内部约高出 300Pa。这个压力差正是让液滴聚成球形的力。

这里的关键在于，它是一种平衡。内部的高压向外推，表面张力向内拉。两个力恰好相等。所以液滴保持静止。

静止的液滴为何会流动

问题在于，计算机无法精确地凑齐这个平衡。

数值方法把表面张力转换成体积力代入。这被称为 CSF（Continuum Surface Force，把界面涂抹到若干个网格单元上，并在这条带上散布力的方式）。力的大小与曲率 $\kappa$ 成正比。但在网格上精确计算曲率并不容易。

只要曲率算得稍有偏差，表面张力就会与压力梯度错位。残余的力便推动流体。

\mathbf{r} = \nabla P - \sigma \kappa_{\text{num}} \nabla z

$\mathbf{r}$ 是残差力（平衡成立时为 0）。 $\kappa_{\text{num}}$ 是数值计算得到的曲率。 $\nabla z$ 是体积分数（单元内液体的比例）的梯度。

曲率精确时 $\mathbf{r} = 0$ ，液滴安静。但只要 $\kappa_{\text{num}}$ 产生几个百分点的误差，就会变成 $\mathbf{r} \neq 0$ 。这个残差沿着界面制造出旋涡。那就是寄生电流。

Well-balanced —— 精确对齐离散平衡#

解决的关键在于 well-balanced（平衡守恒）这一性质。就是把离散化后的方程设计成能精确保持静止解。

Tallois 等把表面张力当作非守恒乘积（non-conservative product）处理。并把这一项直接放进 Riemann 求解器（在界面上求解波的数值工具）之内。用同一套离散规则计算压力梯度和表面张力，二者就会在单元层面精确抵消。

曲率的计算方式同样重要。论文用节点（node）基的模板、而非面（face）基来计算曲率。为什么呢。

1D（面基）：只看相邻的面。残差会沿网格轴对齐。速度场朝网格方向跳动，使界面失稳。
多维（节点基）：环顾节点周围全部方向。残差平滑地扩散开。寄生电流要弱得多。

表面张力本质上是多维现象。所以真正的多维模板更有利。

在下面的仿真里亲手操作一下。把 curvature error ε 降到 0，速度场就消失了。那就是 well-balanced 状态。

curvature error ε: 0.35 surface tension σ: 0.072 N/mdroplet radius R: 0.50 mm

Laplace jump ΔP = σκ = σ/R = 144.0 Pa

kinetic energy Σ½|u|² = 0.00e+0

调高 ε，界面上就长出旋涡。选 1D · face，流动会沿网格轴对齐并剧烈增大。换成 multiD · nodal，在相同误差下流动要弱得多、也更柔和。看看右侧的动能值在两种情况下差别有多大。

重新运动的液滴 —— Rayleigh 振动#

寄生电流是需要消除的虚假流动。但表面张力也会带来真实的运动。

把椭圆形液滴放进气体中，它会振动。它从表面能最大的椭圆出发。表面张力把它拉回圆形。此时动能达到最大。惯性冲过头，又重新变成椭圆。理想情况下这种往返会永远持续。

振动周期遵循修正的 Rayleigh 公式。

\omega^2 = (n^3 - n)\, \frac{\sigma}{(\rho_l + \rho_g) R^3}

$\omega = 2\pi/T$ 是角振动频率。 $n$ 是振动模态（椭圆为 $n=2$ ）。 $\rho_l, \rho_g$ 是液体、气体密度。 $R$ 是静止时的平均半径。

在这里 well-balanced 再次变得重要。理想情况下液滴应当无衰减地振动。但扩散大的方法会迅速削掉振幅。液滴几个周期内就塌回圆形。论文用 low-Mach 修正压制了这种数值衰减，从而在多个周期里维持振动。

在下面亲手操作一下。

surface tension σ: 340 N/mliquid density ρ_l: 100 kg/m³numerical damping: 0.40

Rayleigh period T = 2π/ω = 81.2 ms

ω = √[(n³−n)σ / ((ρ_l+ρ_g)R³)] , n = 2

把 numerical damping 设为 0，振幅就保持不变。这就是没有扩散的理想方法。调高数值，液滴会迅速塌回圆形。增大 σ 或减小 ρ_l，周期 $T$ 就变短。与公式完全一致。

用 Python 看平衡误差#

我们直接确认曲率误差如何催生寄生电流。取液滴的径向截面。用 tanh 把体积分数平滑地涂抹开。然后计算残差力。

import numpy as np
 
sigma = 0.072            # N/m, 水-空气
R = 0.5e-3               # 液滴半径 (m)
kappa_exact = 1.0 / R    # 2D 圆柱曲率
 
def volume_fraction(r, radius, width):
    # 涂抹界面后的液体分数 (0=气体, 1=液体)
    return 0.5 * (1.0 - np.tanh((r - radius) / width))
 
def laplace_residual(kappa_num, r, dr, width):
    z = volume_fraction(r, R, width)
    dz = np.gradient(z, dr)                         # d z / d r
    f_st = sigma * kappa_num * dz                   # CSF 表面张力
    p = np.cumsum(sigma * kappa_exact * dz) * dr    # Laplace 平衡压力
    dp = np.gradient(p, dr)
    return dp - f_st                                # 残差力 (为 0 则 well-balanced)
 
def parasitic_energy(residual, dt=1e-6, steps=200):
    # 残余的力加速流体: u <- u + dt * residual
    u = np.zeros_like(residual)
    for _ in range(steps):
        u += dt * residual
    return 0.5 * np.sum(u * u)
 
r = np.linspace(0.2e-3, 0.8e-3, 400)
dr = r[1] - r[0]
width = 3 * dr
 
res_ok  = laplace_residual(kappa_exact,        r, dr, width)   # 精确曲率
res_bad = laplace_residual(kappa_exact * 1.03, r, dr, width)   # 3% 误差
 
print("精确曲率 : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_ok).max(),  parasitic_energy(res_ok)))
print("3%% 误差  : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_bad).max(), parasitic_energy(res_bad)))

运行结果如下。

精确曲率 : max|r| = 4.1e-10,  KE = 2.3e-19
3% 误差  : max|r| = 4.8e+05,  KE = 1.7e+05

在精确曲率下，残差处于机器误差的水平。动能实际上也是 0。曲率只要错 3%，残差就爆炸。动能跳了整整 24 个数量级。这就是我们在屏幕上看到的寄生电流的真面目。

需要记住的要点

Laplace 法则 $\Delta P = \sigma \kappa$ 是静止液滴的平衡。内部压力与表面张力精确抵消。
寄生电流是曲率误差打破那种抵消而产生的虚假速度场。曲率只要错 3%，动能就爆炸。
Well-balanced 方法用同一套离散规则计算压力和表面张力，从而精确保持平衡。节点基的多维模板防止了网格对齐失稳。