無反射境界条件

By M.Sato

September 8, 2026

1. 概要

圧縮性流体を計算する場合、流出境界条件を適切に設定することが重要である。非定常流れ場では解析領域内で発生した圧力波が流出境界に到達した場合、境界条件を適切に設定しないと、圧力波の反射が発生する。

特に近年は、DNS、LESなどの高精度の解析手法を使う機会が増えており、これらの解析手法は数値粘性が極めて小さくなっているため、流出境界で数値的に発生した反射波が解析結果に重大な誤差をもたらす可能性がある。

ここでは圧縮性流体における無反射境界条件の理論式を紹介し、 OpenFOAMにおける無反射境界条件の実装方法について説明を行った。また、いくつかの流出境界条件の設定に対して解析を実施し、それらが解析結果に与える影響についても紹介した。

2. 無反射境界の理論

2.1. LODI手法

圧縮性流体を記述する方程式はナビエ・ストークス方程式であるが、粘性の効果を省略した一次元のオイラー方程式を元にして無反射境界条件を導く。ここで紹介する手法は、一般的にLODI（Local One Dimensional Inviscid）手法として言及されている。

なお、ナビエ・ストークス方程式に対して適用する無反射境界条件は NSCBC（Navier-Stokes Characteristic Boundary Condition）と呼ばれており、 LODI法に粘性の境界条件を追加したものになっており様々な手法が提案されている。

圧縮性流体の基礎方程式を保存形式で表すと以下のようになる。

$$ \begin{align} \frac{\partial \mathbf{Q}}{\partial t} + \frac{\partial \mathbf{F}^k}{\partial x_k} =0 \notag \end{align} $$

$$ \begin{align} \mathbf{Q}= \begin{bmatrix} \rho \\ \rho u_1 \\ \rho u_2 \\ \rho u_3 \\ \rho e \end{bmatrix}, \quad \mathbf{F}^k= \begin{bmatrix} \rho u_k \\ \rho u_1 u_k + p \delta_{1k} \\ \rho u_2 u_k + p \delta_{2k} \\ \rho u_3 u_k + p \delta_{3k} \\ (\rho e+ p)u_k \end{bmatrix} \notag \end{align} $$

ここで$k=1,2,3$であり座標方向のインデックスを表している。 $\delta_{ij}$はクロネッカーのデルタ、$\rho$は密度、$u_1,u_2,u_3$は1,2,3方向の流速成分、 $p$は圧力、$\rho e$は流体の全エネルギであり次式で表される。

$$ \begin{align} \rho e = \frac{p}{\gamma-1} + \frac{1}{2} \rho u_j u_j \notag \end{align} $$

一方、基礎変数$\mathbf{U}$を使い非保存形式で表すと以下のようになる。なお以下のLODI式の誘導では座標軸方向$k=1$が流出境界に直交する方向、 $k=2,3$が接線方向と仮定している。

$$ \begin{align} \frac{\partial \mathbf{U}}{\partial t} + \mathbf{A}^k \frac{\partial \mathbf{U}}{\partial x_k} =0 \notag \end{align} $$

$$ \begin{align} \mathbf{U}= \begin{bmatrix} \rho \\ u_1 \\ u_2 \\ u_3 \\ p \end{bmatrix}, \quad \mathbf{A}^k = \begin{bmatrix} u_k & \delta_{1k}\rho & \delta_{2k}\rho & \delta_{3k}\rho & 0 \\ 0 & u_k & 0 & 0 & \delta_{1k}/\rho \\ 0 & 0 & u_k & 0 & \delta_{2k}/\rho \\ 0 & 0 & 0 & u_k & \delta_{3k}/\rho \\ 0 & \delta_{1k}\rho a^2 & \delta_{2k}\rho a^2 & \delta_{3k}\rho a^2 & u_k \end{bmatrix} \notag \end{align} $$

ここで$\mathbf{A}^k$はヤコビアン行列と呼ばれるものである。 $a$は音速を表しており、$a=\sqrt{\frac{\gamma p}{\rho}}$と表される。$\gamma$は比熱比である。

また保存量変数$\mathbf{Q}$と非保存量変数$\mathbf{U}$の関係は以下の変換行列$\mathbf{P}$で表される。なおここでは簡単のため、$\tilde{\gamma}=\gamma-1$とした。

