― 4 ―SENAC Vol. 53, No. （3 2020. 7）［共同研究成果］

不均一反応により変化する固体燃料 の大規模シミュレーション ―沼澤松川松下青木結：東北大学工学研究科
嘉也：東北大学工学研究科
洋介：東北大学工学研究科
秀之：東北大学工学研究科

本研究では，X 線 CT を用いてコークスの構造を直接反映させたおおよそ 2 億ボクセルのコークスモデルを作成し，反応温度を 1373 K および 1573 K としてそのコークスモデル内部においてガスの物質移動をともなうガス化反応速度解析を行った．さらに，コークスの不均質な構造を直接反映させた有意性を検討するため，比較対象として均質なモデルに対しても同様の解析を行った．
その結果として，X 線 CT を用いて不均質な構造を直接反映させた場合はコークスが不均一に脆化し，反応後期では均質なモデルを用いた場合には表現できない炭素基質ボクセルが集中する領域のみが残存した．さらに，全体の反応率および見かけの反応速度は均質なモデルを用いた場合
ではほとんど単調に変化したものの，不均質な構造を直接反映させた場合では非線形な変化を示した．

1.緒緒言
多孔質炭素材料であるコークスは高炉内において還元材やスペーサーなど様々な役割を担う．
特に，高炉内の通気・通液性は銑鉄の生産性に大きな影響を及ぼすため，スペーサーとしての役割は重要である．そのため，コークスは高強度であることが求められ，これまでコークスの強度に関して多くの研究が行われてきた[1–7]．これらの研究により，気孔構造がコークスの強度に大きく影響を及ぼすことが明らかとなった．また，高炉内においてコークスはガス化反応によって脆化し，粉化しやすくなる．粉化したコークスは，高炉内の通気性を阻害し，高炉操業に悪影響を及ぼす．そのため，コークスのガス化反応に関して多くの実験がなされてきた[8–10]．
近年では，コークスを含む多孔質炭素材料を対象としたガス化反応の数値シミュレーションも行われている[11–14]．Richter et al.は複数の粒子で構成された多孔質モデルを対象に，気相反応および固気反応を考慮した対流拡散問題を解き，その多孔質体モデル内部および周辺部におけるガスの濃度分布について検討した[12]．Xue et al.は，複数の円錐状の気孔を持つ単一の石炭チャーモデルに対して Richter et al.と同様の解析を行い，気孔であるモデル粒子の粒径が内部の反応速度の分布に及ぼす影響を評価した．さらに，X 線 CT を用いて作成した石炭チャーの多孔質構造を計算に直接反映した研究もある．Fong et al.は，X 線 CT 像から作成した石炭チャーモデルに対して物質移動を考慮したガス化反応の数値シミュレーション手法を提案し，石炭チャーの構造を均質とするモデルと比較した[14]．我々の研究では，コークスの X 線 CT 像から多孔質モデルを作成し，そのモデルの一部を対象として有限要素法に基づいた物質移動を伴う CO2 ガス化反応解析を行った[15]．しかしながら，この手法では計算負荷の観点からより大規模の計算は現実的ではないため，コークス全体の構造変化を詳細に捉えることは困難である．
そこで，本研究では離散化手法を有限要素法から有限体積法に変更し，さらに MPI による並列計算を可能とした計算コードを用いて，X 線 CT 像から作成したおおよそ 2 億ボクセルのコークスモデルを対象に反応温度 1373 K および 1573 K において物質移動を伴うガス化反応速度解析を行った．比較対象として均質なコークスモデルを対象とした解析も行い，実コークスの構造を考慮することの有意性について検討した．

不均一反応により変化する固体燃料の大規模シミュレーション

― 5 ―

2.手手法
2.1

解析対象

