出口は開いているのに波が戻ってくる — NSCBCとσが決める反射

出口は開いています。それでも波は戻ってきます。圧縮性コードの出口境界を外挿で雑に処理して火炎や乱流を計算すると、領域の真ん中に原因不明の振動が育ちます。その周期を測ると、たいてい領域長さを音速で割った値と一致します。計算領域そのものが共鳴管になっているということです。共鳴を生む境界で、計算できる値とでっち上げる値を分けます。その値を動かす係数一つの代償も、同じ表に残します。

Mach スライダーを正の亜音速から音速の向こうへ動かすと、入っていた矢印が出ていきます。符号を負にすると、同じ面の条件数が跳ねます。

亜音速出口の条件は一つです。音速を超えると 0、流入にひっくり返すと 4 になります。矢印が領域の内側を向いている間だけ、その波はでっち上げた値です。

計算領域の端は物理ではない

物理的な領域に端はありません。燃焼器出口の先にも空間は続きます。しかし格子はどこかで終わらなければならず、最後のセルには隣がありません。隣がなければ微分ができず、微分ができなければ支配方程式を進められません。

RANSコードではこの問題が長く隠れていました。乱流粘性と人工粘性が大きいため、境界で誤って作られた波は数セル進むだけで消えます。LESやDNSになると事情が変わります。人工粘性はほぼゼロ、乱流粘性も最小です。境界が作った誤差は消えずに領域を横切って戻ってきます。

PoinsotとLeleが1992年にまとめた処方は、発想を反転させました。境界で変数を外挿する代わりに、境界を横切る波を数え、それぞれの振幅を決めるのです。Euler方程式に対する特性境界条件(ECBC)を粘性項のあるNavier–Stokesへ拡張したものなので、NSCBC(Navier–Stokes Characteristic Boundary Conditions)と呼ばれます。外挿の手続きは一行も使いません。

境界で数えられるものと、でっち上げるもの

境界が $x_1 = L$ にあるとします。$x_1$ 方向の項を波の形にまとめ直すと、連続の式はこうなります。

$\rho$ は密度、$u_i$ は速度成分、$d_1$ は境界に垂直な方向の寄与をまとめた項です。この $d$ ベクトルが特性解析の産物であり、その中に波の振幅 $\mathcal{L}_i$ が入っています。

$c$ は局所音速($c^2 = \gamma p / \rho$)、$p$ は圧力です。$\mathcal{L}_1$ は $x_1$ の負方向へ走る音波、$\mathcal{L}_5$ は正方向へ走る音波の振幅変化率です。残る三つは流体とともに運ばれます。$\mathcal{L}_2$ はエントロピー、$\mathcal{L}_3$ と $\mathcal{L}_4$ は接線方向速度 $u_2$、$u_3$ で、いずれも速度 $u_1$ で移動します。

分かれるのは速度の符号です。波が領域の外へ出るなら、その振幅は内部点の値から計算できます。中へ入ってくるなら、その情報は解の中に存在しません。でっち上げるしかありません。入ってくる波の本数が、そのままその境界に与えられる物理的境界条件の本数になります。

特性線が音波を領域全体でどう分けるかは、超音速流が上流を知らない場所にあります。ここでは、その線が出口の一面で切れる本数だけを数えます。

$M$ を $+0.3$ から $+1.4$ へ上げると、前で向きを変えた矢印が $\mathcal{L}_1$ です。条件数は 1 から 0。符号を負にすると流入口になり、条件は 4 です。亜音速出口で回っていたコードが超音速で発散するなら、この本数を超える条件をもう一つ掛けていることが多いです。

LODI — 入ってくる波の振幅を作る規則#

入ってくる振幅をでっち上げるには根拠が要ります。NSCBCは境界の各点で、接線方向項・粘性項・反応項をすべて落とした局所一次元非粘性系を立てます。これがLODI(Local One Dimensional Inviscid)関係です。

LODI関係は物理的な条件ではありません。実際に解く方程式でもありません。入ってくる $\mathcal{L}_i$ を推定するためだけに使います。手順は三段階です。物理条件が課された保存方程式を系から削り、削った方程式に対応するLODI関係で未知の $\mathcal{L}_i$ を既知の $\mathcal{L}_i$ で表し、残った方程式で他の変数を時間進行させます。

完全無反射出口はここで最も単純な選択をします。外から入ってくる音波はないと宣言するのです。

