圧力をそのまま入れたら静止した界面が震えた

圧力形と自由エネルギー形のforcingはギブス・デュエムにより連続では等しい。しかし格子上では圧力形が界面に幽霊のような力を残す。

静止した液滴がひとりでに流れ始めた

ファンデルワールス状態方程式を格子ボルツマン(LBM, Lattice Boltzmann Method)に載せて二相流体を解いた。初期条件は静止した液滴が一つ。密度場は滑らかで、速度はすべて0です。

数ステップ回すと界面付近で小さな速度が芽生えた。誰も押していないのに流体が流れる。このような人工的な速度を寄生流(parasitic current、物理的原因なしに界面で生じる幽霊のような流れ)と呼びます。

原因を絞り込んでいくとコードの一箇所に行き着いた。力の項に何を入れるかです。圧力 $p$ をそのまま入れるのか、化学ポテンシャル $\mu$ を入れるのか。教科書は両者が同じだと言う。この記事はその「同じ」がどこまで真なのかを帳簿に記す。答えは一行です — 連続では同じで、格子上では分かれる。

temperature を1より下げると等温線が折れます。一つの圧力に密度が三つ付く区間が開きます。

temperature0.900 Tc

Drag temperature down from 1.0: the isotherm folds into an S, and the amber tie line drops in so the two green lobes have equal area. The blue and pink dots are the vapor and liquid densities the LBM must hold apart — their gap is what the forcing term exists to sustain.

横線は二つの緑のローブの面積が等しくなるマクスウェル圧力です。青い点とピンクの点の間隔が、力の項が後で支えなければならない密度差です。

力はどこから入ってくるのか

LBMは分布関数を移流・衝突させて巨視方程式を復元します。理想気体なら自然に出てくる圧力は $c_s^2 \rho$ 一つだけです( $c_s$ は格子音速)。ファンデルワールスのような非理想流体の本当の圧力はそれとは違う。その差額を埋めるのが力の項 $\mathbf{F}$ の役割です。

目標は次の運動量方程式を復元することです。

\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u}\mathbf{u}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \kappa \rho \, \nabla (\nabla^2 \rho)

$p$ はファンデルワールス圧力、 $\boldsymbol{\tau}$ は粘性応力、最後の項は界面を立てるコルテヴェグ(Korteweg) 応力で、 $\kappa$ はその強さです。ストリーミングが $c_s^2\rho$ を与えるので、力は残りの差額だけを埋めればよい。

ここで二つの分岐が生まれます。圧力をそのまま使う形と、自由エネルギーから来た化学ポテンシャルを使う形です。

\mathbf{F}_p = -\nabla\!\left(p - c_s^2 \rho\right) + \kappa \rho \, \nabla(\nabla^2 \rho)

\mathbf{F}_\mu = -\rho \, \nabla \mu + c_s^2 \nabla \rho + \kappa \rho \, \nabla(\nabla^2 \rho)

前者が圧力形(今回の題材「 $p$ をそのまま」)、後者が自由エネルギー形です。等温線が折れる温度は前のダイヤルで既に下げました。ここからは、二つの力が連続で一致する理由だけを帳簿に残します。

ギブス・デュエム — 圧力と化学ポテンシャルは同じことを二度言う

二つの形が同じかどうかは一つの関係式で決着がつく。等温で圧力と化学ポテンシャルを結ぶギブス・デュエム(Gibbs–Duhem) 関係です。

\mathrm{d}p = \rho \, \mathrm{d}\mu \qquad (T = \text{const})

この式を $\mathbf{F}_\mu$ に代入してみましょう。 $-\rho\nabla\mu = -\nabla p$ なので化学ポテンシャル項がただちに圧力勾配に変わる。残る $c_s^2\nabla\rho$ は圧力形が $-\nabla(-c_s^2\rho)$ として持っていた項と正確に噛み合う。結局 $\mathbf{F}_p = \mathbf{F}_\mu$ です。

つまり「 $p$ をそのまま入れてよい」という主張の唯一の根拠はギブス・デュエムです。ファンデルワールスはこの関係を正確に満たす。状態方程式と自由エネルギーが同じ熱力学から出てきたからです。Chapman–Enskog展開で確認した LBM forcingスキームたちがそれぞれ異なる顔をしていても結局同じ巨視方程式を復元するのと同じ種類の等価です。

