界面只有两格宽时，表面张力少了39% — Korteweg应力的各向异性

表面张力不是跨越界面的压力差，而是法向与切向压力之差的积分；格子上的界面一旦薄于约四格，这个积分就会明显变小。

让压力计横穿界面

设想一个静止的平直气液界面。把压力计从气相一侧慢慢移到液相一侧。指针纹丝不动。平直界面没有Laplace压力差(曲率造成的压力跳跃)。可表面张力明明存在。那么它藏在哪里?

把压力计转90度，答案就出现了。沿着与界面平行的方向测量，指针只在界面内部下降。本文用Korteweg应力(由密度梯度产生的应力)推导这个差值从何而来。然后在格子间距为1的LBM格子上，实际测量这个积分随界面宽度丢失了多少。界面宽两格时，丢掉39%。

静止流体的应力不止一个压力

写多相LBM的控制方程时，动量方程中的应力分成三部分。

\boldsymbol{\sigma} = -p\,\mathbf{I} + \boldsymbol{\tau}_K + \boldsymbol{\tau}_{NS}, \qquad \boldsymbol{\tau}_{NS} = 2\mu\,\mathbf{S}^{d}

$p$ 是热力学压力， $\boldsymbol{\tau}_K$ 是Korteweg应力， $\boldsymbol{\tau}_{NS}$ 是Navier–Stokes黏性应力。 $\mathbf{S}^{d}$ 是无迹应变率张量(trace-free rate-of-strain，去掉体积变化后的纯剪切应变率)。

静止时 $\mathbf{S}^{d}=0$ ，黏性应力消失。只剩下 $p$ 和 $\boldsymbol{\tau}_K$ 。在均匀流体中 $\boldsymbol{\tau}_K$ 也为零。它只在有密度梯度的界面内部出现。而且这一项分方向。它在垂直于界面和平行于界面的方向上给出不同的值。

下面的模拟可以直接动手操作。黄色虚线横扫界面，显示该位置流体微元所受的应力。

T/T_c = 0.85 κ = 1.00 probe

z = 0.00 · ρ = 0.000 · p_N−p_b = 0.0e+0 · p_T−p_b = 0.0000 · σ = 0.5222 · σ/√κ = 0.5222 · width = 2.87

绿色箭头(法向分量)在界面各处都与虚线圆一样长。橙色箭头(切向分量)在界面中心变短。增大κ，界面变宽，橙色面积也就是σ按 $\sqrt{\kappa}$ 增长。

从自由能到Korteweg张量#

扩散界面模型从一个自由能泛函出发。

\Psi[\rho] = \int \left( \psi(\rho) + \frac{\kappa}{2}\,|\nabla\rho|^2 \right) dV

$\psi(\rho)$ 是均匀流体的自由能密度， $\kappa$ 是梯度能系数。对这个泛函变分，得到压力张量。

\mathbf{P} = \left( p_0 - \kappa\rho\nabla^2\rho - \frac{\kappa}{2}|\nabla\rho|^2 \right)\mathbf{I} + \kappa\,\nabla\rho\otimes\nabla\rho

$p_0=\rho\psi'-\psi$ 是状态方程压力。最后一项 $\kappa\nabla\rho\otimes\nabla\rho$ 是各向异性的来源。这个张量的散度满足 $\nabla\cdot\mathbf{P}=\rho\nabla\mu$ ，其中 $\mu=\psi'(\rho)-\kappa\nabla^2\rho$ 是化学势。在LBM里把这两种形式作为力加入时的区别，见非理想LBM forcing的两种形式。

本文使用约化单位下的van der Waals状态方程。

p_0(\rho) = \frac{8T\rho}{3-\rho} - 3\rho^2

$\rho$ 、 $T$ 、 $p$ 都除以了各自的临界值。 $T<1$ 时，气相密度 $\rho_g$ 与液相密度 $\rho_l$ 共存。二者由Maxwell条件(两相的 $\mu$ 和 $p$ 相等)确定。

法向分量是平的，切向分量下凹

取界面法向为 $z$ 轴。密度只依赖于 $z$ 。张量的两个分量为:

p_N = p_0 - \kappa\rho\rho'' + \frac{\kappa}{2}\rho'^2, \qquad p_T = p_0 - \kappa\rho\rho'' - \frac{\kappa}{2}\rho'^2

$p_N=P_{zz}$ 是法向分量， $p_T=P_{xx}$ 是切向分量。静力平衡要求 $dP_{zz}/dz=0$ 。所以 $p_N$ 在整个界面上都等于体相压力 $p_b$ 。这就是压力计不动的原因。

两者之差为 $p_N-p_T=\kappa\rho'^2\ge 0$ 。切向压力只在界面内部更低。压力更低，意味着这一层在横向拉扯。这就是表面张力的力学定义。

