😋

catalyticFoamのチュートリアル

に公開

この記事について

本記事では、実務で活用できそうなOSSツールを実際に試して、その使用感を備忘録として記録していきます。
今回取り上げるのはcatalyticFoamというOpenFOAMベースのソルバーです。OpenFOAM上で不均一触媒反応(壁面・表面反応)を「詳細な表面マイクロキネティクス」と連成して解くためのソルバーで、触媒反応の解析に利用可能です。

この記事で分かること
 - 前回:catalyticFoamの概要と特徴
 - 前回:インストール手順
 - 今回:チュートリアルケースの実行
 - 次回以降:化学プロセス解析への応用可能性(余力があれば...)

この記事は実際に手を動かしながら検証した内容を記載しています。生成AIも利用してコード解析を行っているので、誤りがあるかもしれません。間違いを見つけた方は、ぜひ優しく教えてください。

チュートリアル調査概要

触媒表面反応を伴う流体解析ソルバー catalyticFoam の使用例として、主に4つのケースが含まれています。それぞれの特徴と設定内容を以下に詳述します。

catalyticFoam チュートリアルケース一覧

ケース名 計算領域
メッシュ
気相反応 表面反応 ISAT
(高速化手法)
触媒負荷
\alpha_{\text{cat}}
反応モデル 主な解析目的
case 2D 平板チャネル
(blockMesh)
OFF ON OFF 10.0 Rh 詳細反応機構 基本的な触媒表面反応の解析
case
_homohet
2D 平板チャネル
(blockMesh)
ON ON OFF 10.0 Rh 詳細反応機構 気相反応と表面反応の連成解析
case_isat 2D 平板チャネル
(blockMesh)
OFF ON ON 10.0 Rh 詳細反応機構 ISAT法による計算高速化の検証
sphere 3D 球体周り
(snappyHexMesh)
OFF ON OFF 5.0 WGS (水性ガスシフト)
UBIモデル
複雑形状への適用とメッシュ生成

各項目の解説

  • 気相反応 (Homogeneous): 流体内部で起こる化学反応(燃焼など)。case_homohetのみ有効で、流体中での化学種変化も計算されます。
  • 表面反応 (Heterogeneous): 触媒壁面上 (reactingWall) での吸着・反応・脱離プロセス。全てのケースで有効になっています。
  • ISAT (In Situ Adaptive Tabulation): 化学反応計算の結果を動的に表化して再利用し、計算時間を短縮する手法です。case_isatでは、不均一反応に対して許容誤差 5e-3 で有効化されています。
  • 触媒負荷 (\alpha_{\text{cat}}): constant/solverOptions内の alfaCatalyst で設定されるパラメータで、幾何学的表面積に対する有効触媒面積の比率(触媒の担持量や粗さに相当)を表します。
  • 反応モデル: kineticディレクトリ内のファイルで定義された OpenSMOKE++ 形式の反応機構です。
    • Rh 詳細反応機構: 水素、酸素、CO、メタンなどの酸化・改質反応を含む詳細メカニズム。
    • WGS_UBI: 水性ガスシフト反応 (Water Gas Shift) を対象としたモデル。

チュートリアルケース:sphere

1. 解析の概要

ロジウム(Rh)触媒粒子周りにおける水性ガスシフト反応(WGS)解析をシミュレートします。

  • 物理現象:
    • 高温(850 K)かつ常圧(推定)のガス流れ。
    • 不均一反応(Heterogeneous Reaction): 壁面(reactingWall)において、ガス種が触媒表面に吸着・反応・脱離する過程を解いています。
    • 対象反応: 水性ガスシフト反応(\text{CO} + \text{H}_2\text{O} \rightarrow \text{CO}_2 + \text{H}_2)。
  • 解析領域:
    • 球上触媒(sphere)周りの流れ場。
    • 境界条件として「反応する壁(reactingWall)」と「不活性な壁(inertWall)」が定義されています。

2. 数理モデルとソースコード構造

