何も起きないはずの問題で圧力が21%跳ねた

非保存項の離散化は自由に選べるものではありません。保存フラックスを先に決めた時点で、一様流れの条件が残りの項の係数を一つに固定します。

何も起きないはずの問題が最初のテストでした

多成分圧縮性ソルバーの最初の検証ケースを選びました。二つの気体が接触面一つを挟んで向かい合っています。圧力はどこでも100 kPa、速度はどこでも100 m/s。異なるのは二つの気体の密度と比熱比だけです。

厳密解は退屈です。接触面が右へ流れていくだけで、圧力と速度は最後まで一様のままです。格子が粗くても、時間刻みを大きく取っても同じです。この問題には解くべき波がありません。

20ステップ後、界面セルの圧力が78.6 kPaまで下がりました。一様値から21%外れた値です。速度はきっちり100 m/sを保ち、質量とエネルギーも機械精度まで保存されました。フラックスは正常でした。狂っていたのはフラックスの隣に立っていた項です。

下のシミュレーションでその振動を自分で作ってみましょう。

gamma_21.667rho_1/rho_28.0CFL0.50

step 0A  max|P-P0| 0.00 kPa (0.0%)B  max|P-P0| 0.0e+0 Pa

Drag gamma_2 toward 1.400: the red curve flattens onto the green one, because with equal gamma the mixture rule is no longer nonlinear. Push it away, or raise the density ratio, and the red curve dips at the two interfaces while the green line never leaves 100 kPa. CFL changes how fast the damage accumulates, not whether it appears.

二本の曲線は同じ保存変数を同じ風上演算子で運びます。違うのは熱力学変数の運び方だけです。 gamma_2を1.400側へ引くと、赤い曲線が緑の曲線の上に沈み込みます。逆へ押すか密度比を上げると、界面二か所で赤い曲線が下へ落ち込みます。圧力が一様に保たれるかどうかは、スキームの精度ではなくこの選択一つで決まります。

この計算で正確に何が狂ったのか

圧力は保存変数ではありません。状態方程式から取り戻す値です。比熱比 $\gamma$ が空間的に変わると

P = \frac{\rho E - \tfrac{1}{2}\rho u^{2}}{\Gamma}, \qquad \Gamma \equiv \frac{1}{\gamma - 1}

ここで $\rho E$ は単位体積あたりの全エネルギー、 $\Gamma$ は比熱比を担う熱力学変数です。

一様圧力・速度のもとで全エネルギーは $\rho E = P\Gamma + \tfrac{1}{2}\rho u^{2}$ です。風上移流は線形なので、 $\rho E$ の更新値は $P \cdot (\Gamma$ の風上移流 $)$ に運動エネルギー項を足したものになります。すると新しい圧力はこう整理できます。

P^{n+1}_j = P^{n} \, \frac{\mathcal{U}(\Gamma)_j}{\Gamma^{n+1}_j}

$\mathcal{U}$ は $\rho E$ に実際に適用したその風上演算子です。

条件は分数一つに帰着します。 $\Gamma^{n+1}$ を $\mathcal{U}(\Gamma)$ と同じ値にして初めて圧力がその場に留まります。つまり $\Gamma$ の運び方が、エネルギーの運び方と代数的に噛み合う必要があります。

私が最初に書いたコードは質量分率 $Y$ を保存形で運び、混合則から $\Gamma$ を復元していました。 $\Gamma(Y)$ は $Y$ について線形ですが、 $Y^{n+1} = \mathcal{U}(\rho Y)/\mathcal{U}(\rho)$ は密度が変わる場所では $\mathcal{U}(Y)$ になりません。二つの分数がずれると、その差がそのまま圧力誤差です。

Abgrallが1996年に立てた要求は一行です#

Reと Abgrall が弱圧縮性多成分モデルを構築する際に引用した基準は、こういう一文です。「圧力と速度が一様な二相流れは、時間が経っても同じ変数について一様に保たれなければならない。」論文はこれを圧力非擾乱条件(pressure non-disturbance condition)、あるいはAbgrallの基準と呼びます。

