界面幅が2セルなら表面張力の39%が消える

表面張力は界面をまたぐ圧力差ではなく、法線方向と接線方向の圧力差を積分した値です。格子上の界面が4セルより薄くなると、その積分ははっきり小さくなります。

界面を横切って圧力計を動かしてみた

静止した平らな気液界面があります。気相側から液相側へ、圧力計をゆっくり動かしてみます。針は動きません。平らな界面にはLaplace圧力差(曲率が生む圧力の跳び)がないからです。それでも表面張力は確かに存在します。では、表面張力はどこに隠れているのでしょうか。

圧力計を90度回すと答えが見えます。界面に平行な向きで測ると、界面の中でだけ針が下がります。この記事では、その差がどこから来るのかをKorteweg応力(密度勾配が生む応力)から導きます。さらに、格子間隔1のLBM格子で、この積分が界面幅に応じてどれだけ失われるかを実際に測ります。界面幅が2セルなら39%が消えます。

静止流体の応力は圧力ひとつではない

多相LBMの支配方程式を書くと、運動量式の応力は3つに分かれます。

\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$ も0です。密度勾配のある界面の中でだけ現れます。そしてこの項は方向を区別します。界面に垂直な方向と平行な方向に、異なる値を与えます。

下のシミュレーションで実際に操作してみましょう。黄色の点線が界面を横切って動き、その位置の流体要素にかかる応力を表示します。

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でこの2つの形を力として入れたときの違いは、非理想LBM forcingの2つの形で扱いました。

この記事では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$ だけに依存します。テンソルの2成分は次のとおりです。

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$ に等しくなります。圧力計が動かなかった理由です。

2成分の差は $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)$ は2つのバルク状態の間にある自由エネルギーの山です。この式から $\sigma\propto\sqrt{\kappa}$ がすぐに分かります。格子とは無関係な基準値です。

コード:dz=1格子上のvan der Waals界面#

LBMは格子間隔を1とします。すると界面幅はκと $T$ だけで決まります。下のコードは、化学ポテンシャルの残差を0に近づける勾配流(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

上のコードをブラウザで直接動かしてみましょう。初期値は出力の2行目( $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に張り付きます。実行が終わるたびに点が1つ残るので、 $T$ を変えて出力表にない組み合わせを埋めてみてください。

2セルの界面ではσの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|列は、法線圧力がどれだけ平らでないかを示します。連続体なら0です。幅2.08セルでは $2.1\times10^{-2}$ で、接線のくぼみの3分の2にもなります。この残った勾配は、静止界面に偽の力として働きます。スプリアス流れ(spurious current、静止すべき界面の周りに生じる偽の流れ)の原因のひとつです。

itersが30000の2行は収束の失敗ではありません。格子上の界面は平行移動に対して対称ではありません。そのためプロファイルが格子点にかみ合う位置へ、ごくゆっくり滑っていきます。残差が $10^{-7}$ 付近で止まっている間も、σは5桁目まで変わりません。

√κと(1−T)^1.5#

sigma_exact列を見ると、κが4倍になるたびに0.165 → 0.370 → 0.740 → 1.480と2倍ずつ増えます。第一積分公式の $\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$ の側の比が良くなります。

実務ではこの2つがぶつかります。目標のσと密度比を決めるとκと $T$ が決まり、すると幅も決まります。幅が4セルを下回るとσは5%以上小さくなります。保存型フェーズフィールドモデルにも同じ幅の制約があります。保存型Allen–Cahnの曲率項で見た界面厚さWも、同じ理由で4~5セル以上が推奨されます。

κを決める前に一度測っておきたい幅

多相LBMでκを決めたら、本計算の前に平らな界面をひとつ1Dで緩和してみましょう。確認する数字は2つです。ひとつは格子上の $\kappa\sum\rho'^2$ と $\int\sqrt{2\kappa W}\,d\rho$ の比です。もうひとつは $\max|p_N-p_b|$ です。比が0.95を下回るなら、Laplaceテストのσもその分だけ小さく出ます。その場合はκを大きくするか、 $T$ を臨界点側に上げるほうがよいでしょう。法線圧力が平らでないなら、静止液滴テストではまずスプリアス流れを疑いましょう。

界面幅が2セルなら表面張力の39%が消える — Korteweg応力の異方性