2.1.1 X 線 CT から作成したコークスモデル(Case N)
X 線 CT を用いることでコークスの構造を直接反映させた解析対象を作成した．具体的には電気加熱式コークス炉を用いて Goonyella 炭を乾留し，コークスを作製した．コークスを中心部から円柱状に切り出し，コークス試料とした．そのコークス試料の大きさは⌀20×14 mm である．マ
イクロ X 線 CT スキャナー(TOSCANER-32250 hd, Toshiba IT & Control Systems Corporation, Japan)
を用いて，その円柱状のコークスのマイクロ X 線 CT 像を撮像した．なお，画面分解能は 20.6µm/pixel であり，スライス厚は 32 µm である．本研究では，コークスが気孔および炭素基質から構成されると仮定したうえで，Otsu の手法[16]に基づいて 2 値化し、コークスの X 線 CT 像におけるボクセルを気孔と炭素基質に分類した．2 値化した画像を積層させることで，図 1 に示すコークスモデルを作成した．このコークスモデルの要素数は約 4 億ボクセルである．本研究では，高さ方向の中心の面を対称面とし，コークスモデルの半分の領域を解析対象とした．この解析対象の分割数は 950 × 950 × 224 (approx. 0.2 billion)である．また，炭素基質ボクセルの密度は 1400kg/m3 とした．

図 1 およそ 4 億ボクセルのコークスモデル．解析対象はコークスモデルの半分の領域(黒色の領域)である．This figure is reprinted with permission from Elsevier [21].

2.1.2 均質なコークスモデル(Case U)比較対象として均質な解析対象を作成した．この解析対象の大きさは Case N と同様に⌀20×14mm とした．また，炭素基質ボクセルの密度は Case N と同じである．解析対象内のすべてのボクセルに空隙率を与えることで，解析対象内部に均質な気孔構造を形成させた．なお，Case N の解析対象全体における平均の気孔率が 0.49 であるため，初期のタイムステップにおけるボクセルの空隙率をとした．また，Case U の各炭素基質ボクセルの表面積は，コークスモデル全体の表面積が Case N と等しくなるように設定した．

― 6 ―

2.2

SENAC Vol. 53, No. （3 2020. 7）

コークスのガス化反応の大規模シミュレーション

2.2.1 支配方程式
コークスの X 線 CT から作成した解析対象において，拡散による CO, CO2 および N2 ガスの物質移動を解析した．Case N では，X 線 CT 像により解像された気孔ボクセルおよび気相ボクセル内においてガスの拡散を考慮し，その拡散方程式は次式のように表せる:(1)
ここで，Dm,mix は多成分系における化学種 m の拡散係数を，Sm は化学種 m の生成項を，Ym は化学種 m の質量分率を，はガス密度を表す．生成項には，気孔と炭素基質との界面において CO2 ガス化反応を考慮した．また，ガスが炭素基質ボクセル内部に拡散しないと仮定した．
Case U では，解析対象内のすべてのボクセルに空隙率を与えたため，そのボクセル内においてガスの拡散方程式は次式のように表せる:(2)
なお，気相では，空隙率をとした．

2.2.2 ガスの拡散
解析対象のメッシュサイズは数十m オーダーであり，解析対象内では分子拡散が支配的であるとした．そこで，ガスの拡散は分子拡散とし仮定し，次式により化学種 m の拡散係数 Dm,mix を算出した：(3)
ここで，Dlm は二成分系の拡散係数であり，Xm は化学種 m の体積分率を表す．この式は二成分系の拡散係数を多成分系に拡張した式である．2 成分系分子拡散係数は次式に示す Chapman-Enskog
理論により算出した：(4)
ここで，M はガスの分子量，T はガス温度，P は全圧，σave は衝突する分子の衝突直径，ΩD-ave はLennard-Jones ポテンシャルによる衝突積分をそれぞれ表す．

2.2.3 不均一反応
炭素基質と気相および気孔の界面における不均一反応として，以下に示す総括反応を考慮した：(R1)
C + CO2 → 2CO
炭素基質ボクセルでは，CO2 ガス化反応によって炭素の質量を減少させた．一方で，気孔ボクセ
ルおよび気相ボクセルでは，ガス化反応によって CO2 を消費し，CO を生成させた．なお，反応の進行よらず炭素基質ボクセルの表面積は一定とした．また，CO2 ガス化反応の反応速度定数 kは，Kashiwaya et al.が提案した反応速度式[17]に基づいて次式より算出した：,(5)