Pythonで確認した共存密度とギブス・デュエム#

言葉だけでは心もとない。ファンデルワールス状態方程式を換算単位で組み、温度別の共存密度をニュートン法で求めたうえで、ギブス・デュエム欠陥 $\lvert \mathrm{d}p/\mathrm{d}\rho - \rho\,\mathrm{d}\mu/\mathrm{d}\rho \rvert$ を直接測ってみましょう。

import numpy as np
 
A, B, R = 9.0 / 8.0, 1.0 / 3.0, 1.0   # ファンデルワールス、換算単位 (rho_c=1, T_c=1)
 
def p_eos(rho, T):                     # ファンデルワールス圧力
    return rho * R * T / (1.0 - B * rho) - A * rho * rho
 
def mu_eos(rho, T):                    # 化学ポテンシャル mu = df/drho
    return R * T * (np.log(rho / (1.0 - B * rho)) + B * rho / (1.0 - B * rho)) - 2.0 * A * rho
 
def maxwell(T):                        # 等面積則: p と mu が二相で等しくなる密度
    x = np.array([0.30, 1.90]); h = 1e-8
    for _ in range(80):
        f = np.array([p_eos(x[0], T) - p_eos(x[1], T), mu_eos(x[0], T) - mu_eos(x[1], T)])
        J = np.empty((2, 2))
        for k in range(2):
            y = x.copy(); y[k] += h
            g = np.array([p_eos(y[0], T) - p_eos(y[1], T), mu_eos(y[0], T) - mu_eos(y[1], T)])
            J[:, k] = (g - f) / h
        x -= np.linalg.solve(J, f)
    return x[0], x[1]
 
# ギブス・デュエム: dp = rho d(mu). 圧力をそのまま使ってよい唯一の根拠。
print("T/Tc  rho_vap  rho_liq |  max| dp/drho - rho*dmu/drho |")
for T in (0.95, 0.90, 0.85):
    rv, rl = maxwell(T)
    r = np.linspace(rv, rl, 400)
    dp = np.gradient(p_eos(r, T), r)
    dmu = np.gradient(mu_eos(r, T), r)
    err = np.abs(dp - r * dmu).max()
    print(f"{T:4.2f}  {rv:7.4f}  {rl:7.4f} |  {err:.2e}")

T/Tc  rho_vap  rho_liq |  max| dp/drho - rho*dmu/drho |
0.95   0.5790   1.4617 |  2.94e-04
0.90   0.4257   1.6573 |  9.44e-04
0.85   0.3197   1.8071 |  1.98e-03

欠陥は $10^{-3}$ 程度で、それすら有限差分で微分を測ったせいです。温度を下げるほど共存密度の間隔が広がる。ギブス・デュエムは連続では正確に成立します。ここまでは圧力形と自由エネルギー形が完全に同じです。

共存密度を列ごとに引いてみます。 $T/T_c=0.95$ は気体 $0.5790$ 、液体 $1.4617$ 、差 $0.8827$ です。 $0.90$ は $1.6573-0.4257=1.2316$ 、 $0.85$ は $1.8071-0.3197=1.4874$ です。温度を $0.10$ 下げる間に密度差は $0.8827$ から $1.4874$ へ、約 $1.7$ 倍になります。同じ区間で欠陥は $2.94\times10^{-4}$ 、 $9.44\times10^{-4}$ 、 $1.98\times10^{-3}$ と増えます。恒等式 $\mathrm{d}p=\rho\,\mathrm{d}\mu$ が壊れたのではありません。400点の np.gradient が、より急な等温線で差分誤差を多く出しただけです。 $T=0.90$ の液体/気体比は $1.6573/0.4257\approx 3.89$ です。この比が forcing の支える密度ジャンプの大きさです。