$$ \begin{align} \mathbf{P} = \frac{\partial \mathbf{Q}}{\partial \mathbf{U}} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ u_1 & \rho & 0 & 0 & 0 \\ u_2 & 0 & \rho & 0 & 0 \\ u_3 & 0 & 0 & \rho & 0 \\ \frac{1}{2}u_j u_j & \rho u_1 & \rho u_2 & \rho u_3 & \tilde{\gamma} \end{bmatrix} \notag \end{align} $$

$$ \begin{align} \mathbf{P}^{-1} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ -u_1/\rho & 1/\rho & 0 & 0 & 0 \\ -u_2/\rho & 0 & 1/\rho & 0 & 0 \\ -u_3/\rho & 0 & 0 & 1/\rho & 0 \\ \frac{1}{2}u_j u_j \tilde{\gamma} & -\tilde{\gamma} u_1 & -\tilde{\gamma} u_2 & -\tilde{\gamma} u_3 & \tilde{\gamma} \end{bmatrix} \notag \end{align} $$

ヤコビアン$\mathbf{A}^1$の固有値は$\boldsymbol{\Lambda}=\mathrm{diag}(u_1-a,u_1,u_1,u_1,u_1+a)$であり、左固有ベクトルを求め、行列として表すと$\mathbf{S}_1$のようになる。

$$ \begin{align} \boldsymbol{\Lambda} = \mathbf{S}_1 \mathbf{A}^1 \mathbf{S}_1^{-1} \notag \end{align} $$

$$ \begin{align} \mathbf{S}_1 = \begin{bmatrix} 0 & -\rho a & 0 & 0 & 1 \\ a^2 & 0 & 0 & 0 & -1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & \rho a & 0 & 0 & 1 \end{bmatrix} \notag \end{align} $$

$$ \begin{align} \mathbf{S}_1^{-1} = \begin{bmatrix} 1/(2a^2) & 1/a^2 & 0 & 0 & 1/(2a^2) \\ -1/(2\rho a) & 0 & 0 & 0 & 1/(2\rho a) \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 1/2 & 0 & 0 & 0 & 1/2 \end{bmatrix} \notag \end{align} $$

これを元の式に代入すると次式が得られる。（ただし$k=2,3$）

$$ \begin{align} \frac{\partial \mathbf{U}}{\partial t} + \mathbf{S}_1^{-1} \mathcal{L} + \mathbf{A}^k \frac{\partial \mathbf{U}}{\partial x_k} =0 \notag \end{align} $$

流れが一次元的であると仮定し、$\partial/\partial x_2=\partial/\partial x_3=0$とすると次式が得られる。

$$ \begin{align} \frac{\partial \mathbf{W}}{\partial t} + \mathcal{L} =0 \notag \end{align} $$

ここで、特性変数$\mathbf{W}$の微小増分$\delta \mathbf{W}$は次式で表される。

$$ \begin{align} \delta \mathbf{W} = \mathbf{S}_1 \delta \mathbf{U} = \begin{bmatrix} \delta p - \rho a \delta u_1 \\ a^2 \delta \rho - \delta p \\ \delta u_2 \\ \delta u_3 \\ \delta p + \rho a \delta u_1 \end{bmatrix} \notag \end{align} $$

$$ \begin{align} \mathcal{L} = \mathbf{S}_1 \mathbf{A}^1 \frac{\partial \mathbf{U}}{\partial x_1} = \begin{bmatrix} \mathcal{L}_1 \\ \mathcal{L}_2 \\ \mathcal{L}_3 \\ \mathcal{L}_4 \\ \mathcal{L}_5 \end{bmatrix}= \begin{bmatrix} (u_1-a)\left( \frac{\partial p}{\partial x_1} - \rho a \frac{\partial u_1}{\partial x_1}\right) \\ u_1\left( a^2\frac{\partial \rho}{\partial x_1} - \frac{\partial p}{\partial x_1}\right) \\ u_1 \frac{\partial u_2}{\partial x_1} \\ u_1 \frac{\partial u_3}{\partial x_1} \\ (u_1+a)\left( \frac{\partial p}{\partial x_1} + \rho a \frac{\partial u_1}{\partial x_1}\right) \end{bmatrix} = \boldsymbol{\Lambda}\frac{\partial \mathbf{W}}{\partial x_1} \notag \end{align} $$

これより、基礎変数$\mathbf{U}$の時間微分$\partial \mathbf{U}/\partial t$は次式で表される。