不均一反応により変化する固体燃料の大規模シミュレーション― 7 ―ここで，kc および ka は結晶性炭素および非結晶性炭素の反応速度定数を，Nc, Na は全質量に対して結晶および非結晶性炭素が占める質量の割合を表す．炭素の反応速度定数 kc, ka は，Langmuir-Hinshelwood 機構に基づき，式(6,7)のように定義されている:
,(6)

,(7)

ここで，k1,c, k1,a, k2, k3 は素反応の反応速度定数，P は圧力である．なお，Kashiwaya et al.の研究を参考に Nc を 0.25，Na を 0.75 とした．また，炭素基質ボクセルの反応率 X は次式のように定義した：,(8)

ここで，m0 は初期の炭素基質ボクセルの質量，mt は経過時間 t の炭素基質ボクセルの質量を示す．
コークス全体の反応率は，次式で表せる：(9)2.2.4 解析条件
有限体積法に基づいて支配方程式を離散化し，拡散項には中心差分法を適用した．反応による構造変化を観察するため，時間積分法には完全陰解法を用いて解析時間 18000 s まで解析した．
な お ， 時 間 刻 み を 30 s と し た ． 行 列 ソ ル バ ー に AGMG (AGgregation-based algebraic
MultiGrid)[18–20]を用いることで，離散化方程式の相対残差が 10-6 になるまで反復計算を行った．
さらに，生成項が質量分率に依存するため，一タイムステップ中にそれらの項を少なくとも 4 回更新した．また，境界条件に関しては，図 2 に示すように高さ方向の中心の境界面において濃度勾配をゼロとし，この境界面を対称面とした．その他の境界は，N2: CO2: CO = 0.5: 0.5: 0 とし，ガスの質量分率を固定した．また，反応温度が反応によるコークスの構造変化に及ぼす影響を検討するため，反応温度を 1373 K あるいは 1573 K とした．この計算は in-house コードを使用しており，計算時間は東北大学サイバーサイエンスセンター所有の並列コンピュータ LX406Re-2 において 6 ノード 72 コアを用いておおよそ 2 日であった．

図 2. 境界条件．点線で示す対称面において質量分率の勾配をゼロとし，実線で示すその他の面では質量分率を固定した．This figure is reprinted with permission from Elsevier [21].

― 8 ―

SENAC Vol. 53, No. （3 2020. 7）

3.結果と考察3 . 1 構 造 変化 の 比 較
X 線 CT を用いてコークスの構造を直接反映することの有意性を示すため，2 つの Case で計算結果を比較した．図 3 に，炭素基質ボクセルの局所の反応率の空間分布の経時変化を示す．図 3 (a)および (b)に着目すると，Case N ではコークスが不均一に脆化したことがわかる．1373 K ではコークスの構造はほとんど変化しなかったものの，1573 K では反応の進行に伴いコークスの構造が大きく変化し，炭素基質ボクセルが集中する領域のみが残存した．特に，経過時間 18000 s ではイナーチニット(IMDC)のような特徴的な領域のみが残存した．この理由として，これらの領域では比表面積が小さく，見かけの反応性が低いことが考えられる．一方で，Case U では，炭素基質ボクセルはコークス内に一様に存在し，反応進行にともないコークスが縮小した．いずれの反応温度においても，局所の反応率は径方向の中心に近づくにつれ減少した．これらより，図 3 (b)に示す炭素基質ボクセルが集中した領域は X 線 CT 像に由来し，Case N ではこのような特徴的な構造を考慮可能であることが示された．
3.2 コークス全体の反応率および見かけの反応速度の比較
2 つの Case を定量的に比較するため，図 4 にコークス全体の反応率および見かけの反応速度の経時変化を示す．いずれの Case においてもみかけの反応率は，全体の傾向として時間の経過とともに単調に増加し，Case N の 1573 K では反応時間 200 min 以降において傾きが減少する非線形な挙動が見られた．これは図 4 (b)に示す見かけの反応速度からもわかり，反応後期において反応表面積が減少したためである．また，いずれの温度においても Case N の全体の反応率および見かけの反応速度は Case U よりも大きい値を示した．反応温度 1373 K では，Case N の見かけの反応速度が Case U よりもおおよそ 2 倍大きかった．一方，反応温度 1573 K では，Case N の見かけの反応速度は全体の反応率が 0 から 0.5 の範囲では Case U の 2 倍程度大きかった．この差は 2 つのCase で表面積と CO2 の拡散性が異なるためである．また，全体の反応率が 0.5 以上の範囲では両者の差が小さくなる傾向が見られた．これは両者の全体の表面積が同程度となったためである．
図 5 に反応によるコークス全体の表面積の変化を示す．ここで，コークス全体の表面積は炭素基質ボクセルと気孔ボクセルあるいは気相ボクセルの界面の表面積の総和である．1573 K の全体の反応率 0.5 以上の範囲に着目すると，Case N と Case U の全体の表面積の差が減少したことがわかる．したがって，全体の反応率 0.8 付近において見かけの反応速度が概ね良好に一致した．以上
より，Case N と Case U の反応率および反応速度の差は反応による表面積の変化挙動が大きく異なることに起因することが示された．

