にじんだ界面を毎ステップ元に戻す — 反拡散シャープニング(IST)を実装する

薄くするだけではセル数がまた揺れる。各ステップの後にε/aを固定しないと、衝撃の横で幅が戻ってしまう。

衝撃波–気泡問題を5方程式モデルで回してみました。200ステップ後、ヘリウム気泡の縁は6セルから20セルへ膨らんでいました。物理量は合っているのに、界面は粥のようにつぶれています。これが拡散界面法(界面を有限の幅でにじませておく手法)の宿命です。Nguyenら(2021)は、このにじみを毎タイムステップ後の後処理一回で元に戻します。今日はその反拡散手法を自分で実装し、なぜ界面の厚さが一定に固定されるのかを目で確かめます。

numerical diffusion を入れると6セルの縁がセルを飲み込み、sharpening (IST) を入れるとその広がりが ε/a の近くで止まります。

compression a: 1.00regularization ε: 0.0040

measured band: 0 cells

target ≈ ε/a: 2 cells

steady tanh thickness scales like ε/a

黄色い帯が目標厚さです。compression a を上げると帯は細くなり、regularization ε を上げると再び厚くなります。測定セル数がその比を追うかどうかだけ見れば十分です。

論文情報

タイトル: Numerical modeling of multiphase compressible flows with the presence of shock waves using an interface-sharpening five-equation model
著者: Van-Tu Nguyen, Thanh-Hoang Phan, Warn-Gyu Park (Pusan National University)
出典: International Journal of Multiphase Flow, 2021, 103542
DOI: 10.1016/j.ijmultiphaseflow.2020.103542
一行要約: 5方程式圧縮性二相モデルに界面シャープニング(IST)を後処理として付け、衝撃波の隣でも界面を一定の厚さに保つ。

なぜ界面はひとりでに太るのか

5方程式モデルは、保存則に色関数を一つ加えた準保存系です。体積分率 $\alpha_1$ (セル内で流体1が占める割合)の移流方程式はこうなります。

\frac{\partial \alpha_1}{\partial t} + \mathbf{u}\cdot\nabla \alpha_1 = \alpha_1 K\, \nabla\cdot\mathbf{u}

$\alpha_1$ は流体1の体積分率、 $\mathbf{u}$ は混合速度、右辺の Kapila 項は二相の圧縮率の差を補正します。

問題は、これを衝撃波捕捉スキームで解くときに起きます。セル境界で値を滑らかに補間した瞬間、鋭くあるべき $\alpha_1$ の階段が毎ステップ少しずつにじみます。この数値拡散は蓄積します。ステップが重なると界面は数セルから数十セルへ広がり、気泡の形や衝撃波の反射位置といった核心の特徴が丸ごと消えます。

WENO のような高次スキームは拡散を減らしますが、振動を生み、多次元では高価です。そこで著者らは方向を変えます。スキーム自体には触れず、毎ステップの後で $\alpha$ 場を別に手入れすることにします。

反拡散: シャープニング方程式

核心は Shukla ら(2010)の正則化方程式です。物理時間ではなく仮想時間 $\tau$ に対して $\alpha$ を数回だけ進化させます。

\frac{\partial \alpha}{\partial \tau} = \nabla\cdot\!\left(\varepsilon\, \nabla \alpha\right) - \nabla\cdot\!\left(a\, \alpha(1-\alpha)\, \hat{\mathbf{n}}\right)

$\varepsilon$ は厚さを規制する拡散係数、 $a$ は界面を締める圧縮の強さ、 $\hat{\mathbf{n}} = \nabla\alpha / |\nabla\alpha|$ は界面法線方向です。

二つの項が綱引きをします。第一項は普通の拡散なので $\alpha$ を広げます。第二項は符号が逆の**反拡散(anti-diffusion)**で、界面を法線方向に締めます。 $\alpha(1-\alpha)$ の因子のおかげで、圧縮は界面( $0<\alpha<1$ )でのみ働き、純流体領域( $\alpha=0$ または $1$ )では切れます。

