― 6 ―SENAC Vol. 56, No. （1 2023. 1）[ 共同研究成果 ]熱分解反 応シミュ レーショ ン におけ る温度の 変動
松川 嘉也：東北大学大学院工学研究科
松下 洋介：弘前大学大学院理工学研究科
青木 秀之：東北大学大学院工学研究科

本研究では、化学反応を考慮した熱・物質移動を伴う流体計算(CFD)を実施し、壁からの対流伝熱および反応熱に伴う温度の変動の大きさを検討した。また、温度の変動が化学反応に及ぼす影響を検討した。加熱壁近傍の温度勾配が大きい部分で、温度の変動の絶対値が大きくなったものの、その値は小さいことが明らかになった。平均温度に対する 1 %程度の時間変動は、熱分解反応速度にさほど影響を及ぼさないことが示され、これは反応熱が温度の変動を抑制する方向に働くためであると示唆された。

1.

はじめに

熱分解反応は、さまざまな化学工業プロセスにおいて、目的の反応あるいは不本意な副反応として起こっている。プロセスを最適化するために、CFD と化学反応を連成させた解析が多く行われている。化学工業プロセスの多くは乱流場であり、産業界では乱流を表現するモデルとして計算コストの兼ね合いから Reynolds-averaged Navier-Stokes (RANS)による解析が多く行われているが、RANS から得られる温度場は時空間的に平均化されたものである。化学反応速度は Arrheniusの式に代表されるように、温度の上昇に伴い指数関数的に増加するため、温度に変動がある場合に平均温度を絶対反応速度の算出に用いてしまうと実際の化学反応速度とかい離する可能性がある。本報では、熱分解シミュレーションを実施し、温度の変動の大きさを見積もった。

2.

シミュレーションに用いるソフトウエア構成の検討

2.1 シミュレーションに用いたライブラリ・コンパイラ
本研究では FORTRAN で記述された in-house コードにより流体計算を実施した。ライブラリと表 1 に示す行列計算ライブラリを用いた。Open MPI は AOBA-B にインストールされていして、表るものを用いた。AOCL、OpenBLAS および ScaLAPACK について、本研究を開始した 2020 年 11
月当初において、GCC 環境用のライブラリが提供されていなかったか、著者の力不足により使用できなかったため、GCC 環境においてはソースコードからコンパイルしたライブラリを用いた。
それ以外のライブラリはソースコードからコンパイルしたものを用いた。
表 1 シミュレーションに用いたライブラリ
ライブラリ名バージョン
Message Passing Interface
Open MPI
4.0
線形代数演算AOCL
2.2
OpenBLAS
0.3.10並列行列計算ScaLAPACK
2.1.0
疎行列直接解法
MUMPS
5.3.5
メッシュ分割ParMETIS
4.0.3
疎行列反復解法ソルバー
AMGS 並列版[1]1.10
AGMG [2-4]3.3.5
研究開始当初において、コンパイラとして AOCC 2.2、GCC 10 および Intel Compiler が利用可能

熱分解反応シミュレーションにおける温度の変動― 7 ―であった。表 1 のライブラリを本研究の in-house コードにて使用する場合、Intel Compiler では計算速度が GCC と比較して著しく遅いことを確認済みであった。そこで熱分解反応シミュレーショ
ンに先立ち、AOCC 2.2 と GCC 10 との速度比較を実施した。逐次計算で長時間 AOBA-B の計算リソースを使用するのを避けるため、比較的小さな 84 万格子からなる計算対象とした。

2.2 ベンチマーク解析手法・解析対象
コンパイラのベンチマークを、化学反応と熱・物質移動を考慮しない、Large Eddy Simulation(LES)によって比較的レイノルズ数が小さい非定常乱流場の急拡大流れを解析することで実施し
図 1 に示す。解析領域を約 84 万分割した。支配方程式にた。急拡大流れ[5]を模擬した解析対象を図
はフィルター操作を施した非圧縮性 Navier-Stokes 方程式を用いた。有限体積法に基づき離散化した。対流項および拡散項の離散化スキームには二次中心差分法を用いた。時間項の離散化スキームには 3 次の Adams-Bashforth 法を用いた。圧力と速度のカップリングには SMAC 法[6]を適用し、圧力補正値に関する行列方程式の解法には AMGS [1]を用いた。収束判定条件には流入質量流量で規格化した離散化した連続の式の誤差の絶対値の総和を用い、10-6 以下となった時点で収束解が得られたものとした。

図1

ベンチマークの解析対象