\sigma = \int_{-\infty}^{\infty} (p_N - p_T)\,dz = \kappa\int \rho'^2\,dz

把平衡方程 $\kappa\rho''=\psi'(\rho)-\mu_b$ 乘以 $\rho'$ 再积分，得到首次积分: $\frac{\kappa}{2}\rho'^2 = W(\rho) \equiv \psi(\rho)-\mu_b\rho+p_b$ 。代入后，空间积分变成密度积分。

\sigma = \int_{\rho_g}^{\rho_l} \sqrt{2\kappa\,W(\rho)}\;d\rho

$W(\rho)$ 是两个体相状态之间的自由能山丘。由此式立刻看出 $\sigma\propto\sqrt{\kappa}$ 。这是与格子无关的参考值。

代码:dz=1格子上的van der Waals界面#

LBM把格子间距取为1。于是界面宽度只由κ和 $T$ 决定。下面的代码用梯度流(gradient flow，每次沿残差方向下降一小步的迭代)把化学势残差压到零，得到平衡界面。然后在格子上计算 $p_N$ 和 $p_T$ ，比较应力积分得到的σ和首次积分公式给出的σ。

import math
 
# 约化van der Waals: rho_c = T_c = p_c = 1
def p0(r, T):  return 8*T*r/(3 - r) - 3*r*r
def mu0(r, T): return 8*T/3*(math.log(r/(3 - r)) + r/(3 - r)) - 6*r
def psi(r, T): return r*mu0(r, T) - p0(r, T)
 
def maxwell(T):
    # 化学势与压力相等 (Newton法，数值雅可比)
    rg, rl = 1 - 2*math.sqrt(1 - T), 1 + 2*math.sqrt(1 - T)
    h = 1e-7
    for _ in range(50):
        f1 = mu0(rl, T) - mu0(rg, T)
        f2 = p0(rl, T) - p0(rg, T)
        a = (mu0(rl+h, T) - mu0(rl, T))/h; b = -(mu0(rg+h, T) - mu0(rg, T))/h
        c = (p0(rl+h, T) - p0(rl, T))/h;   d = -(p0(rg+h, T) - p0(rg, T))/h
        det = a*d - b*c
        rl -= ( d*f1 - b*f2)/det
        rg -= (-c*f1 + a*f2)/det
    return rg, rl
 
def korteweg_profile(T, kappa, n=100, dz=1.0, tol=1e-7):
    rg, rl = maxwell(T)
    mub, pb = mu0(rl, T), p0(rl, T)
    L = n*dz                                         # 格子单位: 与LB格子一样 dz = 1
    z = [(i + 0.5)*dz - L/2 for i in range(n)]
    r = [0.5*(rl + rg) + 0.5*(rl - rg)*math.tanh(zi*math.sqrt(1 - T)/math.sqrt(kappa)) for zi in z]
    dt = 0.2*min(dz*dz/kappa, 1.0)
    for it in range(1, 30001):
        e = [r[0]] + r + [r[-1]]                     # 两端零梯度
        res = [mu0(r[i], T) - kappa*(e[i+2] - 2*r[i] + e[i])/dz**2 - mub for i in range(n)]
        r = [r[i] - dt*res[i] for i in range(n)]
        if max(abs(x) for x in res) < tol:
            break
    e = [r[0]] + r + [r[-1]]
    pN, pT = [], []
    for i in range(n):
        d1 = (e[i+2] - e[i])/(2*dz)
        lap = (e[i+2] - 2*r[i] + e[i])/dz**2
        iso = p0(r[i], T) - kappa*r[i]*lap
        pN.append(iso + 0.5*kappa*d1*d1)             # Korteweg压力张量的法向分量
        pT.append(iso - 0.5*kappa*d1*d1)             # 切向分量
    sigma_stress = sum(a - b for a, b in zip(pN, pT))*dz
    # 首次积分 kappa/2 rho'^2 = W(rho)  ->  sigma = int sqrt(2 kappa W) d rho
    m = 4000
    sigma_exact = 0.0
    for k in range(m):
        rr = rg + (k + 0.5)*(rl - rg)/m
        W = max(psi(rr, T) - mub*rr + pb, 0.0)
        sigma_exact += math.sqrt(2*kappa*W)*(rl - rg)/m
    w = (rl - rg)/max(abs(e[i+2] - e[i])/(2*dz) for i in range(n))
    return dict(w=w, it=it, rg=rg, rl=rl, spread=max(abs(p - pb) for p in pN),
                dip=min(pT) - pb, s_st=sigma_stress, s_ex=sigma_exact)
 