$$ \begin{align} \frac{\partial }{\partial t} \begin{bmatrix} \rho \\ u_1 \\ u_2 \\ u_3 \\ p \end{bmatrix} + \begin{bmatrix} \frac{1}{a^2}\left( \mathcal{L}_2 + \frac{1}{2}(\mathcal{L}_5 + \mathcal{L}_1 )\right) \\ \frac{1}{2 \rho a}(\mathcal{L}_5 - \mathcal{L}_1 ) \\ \mathcal{L}_3 \\ \mathcal{L}_4 \\ \frac{1}{2} (\mathcal{L}_5 + \mathcal{L}_1 ) \end{bmatrix} = 0 \notag \end{align} $$

あるいは、次式のように空間微分の関係式として記述することもできる。

$$ \begin{align} \frac{\partial }{\partial x} \begin{bmatrix} \rho \\ u_1 \\ u_2 \\ u_3 \\ p \end{bmatrix}= \begin{bmatrix} \frac{1}{a^2}\left( \frac{\mathcal{L}_2}{u_1} + \frac{1}{2}\left(\frac{\mathcal{L}_5}{u_1+a} + \frac{\mathcal{L}_1}{u_1-a} \right)\right) \\ \frac{1}{2 \rho a}\left(\frac{\mathcal{L}_5}{u_1+a} - \frac{\mathcal{L}_1}{u_1-a} \right) \\ \frac{\mathcal{L}_3}{u_1} \\ \frac{\mathcal{L}_4}{u_1} \\ \frac{1}{2} \left(\frac{\mathcal{L}_5}{u_1+a} + \frac{\mathcal{L}_1}{u_1-a} \right) \end{bmatrix} \notag \end{align} $$

これらの関係式を利用することで適切な境界条件を設定することが可能となる。

2.2. 無反射のための条件

特性量$\mathbf{W}$は特性速度$u_1-a$、$u_1+a$、$u_1$で輸送される。亜音速流れの場合、$u_1-a<0$であることから$\mathcal{L}_1$は解析領域外から解析領域内部に伝播する波動となっていることがわかる。

図1 — 流出境界における特性波の伝播方向

特性量$\mathbf{W}$の輸送方程式より、特性波の振幅$\mathcal{L}_i$に様々な条件を設定することで、適切な境界条件を設定できることがわかる。例えば、

流出境界において圧力$p=0$と指定する場合、$\partial p/\partial t=0$より、$\mathcal{L}_5=-\mathcal{L}_1$とする
流出境界において完全無反射条件とするには$\mathcal{L}_1=0$とする

$\mathcal{L}_1=0$とした場合、前節の式に基づいて次式が得られる。これらの式は基礎変数$p,u_1,u_2,u_3$に対する移流方程式になっている。

$$ \begin{align} &\frac{\partial p}{\partial t} + (u_1+a) \frac{\partial p}{\partial x_1} =0 \notag \\ &\frac{\partial u_1}{\partial t} + (u_1+a) \frac{\partial u_1}{\partial x_1} =0 \notag \\ &\frac{\partial u_2}{\partial t} + u_1 \frac{\partial u_2}{\partial x_1} =0 \notag \\ &\frac{\partial u_3}{\partial t} + u_1 \frac{\partial u_3}{\partial x_1} =0 \notag \end{align} $$

なお、圧力を完全無反射として境界条件指定した場合（$\mathcal{L}_1=0$）、圧力の勾配しか指定していないため、圧力がドリフトにより変動する可能性がある。 $\mathcal{L}_1=K(p-p_\infty)$とすることで圧力値の不定性を除去することが可能であるが、この場合、圧力は部分的に反射する条件となる。

💡 輸送量$u_1,u_2,u_3$（流速成分）に関して見ると、$u_1$は移流速度が特性速度$u_1+a$になっているのに対して、$u_2,u_3$は移流速度を特性速度$u_1$とする必要があることがわかる。

3. OpenFOAMにおける無反射境界

OpenFOAMでは、流出境界における圧力波の反射を抑制するために、 advective境界条件、およびwaveTransmissive境界条件が用意されている。これらは基本的にLODI関係式を使った境界条件である。

advective境界は特性速度として境界面の「法線方向速度$u$」を用いている。 waveTransmissive境界はadvective境界を継承したクラスで定義されており、特性速度として「法線方向速度＋音速$u+a$」を用いている。

3.1. advective境界

これは以下の条件式を境界に適用する。

$$ \begin{align} \frac{\partial \phi}{\partial t} + V_n \left( \frac{\partial \phi}{\partial n} + \frac{\phi_P- \phi_\infty}{L_\infty}\right)=0 \notag \end{align} $$