この一文が特異なのは、精度への要求ではないという点です。1次でも5次でも関係ありません。安定性条件でもありません。CFLを下げても振動はそのまま残ります。これは離散化どうしの代数的整合への要求です。ある方程式でどの演算子を使ったかが、別の方程式で許される演算子を決めてしまいます。

同じ事情は二流体モデルの双曲性と界面圧力で界面圧力項を選ぶときにも現れます。あちらでは固有値が実数になるように項を選び、ここでは一様流れが保たれるように項を選びます。どちらも「物理的にもっともらしい離散化」が複数あり、そのうち一つだけが生き残ります。

Pythonで同じ格子に二つの離散化を並べました#

100セルの周期格子に接触面を一つ置き、二つの方式を同じステップ数だけ回しました。保存変数 $\rho$ 、 $\rho u$ 、 $\rho E$ は両方式とも同一の1次風上で更新します。違いは熱力学変数一つだけです。

G1, G2 = 1.4, 1.667      # 二つの気体の比熱比
P0, U0 = 1.0e5, 100.0    # 一様圧力 [Pa], 一様速度 [m/s]
R1, R2 = 1.0, 0.125      # 二つの気体の密度 [kg/m^3]
 
 
def gamma_var(y):
    """質量分率 y -> 1/(gamma-1)。y について線形な混合則。"""
    return y / (G1 - 1.0) + (1.0 - y) / (G2 - 1.0)
 
 
def advect_upwind(q, lam):
    """u>0 の1次風上移流。左側の流入セルは固定。"""
    return [q[0]] + [q[j] - lam * (q[j] - q[j - 1]) for j in range(1, len(q))]
 
 
def initial_state(n):
    x = [(j + 0.5) / n for j in range(n)]
    y = [1.0 if xi < 0.3 else 0.0 for xi in x]
    rho = [R1 if xi < 0.3 else R2 for xi in x]
    return x, y, rho
 
 
def step_massfraction_closure(rho, mom, ene, ry, lam):
    """rho*Y を保存形で運び、混合則から gamma を復元する。"""
    rho_n = advect_upwind(rho, lam)
    mom_n = advect_upwind(mom, lam)
    ene_n = advect_upwind(ene, lam)
    ry_n = advect_upwind(ry, lam)
    y_n = [ry_n[j] / rho_n[j] for j in range(len(rho_n))]
    p_n = [(ene_n[j] - 0.5 * mom_n[j] ** 2 / rho_n[j]) / gamma_var(y_n[j])
           for j in range(len(rho_n))]
    return rho_n, mom_n, ene_n, ry_n, p_n
 
 
def step_gammavar_transport(rho, mom, ene, gv, lam):
    """1/(gamma-1) を非保存(移流)形で同じ風上演算子に乗せる。"""
    rho_n = advect_upwind(rho, lam)
    mom_n = advect_upwind(mom, lam)
    ene_n = advect_upwind(ene, lam)
    gv_n = advect_upwind(gv, lam)
    p_n = [(ene_n[j] - 0.5 * mom_n[j] ** 2 / rho_n[j]) / gv_n[j]
           for j in range(len(rho_n))]
    return rho_n, mom_n, ene_n, gv_n, p_n
 
 
def run_interface_advection(n=100, steps=60, cfl=0.5):
    x, y0, rho0 = initial_state(n)
    gv0 = [gamma_var(v) for v in y0]
 
    rho_a = list(rho0)
    mom_a = [r * U0 for r in rho0]
    ene_a = [P0 * gv0[j] + 0.5 * rho0[j] * U0 ** 2 for j in range(n)]
    ry_a = [rho0[j] * y0[j] for j in range(n)]
 
    rho_b, mom_b, ene_b = list(rho_a), list(mom_a), list(ene_a)
    gv_b = list(gv0)
 
    hist, p_a, p_b = [], None, None
    for k in range(1, steps + 1):
        rho_a, mom_a, ene_a, ry_a, p_a = step_massfraction_closure(
            rho_a, mom_a, ene_a, ry_a, cfl)
        rho_b, mom_b, ene_b, gv_b, p_b = step_gammavar_transport(
            rho_b, mom_b, ene_b, gv_b, cfl)
        if k % 20 == 0:
            ea = max(abs(v - P0) for v in p_a)
            eb = max(abs(v - P0) for v in p_b)
            eu = max(abs(mom_a[j] / rho_a[j] - U0) for j in range(n))
            hist.append((k, ea, eb, eu))
    return x, p_a, p_b, hist
 
 