力の二行を重ねると、差が出る場所は一つです。 $\mathbf{F}_p$ の $-\nabla(p-c_s^2\rho)$ と $\mathbf{F}_\mu$ の $-\rho\nabla\mu+c_s^2\nabla\rho$ は、 $\mathrm{d}p=\rho\,\mathrm{d}\mu$ の下で一致します。コルテヴェグ項 $\kappa\rho\nabla(\nabla^2\rho)$ は両側に同じ形なので差を作りません。 $c_s^2=1/3$ は格子音速が既に入れてある補正であり、界面を押す物理的な力ではありません。連続の帳簿で二つの形が同じだという文はここで終わります。格子上でその文が崩れる大きさは、後の表の $1.60\times10^{-2}$ です。

離散格子で両者が分かれる

問題は格子です。連続の $\mathrm{d}p = \rho\,\mathrm{d}\mu$ が離散微分でも成立する保証はありません。中心差分した $\nabla p$ と $\rho\,\nabla\mu$ は界面のように密度が急に折れる場所で互いにずれる。

静止した平面界面をオイラー・ラグランジュ条件で解き、その上で二つの力の形を直接計算しました。

import numpy as np
 
A, B, R = 9.0 / 8.0, 1.0 / 3.0, 1.0
CS2, KAPPA, T = 1.0 / 3.0, 0.02, 0.90
 
def p_eos(rho):  return rho * R * T / (1.0 - B * rho) - A * rho * rho
def mu_eos(rho): return R * T * (np.log(rho / (1.0 - B * rho)) + B * rho / (1.0 - B * rho)) - 2.0 * A * rho
def dmu(rho):    return R * T * (1.0 / (rho * (1.0 - B * rho)) + B / (1.0 - B * rho) ** 2) - 2.0 * A
def diff1(a, dx): return (np.roll(a, -1) - np.roll(a, 1)) / (2.0 * dx)
def lap(a, dx):   return (np.roll(a, -1) - 2.0 * a + np.roll(a, 1)) / dx ** 2
 
def maxwell():
    x = np.array([0.30, 1.90]); h = 1e-8
    for _ in range(80):
        f = np.array([p_eos(x[0]) - p_eos(x[1]), mu_eos(x[0]) - mu_eos(x[1])])
        J = np.empty((2, 2))
        for k in range(2):
            y = x.copy(); y[k] += h
            g = np.array([p_eos(y[0]) - p_eos(y[1]), mu_eos(y[0]) - mu_eos(y[1])])
            J[:, k] = (g - f) / h
        x -= np.linalg.solve(J, f)
    return x[0], x[1]
 
rv, rl = maxwell()
mu_co = mu_eos(np.array([rv]))[0]
 
# (1) 静止平面界面を解き、二つのforcingの形を比較する
NX = 240
xs = np.arange(NX)
rho = 0.5 * (rl + rv) + 0.5 * (rl - rv) * (np.tanh((xs - NX / 4) / 6.0) - np.tanh((xs - 3 * NX / 4) / 6.0) - 1.0)
for _ in range(6000):                                    # オイラー・ラグランジュ残差の緩和
    rho -= 0.15 * (mu_eos(rho) - KAPPA * lap(rho, 1.0) - mu_co) / dmu(rho)
 
Gp  = -diff1(p_eos(rho) - CS2 * rho, 1.0) + KAPPA * rho * diff1(lap(rho, 1.0), 1.0) - diff1(CS2 * rho, 1.0)
Gmu = -rho * diff1(mu_eos(rho) - KAPPA * lap(rho, 1.0), 1.0) + CS2 * diff1(rho, 1.0) - diff1(CS2 * rho, 1.0)
gd  = np.abs(diff1(p_eos(rho), 1.0) - rho * diff1(mu_eos(rho), 1.0)).max()
 
print(f"coexistence         rho_vap = {rv:.4f}   rho_liq = {rl:.4f}")
print(f"free-energy form    max|G_mu|        = {np.abs(Gmu).max():.2e}   (well-balanced)")
print(f"pressure form       max|G_p|         = {np.abs(Gp).max():.2e}   (spurious force)")
print(f"gap between forms   max|G_p - G_mu|  = {np.abs(Gp - Gmu).max():.2e}")
print(f"discrete Gibbs-Duhem defect          = {gd:.2e}   <- the gap, exactly")
 