定常状態( $\partial_\tau \alpha = 0$ )を一次元で解くと、二つのフラックスが釣り合います。

\varepsilon\, \frac{d\alpha}{dx} = a\, \alpha(1-\alpha) \quad\Longrightarrow\quad \alpha(x) = \tfrac{1}{2}\!\left[1 + \tanh\!\left(\frac{a}{2\varepsilon}\,x\right)\right]

界面は双曲線正接の形へ収束し、その厚さは $\varepsilon/a$ に比例します。つまり厚さを利用者が決めるのです。これが著者らの強調する「一定厚さ」の性質です。界面が常に同じセル数を保てば、衝撃・接触不連続を安定に捕捉できます。

1Dで直接: にじみとシャープニングの綱引き#

上の画面の綱引きを同じ形で書くと下のスクリプトです。移流が残したにじみを真似てから、シャープニングを仮想時間で数回回します。

import numpy as np
 
N, dx = 201, 1.0 / 200
x = np.linspace(0, 1, N)
 
def band_width(a, lo=0.05, hi=0.95):
    # 界面とみなすセル数(厚さの測定)
    return int(np.sum((a > lo) & (a < hi)))
 
def smear_once(a, D=0.16):
    # 移流スキームが毎ステップ界面を太らせる効果
    lap = np.zeros_like(a)
    lap[1:-1] = a[2:] - 2 * a[1:-1] + a[:-2]
    return a + D * lap
 
def sharpen_sweep(a, comp_a, eps, iters=200):
    # フラックス形の反拡散正則化 (Shukla 2010, 式 (38)-(39))
    dtau = 0.9 * min(dx * dx / (2 * eps), dx / comp_a)
    for _ in range(iters):
        af = 0.5 * (a[:-1] + a[1:])          # 面中心の alpha
        dA = a[1:] - a[:-1]
        s = np.sign(dA)
        J = comp_a * af * (1 - af) * s - eps * dA / dx   # 面フラックス
        a[1:-1] = np.clip(a[1:-1] - dtau / dx * (J[1:] - J[:-1]), 0, 1)
    return a
 
# 鋭い階段から始め -> 20回にじませる
a = 0.5 * (1 + np.tanh((x - 0.5) / 0.012))
for _ in range(20):
    a = smear_once(a)
print("にじんだ後のバンド:", band_width(a), "セル")     # -> 広がる
 
a = sharpen_sweep(a, comp_a=1.0, eps=0.004)
print("シャープニング後のバンド:", band_width(a), "セル") # -> 目標厚さへ収縮

にじんだ後のバンドは広く、シャープニング後には $\varepsilon/a$ が決める数セルへ縮みます。スキームはそのままです。sharpen_sweep は結果配列だけを受け取り、どのリーマンソルバが埋めたかは知りません。

既定値は comp_a=1.0、eps=0.004 です。厚さの尺度 $\varepsilon/a$ は 0.004 です。定常プロファイルの傾き係数は $a/(2\varepsilon)=125$ です。コードが始めた tanh の分母は 0.012 なので、目標は初期の階段より急です。 $dx=1/200=0.005$ です。仮想時間の上限は $dx^2/(2\varepsilon)=0.000025/0.008=0.003125$ と $dx/a=0.005$ の小さい方です。拡散側が勝ちます。0.9 を掛けた $d\tau$ は 0.0028125 です。この上限を超えると反拡散が符号を反転し、界面が振動します。

