流出境界に音波を跳ね返させない方法 — Navier

入射波をゼロにすると音波は出ますが、平均圧は歩きます。σ≈0.25の緩和が、その漂流を引き戻します。

圧縮性 LES が 6 時間後に発散しました。ログに NaN もなく、衝撃波もありません。平均圧力が 0.3 bar ほど持ち上がり、音波が出口から領域内へ跳ね返っていました。原因は出口での zero-gradient 外挿でした。今日は Poinsot–Lele (1992) の NSCBC を、1D の音波パルス一本だけで最後まで追いかけます。出口は「値」ではなく「入ってくる波の振幅」を指定すべきだ、という話です。

Mach スライダーは 0 から 1.6 まで動きます。数字を上げると、内側へ戻る矢印が減ります。

Mach number M0.30

亜音速では太い赤い矢印が一つ残ります。1 を超えるとその矢印は境界の外へ出ます。出口が指定すべき信号の数が、その画面にあります。

領域が 6 時間で発散した#

圧縮性ソルバーの外部境界は、衝撃波よりも多くシミュレーションを殺します。zero-gradient 出口は運動量とエネルギーをゴーストセルへコピーしますが、そのコピーが音響インピーダンスの不整合を作ります。出ていくはずの音波がインピーダンスの跳びに当たると、その一部は必ず内側へ戻ります。その跳びでどれだけ透過し、どれだけ戻るかは、界面を通る音波がスネル側から書いています。

LES や DNS、燃焼計算ではこの反射が平均場に蓄積します。レイノルズ平均が落ち着く前に出口が「共鳴器」として働き始め、圧力がゆっくり漂流します。NSCBC はこの問題群を、「入ってくる特性波の振幅を直接指定する」という一つの原理で解決します。

境界から出てくる 5 つの信号#

保存変数 $U = (\rho, \rho u, E)^T$ について 1D Euler 方程式を書くと

\frac{\partial U}{\partial t} + A(U)\frac{\partial U}{\partial x} = 0

ここで $A = \partial F/\partial U$ はフラックスヤコビアンです。右固有分解により 5 個の固有値 (3D の平面境界での 1D 投影で) $\lambda_i \in \{u-c,\, u,\, u,\, u,\, u+c\}$ が得られます。 $u$ は流速、 $c$ は音速です。

それぞれが一つの特性波に対応します。 $u-c$ は左向きの音波、 $u$ はエントロピー・横速度情報 (流体と一緒に運ばれる)、 $u+c$ は右向きの音波です。

上の Mach スライダーで太い線になっていたのが $u-c$ です。亜音速の出口ではそれだけが内側に入ります。必要な境界条件は一つです。圧力・速度・温度を同時に固定すると over-specified になり、ソルバーが死にます。

LODI — 平面波だけ通す近似#

NSCBC の頭脳は Local One-Dimensional Inviscid (LODI) relations です。境界で 1 ステップだけ粘性項と横方向微分を無視できると仮定すると、保存方程式は特性振幅 $\mathcal{L}_i$ の ODE になります。

\frac{\partial \rho}{\partial t} + \frac{1}{c^2}\left[\mathcal{L}_2 + \tfrac{1}{2}(\mathcal{L}_5 + \mathcal{L}_1)\right] = 0

\frac{\partial u}{\partial t} + \frac{1}{2\rho c}(\mathcal{L}_5 - \mathcal{L}_1) = 0

\frac{\partial p}{\partial t} + \frac{1}{2}(\mathcal{L}_5 + \mathcal{L}_1) = 0

ここで $\mathcal{L}_1 = (u-c)\left(\partial p/\partial x - \rho c\, \partial u/\partial x\right)$ 、 $\mathcal{L}_5 = (u+c)\left(\partial p/\partial x + \rho c\, \partial u/\partial x\right)$ 、 $\mathcal{L}_2$ はエントロピー波です。

