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

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

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

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

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$ で移動します。

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

下の図でMach数を直接動かしてみてください。五本の特性線が境界をどちら向きに通過するかがリアルタイムで変わります。

$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$ の側へ引き寄せます。

下のシミュレーションで直接操作してみましょう。左が閉じた管、右が出口です。fire pulse で圧力波を撃ち、sigma スライダーを動かします。

見るべきは二つです。第一に、パルスが出口に届いた瞬間に赤い曲線(入ってくる波 $A^-$)が立ち上がるかどうか。これが反射です。第二に、source q を上げたまま下の圧力履歴を見ると、$\sigma = 0$ では白い曲線が $p_\infty$ の線に戻らず上側に留まります。反射は消したのに圧力を失った状態です。

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

$\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

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

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

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

境界で反射するのは物理的な音波だけではありません。波長が格子間隔の四倍より短い成分は物理的な解ではなく、離散化が作ったものです。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では消えずに残ります。

出口条件を選ぶときの一行

$\sigma$ はチューニングノブではなく、二つの失敗の間に置かれた座標です。左端には何にも縛られない圧力が、右端には共鳴管になった計算領域があります。0.25付近が勧められるのは、その位置で低周波だけを捕まえ、残りを通すからです。

出口で原因不明の振動に出会ったら、順序はこうです。まずその境界で入ってくる特性波の本数を数え、いま強制している条件数がその値と一致するか確認します。次に振動周期が $2L/c$ の倍数かどうかを見ます。一致すれば境界反射であって物理ではありません。最後に $\sigma$ を下げてみます。振動が減れば原因は境界であり、平均圧力が流され始めたなら反対側の失敗へ踏み込んだということです。