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

数値拡散でぼやけた界面を一定の厚さに固定する後処理手法

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

論文情報

タイトル: 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で直接: にじみとシャープニングの綱引き#

言葉より目が速いです。下のシミュレーションで直接操作してみましょう。numerical diffusion を入れると界面がひとりでに広がり始めます。その状態で sharpening (IST) を入れると、広がりが止まり、厚さが目標値の近くに固定されます。

compression a: 1.00regularization ε: 0.0040

measured band: 0 cells

target ≈ ε/a: 2 cells

steady tanh thickness scales like ε/a

compression a を上げると界面が薄くなり、regularization ε を上げると厚くなります。二つのスライダーの比 $\varepsilon/a$ が目標厚さ(黄色い帯)を決めることは、測定バンドのセル数がその値を追うことで確認できます。

同じ論理を numpy に移すとこうなります。移流スキームの拡散を真似た後、シャープニング方程式を数回回します。

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 は結果配列を受け取って手入れするだけで、どのリーマンソルバを使ったかを知りません。

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

界面を薄くすることと、厚さを一定に固定することは違います。著者らが以前の手法(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}}$ に沿って圧縮するため、円の曲率がつぶれずに保たれます。

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

三つが引っかかりました。第一に、 $\hat{\mathbf{n}} = \nabla\alpha/|\nabla\alpha|$ は純流体領域で $|\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)$ )も同じ系統の考えなので、この手法を圧縮性ソルバへ移植するときの比較基準になります。

この論文が変えたこと

スキームと界面手入れの分離。 リーマンソルバをそのままにし、界面厚さだけを後処理で制御する。移植性が高い。
厚さ = $\varepsilon/a$ 。 界面を薄くではなく一定に保つ。衝撃・接触不連続の捕捉に必要な性質を利用者が直接決める。
圧縮性での熱力学一貫性。 $\alpha$ だけでなく保存変数全体を混合則に合わせて再分配し、界面での速度・圧力・温度平衡を守る。