要点はこうです: 出ていく波 $\mathcal{L}_5,\,\mathcal{L}_2$ は内部セルからの風上差分で直接計算できます。入ってくる波 $\mathcal{L}_1$ だけが BC で決められます。この一行が NSCBC のすべてです。

完全非反射の罠 — drifting mean pressure#

最も単純な非反射 BC は入ってくる波をゼロにします。

\mathcal{L}_1 = 0

音波はきれいに出ていきます。しかし LODI を積分すると $\partial p/\partial t = -\mathcal{L}_5/2$ になり、 $\mathcal{L}_5$ の時間平均にわずかな正のバイアスがあれば $p$ が際限なく漂流します。コードは死にませんが、6 時間後には平均圧力がしれっと 1 atm を超えています。

Rudy & Strikwerda (1980) と Poinsot & Lele (1992) は、この漂流を抑える 1 次緩和項を提案しました。

\mathcal{L}_1 = K\,(p - p_\infty), \qquad K = \frac{\sigma\,(1 - M_{\max}^2)\,c}{L}

$p_\infty$ はターゲット圧力、 $L$ は領域の代表長さ、 $M_{\max}$ は境界での最大マッハ数です。 $\sigma$ は緩和強度で、推奨値は $\sigma \approx 0.25$ です。

Poinsot–Lele 緩和 — σ という一つの摘み#

$\sigma$ には二つの顔があります。

$\sigma$ が小さすぎる ( $\sim 0$ ): ほぼ完全非反射。平均の漂流は止まりません。
$\sigma$ が大きすぎる ( $\sim 5$ ): 境界が再び stiff になり、音波がまた反射します。

音響反射係数は主オーダーで $|R| \approx \sigma L/(c\tau)$ となります ( $\tau$ は波の通過時間)。LES では領域長さと平均マッハが決まった後、 $\sigma$ を 0.1–0.5 の範囲でチューニングします。

通過時間を領域横断 $\tau=L/c$ と置くと $L/(c\tau)=1$ になり、 $|R|\approx\sigma$ まで落ちます。推奨の $\sigma\approx 0.25$ は、この読みでは振幅反射約 0.25 です。 $\sigma$ を 0.1 から 0.5 へ動かすと $|R|$ も同じ区間を追います。 $K$ は $\sigma$ に比例するので、0.25 を 0.5 に倍にすると $K$ も倍です。 $(1-M_{\max}^{2})$ はマッハ補正です。 $M_{\max}=0$ なら因子は 1、 $M_{\max}=0.5$ なら $1-0.25=0.75$ で $K$ は 25% 小さくなります。 $M_{\max}$ が 1 に寄ると $K$ は 0 に消え、スライダーで入ってくる矢印が消える向きと同じです。

モードを Zero-gradient と NSCBC で切り替えます。 $\sigma$ のスライダーは、入ってくる波に掛かる係数です。

sigma (reflectivity)0.05

Zero-gradient ではガウシアンパルスが右境界でほぼそのまま跳ね返ります。反射係数は 1 の近くです。NSCBC にするとパルスは出ていきますが、 $\sigma$ を 0.5 より上げると残響がまた積もります。下の曲線は境界の $p'$ の時間履歴です。

Python — 音波パルスが出口に当たるとき#

二つのリーマン不変量 $\mathcal{A}^+ = u + p/(\rho c)$ 、 $\mathcal{A}^- = u - p/(\rho c)$ で書き直した 1D 線形化 Euler モデルです。トイ問題は、静止媒質に右向きガウシアンパルスを 1 発投げ、二つの BC を比較するというものです。

import numpy as np
 