if __name__ == "__main__":
    x, p_a, p_b, hist = run_interface_advection()
    print("step |  max|P-P0| mixrule |  max|P-P0| Gamma-adv |  max|u-U0| mixrule")
    for k, ea, eb, eu in hist:
        print("%4d | %16.2f | %19.2e | %16.3e" % (k, ea, eb, eu))
 
    j = max(range(len(p_a)), key=lambda i: abs(p_a[i] - P0))
    print("\nworst cell x=%.3f  P=%.1f Pa  (uniform value %.0f Pa)" % (x[j], p_a[j], P0))
    print("relative error:  mixrule %.2f%%   Gamma-adv %.1e%%"
          % (max(abs(v - P0) for v in p_a) / P0 * 100,
             max(abs(v - P0) for v in p_b) / P0 * 100))

step |  max|P-P0| mixrule |  max|P-P0| Gamma-adv |  max|u-U0| mixrule
  20 |         21433.25 |            2.91e-11 |        0.000e+00
  40 |         21587.18 |            4.37e-11 |        0.000e+00
  60 |         21442.28 |            4.37e-11 |        2.842e-14
 
worst cell x=0.635  P=78557.7 Pa  (uniform value 100000 Pa)
relative error:  mixrule 21.44%   Gamma-adv 4.4e-14%

読むべきは三行だけです。混合則の側は21.4 kPaから始まり、ステップを増やしても減りません。 $\Gamma$ 移流の側は $10^{-11}$ Pa、相対誤差で $10^{-14}$ % です。倍精度の底です。

速度の列がゼロであることも重要です。運動量方程式は最初から最後まで正しく解けていました。格子を二倍細かくしても21%は21%のまま残ります。収束しない誤差だからです。

論文が $H_u$ を導いた道筋 — スキームを先に決め、項を後から合わせる#

Reと Abgrall の Baer–Nunziato 型(二相がそれぞれ速度・圧力を持つ7方程式系)モデルには、体積分率方程式が別にあります。この方程式は保存形ではありません。

\frac{\partial \alpha_i}{\partial t} + u_I \frac{\partial \alpha_i}{\partial x} = 0

$\alpha_i$ は相 $i$ の体積分率、 $u_I$ は界面速度です。論文はこの項を $H_u(\alpha_i, u_I)$ という名前の離散演算子として置き、その形を仮定しません。代わりに条件から引き出します。

質量方程式 $(\alpha_i\rho_i)$ はRusanovフラックスで既に決まっています。密度と速度が一様な状態を入れると、そのフラックスは $\rho_i \times (\alpha_i$ のRusanovフラックス $)$ に因数分解されます。更新後も $\rho_i = (\alpha_i\rho_i)^{n+1}/\alpha_i^{n+1}$ がそのまま残るには、分母 $\alpha_i^{n+1}$ が分子とまったく同じフラックス差分で更新されなければなりません。だから $H_u$ は選択ではなく結果です。

H_u(\alpha_i, u_I)_j = \frac{1}{2}\Big[ \big(\alpha_{j+1} - \alpha_{j-1}\big) u_{I,j} - \big|u_{I,j}\big| \big(\alpha_{j+1} - 2\alpha_j + \alpha_{j-1}\big) \Big]

前の項は中心差分、後ろの項は $|u_I|$ の付いた拡散です。二つの和はちょうど $\alpha$ に対する Rusanovフラックス差分になります。 $u_I > 0$ なら $\alpha_j - \alpha_{j-1}$ だけが残ります。純粋な風上です。