ここで、$\phi_P$は境界面上における変数値、$V_n$は境界面の法線方向速度、 $L_\infty$は境界面から遠方場までの距離を表している。時間微分項、空間微分項を離散化して表すと次式が得られる。

$$ \begin{align} \frac{\phi_P^{n+1}-\phi_P^{n}}{\Delta t} + V_n \left( \frac{\phi_P^{n+1}-\phi_C^{n+1}}{\Delta x} + \frac{\phi_P^{n+1} - \phi_\infty}{L_\infty}\right)=0 \notag \end{align} $$

ここで上付きの添字$n,n+1$は時間ステップを表している。また$\phi_C$は境界面に接するセルの変数値を表している。 $\Delta t$は時間刻み幅、$\Delta x$は境界面に接するセルの中心から境界面までの距離を表している。上式を整理すると次式が得られる。

$$ \begin{align} \phi_P^{n+1} &=\frac{ \phi_P^{n} + \alpha \phi_C^{n+1}+ k \phi_\infty }{ 1 + \alpha +k } \notag \\ &= \frac{ \phi_P^{n} + k \phi_\infty }{ 1 + \alpha +k } + \frac{ \alpha \phi_C^{n+1}}{ 1 + \alpha +k } \notag \end{align} $$

ここで、$\alpha=\frac{V_n \Delta t}{\Delta x}$、$k=\frac{V_n \Delta t}{L_\infty}$である。

advective境界条件は、mixed境界条件に属しており、境界条件は次の一般式で表される。これは固定値を指定するDirichlet境界条件と、勾配値を指定するNeumann境界条件の両方を組み合わせたものである。

$$ \begin{align} \phi_P = w \phi_{F} + (1-w) ( \phi_C + \mathbf{d} \cdot \nabla \phi) \notag \end{align} $$

ここで、$w$は重み係数、$\phi_{F}$はDirichlet境界条件で指定する固定値、 $\mathbf{d}$は境界面から接するセル中心までの距離ベクトルを表している。 $\nabla \phi$は境界面に接するセルの勾配値を表している。 $w,\phi_F,\nabla \phi$はOpenFOAMのソースコード内ではそれぞれ valueFraction、refValue、refGradという変数名で参照されており、これらの値を適切に設定することで境界条件を設定している。

advective境界条件は、以下のように変形することができる。

$$ \begin{align} \phi_P = \frac{1+k}{1+\alpha+k} \frac{\phi_P^n +k \phi_\infty}{1+k} + \left( 1- \frac{1+k}{1+\alpha+k} \right) \left(\phi_C^{n+1} + \mathbf{d} \cdot 0 \right) \notag \end{align} $$

ここで

$$ \begin{align} &\mathrm{valueFraction} = \frac{1+k}{1+\alpha+k} \notag \\ &\mathrm{refValue} = \frac{\phi_P^n +k \phi_\infty}{1+k} \notag \\ &\mathrm{refGrad} = 0 \notag \end{align} $$

3.2. waveTransmissive境界

advective境界では特性速度が境界における法線速度$V_n$であるが、 waveTransmissive境界では特性速度が$V_n+a$となっている。ここで$a$は音速であり、$a=\sqrt{\frac{\gamma}{\psi}}$、$\psi=\frac{d \rho}{d p}$で表される。それ以外はadvective境界と同じである。

3.3. LODI2D境界

これは、こちらの資料に紹介されている境界条件である。

前述のようにLODI関係式を使った反射波抑制のためには、流出の法線方向成分にはwaveTransmissive境界条件、接線方向成分にはadvective境界条件を使う必要がある。 LODI2D境界条件は、境界における局所座標系に基づき、局所法線方向、局所接線方向の流速成分を自動的に計算し、それぞれの成分に応じて境界条件を切り替え、その後自動的に全体座標系に戻すようにしている。

このLODI2D境界を発展させて3次元モデルに対しても適用できるようにした研究結果はこちらに紹介されている。ただし任意の3次元形状には対応しておらず、流出境界がx軸に直交する平面に限定されている。

4. 解析例

2次元モデルを使い、流出境界条件を変更した場合の計算結果を比較してみた。

4m×4mの矩形領域の中央に直径0.2mの円形領域（spark）を設置した。初期条件として、時刻ゼロにおいて全領域で流速$U=0$、$p=101325$[Pa]とした。温度については、円形領域で初期温度3000[K]、それ以外は300[K]とした。

