[論文レビュー] 振幅の小さい方を選べ — Deng (2018) MUSCL-THINC-BVD によるインターフェース再構成

1セルの純粋な跳躍は片側の差分が0になり、THINCは境界変動の比較の前に外れます。

東京工業大学のXiaoグループは2018年に、あまり公表されない厚さを測りました。圧縮性の二相流、1000ステップ。1セルだったジャンプが8セルになりました。WENOでも止まりません。1ステップの散逸は小さく見えても、時間が積もると界面は崩れます。

持ち込まれた再構成が MUSCL-THINC-BVD です。セルごとに候補は二つです。境界ジャンプが小さい方だけが残ります。

βのスライダーを1.4から2.2へ上げると、ジャンプセルの tanh は狭くなり、ガウス分布の肩はギザギザに反応します。

N = 80β = 1.60CFL = 0.40

MUSCL THINC BVDt = 0.00THINC picked: 0% of cells

矩形波の端だけに青緑のバーが点くかを見ます。1.6付近でそのバーが増えなくなるかも、一緒に見ます。

一ページ要約

著者・誌名: Deng, Inaba, Xie, Shyue, Xiao. Journal of Computational Physics 371 (2018) 945–966.
問題: 五方程式モデルで解く圧縮性二相流では、物質界面が時間とともににじみます。
提案: 各セルで MUSCL(滑らかな区間)と THINC(ジャンプ)を同時に作り、境界ジャンプが小さい方を残す BVD(Boundary Variation Diminishing) です。
違う点: 後処理の anti-diffusion も人工圧縮もありません。体積分率と残りの保存量に同じ規則を掛けるので、変数の間の整合は後から付いてきます。

衝突する二つの要求

圧縮性二相流は、二つのことを同時に求めます。滑らかな場所は正確で、散逸が少ない必要があります。ジャンプ、つまり界面と衝撃では、単調性と薄い厚さが必要です。一つの関数が両方をこなすのは難しいです。

MUSCL はセルの中で直線です。単調性は守りますが、勾配は1次までです。界面を通るたびに少し削れます。1000ステップで1セルが8セルになる蓄積です。THINC はセル内に tanh を当て、段差を1、2セルに閉じ込めます。滑らかな場所にそのまま使うと、偽の段差が出ます。

論文の出発点は決断です。二つを一つの関数に混ぜません。セルごとに一つを選びます。

MUSCL — 信頼できるディフューザー#

基準の候補は minmod の勾配制限を付けた MUSCL です。

\tilde{q}_i^{\,\text{MUSCL}}(x) = \bar{q}_i + \sigma_i (x - x_i), \quad \sigma_i = \text{minmod}(\bar{q}_i - \bar{q}_{i-1},\; \bar{q}_{i+1} - \bar{q}_i)

$\bar q_i$ はセル平均です。 $\sigma_i$ は中心での勾配です。minmod は二つの差分の符号が同じなら絶対値の小さい方を、違うなら 0 を返します。

境界値は次のように読めます。

q^{L,\text{MUSCL}}_{i+1/2} = \bar{q}_i + \tfrac{1}{2}\sigma_i, \quad q^{R,\text{MUSCL}}_{i-1/2} = \bar{q}_i - \tfrac{1}{2}\sigma_i

振動はありません。ジャンプ区間では厚さだけがゆっくり広がります。

THINC — ジャンプを真似る単調関数#

THINC(Tangent of Hyperbola for INterface Capturing)は、セル内のジャンプを双曲線正接で当てはめます。

\tilde{q}_i^{\,\text{THINC}}(x) = \bar{q}_{\min} + \tfrac{\bar{q}_{\max}}{2}\left[1 + \theta \tanh\!\left(\beta\!\left(\tfrac{x - x_{i-1/2}}{\Delta x} - \tilde{x}_i\right)\right)\right]

$\bar q_{\min}$ と $\bar q_{\max}$ は、隣接平均から取った最小値と振幅です。 $\theta$ はジャンプの向きです。 $\beta$ は厚さです。 $\tilde x_i$ は、セル平均が保たれるように解くジャンプ中心です。

$\beta$ は 1.4 から 2.0 で安定です。標準は 1.6 です。それより大きいとジャンプは1セルに閉じ、見た目はきれいです。滑らかな場所に入るとギザギザになります。tanh を衝撃側に付けた記録は TENO と THINC を一つのステンシルに重ねた再構成に別にあります。あちらは候補を選ばず、重みで混ぜます。

BVD — 二つの候補間の変動を測る#