論文は圧力方程式でも同じ $H_u$ を再利用します。 $u^{*}\partial P/\partial x$ を $\partial(Pu^{*})/\partial x - P\partial u^{*}/\partial x$ に分けたあと、残った非保存項を質量方程式と同じ演算子に乗せます。方程式ごとに違う離散化を使えば、今合わせたばかりの整合が壊れるからです。

三点ステンシルの係数を手で回してみます

$H_u$ の二項のうち拡散項に重み $\theta$ を付けると、ダイヤル一つで二つの極端を行き来できます。 $\theta = 1$ なら論文の $H_u$ 、 $\theta = 0$ なら $u_I \partial\alpha/\partial x$ の純粋な中心差分です。

theta1.00CFL0.40

step 0weights [-1.00, 1.00, 0.00]max|rho - 850| 0.0e+0 kg/m^3

Slide theta down from 1 and watch the downwind weight come back to life: the moment the third cell rejoins the stencil, the recovered density leaves 850 and never returns. Flip u_I and theta = 1 still holds, because the flux difference follows the sign of the interface velocity while the centred stencil does not.

theta を1から下げると、下流セルの重みがゼロから復活します。その瞬間、復元した密度が850 kg/m³を離れて戻ってきません。u_I の符号を反転させても $\theta = 1$ は耐えます。フラックス差分は界面速度の符号に従いますが、中心ステンシルはそうならないからです。

数字でも確認しました。一様密度850 kg/m³、界面一つ、同じ格子で $\alpha$ だけを二つの演算子で更新しました。

RHO, UI, N, LAM = 850.0, 1.0, 80, 0.4   # 一様密度 [kg/m^3], 界面速度, セル数, u*dt/dx
 
 
def alpha_profile():
    """界面を挟んだ体積分率。0.02 <-> 0.98 を三セルにわたってつなぐ。"""
    a = []
    for j in range(N):
        if j < 30:
            a.append(0.98)
        elif j < 33:
            a.append(0.98 - 0.32 * (j - 29))
        else:
            a.append(0.02)
    return a
 
 
def rusanov_flux(q, j, vel):
    """セル j と j+1 の間のRusanov数値フラックス(周期境界)。"""
    ql, qr = q[j % N], q[(j + 1) % N]
    return 0.5 * (qr + ql) * vel - 0.5 * abs(vel) * (qr - ql)
 
 
def hu_upwind(a, j):
    """論文の式 (10) の非保存演算子: alpha に対するRusanovフラックス差分。"""
    return rusanov_flux(a, j, UI) - rusanov_flux(a, j - 1, UI)
 
 
def hu_central(a, j):
    """中心差分で u_I * d(alpha)/dx をそのまま離散化した版。"""
    return 0.5 * UI * (a[(j + 1) % N] - a[(j - 1) % N])
 
 
def march(op, steps):
    """alpha*rho はRusanov、alpha は op で前進させ、復元した rho を見る。"""
    a = alpha_profile()
    ar = [RHO * v for v in a]
    for _ in range(steps):
        ar = [ar[j] - LAM * (rusanov_flux(ar, j, UI) - rusanov_flux(ar, j - 1, UI))
              for j in range(N)]
        a = [a[j] - LAM * op(a, j) for j in range(N)]
    return max(abs(ar[j] / a[j] - RHO) for j in range(N))
 
 
if __name__ == "__main__":
    print("steps |  upwind H_u [kg/m^3] |  centred [kg/m^3]")
    for s in (10, 40, 120):
        print("%5d | %20.2e | %17.4f" % (s, march(hu_upwind, s), march(hu_central, s)))
 
    a = alpha_profile()
    lhs = hu_upwind(a, 31)
    rhs = 0.5 * ((a[32] - a[30]) * UI - abs(UI) * (a[32] - 2 * a[31] + a[30]))
    print("\ncell 31:  flux difference %.6f   paper Eq.(10) %.6f   gap %.1e"
          % (lhs, rhs, abs(lhs - rhs)))

steps |  upwind H_u [kg/m^3] |  centred [kg/m^3]
   10 |             2.27e-13 |         2553.2027
   40 |             4.55e-13 |         6697.6224
  120 |             5.68e-13 |         1206.6015
 