これは中心の円形領域において急激に温度が上昇した場合に相当し、圧力波が周囲に伝播することになる。この際に流出境界における圧力波の反射を調べた。流体は空気とし、密度は理想気体の式に従うとした。解析は時刻$t=0.01$[s]まで実施した。

図2 — 解析モデルと初期条件

対称性を考慮し、1/4の領域を解析範囲とした。解析に使用したメッシュを以下に示す。メッシュは四角形セルで構成されており、xy方向にそれぞれ100分割とした。

図3 — メッシュ分割

4.1. 解析ケース

解析に使用した流出境界条件（境界outletに適用）の一覧を示す。

dof	case1	case2	case3	case4
pressure	WT	WT	WT	WT
velocity	PIOV	WT	adv	LODI2D
temperature	adv	adv	adv	adv

ここで、WTはwaveTransmissive境界、PIOVはpressureInletOutletVelocity境界、 advはadvective境界、LODI2DはLODI2D境界を表している。

圧力については全ケースにおいてwaveTransmissive境界、また温度については全ケースにおいてadvective境界を使用している。

流速場については、前述のように本来、境界法線成分、境界接線成分について異なる特性速度を指定する必要があるが、標準のOpenFOAMにはそのような機能は実装されていない。従ってWTおよびadvは全流速成分について同一の「無反射条件」が適用されることになる。 LODI2Dは、カスタマイズされた境界条件であり、任意の流出境界に対して自動的に法線方向、接線方向を区別して境界条件を設定するようになっている。

なお、waveTransmissive境界、advective境界、LODI2D境界では遠方場を指定した。いずれのケースにおいても$L_\infty=10$、$U_\infty=0$、$p_\infty=101325$と指定した。

4.2. 解析結果

時刻$t=0.005$[s]における流速ベクトルおよび圧力分布を以下に示す。波動が流出境界に到達するまでは、いずれのケースにおいても、ほぼ同心円状に波動が伝播している。

図4 — 時刻$t=0.005$[s]における流速ベクトルと圧力分布

次にcase1～case4の最終時刻$t=0.01$[s]における圧力分布、流速強度分布を示す（上段が圧力、下段が流速強度）。

図5 — 時刻$t=0.01$[s]における圧力分布（上段）と流速強度分布（下段）

case1は通常の流体解析で一般的に使われるpressureInletOutletVelocity境界条件であり、反射の効果を考慮していない。最終時刻において若干の反射波による変動が発生していることがわかる。

case2は流速場においてwaveTransmissive境界条件を使った場合であるが、ほとんど反射波が発生していないことがわかる。

case3は流速場においてadvective境界条件を使用している。このケースは最も反射波の影響が大きくなった。

case4は流速場においてLODI2D境界条件を使ったものであり、流出境界における法線、接線流速成分に応じて、それぞれwaveTransmissive境界条件、advective境界条件を切り替えている。これは理想の境界条件と考えられるが、case4はcase2と同様に反射波の影響が小さいことがわかる。

各ケースの圧力コンター、流速強度コンターのアニメーションを以下に示す。

動画1 — case1（PIOV）圧力コンター

動画2 — case1（PIOV）流速強度コンター

動画3 — case2（waveTransmissive）圧力コンター

動画4 — case2（waveTransmissive）流速強度コンター

動画5 — case3（advective）圧力コンター

動画6 — case3（advective）流速強度コンター

動画7 — case4（LODI2D）圧力コンター

動画8 — case4（LODI2D）流速強度コンター

5. まとめ

無反射境界の理論的背景について説明を行った。またOpenFOAMに実装されている無反射境界条件を使いテスト計算を行った。

今回のテスト計算の条件では、流速$U$の境界条件として waveTransmissive、あるいはLODI2Dを使うのが望ましいことが分かった。

今回の解析で使用したソースコードおよび例題データは以下からダウンロードできる。

⬇️ ダウンロード

ソースコード＆例題データ

参考文献

T. J. Poinsot, S. K. Lele, “Boundary conditions for direct simulations of compressible viscous flows”, Journal of Computational Physics, Vol.101, No.1, pp.104-129, 1992.
L. Lucchese, “Implementation of a non-reflecting boundary condition”, CFD with OpenSource Software, Chalmers University of Technology, 2022.
B. Jarfors, “Non-reflecting boundary conditions in three dimensions”, CFD with OpenSource Software, Chalmers University of Technology, 2024.

弊社の解析事例

弊社の流体解析事例については、下記のリンクからご覧ください。

👉 解析事例のご紹介 OpenFOAM

👉 流体解析実績及び事例