入ってくる波がゼロなら反射もゼロ。きれいに見えます。ところがこの式のどこにも外部圧力 $p_\infty$ がありません。

σ = 0 の代償: 圧力が p∞ に戻らない#

$p_\infty$ がないということは、境界が「今の圧力がいくらであるべきか」を知らないということです。領域内で熱が発生して圧力が上がっても、そのオフセットを戻す復元力がどこにもありません。問題が適切に定義されなくなります。

RudyとStrikwerdaの処方は、入ってくる波を完全にゼロに置く代わりに、圧力差に比例させて縛ることです。

$L$ は領域の代表長さ、$M$ は最大Mach数、$\sigma$ はこの処方唯一の自由パラメータです。$\sigma = 0$ なら完全無反射に戻ります。$\sigma$ を上げると境界が圧力を $p_\infty$ の側へ引き寄せます。

閉じた管の sigma と source q は、出力の表の直後に置きます。

二つの失敗にはさまれた狭い窓

$\sigma$ は互いに逆向きへ引っ張る二つの失敗の間にあります。

$\sigma$ が小さいと平均圧力が流されていき、大きいと境界が硬くなって音響エネルギーを投げ返します。PoinsotとLeleが実際に使った値は $\sigma \approx 0.25$、外挿ベースのB1に対応する係数では $\sigma = 0.58$ でした。これらは理論から導かれたものではなく、二つの失敗の間から選び出した妥協点です。

周波数依存性を見れば、なぜ妥協なのかがはっきりします。角振動数 $\omega$ の音波に対して、この境界の反射係数は

となります。低周波ほど反射が大きい。つまり $\sigma$ は「低い周波数は捕まえ、高い周波数は通す」フィルタです。平均圧力は $\omega \to 0$ の成分なので捕まり、外へ出したい音波は通ります。この分離が効く区間が狭い窓なのです。

コードで測る反射率と圧力オフセット

一次元線形音響なら、状態は二つの特性振幅にちょうど分かれます。$A^\pm = p' \pm \rho c u'$ がそれぞれ $\pm c$ で移動します。$\Delta t = \Delta x / c$ に取れば毎ステップ正確に一セルずらすだけなので、画面に見える揺れはすべて境界条件が作ったものです。

import numpy as np
 
N, C, L = 240, 1.0, 1.0
DX = L / N
DT = DX / C                       # ちょうど一セル移動 — スキーム拡散なし
 
 
def outlet_relax_k(sigma, mach=0.0):
    """NSCBC緩和係数 K = sigma (1 - M^2) c / L"""
    return sigma * (1.0 - mach ** 2) * C / L
 
 
def duct_step(ap, am, am_b, k_relax, q):
    """A+ は右へ一セル、A- は左へ一セル。出口では A- をでっち上げる。"""
    ap[1:] = ap[:-1].copy()
    ap[0] = 0.0
    am[:-1] = am[1:].copy()
 
    p_b = 0.5 * (ap[-1] + am[-1])
    am_b -= DT * k_relax * p_b     # L1 = K (p - p_inf)
    am[-1] = am_b
 
    ap[0] = am[0]                  # 左は閉じた端 (u = 0)
    ap += q * DT                   # 弱い一様な発熱
    am += q * DT
    return am_b
 
 
def measure_outlet(sigma, q, steps, pulse):
    ap, am, am_b = np.zeros(N), np.zeros(N), 0.0
    x = (np.arange(N) + 0.5) * DX
    if pulse:
        ap += np.exp(-((x - 0.30) / 0.09) ** 2)
    k = outlet_relax_k(sigma)
    refl = 0.0
    for n in range(steps):
        am_b = duct_step(ap, am, am_b, k, q)
        if pulse and n > 0.85 * N:
            refl = max(refl, np.abs(am[:-8]).max())
    return refl, float(np.mean(0.5 * (ap + am)))
 
 
for sigma in (0.0, 0.25, 1.0, 4.0, 10.0):
    r, _ = measure_outlet(sigma, q=0.0, steps=650, pulse=True)
    _, p = measure_outlet(sigma, q=0.3, steps=6000, pulse=False)
    print(f"sigma={sigma:5.2f}   反射={r * 100:5.1f}%   平均圧力-p_inf={p:+.4f}")

実行結果です。