print("  T   kappa  width  iters  max|pN-pb|  min(pT-pb)  sigma_stress  sigma_exact  ratio")
for T in (0.85, 0.95):
    for kappa in (0.1, 0.5, 2.0, 8.0):
        d = korteweg_profile(T, kappa)
        print(f"{T:.2f}  {kappa:4.1f}  {d['w']:5.2f}  {d['it']:5d}   {d['spread']:.2e}   {d['dip']:+.4f}"
              f"     {d['s_st']:.5f}       {d['s_ex']:.5f}    {d['s_st']/d['s_ex']:.4f}")

输出如下。

  T   kappa  width  iters  max|pN-pb|  min(pT-pb)  sigma_stress  sigma_exact  ratio
0.85   0.1   2.08     47   2.09e-02   -0.0309     0.10067       0.16545    0.6085
0.85   0.5   2.63    275   1.08e-02   -0.1550     0.32864       0.36995    0.8883
0.85   2.0   4.35  30000   9.00e-03   -0.2429     0.71074       0.73990    0.9606
0.85   8.0   8.21    414   2.52e-03   -0.2650     1.46491       1.47980    0.9899
0.95   0.1   2.39    254   2.67e-03   -0.0112     0.02759       0.03239    0.8518
0.95   0.5   3.99  30000   9.96e-04   -0.0255     0.06887       0.07243    0.9508
0.95   2.0   7.38    290   3.61e-04   -0.0290     0.14301       0.14485    0.9873
0.95   8.0  14.48   1168   9.36e-05   -0.0298     0.28878       0.28971    0.9968

在浏览器里直接运行上面的代码。默认值对应输出第二行( $T=0.85$ ， $\kappa=0.5$ )。

T/T_c = 0.85 κ = 0.5 iters/frame 20

iter 0 · residual 1.00e+0 · width 0.00 cells · σ_stress 0.00000 · σ_exact 0.00000 · ratio 0.0000 · max|p_N−p_b| 0.00e+0

把κ滑块降到0.1，右下角的点会落到比值0.61附近。升到8，点贴在0.99。每跑完一次就留下一个点，可以改变 $T$ ，把输出表里没有的组合补齐。

两格宽的界面丢掉了39%的σ#

关键是ratio一列。界面宽度(密度差除以最大斜率)为2.08格时，比值是0.6085。格子上丢掉了39%的表面张力。2.63格时为0.888，4.35格时为0.961，8.21格时为0.990。比值几乎只由宽度决定，与温度无关。 $T=0.95$ 下3.99格宽给出0.951，与 $T=0.85$ 下4.35格的结果接近。

min(pT-pb)一列说的是同一件事。在 $T=0.85$ 下，κ越大，切向压力的凹坑从-0.031加深到-0.265。它在逼近连续介质的值。界面太薄时，坑底藏在了格点之间。

max|pN-pb|一列显示法向压力离"平"有多远。在连续介质中它为零。宽度2.08格时是 $2.1\times10^{-2}$ ，达到切向凹坑本身的三分之二。这个残留梯度会作为虚假的力作用在静止界面上。它是伪速度(spurious current，本应静止的界面周围出现的虚假流动)的来源之一。

iters为30000的两行并不是没有收敛。格子上的界面在平移下不对称。因此剖面会非常缓慢地滑向与格点咬合的位置。残差停在 $10^{-7}$ 附近的这段时间里，σ前五位有效数字都没有变化。

√κ与(1−T)^1.5#

看sigma_exact一列:κ每增大四倍，σ就翻一倍，0.165 → 0.370 → 0.740 → 1.480。这正是首次积分公式中的 $\sqrt{\kappa}$ 。宽度同样按 $\sqrt{\kappa}$ 缩放。增大κ换来更厚的界面，同时表面张力也随之变大。

温度依赖遵循平均场临界指数。κ=0.5时， $T=0.85$ 与 $T=0.95$ 的σ之比为0.36995/0.07243 = 5.11。对比 $(0.15/0.05)^{3/2}=5.20$ ，二者接近。越靠近临界点，两相密度差越小，界面越宽。所以同样的κ下， $T=0.95$ 那几行的比值更好。

实际工作中两者相互牵制。确定目标σ和密度比，就确定了κ和 $T$ ，宽度也随之确定。宽度低于四格时，σ会少5%以上。守恒型相场模型也有同样的宽度限制。守恒型Allen–Cahn的曲率项中的界面厚度W，也出于同样原因建议取四到五格以上。

定下κ之前先量一次宽度

在多相LBM中选好κ之后，正式计算之前先用1D松弛一个平直界面。检查两个数。一个是格子上的 $\kappa\sum\rho'^2$ 与 $\int\sqrt{2\kappa W}\,d\rho$ 之比。另一个是 $\max|p_N-p_b|$ 。比值低于0.95时，Laplace测试得到的σ也会偏低同样的比例。这时应增大κ，或把 $T$ 往临界点方向调高。如果法向压力不平，在静止液滴测试里应首先怀疑伪速度。