def linearized_acoustic_pulse(bc='nscbc', sigma=0.05, N=240, n_steps=400):
    """1D 線形化 Euler。右端は音響流出境界。
    各 BC モードでの終端時刻の p'(x) を返す。"""
    # 格子 / 音速 / CFL
    c, cfl = 1.0, 0.5
    dx = 1.0 / N
 
    Ap = np.exp(-((np.linspace(0, 1, N) - 0.22) ** 2) / 0.003)  # 右向きパルス
    Am = np.zeros(N)                                             # 左向きは 0
 
    for _ in range(n_steps):
        Ap_new = Ap.copy()
        Am_new = Am.copy()
 
        # Ap は +c で進む → 左風上
        Ap_new[1:] = Ap[1:] - cfl * (Ap[1:] - Ap[:-1])
        Ap_new[0] = 0.0  # 静止流入
 
        # Am は -c で進む → 右風上
        Am_new[:-1] = Am[:-1] - cfl * (Am[:-1] - Am[1:])
 
        # 右境界 — NSCBC の心臓
        if bc == 'wall':
            # zero-gradient: 入ってくる Am は Ap と同値 → 完全反射
            Am_new[-1] = Ap_new[-1]
        else:
            # NSCBC: 入ってくる Am を反射率 sigma で固定
            Am_new[-1] = sigma * Ap_new[-1]
 
        Ap, Am = Ap_new, Am_new
 
    p_prime = 0.5 * (Ap - Am)   # p' (正規化単位)
    return p_prime
 
p_reflect = linearized_acoustic_pulse(bc='wall')
p_nscbc   = linearized_acoustic_pulse(bc='nscbc', sigma=0.05)
 
print(f"max |p'| (reflect): {np.max(np.abs(p_reflect)):.3e}")
print(f"max |p'| (nscbc):   {np.max(np.abs(p_nscbc)):.3e}")
# reflect: 0.998  /  nscbc: 0.025  → 2 桁の差

同じパルスに対して、zero-gradient はほぼそのまま保存 (反射) しますが、NSCBC は 2 桁減衰させます。 $\sigma$ を掃くと反射係数はほぼ線形に追従します。

コメントのピークは 0.998 と 0.025 です。比は約 40 です。コードは出ていく波に $\sigma=0.05$ を掛けますが、ピーク比は 0.05 より小さい 0.025 側に寄ります。パルスが領域を出て広がりながら低くなるからです。0.3 bar を 6 時間で割ると 0.05 bar/h です。 $\mathcal{L}_1=0$ の $\partial p/\partial t=-\mathcal{L}_5/2$ が、その傾きの符号を決めます。

実装に移す前のチェックリスト

まず境界の種類を数える。 亜音速流出 → BC 1 個、超音速流出 → 0 個、亜音速流入 → 4 個、超音速流入 → 5 個。上の図で確認。
σ のチューニングは領域長さと連動。 領域を 2 倍にすれば $K$ も変わる。絶対値で打ち込まないでください。
横方向項を漫然と落とさない。 3D LES では Yoo & Im (2007) の "relaxed" 形がより安定です。
衝撃が出口に当たる場合は NSCBC 単独では不十分。 バッファゾーン (sponge layer) と併用して非線形反射を吸収します。
粘性フラックスは別途処理。 NSCBC が自動でやってくれるわけではありません。

σ を 0 に置くと平均が歩き出す#

$\mathcal{L}_1=0$ は音波をきれいに出しますが、 $\partial p/\partial t=-\mathcal{L}_5/2$ が残ります。 $\mathcal{L}_5$ の時間平均が小さな正なら、圧力に上限はありません。6 時間後の 0.3 bar は、その積分がログに残した跡です。 $\sigma\approx 0.25$ を $K=\sigma(1-M_{\max}^{2})c/L$ に入れると漂流は引き戻されます。 $\tau=L/c$ の見積もりでは、振幅の約 1/4 も戻ります。

$\sigma$ を 5 付近まで上げると境界が再び硬くなり、パルスが跳ね返ります。0 と 5 のあいだで、コードの 0.05 と推奨の 0.25 は同じ直線上の二点です。ここでは、入ってくる波が一つだというところまで来ました。その係数をもう一段しめる話は、出口の σ が反射をどう固定するかにあります。

流出境界に音波を跳ね返させない方法 — Navier–Stokes 特性境界条件 (NSCBC)