OpenFOAMを使った溶岩流の解析

By M.Sato

January 13, 2026

概要

今回は、OpenFOAMを使った溶岩流の解析例を紹介する。

📖 参考にした元の論文はここにある。

溶岩の流れは、温度を考慮して流動場と自由表面を計算する必要がある。

OpenFOAMには、等温場のニ相流を計算するためにinterFoamというソルバーが用意されているが、温度場の計算は含まれていない。ここでは、元論文の著者に習い、 interFoamを拡張し温度場の計算を追加したinterThermRadConvFoamというソルバーを作成した。

🔧 元論文の著者が作成したソルバーはv1912用であり、最新版のOpenFOAMでは動作しないため、 最新版(v2506)で動作するように修正を行った。

解析のモデル化について紹介を行った。いくつかの解析を実施し、新しく作成したソルバーはv1912用のソルバーと同等の計算ができることを確認した。

図：溶岩流の解析モデル概念図（元論文のFig.3より引用）

基礎方程式

解析に使用している基礎方程式を示す。

$$ \begin{align} \frac{\partial \alpha}{\partial t} + \nabla \cdot (\mathbf{u} \alpha) &= 0 \\ \nabla \cdot \mathbf{u} &= 0 \\ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \mathbf{u}^T) &= -\nabla p + \nabla\cdot \mathbf{\tau} + \rho \mathbf{g} + \mathbf{f}_\Sigma \\ \frac{\partial \rho c_p T}{\partial t} + \nabla \cdot (\rho c_p T \mathbf{u} ) - \chi_\Sigma \nabla^2 (kT) &= -\frac{\epsilon \sigma_{SB} f A_{fs}}{V} (T^4 - T_{env}^4)-\frac{h f A_{fs}}{V} (T - T_{env}) \end{align} $$

ここで、

$\alpha$：溶岩の体積分率
$\mathbf{u}$：速度ベクトル
$t$：時間
$\rho$：密度
$p$：圧力、
$\mathbf{\tau}=\mu\left[\nabla\mathbf{u}+(\nabla\mathbf{u})^T\right]$：粘性応力テンソル
$\mathbf{g}$：重力加速度ベクトル

である。 $\mathbf{f}_\Sigma$は表面張力に起因する体積力であり、以下の式で表される。

$$ \mathbf{f}_\Sigma = \sigma_\Sigma \kappa_\Sigma \nabla \alpha,\quad \kappa_\Sigma = - \nabla \cdot \left( \frac{\nabla \alpha}{|\nabla \alpha|} \right) $$

$\kappa_\Sigma$：界面の曲率
$\sigma_\Sigma$：表面張力係数

を表す。なお以下の解析では表面張力の効果は無視している。

$c_p$：比熱
$T$：温度
$k$：熱伝導率
$\epsilon$：放射率、$\sigma_{SB}$はステファン・ボルツマン定数

である。

$f$：界面においてlava-inner-coreが露出している割合を表すための変数であり以下の解析では1.0と仮定している。\
$T_{env}$：周囲温度
$h$：対流熱伝達率
$\chi_\Sigma$：界面位置を示す指示関数であり、界面での熱拡散を抑制するために使われる。その値は、界面上のセルでゼロ、それ以外のセルで１の値を持つ関数である。

VOF法を使って計算することになるため、流体の物性値は各相の体積分率に基づいて以下のように計算される。

$$ \begin{align} \rho &= \alpha \rho_{lava} + (1-\alpha) \rho_{air} \notag \\ \mu &= \alpha \mu_{lava} + (1-\alpha) \mu_{air} \notag \\ \rho{c_p} &= \alpha \rho_{lava}{c_{p,lava}} + (1-\alpha) \rho_{air}{c_{p,air}} \notag \\ k &= \alpha k_{lava} + (1-\alpha) k_{air} \notag \end{align} $$

粘度式

溶岩の粘度には非アレニウス型のVogel-Fulcher-Tammann (VFT)式を用いている。これは以下の式で表される。

$$ \begin{align} log(\mu) &= A + \frac{B}{T - C} \notag \end{align} $$

ここで、$\mu$は粘度[Pa.s]、$T$は絶対温度[K]である。 $A,B,C$は溶岩の組成に依存しており実験データに基づいて決定される定数である。

溶岩表面からの熱移動のモデル化

溶岩流が空気と接する部分では、溶岩表面から空気中に熱が放出される。

この熱放出の原因は、熱輻射による効果と対流熱伝達による効果がある。

また溶岩流が地面と接触する部分では、地面に熱が移動することになる。

