A Two-Cell Interface Loses 39% of Its Surface Tension — The Anisotropy of the Korteweg Stress

Surface tension is not a pressure jump across the interface but the integral of the normal-minus-tangential pressure, and on a lattice that integral shrinks visibly once the interface is thinner than about four cells.

Sliding a pressure gauge across the interface#

Take a flat liquid–vapor interface at rest. Move a pressure gauge slowly from the vapor side into the liquid. The needle does not move. A flat interface has no Laplace pressure jump (the jump that curvature creates). Yet the surface tension is clearly there. So where is it hiding?

Rotate the gauge by 90 degrees and you see it. Measured parallel to the interface, the needle drops, but only inside the interface. This post derives that difference from the Korteweg stress (the stress produced by density gradients). Then it measures, on an LBM lattice with unit spacing, how much of that integral disappears as the interface gets thinner. With a two-cell interface, 39% is gone.

The stress in a fluid at rest is more than one pressure#

When you write the momentum equation of a multiphase LBM, the stress splits into three pieces.

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

$p$ is the thermodynamic pressure, $\boldsymbol{\tau}_K$ the Korteweg stress, and $\boldsymbol{\tau}_{NS}$ the Navier–Stokes viscous stress. $\mathbf{S}^{d}$ is the trace-free rate-of-strain tensor (pure shear strain with volume change removed).

At rest $\mathbf{S}^{d}=0$ , so the viscous stress vanishes. Only $p$ and $\boldsymbol{\tau}_K$ remain. In a uniform fluid $\boldsymbol{\tau}_K$ is zero too. It comes alive only inside the interface, where the density has a gradient. And it is directional. It gives one value normal to the interface and another value parallel to it.

Try it yourself in the simulation below. The yellow dashed line sweeps across the interface and shows the stress acting on a fluid element at that position.

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

The green arrows (normal component) match the dashed circle everywhere in the interface. The orange arrows (tangential component) shrink at the center. Raise κ and the interface widens, and the orange area, which is σ, grows like $\sqrt{\kappa}$ .

From the free energy to the Korteweg tensor#

A diffuse-interface model starts from a single free-energy functional.

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

$\psi(\rho)$ is the free-energy density of a uniform fluid, and $\kappa$ is the gradient-energy coefficient. Varying this functional yields the pressure tensor.

\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$ is the equation-of-state pressure. The last term, $\kappa\nabla\rho\otimes\nabla\rho$ , is the source of the anisotropy. The divergence of this tensor satisfies $\nabla\cdot\mathbf{P}=\rho\nabla\mu$ , where $\mu=\psi'(\rho)-\kappa\nabla^2\rho$ is the chemical potential. How the two forms differ when you feed them into LBM as a force is covered in the two forms of non-ideal LBM forcing.

This post uses the van der Waals equation of state in reduced units.

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

$\rho$ , $T$ , and $p$ are all divided by their critical values. For $T<1$ , a vapor density $\rho_g$ and a liquid density $\rho_l$ coexist. The Maxwell conditions (equal $\mu$ and equal $p$ in both phases) fix them.

The normal component is flat, the tangential one dips#

Put the interface normal along $z$ . The density depends on $z$ only. The two components of the tensor are:

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}$ is the normal component and $p_T=P_{xx}$ the tangential one. Static equilibrium requires $dP_{zz}/dz=0$ . So $p_N$ equals the bulk pressure $p_b$ through the whole interface. That is why the gauge did not move.

Their difference is $p_N-p_T=\kappa\rho'^2\ge 0$ . The tangential pressure is lower, and only inside the interface. A lower pressure means the layer is pulling sideways. This is the mechanical definition of surface tension.

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

Multiply the equilibrium equation $\kappa\rho''=\psi'(\rho)-\mu_b$ by $\rho'$ and integrate to get a first integral: $\frac{\kappa}{2}\rho'^2 = W(\rho) \equiv \psi(\rho)-\mu_b\rho+p_b$ . Substituting turns the spatial integral into a density integral.

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