不均一反応により変化する固体燃料の大規模シミュレーション

― 9 ―

600 s

12000 s(a) Case N at 1373 K

18000 s600 s

12000 s(b) Case N at 1573 K

18000 s600 s

12000 s(c) Case U at 1373 K

18000 s600 s

12000 s18000 s(d) Case U at 1573 K
図 3 1373 K および 1573 K における Case U および Case N の反応率の空間分布の時間変化．
This figure is reprinted with permission from Elsevier [21].

― 10 ―

SENAC Vol. 53, No. （3 2020. 7）(a)(b)
図 4 ガス化反応によるコークスモデル全体の反応率および見かけの反応速度の変化．
This figure is reprinted with permission from Elsevier [21].

図 5 ガス化反応によるコークスモデル全体の表面積の変化．
This figure is reprinted with permission from Elsevier [21].

不均一反応により変化する固体燃料の大規模シミュレーション― 11 ―

4.結言
本研究では，実コークスの構造の複雑さを考慮するため，反応温度 1373 K および 1573 K において X 線 CT 像から作成したおおよそ 2 億ボクセルのコークスモデルを対象に物質移動を伴うガス化反応速度解析を実施した．また，コークスの複雑な構造を考慮することの有意性を示すため，比較対象として均質なモデルに対しても同様の解析を行った．X 線 CT 像を反映させた場合ではコークスが反応の進行にともない不均一に脆化し，炭素基質ボクセルが集中する領域のみが残存した．一方で，均質なモデルを用いた場合では炭素基質ボクセルはコークス内に一様に存在し，反応進行にともないコークスが縮小した．全体の反応率は，全体の傾向として時間の経過とともに単調に増加し，X 線 CT 像を反映させた場合の 1573 K ではその傾きが反応の進行にともない減少する非線形な挙動が見られた．これは反応後期において全体の表面積が減少したためである．
また，X 線 CT 像を反映させた場合の見かけの反応速度は，1373 K および 1573 K の反応初期において均質なモデルを用いた場合よりも大きい値を示した．これは両者の場合で表面積と CO2 の拡散性が異なるためである．1573 K の反応後期では，2 つ場合における見かけの反応速度の差が小さくなった．これは両者の全体の表面積が同程度となったためである．これらより，X 線 CT 像を反映させた場合は，均質なモデルを用いた場合では見られない不均一な脆化を表現可能であり，さらに反応中における全体の表面積の違いによって全体の反応率および見かけの反応速度が均質なモデルを用いた場合と大きく異なることが示された．

謝辞
本研究は，東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができた．また，研究にあたっては同センター関係各位に有益なご指導とご協力をいただいた．本研究は科研費 19J20961 の助成を受けたものです．

― 12 ―

SENAC Vol. 53, No. （3 2020. 7）

参考文献
[1]
[2]
[3]

[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]

[13]

[14]
[15]

[16]
[17]
[18]
[19]
[20]
[21]