2.3 ベンチマーク結果
並列数ごとの AOCC 2.2 と GCC 10 の実行速度を図
図 2 に示す。並列数 20 以上で並列化効率が悪化しているのは、計算対象が 84 万格子と小さいためである。並列数 2 の時には AOCC 2.2 の方が高速であったが、並列数 3 以上では GCC 10 が上回り、並列数の増加とともに計算速度差が広がっていることがわかる。AOCC 2.2 を本ソフトウエア構成において利用した場合にはメモリアクセスが非効率であることが示唆される。並列数を 128 とした場合には、84 万格子と小さな計算対象ではコア間の通信が律速になるためか AOCC 2.2 と GCC 10 との速度差が小さくなった。in-house
コードが FORTRAN で記述されているため、C 言語や C++で記述されたコードにおいても同様となるかは不明であり、また、AOCC のバージョンアップによって状況が改善している可能性もある。線形代数演算ライブラリおよび疎行列反復解法ソルバーについても同様に速度比較を実施してそれぞれ OpenBLAS および AMGS 並列版を用いることとした。

― 8 ―

SENAC Vol. 56, No. （1 2023. 1）61/計算時間 [1000/s]

543

GCC 102AOCC 2.2

10020406080100120

並列数

図23.

AOCC 2.2 および GCC 10 の実行速度比較

熱分解シミュレーション

熱分解反応中における温度の変動の大きさを検討するために、in-house コードによりシミュレーションを実施した。

3 . 1 熱 分 解反 応 モ デ ル
本研究において化学反応モデルに要求されることは化学反応に伴う吸熱を表現することである。
そのため、反応の予測精度の議論は行わず、Appel et al.による 99 化学種からなる詳細化学反応機
構[7]に対し、反応熱の予測精度が低下を最小限にとどめるように簡略化を施した 20 化学種からなる簡略化反応機構を用いた。反応速度の時間積分には VODE [8]を用いた。
3.2 熱分解シミュレーションの解析手法
簡略化反応機構であっても計算負荷は大きく、LES 計算には多大な計算コストを要するため、研究の第一段階として、RANS による定常計算を行った。流体に関する支配方程式として、Reynolds平均および Favre 平均を施した Navier-Stokes 方程式、乱流エネルギーk の輸送方程式、比消散率 ωの輸送方程式および化学種の保存式に加えて、温度および温度の分散を考慮するための式を用いた。温度を考慮するための式をエンタルピーの輸送方程式あるいは温度を変数としたエネルギー"2 の輸送方程保存式とした。温度の分散を考慮するための式を、式式(1)に示すエンタルピーの分散ℎ�式あるいは式式(2)に示す温度の分散𝑇𝑇𝑇𝑇�"2 の輸送方程式とした。

2"2𝜇𝜇𝜇𝜇𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕�𝜇𝜇𝜇𝜇𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕�𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜇𝜇𝜇𝜇��
"2"2"2 𝜔𝜔𝜔𝜔�� +�
�+2�𝜌𝜌𝜌𝜌𝜌 𝜌 � +�𝜌𝜌𝜌𝜌𝜌 𝜌𝜌𝜌𝜌𝜌𝑗𝑗𝑗𝑗 ℎ � =� � − 𝐶𝐶𝐶𝐶ℎ 𝜌𝜌𝜌𝜌𝜌 𝜌�𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗 𝜎𝜎𝜎𝜎ℎ 𝜎𝜎𝜎𝜎ℎ𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗𝜎𝜎𝜎𝜎ℎ𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕

2"2������������������
𝜕𝜕𝜕𝜕
𝜇𝜇𝜇𝜇𝜇𝜇𝜇𝜇𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕�
𝜇𝜇𝜇𝜇𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕�
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
2 ""2 � +"2 � ="2 𝜔𝜔𝜔𝜔 +�𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌�
�𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝑗𝑗𝑗𝑗 𝑇𝑇𝑇𝑇�
�� +�
�+2
� � − 𝐶𝐶𝐶𝐶ℎ 𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌𝜌�
𝑇𝑇𝑇𝑇 � 𝑆𝑆𝑆𝑆𝚤𝚤𝚤𝚤 ℎ𝑓𝑓𝑓𝑓𝚤𝚤𝚤𝚤
�𝑝𝑝𝑝𝑝
𝜎𝜎𝜎𝜎ℎ𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗 𝜎𝜎𝜎𝜎ℎ 𝜎𝜎𝜎𝜎ℎ𝑡𝑡𝑡𝑡 𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑗𝑗𝑗𝑗𝐶𝐶𝐶𝐶(1)(2)𝚤𝚤𝚤𝚤