$W(\rho)$ is the free-energy hill between the two bulk states. The formula shows $\sigma\propto\sqrt{\kappa}$ at a glance. It is a reference value that knows nothing about the lattice.

Code: a van der Waals interface on a dz=1 lattice#

LBM sets the lattice spacing to 1. The interface width is then set only by κ and $T$ . The code below builds the equilibrium interface by gradient flow (stepping a little against the chemical-potential residual each iteration). It then computes $p_N$ and $p_T$ on the lattice and compares the stress integral σ with the first-integral σ.

import math
 
# reduced 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):
    # equal chemical potential and equal pressure (Newton, numerical Jacobian)
    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                                         # lattice units: dz = 1 like an LB grid
    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]]                     # zero-gradient ends
        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)             # normal component of Korteweg pressure tensor
        pT.append(iso - 0.5*kappa*d1*d1)             # tangential component
    sigma_stress = sum(a - b for a, b in zip(pN, pT))*dz
    # first integral 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}")

The output:

  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

Now run the same code in your browser. The default is the second row of the output ( $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

Drop κ to 0.1 and the dot in the lower right falls near a ratio of 0.61. Push it to 8 and the dot sits at 0.99. Every finished run leaves a dot, so change $T$ and fill in combinations the output table does not have.

A two-cell interface loses 39% of σ#

The ratio column is the key. With an interface width (density difference over maximum slope) of 2.08 cells, it is 0.6085. The lattice has lost 39% of the surface tension. At 2.63 cells it is 0.888, at 4.35 cells 0.961, and at 8.21 cells 0.990. The ratio is set almost entirely by the width, regardless of temperature. At $T=0.95$ , a width of 3.99 cells gives 0.951, close to the 4.35-cell case at $T=0.85$ .

The min(pT-pb) column tells the same story. At $T=0.85$ , the tangential dip deepens from -0.031 to -0.265 as κ grows. It is converging to the continuum value. In a thin interface, the bottom of the dip hides between lattice nodes.

The max|pN-pb| column shows how far the normal pressure is from flat. In the continuum it would be zero. At a width of 2.08 cells it is $2.1\times10^{-2}$ , two-thirds of the tangential dip itself. That leftover gradient acts as a fake force on a static interface. It is one source of spurious currents (artificial flow around an interface that should be at rest).

The two rows with iters at 30000 are not failures to converge. A lattice interface is not symmetric under translation. The profile slowly slides toward a position locked to the nodes. While the residual stalls near $10^{-7}$ , σ does not change in the first five digits.

√κ and (1−T)^1.5#

In the sigma_exact column, each fourfold increase of κ doubles σ: 0.165 → 0.370 → 0.740 → 1.480. That is the $\sqrt{\kappa}$ of the first-integral formula. The width scales with $\sqrt{\kappa}$ too. A larger κ buys a thicker interface, but it also raises the surface tension.

The temperature dependence follows the mean-field critical exponent. At κ=0.5, the ratio of σ between $T=0.85$ and $T=0.95$ is 0.36995/0.07243 = 5.11. Compare $(0.15/0.05)^{3/2}=5.20$ . Closer to the critical point, the density difference shrinks and the interface widens. That is why the $T=0.95$ rows score better at the same κ.

In practice these two pull against each other. Pick a target σ and density ratio, and you have picked κ and $T$ , and therefore the width. Below four cells, σ drops by more than 5%. Conservative phase-field models hit the same width limit. The interface thickness W in the curvature term of conservative Allen–Cahn is kept at four to five cells or more for the same reason.

Measure the width before you commit to a κ#

Once you have chosen κ for a multiphase LBM run, relax a single flat interface in 1D before launching anything. Check two numbers. The first is the ratio of the lattice $\kappa\sum\rho'^2$ to $\int\sqrt{2\kappa W}\,d\rho$ . The second is $\max|p_N-p_b|$ . If the ratio is below 0.95, the σ from your Laplace test will come out low by the same amount. Raise κ, or move $T$ toward the critical point. If the normal pressure is not flat, suspect spurious currents first in your static-droplet test.