候補が二つあるなら、隣接の再構成と境界で会ったとき、どちらのジャンプが小さいかを測ればよいのです。セル $i$ の総境界変動(Total Boundary Variation, TBV)は次のように定義されます。

\text{TBV}_i^{P} = |q^{L,\text{MUSCL}}_{i-1/2} - q^{R,P}_{i-1/2}| + |q^{L,P}_{i+1/2} - q^{R,\text{MUSCL}}_{i+1/2}|

$P$ は MUSCL か THINC です。両隣は MUSCL に固定します。中央のセルだけが替わります。

選択規則はこうです。

\tilde{q}_i^{\,\text{BVD}}(x) = \begin{cases} \tilde{q}_i^{\,\text{THINC}}(x) & \text{if } \delta < C_i < 1 - \delta,\ (\bar{q}_{i+1}-\bar{q}_i)(\bar{q}_i-\bar{q}_{i-1}) > 0,\ \text{TBV}^{\text{THINC}}_i < \text{TBV}^{\text{MUSCL}}_i \\ \tilde{q}_i^{\,\text{MUSCL}}(x) & \text{otherwise} \end{cases}

$C_i$ は、ジャンプがセルのどこにあるかです。 $\delta$ はおよそ $10^{-4}$ です。

選択式の最初の関門は差分の積です。隣接する二つの片側差分を掛けて、0より大きい必要があります。片方が0なら積も0で、関門を越えません。0と1が1セルで接する段差が、まさにその形です。ジャンプ直前のセルは右側の差分だけを持ち、左側の差分はありません。総境界変動を測る前に MUSCL へ戻ります。

二つ目の関門は位置の尺度です。 $\delta$ は $10^{-4}$ の付近です。隣接が取る振幅のその割合の中にセル平均が入ると、THINC の資格がありません。平らな頂を切り落とす装置です。

五方程式モデル上の整合性

再構成の比較で終わらない場所がここです。五方程式モデルは体積分率 $\alpha_1$ 、相密度 $\alpha_1\rho_1$ と $\alpha_2\rho_2$ 、運動量 $\rho u$ 、全エネルギー $E$ を一緒に解きます。五つの変数に同じ BVD 規則を掛けます。

以前の方法は体積分率だけを鋭くし、残りは圧力振動が出ないように別々に直していました。BVD は五つの変数が同じセルで同じ決定に従うようにします。ジャンプセルなら五つとも THINC です。滑らかなら五つとも MUSCL です。変数の間の整合は付いてきて、anti-diffusion の段を別に置きません。

NumPy で見る BVD 決定木#

import numpy as np
 
def muscl_minmod_edges(q):
    """ (q[i-1], q[i], q[i+1]) を受け取り、MUSCL の左右端値を返す """
    qm, q0, qp = q
    a, b = q0 - qm, qp - q0
    if a * b <= 0.0:
        slope = 0.0
    else:
        slope = np.sign(a) * min(abs(a), abs(b))
    return q0 - 0.5 * slope, q0 + 0.5 * slope  # qL, qR
 
def thinc_jump_edges(q, beta=1.6, eps=1e-20):
    """ tanh フィットによるジャンプ: THINC の左右端値 """
    qm, q0, qp = q
    qmin = min(qm, qp)
    qmax = max(qm, qp) - qmin
    if qmax < 1e-12:
        return q0, q0
    theta = np.sign(qp - qm)
    C = (q0 - qmin + eps) / (qmax + eps)
    if C <= 1e-6 or C >= 1 - 1e-6:
        return q0, q0
    B = np.exp(theta * beta * (2.0 * C - 1.0))
    A = (B / np.cosh(beta) - 1.0) / np.tanh(beta)
    qR = qmin + 0.5 * qmax * (1.0 + theta * A)
    num = 1.0 + theta * A * np.tanh(beta) + theta * np.tanh(beta)
    den = 1.0 + A * np.tanh(beta)
    qL = qmin + 0.5 * qmax * (num / den)
    return qL, qR
 