面フラックスの圧縮成分は $\alpha(1-\alpha)$ に比例します。その積は $\alpha=0.5$ で 0.25 が最大で、純流体では 0 なので圧縮は切れます。バンド関数は 0.05 と 0.95 の間を数えます。 $\varepsilon$ だけを2倍にすると目標厚さも2倍、 $a$ だけを2倍にすると半分です。導入の200ステップは物理時間です。6セルが20セルになるとバンドは約3.3倍です。smear_once の D=0.16 とシャープニングの iters=200 は別の層です。論文のステップあたり1〜3回は、物理ステップの拡散が小さいときの話です。Kapila の右辺 $\alpha_1 K\nabla\cdot\mathbf{u}$ は、衝撃で発散がスパイクするとき体積分率を直接押します。数値拡散とは別のソースなので、IST はそのソースを消しません。圧力緩和を別に置く Saurel(2009) の手順とは違い、ここでは界面のセル数だけを戻します。

厚さを一定にするということ

界面を薄くすることと、厚さを一定に固定することは違います。著者らが以前の手法(Tiwari ら 2013)と一線を画すのがここです。ソース項にシャープニング関数を加える方式は拡散誤差を減らしますが、厚さが時間・空間でばらつきます。一方、後処理で $\varepsilon/a$ を固定すれば、どこでも同じ厚さになります。

圧縮性流れでは、ここにもう一層加わります。 $\alpha$ を手入れすると、混合密度・エネルギーも一緒に変えないと熱力学が一貫しません。そこで著者らは $\alpha$ だけを直すのではなく、保存変数全体( $\alpha_1\rho_1$ 、 $\alpha_2\rho_2$ 、運動量、エネルギー)を混合則に合わせて再分配します。この「混合一貫正則化」のおかげで、界面での速度・圧力・温度の平衡が崩れません。論文の1D純界面移流テストがこれを示します。水と空気が同じ速度で流れるとき、シャープニングを入れても圧力・温度が保たれます。

2D気泡で: 円形界面を保つ#

1Dで通ったので2D気泡へ進みます。論文の衝撃波–気泡図は、界面が数セル厚の円を最後まで保ちます。numerical diffusion を入れると縁がぼやけ、sharpening と compression a を上げると再びくっきりしたリングへ締まります。

compression a: 1.20

interface band: 0 cells

teal ring = cells with 0.1 < α < 0.9. Fewer = sharper.

青緑のリングが締まっても気泡の中心はそのままです。大きくも小さくもなりません。減るのはリングのセル数だけです。圧縮が $\hat{\mathbf{n}}$ に沿うので、曲率はつぶれません。

再現しながらぶつかったこと

純流体で $|\nabla\alpha|\to 0$ になると法線はゼロ割です。論文は $a\,\alpha(1-\alpha)$ がそこでゼロだからと流しますが、分母に小さな数を足さないと NaN になります。 $d\tau$ が大きいと界面は振動します。CFL 型の上限が要ります。論文のステップあたり1〜3回は、物理ステップの拡散が小さいときの話です。粗い格子や強い衝撃ではその回数では足りないことがあります。

もう一つ。後処理という優雅さには代償があります。シャープニングは支配方程式とは別に $\alpha$ を変えるので、その瞬間の移流は元の PDE の解ではありません。厚さが一定に保たれる代わりに、界面近くの局所保存がわずかにずれる余地があります。著者らは混合一貫再分配でこれを最小化しますが、完全にゼロではありません。OpenFOAM の interfoam が使う圧縮項( $\nabla\cdot(\alpha(1-\alpha)\mathbf{u}_r)$ )も同じ系統の考えなので、この手法を圧縮性ソルバへ移植するときの比較基準になります。

ε/aを固定しないとセル数がまた浮く#

ソース項にシャープニングを足すと拡散誤差は減っても、厚さは時空間で揺れます。ステップの後に $\varepsilon/a$ を固定して、衝撃の横でも同じセル数を残します。 $\alpha$ だけ直して保存変数を混合則で分配し直さないと、そのセルで圧力と温度の平衡が先に壊れます。

$d\tau$ が $dx^2/(2\varepsilon)$ より大きいと、反拡散は界面を締める代わりに振動させます。厚さを定数で買う代償は、その瞬間の移流が元の PDE の解ではないことです。局所保存がわずかにずれる場所は、Abgrall 基準が捕まえる非保存項と同じ種類です。