使用されているソルバーcatalyticPimpleFOAMは、OpenFOAMの流体解析に 詳細化学反応計算ライブラリ(OpenSMOKEなど) をカップリングしたものです。

  • 支配方程式:
    • ナビエ・ストークス方程式: 流れの運動量保存。
    • 化学種輸送方程式: N_2, CO, H_2O, CO_2, H_2 などの移流拡散。
    • 表面被覆率発展方程式: 触媒表面上の吸着種(Rh_s, CO_s, O_s など)の占有率(Site Fraction)の時間発展を解いています。
  • 反応モデル (constant/solverOptions):
    • 気相反応: off(気相中では反応しない)。
    • 表面反応: on(壁面でのみ反応)。
    • 拡散モデル: multi-component(多成分拡散を考慮)。
    • ODEソルバ: OpenSMOKE(硬い常微分方程式系を解くための外部ソルバ)。
  • 反応メカニズム:
    • WGS_UBIフォルダ内のKinetics.kin, Surface.sur, Thermo.tdc,Transport.traで定義されています。これらはCHEMKIN形式に準拠した化学反応定義ファイルです。

3. 計算設定の詳細

過渡計算(PIMPLE法)として設定されています。

  • 時間設定 (system/controlDict):
    • 開始時刻: 0.0秒 / 終了時刻: 0.1秒
    • 時間刻み: 1.0 \times 10^{-8} 秒(非常に小さい刻みで、速い反応スケールを解像しています)
  • 離散化スキーム (system/fvSchemes):
    • 時間項: Euler(1次精度陰解法)。
    • 勾配項: Gauss upwind(1次風上差分)。安定性重視の設定です。
    • 対流項(化学種): Gauss limitedLinear01 1(有界性を保つ高次スキーム)。質量分率が0~1を超えないように制限されています。
  • 解法アルゴリズム (system/fvSolution):
    • PIMPLE法: 各タイムステップ内で反復計算(Outer: 1回, Inner: 2回)を行い、連成を安定化させています。
    • 線形ソルバ:
      • 圧力: GAMG(幾何学的マルチグリッド)。
      • 速度・化学種: PBiCGStab(安定化双共役勾配法) + DILU 前処理。
  • 初期・境界条件 (0/ ディレクトリ):
    • 流入組成: N_2 (96%), CO (2%), H_2O (2%)。生成物である CO_2, H_2 はゼロ。
    • 温度: 850 K 一様。
    • 壁面: reactingWall に対して type catalyticWall が指定されており、ここで表面反応速度が計算され、境界での流束(Flux)として流体場にフィードバックされます。

以上の設定により、触媒表面での反応生成物が境界層を通じて拡散していく様子をシミュレーションしています。

4. 計算実行

ターミナルで以下のコマンドを実行してメッシュを作成します。ケースファイルの条件を変更して再度計算を行う場合は、foamCleanTutorialsを実行することでフォルダを初期状態に戻せます。

./makeMesh_of9x

system/decomposeParDictに並列計算用のメッシュ分割設定が記載されているので、以下のコマンドを実行することで並列計算が行われます。

decomposePar -force
mpirun -np 4 catalyticPimpleFOAM -parallel
reconstructPar

計算結果をParaviewにて表示します。

touch foam.foam
paraview

5. 計算結果

この解析は、高温(850 K)のガス流中に置かれた球状触媒粒子表面における化学反応の過渡変化を、0.0~0.1秒の時間範囲でシミュレーションしたものです。
4つの球状触媒表面で水性ガスシフト反応(\text{CO} + \text{H}_2\text{O} \rightarrow \text{CO}_2 + \text{H}_2)が正常に進行しており、
反応物が消費され、生成物が生成される様子が確認できました。

  • 反応物(CO, H_2O)
    右側入口よりCOH_2Oが流入して球状触媒表面に付着していることが分かります。

  • 生成物(CO_2, H_2)
    球状触媒上で反応が進行しており、後段の3つ目以降の球状触媒付近でCO_2H_2が多く存在していることが分かります。

  • 反応ガス組成の推移
    球状触媒付近の気相中のガス組成を確認すると、COH_2Oが減少傾向を示し、CO_2H_2が増加傾向を示した。

Discussion