def bvd_decide(q_window, beta=1.6, delta=1e-4):
    """
    q_window: 長さ5の配列 [q[i-2], q[i-1], q[i], q[i+1], q[i+2]]
    戻り値: (qL_i, qR_i, picked_thinc)
    """
    qm2, qm1, q0, qp1, qp2 = q_window
    monotone = (qp1 - q0) * (q0 - qm1) > 0.0
    qmin = min(qm1, qp1)
    qmax = max(qm1, qp1) - qmin
    C = 0.5 if qmax < 1e-12 else (q0 - qmin) / qmax
    eligible = monotone and (delta < C < 1.0 - delta)
 
    qL_M, qR_M = muscl_minmod_edges([qm1, q0, qp1])
    if not eligible:
        return qL_M, qR_M, False
 
    qL_T, qR_T = thinc_jump_edges([qm1, q0, qp1], beta=beta)
    # 隣接を MUSCL に固定した TBV 比較
    _, qR_left = muscl_minmod_edges([qm2, qm1, q0])
    qL_right, _ = muscl_minmod_edges([q0, qp1, qp2])
    tbv_M = abs(qR_left - qL_M) + abs(qR_M - qL_right)
    tbv_T = abs(qR_left - qL_T) + abs(qR_T - qL_right)
    if tbv_T < tbv_M:
        return qL_T, qR_T, True
    return qL_M, qR_M, False
 
# テスト: 滑らか + ジャンプ混合プロファイル
N = 80
x = (np.arange(N) + 0.5) / N
q = np.where((x >= 0.6) & (x <= 0.8), 1.0, 0.0)
q += 0.5 * np.exp(-((x - 0.25) / 0.04) ** 2)
 
picks = 0
for i in range(2, N - 2):
    _, _, used_thinc = bvd_decide(q[i-2:i+3])
    picks += int(used_thinc)
print(f"THINC 採用率: {picks}/{N-4}")
# 期待値: ジャンプ両側の数セルのみ

スニペット末尾の注は、採用 4/76 と書きます。すぐ上の関門と、その数字は同じ形を指しません。コードの矩形波は条件分岐で 0 か 1 です。1セルの跳躍なので差分の一つが 0 です。注の四つのセルは、勾配が2、3セルにわたって同じ符号で上がるランプを数えて初めて出ます。1000ステップ後の八つのセルは、すでにそのランプです。厚さが1セルのままである間、この分岐は THINC を呼びません。にじんだ後ではじめて呼びます。

直線では逆です。隣接差分が等しければ minmod の勾配はその差分です。左の面値は左隣の右面値と会い、右側も同じです。MUSCL 側の総境界変動は 0 です。THINC の tanh は直線に曲がりを入れます。面に隙間ができ、変動が大きくなります。不等号は MUSCL 側です。ガウス分布がセルより広ければ、局所的にはこの直線に近いです。滑らかな区間が候補の入れ替えを拒む理由です。

滑らかなガウス分布は局所的に直線に近いので、入れ替えは拒まれます。矩形波が1セルで付いたままでも、関門は開きません。

セルごとにどちらが勝つか

上の移流は三つの再構成を一緒に進め、BVD が THINC を選んだセルをバーで塗ります。 $\beta$ を 1.4 から 2.2 へ上げると、そのセルのジャンプだけが狭くなります。大きすぎるとガウス分布の肩にギザギザが見えます。1.6 でバーの数が増えなくなる場所が標準値です。矩形波が少しにじんでから、端に青緑が付きます。厚さが1セルのときはバーが空です。

どこまで信じられるか

収束は2次に届きません。候補の片方が MUSCL なので、滑らかな区間の精度は MUSCL で止まります。WENO のような高次候補に替えれば5次も可能ですが、この論文は単純な側を選びました。

TBV は隣が MUSCL だと仮定します。両側が一緒に THINC へ行く場合は、一回のスイープで無視します。論文4節の図はその単純化でも許容できました。より一般的な形は後に回されています。

多次元では $\beta$ が面法線の方向に分かれます。界面法線は Young のアルゴリズムで推定します。角と折れでは法線がずれ、ジャンプが一方へ傾きます。

時間積分も別です。論文は SSPRK3 を使います。BVD の判定は Runge–Kutta の段ごとにやり直す必要があります。段の間で判定が揺れると界面が震えます。その震えは後続の研究が再び扱いました。

1セルの跳躍で選択が止まる場所#

1セルで付いたままの間、この分岐は THINC を呼びません。厚さが2、3セルににじんでから、差分の積が正になります。界面を1セルに留める仕事と、留めた後でなければ選べない仕事は、互いにずれます。

そのずれが時間刻みの側へ漏れるかが次の問題です。界面を薄く保っても、移流の上限は残ります。毛管項を外したあとも CFL 0.05 が残った界面移流は、その上限がどこにあるかを問います。BVD がジャンプを1セルに閉じるなら、そのセルを運ぶ時間刻みはどこまで大きくできるでしょうか。