P1輻射モデル
By M.Sato
1. 概要
P1輻射モデルは、計算流体力学（CFD）や熱流体解析において、 熱放射（輻射熱伝達）現象を高速かつ安定に計算するための標準的な数値解析モデルの一つである。
ここでは、P1モデルの理論的背景、数値解析における利点と制限、 境界条件の取り扱い、そして簡単な検証例題を紹介する。
2. 熱輻射輸送方程式
熱輻射輸送方程式（RTE: Radiative Transfer Equation）は 単位立体角当たりのエネルギーの収支を表す積分微分方程式である。 輻射のエネルギーは光子の流れとして表されるため、 連続体の温度場を解く方程式とは別の方程式として計算することになる。
$$ \begin{align} \frac{d I(\mathbf{r},\mathbf{s})}{ds} = \hat{\mathbf{s}}\cdot\nabla I = \kappa_\eta I_{b\eta} -\beta_\eta I_\eta + \frac{\sigma_{s\eta}}{4\pi} \int_{4\pi} I_\eta(\mathbf{r},\mathbf{s’}) \Phi(\mathbf{s},\mathbf{s’}) d\Omega' \end{align} $$
ここで$I$は放射強度、$\kappa_\eta$は吸収係数、$\sigma_{s\eta}$は散乱係数、 $\beta_\eta = \kappa_\eta + \sigma_{s\eta}$は減衰係数、$I_{b\eta}$は黒体放射強度、 $\Phi(\mathbf{s},\mathbf{s’})$は散乱相関関数である。 添字$\eta$は波長依存性を表している。 RTEは、空間座標$\mathbf{r}$と方向座標$\mathbf{s}$、 および波長$\eta$に依存する方程式になっている。
熱輻射輸送方程式を計算するために、いくつかの計算モデルが提案されている。 P1モデルは、この放射強度の空間分布を球面調和関数（Spherical Harmonics）の 第1次項（$P_1$項）までで展開（近似）することにより、方向依存性を簡略化し、 無数の光線を追跡することなく、拡散方程式に類似した形式にして計算を行う手法である。
3. P1近似モデル
3.1. 導出
位置$\mathbf{r}$、方向$\hat{\mathbf{s}}$における放射強度$I(\mathbf{r}, \hat{\mathbf{s}})$は、 球面調和関数の展開により次のように表すことができる。 $\hat{\mathbf{s}}$は球座標系で表した 方向余弦$\hat{\mathbf{s}}=(\sin\theta \cos\phi,\sin\theta \sin\phi,\cos\theta)$である。 $\theta$と$\phi$はそれぞれ方向ベクトル$\hat{\mathbf{s}}$を表す極角と方位角である。
なお以下では簡単のため波長依存は考えない（灰色近似）。また散乱効果も考えないこととする。
$$ \begin{align} I(\mathbf{r}, \hat{\mathbf{s}}) = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} A_l^m(\mathbf{r}) Y_l^m(\hat{\mathbf{s}}) \notag \end{align} $$
ここで、$A_l^m(\mathbf{r})$は位置依存の係数、$Y_l^m(\hat{\mathbf{s}})$は 角度依存性を表す係数であり、以下の正規化された球面調和関数で定義される。
$$ \begin{align} Y_l^m(\hat{\mathbf{s}}) = \left[ \frac{2l+1}{4\pi} \frac{(l-|m|)!}{(l+|m|)!} \right]^{1/2} e^{jm\phi}P_l^{|m|}(\cos\theta) \notag \end{align} $$
ここで$j=\sqrt{-1}$であり、$e^{jm\phi}$は調和関数$\cos(m\phi),\sin(m\phi)$を表している。 $P_l^{|m|}(\cos\theta)$は第一種の陪ルジャンドル多項式（degree=$l$、order=$m$）であり、 次のように与えられる。
$$ \begin{align} P_l^{|m|}(\mu) = \frac{(1-\mu^2)^{|m|/2}}{2^l, l!} \frac{d^{l+|m|}}{d\mu^{l+|m|}} \left(\mu^2 - 1\right)^l \notag \end{align} $$
$l$を無限個数まで取れば輻射強度$I$の厳密解が得られるが、 P1モデルでは$l$に関する級数和を$l=1$で打ち切る。 従って高波数成分の$Y_l^m$関数は表現できないため、 光学厚さが小さい場合には適切に放射強度を表すことができない。
$I(\mathbf{r}, \hat{\mathbf{s}})$を決定するためには、係数$A_l^m(\mathbf{r})$を決める必要がある。 そのために、位置$\mathbf{r}$における局所放射強度に方向余弦$\hat{\mathbf{s}}$を掛けて モーメント方程式を作成し、その値を評価する。
$$ \begin{align} I^0(\mathbf{r})&=\int_{4\pi} I(\mathbf{r}, \hat{\mathbf{s}}) d\Omega \notag \\ I^1(\mathbf{r})&=\int_{4\pi} \hat{\mathbf{s}}_{(1)} I(\mathbf{r}, \hat{\mathbf{s}}) d\Omega = \int_{4\pi} \sin\theta \cos\phi , I(\mathbf{r}, \hat{\mathbf{s}}) d\Omega \notag \\ I^2(\mathbf{r})&=\int_{4\pi} \hat{\mathbf{s}}_{(2)} I(\mathbf{r}, \hat{\mathbf{s}}) d\Omega = \int_{4\pi} \sin\theta \sin\phi , I(\mathbf{r}, \hat{\mathbf{s}}) d\Omega \notag \\ I^3(\mathbf{r})&=\int_{4\pi} \hat{\mathbf{s}}_{(3)} I(\mathbf{r}, \hat{\mathbf{s}}) d\Omega = \int_{4\pi} \cos\theta , I(\mathbf{r}, \hat{\mathbf{s}}) d\Omega \notag \end{align} $$
モーメント値を使えば放射強度は次式で表すことができる。
$$ \begin{align} 4 \pi I(\mathbf{r}, \hat{s}) = I^0(\mathbf{r}) + 3I^1(\mathbf{r}) \sin\theta \cos\phi + 3I^2(\mathbf{r}) \sin\theta \sin\phi + 3I^3(\mathbf{r})\cos\theta \notag \end{align} $$
最初の項$I^0$は方向に依存せず、2,3,4番目の項は方向ベクトル $\hat{\mathbf{s}} = \sin\theta\cos\phi,\hat{i} + \sin\theta\sin\phi,\hat{j} + \cos\theta,\hat{k}$の 方向成分に対応している。上式は以下のように書き直すことができる。
$$ \begin{align} I(\mathbf{r}, \hat{\mathbf{s}}) = a(\mathbf{r}) + \mathbf{b}(\mathbf{r}) \cdot \hat{\mathbf{s}} \end{align} $$
ここで$a(\mathbf{r})$はスカラー成分（等方性）を表し、 $\mathbf{b}(\mathbf{r})$は非等方性寄与を記述するベクトルである。
入射輻射$G$は、すべての方向にわたって積分された総エネルギーフラックス密度として定義される。
$$ \begin{align} G(\mathbf{r}) = \int_{4\pi} I(\mathbf{r}, \hat{\mathbf{s}}) , d\Omega \end{align} $$
(2)式を(3)式に代入すると
$$ \begin{align} G(\mathbf{r}) = a(\mathbf{r}) \int_{4\pi} d\Omega + \mathbf{b}(\mathbf{r}) \cdot \int_{4\pi} \hat{\mathbf{s}} , d\Omega = 4\pi a(\mathbf{r}) \notag \end{align} $$
ここで次式が成り立つことを利用した。
$$ \begin{align} \int_{4\pi} \hat{\mathbf{s}} , d\Omega = \int_0^{2\pi} \int_0^{\pi} \begin{pmatrix} \sin\theta \cos\phi \\ \sin\theta \sin\phi \\ \cos\theta \end{pmatrix} \sin\theta , d\theta , d\phi = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} = 0 \end{align} $$
ランベルトの余弦則を使うと輻射熱フラックスは次式で定義される。
$$ \begin{align} \mathbf{q}(\mathbf{r}) = \int_{4\pi} I(\mathbf{r}, \hat{\mathbf{s}}) \hat{\mathbf{s}} , d\Omega \notag \end{align} $$
(2)式を利用すると次式が得られる。
$$ \begin{align} \mathbf{q}(\mathbf{r}) = a(\mathbf{r}) \int_{4\pi} \hat{\mathbf{s}} , d\Omega + \mathbf{b}(\mathbf{r}) \int_{4\pi} \hat{\mathbf{s}}\hat{\mathbf{s}} , d\Omega = \frac{4\pi}{3} \mathbf{b}(\mathbf{r}) \notag \end{align} $$
ここで以下の関係式を利用した。
$$ \begin{align} \int_{4\pi} \hat{\mathbf{s}}\hat{\mathbf{s}} , d\Omega &= \int_0^{2\pi} \int_0^{\pi} \begin{pmatrix} \sin^2\theta \cos^2\phi & \sin^2\theta \sin\phi \cos\phi & \sin\theta \cos\theta \cos\phi \\ \sin^2\theta \sin\phi \cos\phi & \sin^2\theta \sin^2\phi & \sin\theta \cos\theta \sin\phi \\ \sin\theta \cos\theta \cos\phi & \sin\theta \cos\theta \sin\phi & \cos^2\theta \end{pmatrix} \sin\theta , d\theta , d\phi \notag \\ &= \int_0^{\pi} \begin{pmatrix} \pi \sin^2\theta & 0 & 0 \\ 0 & \pi \sin^2\theta & 0 \\ 0 & 0 & 2\pi \cos^2\theta \end{pmatrix} \sin\theta , d\theta \notag \\ &= \frac{4\pi}{3} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} = \frac{4\pi}{3} \delta \end{align} $$
ここで$\delta$は単位テンソルであり、$\mathbf{b} \cdot \delta = \mathbf{b}$である。 従って(2)式は入射放射$G$と輻射熱フラックス$\mathbf{q}$を使って次式のように表すことができる。
$$ \begin{align} I(\mathbf{r}, \hat{\mathbf{s}}) = \frac{1}{4\pi} \left[ G(\mathbf{r}) + 3\mathbf{q}(\mathbf{r}) \cdot \hat{\mathbf{s}} \right] \end{align} $$
放射強度$I$の表現式が$G$と$\mathbf{q}$の関数として得られたので、 散乱項や境界条件に関する条件を導くことができる。 以下では簡単のため、散乱係数$\sigma_s=0$とする。 RTE式に(6)式を代入すると次式が得られる。
$$ \begin{align} \frac{1}{4\pi} \nabla \cdot \left[ \hat{\mathbf{s}} \left( G + 3\mathbf{q} \cdot \hat{\mathbf{s}} \right) \right] + \frac{\kappa}{4\pi} \left( G + 3\mathbf{q} \cdot \hat{\mathbf{s}} \right) = \kappa I_b \notag \end{align} $$
上式を全立体角に渡って積分を行う。
$$ \begin{align} &\frac{1}{4\pi} \left[ \nabla \cdot G \int_{4\pi} \hat{\mathbf{s}} d\Omega + 3 \nabla \cdot \mathbf{q} \cdot \int_{4\pi} \hat{\mathbf{s}}\hat{\mathbf{s}} , d\Omega \right] + \frac{\kappa}{4\pi} G \int_{4\pi} d\Omega + \frac{3\kappa}{4\pi} \mathbf{q} \cdot \int_{4\pi} \hat{\mathbf{s}} , d\Omega = \kappa I_b \int_{4\pi} d\Omega \notag \end{align} $$
整理すると次式が得られる。
$$ \begin{align} \nabla \cdot \mathbf{q} = \kappa(4\pi I_b - G) \end{align} $$
次に、RTE式に方向余弦$\hat{\mathbf{s}}$を掛けて全立体角に渡って積分を行うと次式が得られる。
$$ \begin{align} &\frac{1}{4\pi} \left[ \nabla \cdot G \int_{4\pi} \hat{\mathbf{s}}\hat{\mathbf{s}} , d\Omega + 3 \nabla \cdot \mathbf{q} \int_{4\pi} \hat{\mathbf{s}}\hat{\mathbf{s}}\hat{\mathbf{s}} , d\Omega \right] + \frac{\kappa}{4\pi} G \int_{4\pi} \hat{\mathbf{s}} , d\Omega + \frac{3\kappa}{4\pi} \mathbf{q} \cdot \int_{4\pi} \hat{\mathbf{s}}\hat{\mathbf{s}} , d\Omega = \kappa I_b \int_{4\pi} \hat{\mathbf{s}} , d\Omega \notag \end{align} $$
ここで$\int_{4\pi} \hat{\mathbf{s}}\hat{\mathbf{s}}\hat{\mathbf{s}} , d\Omega = 0$であることを 利用すると最終的に次式が得られる。
$$ \begin{align} \frac{1}{3} \nabla \cdot (G\delta) + \kappa \mathbf{q} \cdot \delta = 0 \notag \end{align} $$
あるいは
$$ \begin{align} \mathbf{q} = -\frac{1}{3\kappa} \nabla G = -\Gamma \nabla G \end{align} $$
ここで$\Gamma=1/(3\kappa)$である。 最終的に未知変数4個に対して4個の方程式が得られる。
$$ \begin{align} \nabla \cdot \mathbf{q} &= \kappa(4\pi I_b - G) \notag \\ \nabla \cdot G &= -3\kappa \mathbf{q} \notag \end{align} $$
あるいは$\mathbf{q}$を消去して$G$に関する方程式にすると以下のようになる。 ただしここでは黒体の放射強度を$I_b=\frac{\sigma T^4}{\pi}$とした （$\sigma$はステファン・ボルツマン定数）。
$$ \begin{align} \nabla \cdot (\Gamma \nabla G) -\kappa G &= -4 \kappa \sigma T^4 \end{align} $$
これは$G$に関する拡散方程式（正確にはヘルムホルツ方程式）になっていることがわかる。
3.2. 特徴
P1モデルの特徴を以下に示す。
-
利点
- 熱輻射輸送方程式（RTE）を1変数（入射輻射$G$）の拡散方程式に変換するため、計算が高速
- 散乱効果、すなわち粒子（液滴、すす粒子など）の影響を考慮できる
- 光学厚さが大きい場合には適切な結果を得ることができる
-
制限
- すべての境界表面は拡散面であると仮定している
- 光学厚さが小さい場合には適切な結果を得ることができない
- 輻射熱フラックスを過大に評価する傾向がある
3.3. Marshak境界条件
P1モデルを使って輻射の計算を行う場合、どのように境界条件を設定すればよいかを考えてみる。 簡単のため、空間内に水平に平板が設置されており、上側が流体、下側が固体であるとする。 法線方向は下向きであり、拡散灰色壁（放射率$\varepsilon$、温度$T_w$）とする。
前項で説明したように$I$を1次までのルジャンドル展開で近似すると次式のようになる。 ここで$q$は$\mathbf{n}$方向の輻射熱フラックスである。
$$ \begin{align} I_{P1}(x, \mu) = \frac{1}{4\pi}\left[ G(x) + 3 q(x)\mu \right] \end{align} $$
P1モデルの放射強度$I$は$\mu$の1次式なので、 壁面での不連続な角度分布（壁から出る方向と、気体から来る方向）を 全方向で満たすことはできない。 そのため、壁から気体へ向かう半球（$\mu < 0$）で、角度分布そのものではなく、 $\mu$を重みにした半球積分（熱流束）を一致させるようにする。
拡散灰色壁から気体側へ出る放射$I_w$（$\mu < 0$の方向）は、壁自身の放射と反射の和で表される。
$$ \begin{align} I_w(\mu < 0) = \frac{\varepsilon \sigma T_w^4}{\pi} + \frac{(1-\varepsilon), q_{in}}{\pi} \notag \end{align} $$
ここで$q_{in}$は壁へ入射する熱流束（$\mu > 0$の半球の入射分）である。 壁から出る熱流束について以下の等式が成り立つとする。
$$ \begin{align} \int_{\mu < 0} I_{P1}, \mu , d\Omega = \int_{\mu < 0} I_w, \mu , d\Omega \notag \end{align} $$
図1 — 壁面における放射強度の角度分布
左辺はP1の分布を代入して以下のように計算できる。
$$ \begin{align} \int_{\mu < 0} I_{P1}, \mu , d\Omega = \frac{1}{2}\int_{-1}^{0} \left( G\mu + 3q\mu^2 \right) d\mu = -\frac{G}{4} + \frac{q}{2} \notag \end{align} $$
右辺は壁から出る熱流束の符号を負にしたものである。
$$ \begin{align} -q_{out} = -\left[ \varepsilon \sigma T_w^4 + (1-\varepsilon), q_{in} \right] \notag \end{align} $$
P1での入射熱流束（$\mu > 0$の半球）も同様に計算でき、$q_{in} = G/4 + q/2$となる。 これを代入すると次のようになる。
$$ \begin{align} -\frac{G}{4} + \frac{q}{2} = -\varepsilon \sigma T_w^4 - (1-\varepsilon)\left( \frac{G}{4} + \frac{q}{2} \right) \notag \end{align} $$
図2 — Marshak境界条件における熱流束の釣り合い
輻射熱フラックス$q$について整理する。
$$ \begin{align} \left(\frac{1}{2}+\frac{1-\varepsilon}{2} \right)q =-\varepsilon \sigma T_w^4 +\left(\frac{1}{4}-\frac{1-\varepsilon}{4} \right)G \notag \end{align} $$
$$ \begin{align} q = \frac{\varepsilon}{2(2-\varepsilon)} \left( G - 4\sigma T_w^4 \right) \notag \end{align} $$
これがMarshak境界条件で指定される関係式である。
3.4. OpenFOAMにおける実装
輻射熱フラックスは次式で表される。
$$ \begin{align} q = -\frac{1}{3\kappa} \frac{\partial G}{\partial n} = -\Gamma \frac{\partial G}{\partial n} \notag \end{align} $$
従ってMarshak境界条件は次式で表される。
$$ \begin{align} -\Gamma \frac{\partial G}{\partial n} = \frac{\varepsilon}{2(2-\varepsilon)}\left( G- 4\sigma T^4 \right) \notag \end{align} $$
これは、$G$と$\partial G/\partial n$を組合せたmixed型の境界条件になっていることがわかる。 face値に添字$w$、cell値に添字$c$をつけて離散化した表現式を記述すると以下のようになる。
$$ \begin{align} -\Gamma_c\frac{G_w - G_c}{\Delta x} = \frac{\varepsilon_w}{2(2-\varepsilon_w)}\left( G_w - 4\sigma T_w^4 \right) \notag \end{align} $$
変形すると次式が得られる。
$$ \begin{align} \left( \frac{\varepsilon_w}{2(2-\varepsilon_w)} + \frac{\Gamma_c}{\Delta x}\right)G_w = \frac{\varepsilon_w}{2(2-\varepsilon_w)} 4\sigma T_w^4 + \frac{\Gamma_c}{\Delta x}G_c \notag \end{align} $$
ここで$E=\varepsilon_w/\left(2(2-\varepsilon_w)\right)$と置くと、次式が得られる。
$$ \begin{align} G_w = \frac{1}{1+\frac{\Gamma_c}{E\Delta x}} 4\sigma T_w^4 + \left(1- \frac{1}{1+\frac{\Gamma_c}{E\Delta x}} \right) G_c \notag \end{align} $$
従って、OpenFOAMでは以下のようにmixed型境界条件のパラメータを設定する。
valueFraction= $\frac{1}{1+\frac{\Gamma_c}{E\Delta x}}$refValue= $4\sigma T_w^4$refGrad= $0$
4. 解析例
簡単な一次元モデルを使い、OpenFOAMのP1モデルの検証計算を行ってみた。
この解析では、温度$T$と入射輻射$G$のみを計算する （実際には温度は全領域で拘束指定するため計算しない）ことになるため、 ソルバー thermoFoam を使った。
4.1. 解析モデル
以下に計算に使用したモデルを示す。 無限の広がりを持つ平面壁が距離$L$だけ離れた位置で対向しているとした。
- 気体温度$T_g$（一様、拘束で固定）、壁温$T_w$、壁は黒体（$\varepsilon_w = 1$）
- 吸収係数$\kappa$（
absorptionCoeffとemissionCoeffは同値）、散乱なし - 座標$x$は層の中心を原点とし、$x \in [-L/2,\ L/2]$
- 光学厚さ$\tau = \kappa L$
- 参考パラメータ: $\kappa = 1\ \mathrm{m^{-1}}$、$L = 0.1,\ 1,\ 10\ \mathrm{m}$（$\tau = 0.1,\ 1,\ 10$）、$T_g = 1000\ \mathrm{K}$、$T_w = 300\ \mathrm{K}$
図3 — 解析モデル
4.2. メッシュ分割
x方向に100分割してメッシュを作成した。下図は$L=1$[m]の場合である。
図4 — メッシュ分割（$L=1$[m]）
4.3. 理論解
4.3.1. P1モデルを使った場合の理論解
P1モデルの支配方程式は、定常・散乱なしならば、次のようになる。
$$ \begin{align} \frac{d^2 G}{dx^2} = 3\kappa^2 \left( G - 4\sigma T_g^4 \right) \notag \end{align} $$
壁でのMarshak境界条件（黒体壁）は、壁が$x = L/2$のとき次の形になる。
$$ \begin{align} \frac{dG}{dx} = -\frac{3\kappa}{2} \left( G - 4\sigma T_w^4 \right) \notag \end{align} $$
$\Delta = T_g^4 - T_w^4$、$\lambda = \sqrt{3},\kappa$として解くと、次が得られる。
$$ \begin{align} G(x) = 4\sigma T_g^4 + C \cosh(\lambda x) \notag \end{align} $$
$$ \begin{align} C = -\frac{6\sigma\Delta}{D}, \qquad D = \sqrt{3},\sinh\!\left(\frac{\sqrt{3},\tau}{2}\right) + 1.5\cosh\!\left(\frac{\sqrt{3},\tau}{2}\right) \notag \end{align} $$
壁面の正味放射熱流束は次のとおりである。
$$ \begin{align} q_w = \frac{G_w - 4\sigma T_w^4}{2} = \sigma\Delta \left( 2 - \frac{3\cosh\!\left(\sqrt{3},\tau/2\right)}{D} \right) \notag \end{align} $$
4.3.2. 輻射輸送方程式の厳密解
$E_2,\ E_3$を指数積分関数として、次のように書ける。
$$ \begin{align} G(x) = 4\sigma T_w^4 + 4\sigma\Delta \left[ 1 - \frac{1}{2}\left( E_2\!\left(\kappa\left(\tfrac{L}{2} - x\right)\right) + E_2\!\left(\kappa\left(\tfrac{L}{2} + x\right)\right) \right) \right] \notag \end{align} $$
$$ \begin{align} q_w = \sigma\Delta \left( 1 - 2E_3(\tau) \right) \notag \end{align} $$
ここで、指数積分関数$E_n(x)$は次のように定義される。
$$ \begin{align} E_n(x) = \int_1^\infty \frac{e^{-xt}}{t^n},dt \notag \end{align} $$
指数積分関数はPythonの scipy.special.expn(2, ·) と expn(3, ·) を使うことで計算できる。
4.4. 計算結果
OpenFOAMの計算で得られた入射放射$G$のコンター図を示す。 これは$L=1$[m]（光学厚さ$\tau=1$）のときの結果である。
図5 — 入射放射$G$のコンター図（$\tau=1$）
次に、光学厚さ$\tau$を変化させたときの入射放射$G$の分布を示す。 左から$\tau = 0.1,\ 1,\ 10$である。 光学厚さが小さい場合、P1モデルは厳密解と大きく異なることがわかる。
図6 — 入射放射$G$の分布（左から$\tau = 0.1,\ 1,\ 10$）
以下に、輻射熱フラックス$q_w$の光学厚さ依存性を示す。 P1モデルは光学厚さが大きい場合には輻射熱フラックスを過大評価する傾向があることがわかる。
図7 — 壁面輻射熱フラックス$q_w$の光学厚さ依存性
5. まとめ
P1輻射モデルの誘導過程を示し、また典型的な境界条件として使われる Marshak境界について説明を行った。 また簡単な一次元モデルを使って検証計算を行い、解の特徴を調べた。
今回の解析で使用した例題データおよび検証用プログラムは以下からダウンロードできる。
⬇️ ダウンロード
計算結果データ: $\tau=0.1$ / $\tau=1$ / $\tau=10$
弊社の解析事例
弊社の流体解析事例については、下記のリンクからご覧ください。