温度の分散の輸送方程式において最終項の寄与は小さく、かつ、温度の変動𝑇𝑇𝑇𝑇 " を見積もることが困難なため、これを無視した。
対流項の離散化スキームに min-mod の流束制限関数[9]に基づく TVD scheme を用い、拡散項の離散化に二次中心差分法を用いた。運動量保存式の時間積分を対流項と拡散項に分離し、対流項に 3 次の Adams-Bashforth 法を、拡散項に Crank-Nicolson 法を用いた。圧力速度結合解法に

熱分解反応シミュレーションにおける温度の変動― 9 ―SMAC 法 [6]を適用し、圧力補正式を AMGS [1]により解いた。全計算格子の温度の変化が 1 K/sを下回った時点で定常状態に達したと判断し、そこからさらに 0.1 s 経過するまで計算を実施した。
離散化した保存式の行列方程式の解法に多項式前処理付き安定化双共役勾配法 (Bi-ConjugateGradient Stabilized Method with polynomial preconditioning, Bi-CGSTAB) [10]を用いた。連続の式の残差を流入量で規格化し、それ以外の支配方程式の残差を左辺の絶対値の総和で規格化した。連続の式の規格化した残差が 10-6、それ以外の支配方程式の残差が 10-4 を下回った時点で収束解が得られたものとした。本研究の流体シミュレーションにおいては温度の変動が反応速度に及ぼす影響を無視できると仮定し、平均温度から反応速度を算出した。本来であれば、平均温度を Arrhenius式に基づく反応速度式に代入するのは不適切である。吸熱反応である熱分解においては、のちに議論するように温度の変動が化学反応速度に及ぼす影響を無視することは、温度変動の大きさをわずかに過小評価する。本検討の目的は、変動の大きさが取りうる最大値を把握することであり、この仮定は大きな問題にならないと考える。

3.3 熱分解シミュレーションの解析対象・解析条件
図 3 に示す解析対象を、約 9 万の計算格子数に分割した。あらかじめ入口流速を 40 m/s とし、Reynolds 平均および Favre 平均を施した Navier-Stokes 方程式、k の輸送方程式および ω の輸送方
程式のみを考慮した定常計算を実施し、十分に発達させた流れ場の結果から主流方向の速度、k および ω の分布を得て流入境界として与えた。入口温度を 600℃とし、反応管の内壁面の温度を900℃で一定とした。入口組成として窒素のモル分率を 0.8 およびエタンのモル分率を 0.2 と与えた。

inflow

outflow

Φ38

No slip1250

図3

熱分解シミュレーションの解析対象

3.4 熱分解シミュレーションの結果と考察
エンタルピーの分散の輸送方程式を解いた場合および温度の分散の輸送方程式を解いた場合を式(3)に比較した。エンタルピーの分散を解いた場合に、確率密度関数から分散を求める式である式より温度の分散を算出した。
"2 = �
𝑇𝑇𝑇𝑇�

𝑇𝑇𝑇𝑇𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚

𝑇𝑇𝑇𝑇𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚(𝑇𝑇𝑇𝑇 − 𝑇𝑇𝑇𝑇�)2 𝑃𝑃𝑃𝑃�(𝑇𝑇𝑇𝑇)𝑑𝑑𝑑𝑑𝑇𝑇𝑇𝑇(3)

図 4 に入口から主流方向に 10 m の位置における温度の分散の半径方向分布を示す。エンタルピ
ーの分散を解いた場合と温度の分散を解いた場合で、温度の分散の値の範囲は概ね等しかったものの分布の形状が異なった。これは、温度の分散を解く場合は最終項の温度の変動項の寄与を 0と考え、エンタルピーの分散を解く場合には化学種の分散がエンタルピーの分散に従うとしたためであると考える。分散の最大値が 140 K2 だったことから、変動の大きさは最大で 12 K、すなわ
ち最大温度の 1 %程度と見積もられる。

― 10 ―

SENAC Vol. 56, No. （1 2023. 1）

図4

入口から 10 m の位置における温度の分散の半径方向分布

熱分解反応中での温度の変動の大きさを推算できたので、温度の変動が熱分解反応に及ぼす影響を検討する。温度の変動には時間的な変動と空間的な変動があるが、本研究では時間的な変動を考える。詳細化学反応機構に基づく 0 次元の反応動力学計算により簡易な検討を実施した。温式(4)に基づいて与え、温度を変動させた。
度変動を模した熱流束 q を式

𝑞𝑞𝑞𝑞 = 𝐴𝐴𝐴𝐴𝑞𝑞𝑞𝑞 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐�𝜔𝜔𝜔𝜔𝑞𝑞𝑞𝑞 𝜕𝜕𝜕𝜕�

