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

液滴本来是静止的。力是齐的。计算一开，表面长出小涡。没有人推。流体自己流起来。这个虚假速度场叫寄生电流（parasitic current）。Tallois 等（2025）用 well-balanced 表面张力把这股流动关掉。

把 curvature error ε 降到 0，速度场就消失。把同一个滑块抬高，界面上的涡又回来。

ε 为 0 的那一帧就是平衡成立的状态。换成 1D · face，同样的误差会把流动贴到网格轴上。multiD · nodal 把同样的误差散成弱得多的场。右侧的动能就是这个差别的数字。

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

弯曲的界面（两种流体相遇的边界）造出压力跳跃。这是 Laplace 法则。

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

二维液滴是 $\kappa=1/R$。$R$ 越小，跳跃越大。直径 1 mm，也就是 $R=0.5\,\mathrm{mm}$，按球来读，$\kappa=2/R$，$\Delta P=\sigma\cdot 2/R$。代入水的 $\sigma=0.072\,\mathrm{N/m}$，得到 $0.072\times 4000=288\,\mathrm{Pa}$，就是常写的约 300 Pa。代码是二维圆柱，$\kappa=1/R=2000\,\mathrm{m^{-1}}$，跳跃是 $0.072/0.0005=144\,\mathrm{Pa}$。球的一半。

内部压力向外推，表面张力向里拉。两者相等，液滴就不动。这个状态是平衡。

静止的液滴为何会流动

计算机在单元上凑不齐这个平衡。

数值方法把表面张力变成体积力。这就是 CSF（Continuum Surface Force）。界面被抹到几个单元上，力撒在这条带里。力的大小正比于曲率 $\kappa$。网格上的曲率很容易算歪。把同一股力推到壁面接触角上，是另一件收口的事。这里先坏的是更前面的一步，曲率和压力不再抵消。

曲率稍有偏差，表面张力就和压力梯度错开。剩下的力推动流体。

$\mathbf{r}$ 是残差力。平衡时它是 0。$\kappa_{\text{num}}$ 是数值曲率。$\nabla z$ 是体积分数（单元里液体所占的比例）的梯度。

$\kappa_{\text{num}}$ 精确时 $\mathbf{r}=0$。偏几个百分点，$\mathbf{r}$ 就不是 0。这个残差沿着界面造出旋涡。上图里的 ε 就是那几个百分点。

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

离散化若把静止解仍当成解，就是 well-balanced（平衡守恒）。

Tallois 等把表面张力当成非守恒乘积（non-conservative product），放进 Riemann 求解器里。压力梯度和表面张力用同一套离散规则来算，就会逐个单元抵消。本该什么都不发生的问题里，这个抵消一破，压力就会跳。那个判据在 非守恒项打破静止解的地方 另有一处。在液滴上，这个项就是曲率。

曲率模板也会把结果分开。基于面（face）的一维模板只看相邻的面。残差贴着网格轴，速度朝那个方向窜，界面跟着晃。基于节点（node）的多维模板看节点四周。残差散开，寄生电流变弱。表面张力是多维现象，模板也得是多维的。上图里 1D · face 和 multiD · nodal 的动能差，就是这个选择。

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

寄生电流是该除掉的虚假流动。表面张力也会造出真实的运动。

椭圆放进气体里会振动。它从表面能最大的椭圆出发。表面张力把它拉回圆形。那时动能最大。惯性冲过头，椭圆又回来。没有衰减，这个来回会继续。

周期遵循修正过的 Rayleigh 公式。

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

$n=2$ 时 $n^3-n=8-2=6$。用代码里的 $\sigma=0.072$、$R=5\times 10^{-4}$，$R^3=1.25\times 10^{-10}$，$6\sigma/R^3=3.456\times 10^{9}$。$\omega^2$ 是这个数除以 $\rho_l+\rho_g$。这篇没有把两个密度写成数字。分母变大，周期变长。$\sigma$ 变大，分子变大，周期变短。

理想的液滴无衰减地振动。扩散大的方法几个周期内就把振幅削成圆。论文用 low-Mach 修正压住这种数值衰减。

用 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

精确曲率的最大残差是 $4.1\times 10^{-10}$。动能 $2.3\times 10^{-19}$ 实际上是 0。3% 误差，也就是把 $\kappa$ 乘上 1.03 的那一侧，最大残差是 $4.8\times 10^{5}$，动能是 $1.7\times 10^{5}$。两者之比是 $1.7\times 10^{5}/2.3\times 10^{-19}\approx 7.4\times 10^{23}$，取对数大约 23.9 位。文中的 24 位就是这个数。

这个 $4.8\times 10^{5}$ 是压力跳跃的 3% 变成梯度之后的大小。$\sigma\cdot 0.03\cdot\kappa=0.072\times 60=4.32\,\mathrm{Pa}$，是留在界面上的压力误差。把 $dt=10^{-6}$ 加 200 次，时间间隔之和是 $2\times 10^{-4}$。再乘最大残差，速度尺度是 $2\times 10^{-4}\times 4.8\times 10^{5}=96$。3% 留下的力，把这个尺度的速度堆到界面上。

把 numerical damping 设为 0，椭圆的振幅保持住。把数值抬高，变成圆的速度就加快。增大 σ 或减小 ρ_l，周期 $T$ 变短。上面算出的 $6\sigma/R^3$ 是分子，所以把 $\sigma$ 加倍，分子加倍，$\omega$ 变成 $\sqrt{2}$ 倍。

半径 0.5 mm 落在采样区间的正中。点从 0.2 mm 排到 0.8 mm，一共 400 个。两端相距 0.6 mm，间隔 399 段，一段大约 1.50 μm。界面宽度是三段，大约 4.51 μm。半径大约是这个宽度的 111 倍，抹开的带相对液滴很薄。

用同一个表面张力配上球的曲率 2/R，压力跳跃是 288 Pa，其中 3% 是 8.64 Pa。代码里的圆柱曲率只有一半，留下的压力也是一半，4.32 Pa。把球的公式塞进这个二维残差，误差力会被读成两倍。

半径加倍时，立方变成 8 倍，角频率降到 8 的平方根，大约 2.83，周期按这个倍数变长。表面张力加倍只把频率乘上 2 的平方根，大约 1.41。半径对周期的推动更大。密度之和仍然没有数字。周期的秒数在这里闭不上。

留下那 3%，只把时间步放大，还是会坏#

静止液滴如果不安静，先把曲率对齐，再去量周期。残差 $4.8\times 10^{5}$ 还活着的时候，把 $\Delta t$ 加大，只是让速度尺度 96 堆得更快。well-balanced 把这个残差降到 $4.1\times 10^{-10}$，才有资格谈时间步。模板如果还留在面上，轴向流动会把这个资格拿回去。

残差死了以后，毛细时间步限制还在。$\Delta t$ 能推到几倍，这个残差式回答不了。把表面张力改成隐式之后，时间步实际能打开多宽 带着下一次失败。