# (2) 同じ物理界面をdxだけ細かくすると欠陥は速やかに0へ収束する
print("\ncells/interface |  Gibbs-Duhem defect  order")
prev = None
for n in (10, 20, 40, 80):
    L = 40.0; N = int(L * n / 10)
    z = np.linspace(-L / 2, L / 2, N, endpoint=False); dx = z[1] - z[0]
    r = 0.5 * (rl + rv) - 0.5 * (rl - rv) * np.tanh(z / (0.1 * n))
    d = np.abs(diff1(p_eos(r), dx) - r * diff1(mu_eos(r), dx))[N // 4:3 * N // 4].max()
    order = "" if prev is None else f"{np.log(prev / d) / np.log(2.0):5.2f}"
    print(f"{n:9d}       |  {d:.3e}          {order}")
    prev = d

coexistence         rho_vap = 0.4257   rho_liq = 1.6573
free-energy form    max|G_mu|        = 2.02e-16   (well-balanced)
pressure form       max|G_p|         = 1.60e-02   (spurious force)
gap between forms   max|G_p - G_mu|  = 1.60e-02
discrete Gibbs-Duhem defect          = 1.60e-02   <- the gap, exactly
 
cells/interface |  Gibbs-Duhem defect  order
       10       |  1.140e-02          
       20       |  1.428e-03           3.00
       40       |  5.115e-05           4.80
       80       |  1.611e-06           4.99

resolution を上げると界面がまたがるセルが増え、赤い峰の高さが一緒に死にます。

resolution (cells / interface)2.5

The red hump is the force the pressure form adds over the free-energy form — the discrete Gibbs–Duhem defect. It lives exactly on the interface and pushes a fluid that should be at rest. Drag resolution up: the hump collapses onto the green zero line. It was never physics, only the price of writing p on too few cells.

赤い峰は圧力形が自由エネルギー形の上に余分に載せる力、つまり $G_p - G_\mu$ です。緑の点線は0です。セルを増やすと峰は沈みます。幽霊のような力は物理ではなく、この差分の余りです。

自由エネルギー形は界面で $2.02\times10^{-16}$ です。圧力形は $1.60\times10^{-2}$ を残します。二つの差は離散ギブス・デュエム欠陥と小数点まで一致します。 $1.60\times10^{-2}/2.02\times10^{-16}\approx 7.9\times10^{13}$ です。機械精度と幽霊のような力の桁が同じ行に並びます。

格子の表も同じ帳簿です。界面あたり10セルで欠陥は $1.140\times10^{-2}$ 、20セルで $1.428\times10^{-3}$ 、観測次数は $3.00$ です。40セルは $5.115\times10^{-5}$ (次数 $4.80$ )、80セルは $1.611\times10^{-6}$ (次数 $4.99$ )です。10セルから80セルへ8倍細かくするだけで、 $1.140\times10^{-2}/1.611\times10^{-6}\approx 7.1\times10^{3}$ 倍減ります。中心差分が界面勾配で $\mathrm{d}p=\rho\,\mathrm{d}\mu$ を取りこぼした大きさが、二つの力の差そのものです。この力が静止した界面を押します。

圧力形を残すと薄い界面が代金を請求する

自由エネルギー形( $-\rho\nabla\mu$ )は定義上、界面で均衡が取れます。粗い格子でも幽霊のような力は $2.02\times10^{-16}$ に張りつきます。圧力をそのまま使うなら、 $\nabla p$ を勝手に差分せず $\rho\,\nabla\mu$ と同じ差分で組む必要があります。そうして初めて離散ギブス・デュエムが生き、表の $1.60\times10^{-2}$ が下がります。

界面を4〜5セル以上に解けると、圧力形の欠陥は次数 $4.99$ まで落ちます。実務の界面は3セル前後です。その厚さでは幽霊のような力が密度差の二乗に比例して残ります。力の項が連続の恒等式だけを満たし離散モーメントを空にすると、同じ残差が残ります。その帳簿は力項の3次モーメント限界にあります。「 $p$ をそのまま」の代償は界面厚さで請求されます。

圧力をそのまま入れたら静止した界面が震えた — 非理想LBM forcingの二つの形