Grant MGK, Chaklader ACD, Price JT. Factors affecting the strength of blast furnace coke. vol. 70.
1991. doi:10.1016/0016-2361(91)90150-9.
Patrick JW, Stacey AE. The strength of industrial cokes: Part 1. Variability of tensile strength in relation
to fissure formation. vol. 51. 1972. doi:10.1016/0016-2361(72)90044-0.
Patrick JW, Stacey AE. The strength of industrial cokes. 7. Further studies of the influence of additives
in a coke-oven charge on the tensile strength of coke. Fuel 1978;57:258–64.
doi:10.1016/0016-2361(78)90001-7.
Patrick JW, Stacey. The strength of industrial cokes Part 6 Further studies of the influence of cokebreeze in a coal charge on tensile strength of coke. Fuel 1975;54 Oct.:256–64.
Jenkins DR, Lomas H, Mahoney M. Uniaxial compression of metallurgical coke samples withprogressive loading. Fuel 2018;226:163–71. doi:10.1016/j.fuel.2018.03.173.
Kubota Y, Nomura S, Arima T, Kato K. Quantitative Evaluation of Relationship between Coke Strengthand Pore Structure. ISIJ Int 2011;51:328–36.
Saito Y, Kanai T, Yamazaki Y, Matsushita Y, Aoki H, Nomura S. Effect of Meso-scale Pore Structure
on Coke Strength. J Japan Inst Energy 2017;96:93–101.
Xing X, Rogers H, Zhang G, Hockings K, Zulli P, Ostrovski O. Coke degradation under simulated blastfurnace conditions. ISIJ Int 2016;56:786–93. doi:10.2355/isijinternational.ISIJINT-2015-704.
Xing X, Rogers H, Zulli P, Hockings K, Ostrovski O. Effect of coal properties on the strength of coke
under simulated blast furnace conditions. Fuel 2019;237:775–85. doi:10.1016/j.fuel.2018.10.069.
Pusz S, Krzesinska M, Smedowski, Majewska J, Pilawa B, Kwiecinska B. Changes in a coke structure
due to reaction with carbon dioxide. Int J Coal Geol 2010;81:287–92. doi:10.1016/j.coal.2009.07.013.
Richter A, Nikrityuk PA, Kestel M. Numerical investigation of a chemically reacting carbon particle
moving in a hot O2/CO2 atmosphere. Ind Eng Chem Res 2013;52:5815–24. doi:10.1021/ie302770j.
Richter A, Nikrityuk PA, Meyer B. Three-dimensional calculation of a chemically reacting porous
particle moving in a hot O2/CO2 atmosphere. Int J Heat Mass Transf 2015;83:244–58.
doi:10.1016/j.ijheatmasstransfer.2014.11.090.
Xue Z, Guo Q, Gong Y, Xu J, Yu G. Numerical study of a reacting single coal char particle withdifferent pore structures moving in a hot O2/CO2 atmosphere. Fuel 2017;206:381–9.
doi:10.1016/j.fuel.2017.06.035.
Fong GH, Jorgensen S, Singer SL. Pore-resolving simulation of char particle gasification using
micro-CT. Fuel 2018;224:752–63. doi:10.1016/j.fuel.2018.03.117.
Numazawa Y, Igawa D, Matsuo S, Saito Y, Matsushita Y, Aoki H, et al. Numerical analysis of reaction
degradation for threedimensional coke pore structure. ISIJ Int 2018;58:1420–6.
doi:10.2355/isijinternational.ISIJINT-2018-144.
Otsu N. A Threshold Selection Method from Gray-Level Histograms. IEEE Trans Syst Man Cybern
1979;9:62–6. doi:10.1109/TSMC.1979.4310076.
Kashiwaya Y, Ishii K. Kinetic analysis of coke gasification based on non-crystal/crystal ratio of carbon.
Tetsu-to-Hagané 1990;31:440–8. doi:10.2355/isijinternational.31.440.
Notay Y. An aggregation-based algebraic multigrid method. Electron Trans Numer Anal
2010;37:123–46.
Notay Y. Aggregation-Based Algebraic Multigrid for Convection-Diffusion Equations. SIAM J Sci
Comput 2012;34:A2288–A2316.
Notay Y, Napov A. An Algebraic Multigrid Method with Guaranteed Convergence Rate. SIAM J Sci
Comput 2012;34:A1079–A1109. doi:10.1137/090750688.
Numazawa Y, Saito Y, Matsushita Y, Aoki H. Large-scale simulation of gasification reaction with masstransfer for metallurgical coke: Model development. Fuel 2020;266:117080.
doi:10.1016/J.FUEL.2020.117080.