(4)

ここで、Aq は振幅、ωq は角速度である。反応による発熱・吸熱がなければ、この熱流束は、温度に対して Aq に比例した変動を与える。本研究では、反応による発熱・吸熱がない場合に初期温度に対して 変動の大きさが 1 %または 5 %となるような熱流束を与えた。
図 5 に示すように、熱流束の振動が温度に与える影響は、100 Hz の振動を与えた場合にはかなり限定される。一方、10 Hz と 50 Hz の場合では温度の変動が見られるようになった。これは、熱図 6 に、熱流束の振分解反応のような吸熱反応では熱流束の短周期振動が減衰することを示す。図動が原料のエタン C2H6 と主要な生成物であるエチレン C2H4 の質量分率に与える影響を示す。100
Hz の振動を与えた場合の消費量/生産量は、変動がない場合とほぼ同じであった。これは、これらのケースで温度変動が見られなかったことと一致する。また、10 Hz と 50 Hz で 1 %振動させた場合の質量分率は、温度のばらつきがあるにもかかわらず、ばらつきのない場合と差が小さかった。
このことは、温度のばらつきが小さい場合には反応にあまり影響を与えないことを示唆している。
10 Hz と 50 Hz の振動を 5 %加えたケースでは、熱流束の振動が反応に影響を及ぼしていることがわかる。長周期の温度変動は LES で解像することができるため、LES ではアレニウス式から反応速度を直接利用できる可能性が高い。一方、RANS モデルは原理的に時間変動を解像することができない。したがって、長周期の温度変動が存在する場合、RANS モデルで得た平均温度をArrhenius 式に基づく反応速度に適用すると誤差を生じることが懸念される。

熱分解反応シミュレーションにおける温度の変動

図5― 11 ―熱流束の振動が温度に及ぼす影響(a) エタン

図64.

(b) エチレン
熱流束の振動が(a)エタンおよび、(b)エチレン質量分率に及ぼす影響

結論

本研究の RANS に基づく熱分解シミュレーションにより熱分解反応中における温度の変動の大きさが小さいことが示唆された。また、温度の変動の大きさが熱分解反応に及ぼす影響を検討し、1 %未満の温度の変動が熱分解反応に及ぼす影響が極めて小さいことを明らかにした。さらなる詳細な検討として LES に基づく非定常解析が必要であり、次報にて報告する予定である。

謝辞
本研究は、東北大学サイバーサイエンスセンターのスーパーコンピュータを利用すること

― 12 ―

SENAC Vol. 56, No. （1 2023. 1）で実現することができた。また、研究にあたっては同センター関係各位に有益なご指導とご協力をいただいた。ここに謝意を示す。
参考文献[1] A. Fujii, T. Tanaka, H. P. C. Laboratory, AMGS library, 2010.
[2] Y. Notay, An aggregation-based algebraic multigrid method, Electronic Transactions on
Numerical Analysis 37 (2010) 123-146.
[3] A. Napov, Y. Notay, An Algebraic Multigrid Method with Guaranteed Convergence Rate, SIAM
Journal on Scientific Computing 34 (2012) A1079-A1109.
[4] Y. Notay, Aggregation-Based Algebraic Multigrid for Convection-Diffusion Equations, SIAM
Journal on Scientific Computing 34 (2012) A2288-A2316.
[5] N. Kasagi, A. Matsunaga, Three-dimensional particle-tracking velocimetry measurement ofturbulence statistics and energy budget in a backward-facing step flow, International Journal ofHeat and Fluid Flow 16 (1995) 477-485.
[6] A.A. Amsden, F.H. Harlow, A simplified MAC technique for incompressible fluid flow
calculations, Journal of Computational Physics 6 (1970) 322-325.
[7] J. Appel, H. Bockhorn, M. Frenklach, Kinetic modeling of soot formation with detailed
chemistry and physics: laminar premixed flames of C2 hydrocarbons, Combust. Flame 121 (2000)
122-136.
[8] P. Brown, G. Byrne, A. Hindmarsh, VODE: A Variable-Coefficient ODE Solver, SIAM Journal
on Scientific and Statistical Computing 10 (1989) 1038-1051.
[9] D.B. Spalding, A novel finite difference formulation for differential expressions involving both
first and second derivatives, International Journal for Numerical Methods in Engineering 4(1972) 551-559.
[10] H.A.v.d. Vorst, Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for theSolution of Nonsymmetric Linear Systems, SIAM Journal on Scientific and Statistical Computing
13 (1992) 631-644.