sigma= 0.00   反射=  0.0%   平均圧力-p_inf=+0.3000
sigma= 0.25   反射=  1.9%   平均圧力-p_inf=+0.0007
sigma= 1.00   反射=  7.3%   平均圧力-p_inf=+0.0006
sigma= 4.00   反射= 24.5%   平均圧力-p_inf=+0.0227
sigma=10.00   反射= 46.8%   平均圧力-p_inf=-0.1142

fire pulse のパルスが出口に届いたとき、赤い曲線 $A^-$ が立ち上がれば反射です。source q を上げたまま $\sigma = 0$ だと、白い曲線は $p_\infty$ の上に留まります。反射はなく、圧力は手放した状態です。

$\sigma = 0$ は反射を完全に消しますが、発熱が作ったオフセット0.3をそのまま残します。$\sigma = 0.25$ では反射が2%未満で、圧力は $p_\infty$ から0.0007以内に抑えられます。ここまでは予想どおりです。

予想外なのは最後の二行です。$\sigma$ を4、10と上げると反射が25%、47%へ跳ね上がるのは当然として、圧力オフセットも再び悪化します。境界が硬すぎると自分で振動を作り、その振動が平均値を揺らすのです。$\sigma$ を上げることは、圧力を捕まえる代わりに反射を買う取引ではありません。ある点を越えると両方を失います。

出力の反射列は 0.0、1.9、7.3、24.5、46.8% です。コードは $C = 1$、$L = 1$、Mach 0 なので、$K = \sigma \left( 1 - M^2 \right) \frac{c}{L}$ は $K = \sigma$ に縮みます。σ を 0.25 から 10 へ、40 倍にすると反射は 1.9 から 46.8 です。約 25 倍。線形ではありません。$|R| = \frac{K}{\sqrt{K^2 + 4\omega^2}}$ は K が大きくなるほど 1 に飽和します。パルスは単一周波数ではなく、幅 0.09 のガウスです。一つの ω からこの百分率を逆算はしません。向きが合えば足ります。K を上げると反射は増え、そのあと鈍ります。

圧力列はその単調に従いません。$q = 0.3$ の実行で $\sigma = 0$ の平均は +0.3000、熱源と同じ大きさです。$\sigma = 0.25$ で +0.0007、$\sigma = 1$ で +0.0006。0.3000 に対する 0.0007 は 0.23% です。ここまでは緩和がオフセットをほぼ消します。$\sigma = 4$ では +0.0227。0.0007 の約 32 倍です。$\sigma = 10$ では −0.1142。符号が反転します。0.0007 に対する絶対値は約 160 倍です。反射 46.8% の境界は平均を捕まえられず、揺らします。表でオフセットが底なのは 0.25 と 1.0 です。Poinsot と Lele が使った 0.25 はその底の左端です。0.58 は B1 側の対応値なので、この列にはそのまま乗りません。

格子が作った波は逆向きに戻ってくる

境界で反射するのは物理的な音波だけではありません。波長が格子間隔の四倍より短い成分は物理的な解ではなく、離散化が作ったものです。PoinsotとLeleはこれを「q波」と呼び、物理的な「p波」と区別しました。

q波の特徴は群速度です。移動速度 $V$ の一次元移流方程式でも、短波長の群速度 $u_g$ は $V$ と符号が逆になります。流れは右へ進むのに、数値誤差は左へ遡ります。しかも $|u_g / V|$ はスキームの次数が高いほど大きくなります。高次スキームほど危ないということです。

ですから境界条件は二つの反射係数で評価しなければなりません。物理波の反射 $A_p / A_1$ と数値波の反射 $A_q / A_1$ です。使える境界条件ならどんな場合でも $A_q / A_1 \ll 1$ でなければならず、無反射をうたうなら $A_p / A_1 \ll 1$ まで満たす必要があります。計算開始時に初期場へ急峻な勾配を入れるだけでq波は生成され、DNSでは消えずに残ります。

周期が 2L/c なら、先に出口の σ を疑います#

$\sigma$ はチューニングのつまみではありません。一方の端は何にも縛られない圧力、もう一方は共鳴管になった領域です。0.25 付近は低周波だけを捕まえ、残りを通す位置です。

原因の見えない振動なら、入ってくる特性波の本数と、強制している条件数が一致するかを数えます。周期が $2L/c$ の倍数なら物理ではなく境界です。$\sigma$ を下げて振動が減れば原因は境界、平均圧力が流れ出せば反対側へ越えています。同じ出口を波の本数から数える手順は、NSCBC を無反射出口として組んだ前の整理にあります。