もし溶岩流表面が解析領域の境界になっている場合、直接熱フラックス境界条件を指定できるが、今回の計算モデルではVOF法を使っているため、溶岩表面は解析領域内部に含まれことになる。

そのため、溶岩流が空気と接している部分では、流体のエネルギー保存式に等価な熱ソース項として指定した。

図：溶岩表面からの熱移動モデル

地面からの熱フラックス

地面からの熱移動を考慮するために熱伝達境界を設定した。

元論文では、下図に示すように地面境界に対してexprMixedタイプの境界条件を使用しているが、境界条件の表現式が正しくないと思われるため、ここではcodedMixedタイプの境界条件に書き直して使用した。

⚠️ 元論文の問題点：

熱伝導率をPr数から逆算する式になっているが粘度が含まれていない

流体の相分率を考慮した熱伝導率になっていない

mag(pos())は原点からの距離であり壁面からの距離ではない

図：元論文で使用されているexprMixed境界条件

mixedタイプの境界条件の設定方法を説明した図を以下に示す。

図：mixedタイプ境界条件の設定方法

対流熱伝達効果の考慮

溶岩表面(空気との界面)からの対流熱フラックスは次式で表される。

$$ q_{conv} = h (T - T_{env}) $$

ここで、

$q_{conv}$：単位面積あたりの対流熱フラックス
$h$：対流熱伝達率
$T$：溶岩の表面温度
$T_{env}$：周囲環境の温度

を表す。熱フラックスにface面積をかけてセル体積で割ることにより、溶岩自由表面セルに対する等価な発熱項を得ることができる。

この発熱項を線形近似してエネルギー保存式に組み込む。

$$ \begin{align} Q_{conv} &= \frac{ h f A_{fs} }{V} (T - T_{env}) \notag \end{align} $$

この発熱項は線形式となっており、そのままエネルギー保存式に組み込むことができる。

熱輻射効果の考慮

溶岩表面(空気との界面)からの輻射熱フラックスは次式で表される。

$$ \begin{align} q_{rad} = \epsilon \sigma_{SB} (T^4 - T_{env}^4) \notag \end{align} $$ ここで、

$q_{rad}$：単位面積あたりの放射熱フラックス
$\epsilon$：溶岩の放射率
$\sigma_{SB}$：ステファン・ボルツマン定数
$T$：溶岩の表面温度
$T_{env}$：周囲環境の温度

を表す。熱フラックスにface面積をかけてセル体積で割ることにより、溶岩自由表面セルに対する等価な発熱項を得ることができる。

この発熱項を線形近似してエネルギー保存式に組み込む。

$$ \begin{align} Q_{rad} &= \frac{ \epsilon \sigma_{SB}f A_{fs} }{V} \left( T^4 - T_{env}^4 \right) \notag \\ &=\beta (T_{env}^4 + 3 T^4)-4 \beta T^3 T \notag \end{align} $$

ここで $\beta=\frac{\epsilon\sigma_{SB}f A_{fs}}{V}$である。

解析手法

溶岩流を解析するためには、２流体（溶岩と空気）の界面を追跡する必要がある。高温の溶岩は流出後は周囲温度によって冷却される。

溶岩の粘度は温度に依存するため、温度場の計算も必要になる。

ここでは、元論文の著者に習い、二流体の等温場の解析追跡ソルバー interFoamを拡張し、温度場の計算を追加したinterThermFoamというソルバーを作成した。

なお、標準版のOpenFOAMには、compressibleInterFoamというソルバーがあり、圧縮性流体として同等の解析を行うことができると考えられるが、今回作成したソルバーは二相の流体をそれぞれ非圧縮性流体として扱っている点が異なる。

ソースプログラムの修正

元論文の著者はOpenFOAM v1912を使っているが、現在のOpenFOAMの最新版で動作するようにソースプログラムを修正した。

主な修正点は以下の通りである。

熱伝達境界の実装： 前述のように地面と接触する部分では、熱伝達境界を設定する。前述のように元論文では、OpenFOAMのexprMixedタイプの境界条件を使っているが、正しく実装されていないと思われる。ここではcodedMixed境界条件に変更し、正しく熱伝達境界条件を実装した。
流体の熱伝導率を入力： 元論文では、流体の熱伝導率をPr数から逆算しているが、粘度が温度関数となっているため熱伝導率は粘度に比例するような表現式になる。実際には溶岩の粘度と熱伝導率は独立のパラメータとして指定すべきであるため、ここでは流体の熱伝導率を直接入力するように修正した。