cell 31:  flux difference -0.320000   paper Eq.(10) -0.320000   gap 0.0e+00

風上 $H_u$ は120ステップ回しても $10^{-13}$ kg/m³ に留まります。中心差分は10ステップで2,553 kg/m³ も外れます。120ステップで値が1,207まで下がったのは回復ではなく発散です。 $\alpha$ がゼロ近傍を行き来し、割り算が何でも吐き出す領域に入ったという意味です。

最後の行は導出が正しかったことの確認です。 $\alpha$ に対するRusanovフラックス差分と論文の式 (10) の閉じた形が、セル31で桁まで一致しています。二つの表現は代数的に同じものです。

この条件が精度と無関係だという言葉の意味

一つはっきりさせておきます。 $H_u$ を使えば解が正確になる、ということではありません。 1次風上は依然として界面をなまします。上の図でも $\alpha$ の階段はステップごとに厚くなります。

Abgrall基準が保証するのは別の種類のことです。外れるとしても物理的に筋の通る方向に外れる。 界面がにじむのは数値拡散で、それは格子を細かくすれば減ります。一様圧力場に21 kPaのスパイクが立つのは数値拡散ではありません。それは格子を細かくしても減らず、状態方程式が硬い(stiff)ほど大きくなり、しばしば負圧で計算を止めます。

陰的表面張力が開く時間刻みの幅を扱ったときにも同じ区別が出てきました。ある制約は精度を買えば解け、ある制約は買っても解けません。 Abgrall基準は後者です。高次スキームに上がっても、この条件は別途あらためて合わせる必要があります。

次に界面で圧力が跳ねたらどこから見るか

私がこの問題に再会したら、順序はこうです。

まず一様圧力・速度の問題を回します。波のない問題で圧力変動がゼロでなければ、フラックスを見る必要はありません。答えは非保存項か状態方程式の復元の側にあります。

次に格子を二倍細かくして同じ値を測ります。振動が半分に減れば数値拡散の問題です。そのままなら整合性の問題です。この一度の実行が二つの原因を切り分けてくれます。

最後に、方程式ごとに非保存項へどのステンシルを使ったかを並べて書き出します。質量はRusanov、体積分率は中心差分、エネルギーはまた別の何か、だったならそこが原因です。論文が $H_u$ 一つを三つの方程式に使い回したのは、コードを短くしたいという趣味ではありませんでした。

参考文献

B. Re, R. Abgrall, Non-equilibrium Model for Weakly Compressible Multi-component Flows: the Hyperbolic Operator, arXiv:1911.00270 — §2.2 The discretization
R. Abgrall, How to Prevent Pressure Oscillations in Multicomponent Flow Calculations: A Quasi Conservative Approach, J. Comput. Phys. 125 (1996)

何も起きないはずの問題で圧力が21%跳ねた — Abgrall条件と非保存項

何も起きないはずの問題が最初のテストでした

この計算で正確に何が狂ったのか

Abgrallが1996年に立てた要求は一行です#

Pythonで同じ格子に二つの離散化を並べました#

論文が $H_u$ を導いた道筋 — スキームを先に決め、項を後から合わせる#

三点ステンシルの係数を手で回してみます

この条件が精度と無関係だという言葉の意味

次に界面で圧力が跳ねたらどこから見るか

関連記事

何も起きないはずの問題で圧力が21%跳ねた — Abgrall条件と非保存項

何も起きないはずの問題が最初のテストでした

この計算で正確に何が狂ったのか

Abgrallが1996年に立てた要求は一行です#

Pythonで同じ格子に二つの離散化を並べました#

論文が HuH_uHu​ を導いた道筋 — スキームを先に決め、項を後から合わせる#

三点ステンシルの係数を手で回してみます

この条件が精度と無関係だという言葉の意味

次に界面で圧力が跳ねたらどこから見るか

関連記事

論文が $H_u$ を導いた道筋 — スキームを先に決め、項を後から合わせる#