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

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

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

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

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

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

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

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$ をそのまま」)、後者が自由エネルギー形です。両者が本当に同じか見るには、まずファンデルワールス状態方程式の姿を知る必要があります。以下で温度を直接下げてみましょう。

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.

temperature を1.0より下げると等温線がS字に折れます。一つの圧力に密度が三つ対応する。このとき物理的な二相圧力は二つの緑のローブの面積が等しくなる位置(マクスウェルの等面積則)で決まります。青い点とピンクの点の間隔がまさにforcing項が支えなければならない密度差です。

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

二つの形が同じかどうかは一つの関係式で決着がつく。等温で圧力と化学ポテンシャルを結ぶギブス・デュエム(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}$ 程度で、それすら有限差分で微分を測ったせいです。温度を下げるほど共存密度の間隔が広がる。ギブス・デュエムは連続では正確に成立します。ここまでは圧力形と自由エネルギー形が完全に同じです。

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

問題は格子です。連続の $\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

三行が核心です。自由エネルギー形は界面で力が $10^{-16}$ 、事実上0です。圧力形は $1.6\times10^{-2}$ の力を残す。そして二つの形の差が離散ギブス・デュエム欠陥と小数点まで正確に一致する。圧力形がこぼした幽霊のような力の正体がまさにこの欠陥です。この力が静止した界面を押して寄生流を作ります。

以下で界面をどれだけ広く広げたか(格子解像度)を変えてみましょう。

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.

resolution を上げて界面をより多くのセルにまたがらせると、赤い圧力形の曲線の峰が0へ沈み込みます。緑の自由エネルギー形は最初から最後まで0に張りついている。幽霊のような力は物理ではなく離散化の副産物でした。

では何を入れるべきか

まとめると選択は二つです。第一に、自由エネルギー形( $-\rho\nabla\mu$ )を使う。この形は定義上界面で均衡が取れるので、粗い格子でも幽霊のような力がありません。第二に、本当に圧力をそのまま使いたいなら $\nabla p$ を自由に離散化せず $\rho\,\nabla\mu$ と一致するように組む。そうすれば離散でもギブス・デュエムが成立して均衡が生きる。

界面を4〜5セル以上に十分解けるなら圧力形の誤差は上の表のように速やかに消えます。しかし実務の界面はたいてい3セル前後で薄い。その領域で圧力形は密度差の二乗に比例する幽霊のような力を出す。同じ症状を寄生流とwell-balancedな界面張力で表面張力の側から見たことがあります。根は一つです — 連続で均衡する項を離散でも均衡に移せたか。

「 $p$ をそのまま使う」という便利さはタダではない。その代償は界面厚さという通貨で請求されます。

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