Prediction of Pulverized Coal Combustion Behavior around Tuyere by Using LES and Extended CPD Model

Synopsis:

Pulverized coal injection into blast furnace is the precious technology to decrease the iron making cost. To inject more pulverized coal, it is necessary to investigate the combustion behavior of pulverized coal around the tuyere and many numerical simulations as well as experiments have been carried out. In this paper, by using Extended CPD model which is new devolatilization model and Large Eddy Simulation model, the combustion behavior of individual pulverized coal particles was investigated. The volume fraction of CO+CO₂ which is the indicator of how much carbon included in pulverized coal changed to gas phase was reasonably predicted by new simulation model. It is suggested the devolatilization processes of pulverized coals under 40 µm varied widely between individual particles.

1. 緒言

高炉への微粉炭吹込みは，老朽化が進行したコークス炉の延命やコークス比削減による製銑コストの低廉化などの観点から製銑工程における基幹技術となっており，微粉炭の多量吹込みが国内外を問わず指向されている^1,2)。しかし，微粉炭比が100 kg/t-pigの操業の場合には炉内圧損への影響が小さい^3,4)微粉炭吹込み操業も，150 kg/t-pig，200 kg/t-pigと微粉炭の吹込み量が増加するに従い，羽口部におけるコークス粉が増加し，高炉内の炉下部の通気性が悪化するといった課題が顕在化している⁵⁾。コークスの粉化による炉内の通気性の悪化には微粉炭の燃焼性の改善が有効なため⁶⁾，高燃焼効率を達成できるランスの開発²⁾や異種の微粉炭を混ぜることによる燃焼性改善の検討⁷⁾が行われている。これらの検討には微粉炭の燃焼実験と同時に数値シミュレーションが微粉炭燃焼の挙動を把握するための有用なツールとして用いられてきた。

高炉への微粉炭吹込みにおける燃焼場の数値シミュレーションへの取り組みの歴史は古く1980年代には1次元(1D)シミュレーションによる検討が行われていた⁸⁾。1990年代には2Dに拡張され羽口，レースウェイ内の流れ，伝熱，微粉炭の燃焼反応の解明に向けた取り組みが行われた^9,10)。最近では3Dシミュレーションによる高炉微粉炭吹込み時の燃焼場のシミュレーションを行った事例も報告されている^11–13)。しかしながら，これらの事例では，流れ場のあらゆる渦スケールをモデル化し，いわゆるレイノルズ平均流れを評価するReynolds-averaged Navier-Stokes(RANS)シミュレーション法が用いられる。RANSシミュレーション法は，計算負荷は低いものの，火炎の揺らぎ等，本質的に非定常な現象である燃焼場の再現性に乏しく，燃焼流の構造や微粉炭粒子個々の詳細な燃焼特性を理解することは難しい。加えて，石炭の反応モデルのうち，燃焼場の特性を決める重要な因子の一つである熱分解に伴う揮発分放出のモデルに1ステップモデル¹⁴⁾や2ステップモデル¹⁵⁾など簡易なモデルが用いられることが多い。これらのモデルは，石炭の発熱量や元素組成に合うよう仮想化学種を用いること，また揮発分放出速度や放出量，揮発分組成といった極めて重要な物理量の粒子昇温速度依存性を無視すること等，高炉への微粉炭吹込みのような急速昇温条件への適用性に乏しいと考えられる。

他方，基礎研究を中心として，近年微粉炭燃焼場のシミュレーションに対して，より高精度に乱流場の再現が可能なLarge Eddy Simulation(LES)法の適用が検討されてきており，温度場や主要化学種の分布に加えて，NO濃度や粒子の分散挙動等，より詳細な燃焼特性の予測技術の開発が進められてきている^16–19)。また，揮発分放出モデルについては，石炭の¹³C-NMR分析値や工業分析値，元素分析値などの情報から構築される化学的ネットワーク構造モデルに基づき，微粉炭粒子の熱履歴に応じた揮発放出速度や組成などを推定するChemical percolation devolatilization(CPD)モデル²⁰⁾やFLACHCHAINモデル²¹⁾などが提案されて，主に火力発電用微粉炭ボイラーの数値シミュレーションへの適用が図られているものの，タール分等の多環芳香族を仮想化学種に置き換えるなど，数値シミュレーションとの親和性の点で不十分である²²⁾。極く最近，Umemotoら²³⁾はCPDモデルにキュリーポイントパイロライザによる急速熱分解実験データを加味することにより，初期熱分解における多環芳香族の放出過程を具体的に記述する拡張CPD(Ex-CPD)モデルを提案したが，高炉への微粉炭吹込みのような高速・高温条件への適用性を検証した報告はない。

ランスより吹き込まれた微粉炭(Pulverized coal: PC)は，羽口からレースウェイへと侵入する。レースウェイの奥行きは1~1.5 m程度²⁵⁾であり，羽口風速を150 m/sとすれば，レースウェイにおける微粉炭の滞留時間は約10 msと推定できる。このような微小なタイムスケールでの燃焼場を正確に把握するためには，ランス出口付近での微粉炭の着火挙動を含む非定常性や粒子の分散性を正確に再現する必要があり，LES法が適している。また，UmemotoらのEx-CPDモデルでは，従来CPDモデルでは不明であった軽質ガスの組成をCO，CO₂，CH₄，H₂O，タールの組成をベンゼン，ナフタレン，フェナトレンとすることで，初期熱分解によって放出される揮発分のガス組成を具体的に考慮できる。そのため微粉炭の初期燃焼挙動の精度が重要である高炉吹込み条件下ではEx-CPDモデルを用いることで従来のCPDモデルよりも高い精度で微粉炭燃焼場の再現が可能になると考えられる。

以上より本報ではLESとEx-CPDモデルの高炉への微粉炭吹込みへの適用性を検証するために，まず最も簡易的な微粉炭吹込み条件でLESとEx-CPDモデルを組み合わせた数値シミュレーションとコークス充填層型小型燃焼炉を用いた燃焼実験結果を比較した。また，シミュレーション結果より微粉炭粒子個々の揮発分放出速度について検討した。

2. 数値解析モデル

2・1 気相の支配方程式と解析手法

本研究では非構造格子系の有限体積法コードFrontFlowRed-Comb(NuFD/FFR extended by Kyushu University, CRIEPI, Kyoto University and NuFD)を用いた。また，LES(Large-eddy simulation)により乱流を考慮した。LESでは，渦スケールによる渦特性の違いを活かし，系に依存するより大きな渦をGS(Grid scale)成分として直接計算し，一方で普遍的な挙動を示す小さな渦(SGS(Subgrid scale)成分)²⁴⁾を乱流モデルで考慮することで非定常な流れ場を精度よく再現することが可能である。これは気相の支配方程式に格子程度のフィルタを掛けることで可能となる。乱流燃焼場のような密度変化を伴う現象を数値シミュレーションで取り扱う場合，Favre平均を用いる。Favre平均は以下のように定義される。

  

$$\tilde{\varphi }\equiv \frac{\overline{\rho \varphi }}{\overline{\rho }}$$(1)

  

$$\varphi =\tilde{\varphi }+\overline{\varphi }$$(2)

ここで， － は空間フィルタ操作を， ～ はFavre平均操作を表す。Favre平均を施した連続の式，ナビエ・ストークス式，エネルギ保存の式および，各化学種の質量分率保存の式から気相のLES方程式は構成される²⁶⁾。

  

$$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial }{\partial x_j}\left(\overline{\rho }{\tilde{u}}_j\right)=S_m$$(3)

  

$$\frac{\partial \left(\overline{\rho }{\tilde{u}}_j\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho {\tilde{u}}_i{\tilde{u}}_j+\overline{p}{\delta }_{ij}-{\tilde{\tau }}_{ij}-{\tau }_{u_i{u}_j}\right)-\overline{\rho }{g}_i={\tilde{S}}_{ui}$$(4)

  

$$\frac{\partial \left(\overline{\rho }\tilde{h}\right)}{\partial t}+\frac{\partial }{\partial x_j}\left[\overline{\rho }{\tilde{u}}_j\tilde{h}-\overline{\rho }a\frac{\partial \tilde{h}}{\partial x_j}+{\displaystyle \sum _{k=1}^n{\tilde{h}}_k\left(\overline{\rho }a-\overline{\rho }{D}_k\right)}\frac{\partial {\tilde{Y}}_k}{\partial x_j}-{\tau }_{u_{j}h}\right]=S_h$$(5)

  

$$\frac{\partial \left(\overline{\rho }{\tilde{Y}}_k\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\overline{\rho }{\tilde{u}}_j{\tilde{Y}}_k-\overline{\rho }{D}_k\frac{\partial {\tilde{Y}}_k}{\partial x_j}-{\tau }_{u_j}{}_{Y_k}\right)=\overline{{\dot{\omega }}_k}+{\overline{S}}_{Y_k}$$(6)

式(4)~(6)の左辺第2項に現れるτ_uiuj，τ_uih，τ_uiYkは，SGSの影響を表しており，別途モデル化が必要となる項である。SGS項のモデル化には，Dynamic Smagorinskyモデル^27,28)を用いた。圧力と速度の補正にはSMAC法²⁹⁾を用いた。式(4)の移流項の解析スキームは2次中心差分と1次風上差分を9:1でブレンドした。時間積分アルゴリズムには1次Euler陰解法を用いた。時間刻みはdt=1.6 μsとし，計算開始後5000ステップから7500ステップの間での平均値をとり，統計量を算出した。

2・2 粒子の支配方程式と解析手法

微粉炭粒子の挙動は，Lagrange法を用いて求めた。気相と粒子の相互作用はPSI-Cell法³⁰⁾を使用した。PSI-Cell法では，気相の質量，運動量，物質およびエンタルピの交換を粒子が存在する質点で行う。これにより，気相の質量，運動量，物質およびエンタルピの保存に関する支配方程式では，セル内に存在する粒子の影響の総和をソース項として与えることが可能となる。粒子の位置(x_p,i)，速度(u_p,i)，温度(T_p)，質量(m_p)の保存式は次式から構成される。

  

$$\frac{d{x}_{p,i}}{dt}=u_{p,i}$$(7)

  

$$m_p\frac{d{u}_{p,i}}{dt}=\frac{1}{2}{C}_d{A}_p{\mathrm{Re}}_p\left|u_i-u_{p,i}\right|\left(u_i-u_{p,i}\right)$$(8)

  

$$C_d=\frac{24\left(1+0.15{\text{Re}}_p{}^{0.687}\right)}{{\text{Re}}_p}$$(9)

  

$${\text{Re}}_p=\frac{D_p\left|u_i-u_{p,i}\right|}{\nu }$$(10)

  

$$\frac{d{T}_p}{dt}=-\frac{1}{c_{p,p}{m}_p}\left[A_p{q}_p{}^{″}-A_p{\epsilon }_p\sigma \left({\Theta }_R^4-T_p^4\right)-Q_p-L_v\frac{d{m}_p}{dt}\right]$$(11)

  

$$\frac{d{m}_p}{dt}=\frac{dV}{dt}+\frac{dC}{dt}$$(12)

ここで式中のΘ⁴_RはΘ⁴_R=(I/4σ)^1/4とした。また，輻射強度Iはdiscrete ordinate法³¹⁾によって求めた。式中のdV/dt，dC/dtは微粉炭中の揮発分の放出速度と固定炭素の反応速度を表す。乱流による粒子運動への影響はDynamic random walk SGSモデル³²⁾により考慮した。Dynamic random walk SGSモデルは式(8)中の気相の速度u_i=${\overline{u}}_i$+u_i'についてSGS成分であるu_i'を表現するためのモデルである。Dynamic random walk SGSモデルではu_i'を等方性仮定，標準偏差の$\sqrt{2/3k_s}$の正規分布の乱数をラグランジュ的にスケーリングすることによってu_i'を与えている。Dynamic random walk SGSモデルはSGS成分を加味したモデルのため，粒子分散の高精度化が期待できる。筆者らは過去文献^16–19)等において，実験との比較を通じて粒子運動モデルの精度検証を行うとともに，粒子分散挙動に関する現象解明に取り組んできた。さらに，Zhangら³³⁾は粒径が粒子分散挙動に及ぼす影響について調査し，本研究と同様の粒子運動モデル用いたシミュレーションは，粒径毎の分散挙動に関しても実験結果と良好に一致することを示している。

2・3 微粉炭燃焼モデル

本節では，微粉炭燃焼の解析モデルについて説明する。

2・3・1 揮発分放出モデル

ランスからブローパイプへ投入された微粉炭粒子は周囲からの熱伝達および輻射により急速に加熱され，炭化水素，CO，H₂を主成分とする揮発分を放出する。本研究ではEx-CPDモデルおよび1ステップモデルを用いた検討を行った。

オリジナルのCPDモデルは，石炭の化学構造を考慮可能なネットワーク型初期熱分解モデルの一つである。石炭の化学構造は，石炭を構成する最小単位である芳香核クラスターとそれをつなぐ結合からなるベーテ構造と仮定し，NMRデータと元素分析値のみから初期石炭構造を規定する。CPDモデルでは，熱分解で切断され減少する結合の数を速度論から計算することで，チャー，軽質ガス，およびタールの収率を得ることができるが，結合については分解性結合と非分解性結合の2種，官能基について1種のみが考慮されるため，生成する軽質ガスやタールの組成は不明であった。Ex-CPDモデルでは，新たに分解性結合2種，非分解性結合1種，および官能基4種を考慮し，更にNMRデータから得られる芳香核クラスターの平均分子量と迅速熱分解試験の結果から，石炭中に存在する1~3環の割合を決めることで，熱分解過程で放出される軽質ガスをCO，CO₂，CH₄，およびH₂Oとし，タールをベンゼン，ナフタレン，およびフェナントレンと具体的に表現することが可能となっており，数値解析への適用性が高い。

一方で1ステップモデルでは揮発分の放出を下記のように単純化し扱う。

  

$$Coal\to VM+Char$$(13)

揮発分放出速度は以下の式で与えられる。

  

$$\frac{dV}{dt}=k_v\left(V^{*}-V\right)$$(14)

  

$$k_v=A_v\text{exp}\left(-\frac{E_v}{R{T}_p}\right)$$(15)

ここでV*は微粉炭中に含まれる揮発分の質量分率を，VはV*のうち現時点において気相中に揮発した揮発分の質量分率を表す。またA_vは頻度因子，E_vは活性化エネルギ，Rは気体定数，T_pは粒子温度である。微粉炭中の揮発分の質量分率V*は加熱条件で変化することが知られている。V*と工業分析で得られた揮発分量V_proはQファクター(Q_fac)を用いると以下の関係式が成り立つ。

  

$$V^{*}=Q_{fac}\times V_{pro}$$(16)

本研究では，揮発分の放出過程における微粉炭粒子の粒径は一定とした。

2・3・2 チャーのガス化モデル

揮発分の放出が終了すると微粉炭粒子は，炭素を主成分とするチャーとなる。チャー中の固定炭素(FC)のガス化反応によってCOやH₂が生じる。チャーのガス化反応は以下の1つの発熱反応(式(17))と2つの吸熱反応(式(18)，式(19))を仮定した。O₂によるガス化反応はCO₂やH₂Oによるガス化反応よりも反応速度は速いが，高炉への微粉炭吹込みでは微粉炭の完全燃焼に必要な酸素量よりも熱風中の酸素量は少ない。このような条件下ではCO₂によるガス化やH₂Oによるガス化も流れ場に大きな影響を与えると考えられる。

  

$$\text{C}+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\text{O}}_2\stackrel{{\text{k}}_1}{\to }\text{CO}$$(17)

  

$$\text{C}+{\text{CO}}_2\stackrel{{\text{k}}_2}{\to }2\text{CO}$$(18)

  

$$\text{C}+{\text{H}}_2\text{O}\stackrel{{\text{k}}_3}{\to }\text{CO}+{\text{H}}_2$$(19)

ガス化反応の速度式はBhatia and Perlmutterが提案している以下のRandom Pore Model³⁴⁾を用いた。

  

$$\frac{d\gamma }{dt}=A_{0i}{P}_{Ai}^n\text{exp}\left(-E/R{T}_p\right)\left(1-\gamma \right)\sqrt{1-{\Psi }_i\text{ln}\left(1-\gamma \right)}$$(20)

  

$$\gamma =1-M_{C_i}/M_{C_0}$$(21)

ここで，Ψ_iは細孔構造係数，P_Aiはガス化剤分圧，γはチャーの反応率，M_Ciはシミュレーション中の固定炭素量，M_C0は初期の固定炭素量である。本研究では組成や元素割合が近いDT炭²³⁾の反応速度パラメータを用いた。

本研究では，チャーのガス化過程における微粉炭粒子の粒径は一定とした。

2・4 気相反応モデル

気相中のおける揮発分と熱風との燃焼反応は，以下の総括反応に従うものとした。本研究で考慮した反応経路を以下に示す。

  

$$C{H}_4+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.O_2\stackrel{k_4}{\to }CO+2H_2$$(22)

  

$$H_2+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.O_2\stackrel{k_5}{\to }2H_{2}O$$(23)

  

$$CO+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.O_2\stackrel{k_6}{\to }C{O}_2$$(24)

  

$$CO+H_{2}O\stackrel{k_7}{\leftrightarrow }C{O}_2+H_2$$(25)

  

$$C_6{H}_6+3O_2\stackrel{k_8}{\to }6CO+3H_2$$(26)

  

$$C_{10}{H}_8+5O_2\stackrel{k_9}{\to }10CO+4H_2$$(27)

  

$$C_{14}{H}_{10}+7O_2\stackrel{k_{10}}{\to }14CO+5H_2$$(28)

  

$$C_{16}{H}_{10}+8O_2\stackrel{k_{11}}{\to }16CO+5H_2$$(29)

乱流燃焼モデルにはLESのSGS成分の影響を考慮できるScale similarity filtered reaction rate model(SSFRRM)³⁴⁾を採用した。SSFRRMモデルは，LESを用いた燃焼シミュレーションのために開発されたモデルである。LESによる燃焼シミュレーションで必要になる空間フィルタを掛けた反応速度は，温度や化学種濃度のGS成分からスケール相似則により直接算出される。SSFRRMにおける${\overline{\dot{\omega }}}_F$は以下のように表される。

  

$${\overline{\dot{\omega }}}_F=\overline{\dot{\omega }\left(\rho ,Y_k,T\right)}=\overline{\dot{\omega }\left(\overline{\rho },{\tilde{Y}}_k,\tilde{T}\right)+K_1\left(\overline{\overline{\dot{\omega }\left(\overline{\rho },{\tilde{Y}}_k,\tilde{T}\right)}}-\overline{\overline{\dot{\omega }\left(\overline{\overline{\rho }},{\tilde{\tilde{Y}}}_k,\tilde{\tilde{T}}\right)}}\right)}$$(30)

ここで，K₁はモデル定数であり推奨値であるK₁=1.0を採用した。また，$\dot{\omega }$(ρ，Y_k，T)は以下のArrhenius式モデルで与えられる。

  

$$\dot{\omega }\left(\rho ,Y_k,T\right)=A_k{T}_g{}^n\text{exp}\left(-\frac{E_k}{RT}\right)\cdot {\left[F\right]}^d{\left[O\right]}^e$$(31)

Table 1に式(31)における反応速度定数を示す。本研究では芳香族の速度定数にはベンゼンの定数を用いた。

Table 1. Kinetic parameters for gas reactions.

	Ak	n	Ek [J/mol]	Reference
k4	3.0 × 108	0	1.26 × 105	32)
k5	6.8 × 1015	–1	1.68 × 105	32)
k6	2.2 × 1012	0	1.67 × 105	33)
k7	Functional Chemistry			34)
k9	2.4 × 1011	0	1.26 × 105	33)
k10	2.4 × 1011	0	1.26 × 105	33)
k11	2.4 × 1011	0	1.26 × 105	33)
k12	2.4 × 1011	0	1.26 × 105	33)

式(25)は水性ガスシフト反応と呼ばれる反応である。微粉炭のガス化反応では式(18)，式(19)に示す通りH₂OやCO₂が重要な役割を果たす。そのため気相中に含まれるH₂OおよびCO₂の量を正確に把握することが重要である。しかし，水性ガスシフト反応の総括反応における汎用的な速度定数は未だ得られていない。そこで，本研究では水性ガスシフト反応の適切な速度定数を求めるためにWatanabeらが提案した，Functional chemistryモデルを採用した³⁸⁾。Functional chemistryモデルでは，想定される流れ場のガス温度・化学種濃度条件における素反応解析を複数実施し反応速度パラメータのテーブルを事前に作成する。微粉炭燃焼場のシミュレーション時には計算格子毎に温度，化学種濃度に応じた反応速度パラメータをテーブルから決定する。素反応解析にはCHEMKIN-PRO³⁹⁾を使用した。CHEMKIN-PROでは，何千もの反応の組み合わせ(反応機構)を解き，反応の進行に伴う組成変化や温度変化などをシミュレーションすることで化学反応間の依存関係を明らかにすることができる。詳細反応メカニズムには，53の化学種および325の素反応から成るGRI-Mech3.0⁴⁰⁾を用いた。反応モデルには完全混合反応器モデルを用いた。本手法を採用することにより低負荷かつ高精度に水性ガスシフト反応を再現することが可能となった。

3. 解析対象系と解析条件

3・1 解析対象系

本実験に用いた実験装置をFig.1に示す。本装置はコークス充填層部，羽口部，ブローパイプ部から構成されている。充填層は幅400 mm，奥行き1000 mm，高さ1400 mm，ブローパイプ部の内径はΦ90 mm，羽口径はΦ65 mmであり，ブローパイプには8つの観察窓が付いている。観察窓にはそれぞれ8 Nm³/hの冷却用N₂ガスを導入している。高温熱風はLPGバーナーで燃焼ガスを作成後，酸素を混合し，Table 2の組成となるように調整してブローパイプ部に送った。ランスからは，キャリアガスであるN₂と微粉炭が吹き込まれ，ブローパイプ中の熱風と混合する。吹き込まれた微粉炭は炉壁，火炎，レースウェイ内の赤熱コークスからの輻射および熱風からの熱伝達によって急速に加熱され，燃焼しながらブローパイプ中を移動し，羽口よりレースウェイ内に進入する。また，ランス先端から37 mmおよび300 mmの位置からプローブを挿入し，燃焼ガスを径方向(水平方向)にサンプリングした。

Fig. 1.

Experimental apparatus.

Table 2. Hot blast conditions.

Volume [Nm3/h]	Temperature [K]	Gas Compositions [vol%]
269.4	1723	O2	N2	CO2	H2O
		32.1	50.5	7.46	9.95

本研究の解析形状および境界条件をFig.2に示す。今後の記述における座標系はFig.2に従うものとする。ブローパイプ部，ランス，観察窓および羽口は実験装置と同形状を再現した。本報での微粉炭燃焼挙動の評価点はブローパイプ内のため，レースウェイは簡略化した。計算格子点数は約320万点であり，ランス出口近傍において格子密度が高くなるように配置した。レースウェイの赤熱コークスの影響を考慮するため下流側の壁には1923 Kの温度境界条件を設定した。

Fig. 2.

Boundary conditions.

3・2 供試炭性状

実験に使用した石炭の性状をTable 3に示す。本研究で対象とした石炭は揮発分量(VM:volatile matter)と固定炭素(FC:fixed carbon)の比(燃料比)は1.61である。

Table 3. Coal properties.

Volatile matter [wt%]	33.5
Fixed carbon [wt%]	51.8
Ash [wt%]	12.4
C [wt%, d.a.f]	82.6
H [wt%, d.a.f]	5.24
N [wt%, d.a.f]	2.19
S [wt%, d.a.f]	0.49
O [wt%, d.a.f]	9.44
LHV [kJ/kg]	28670

石炭の粒径分布をFig.3に示す。本シミュレーションでは実験で測定した累積質量分率(η)を以下のRosin-Rammler分布の式で近似した。解析上は5~200 μmの範囲で式(32)の分布に従い微粉炭粒子がタイムステップ毎に投入される。なお，微粉炭の密度は1000 kg/m³とした。

  

$$\eta =100\times \left[1-\text{exp}\left(-{\left(\frac{D_p}{93.1}\right)}^{1.08}\right)\right]$$(32)

Fig. 3.

Particle diameter distribution.

3・3 揮発分放出モデルパラメータ

本研究ではEx-CPDモデル(Case1)と1ステップモデル(Case2)を利用した数値シミュレーションを実施した。

Ex-CPDモデルの設定に必要な芳香核クラスターの平均分子量は¹³C NMRの結果から求めた。また，タール成分中の1環，2環および3環の割合は組成や元素割合が近いDT炭²³⁾のパラメータを用いた。1ステップモデルで利用した

揮発分放出速度のパラメータをTable 4に示す。Case2は下記の揮発分放出速度の推定式を用いた⁴¹⁾。Qファクターは1.2とした。

  

$$A_v={10}^{22.1-0.067\times C\%}$$(33)

  

$$E_v=23.7\times {10}^{-3}\text{exp}\left(\frac{C\%}{6.35}\right)+232907$$(34)

Table 4. Devolatilization model parameter.

	Model	Av [1/s]	Ev [J/mol]	Qfac	Ref
Case1	Ex-CPD	–	–	–
Case2	1step	1.675 × 1017	2.35 × 105	1.2	37)

4. 結果

ランス先端から300 mmの位置でのz方向のCO+CO₂の体積分率およびN₂の体積分率を実験結果と各Caseでのシミュレーション結果とをあわせてFig.4に示す。r/Rは無次元半径であり，測定位置をブローパイプ半径で除して規格化した。r/R=0がブローパイプ中心部を表し，r/R=1が壁面を表す。Case1，Case2ともに微粉炭中に含まれる炭素分のガス化率を示す指標であるCO+CO₂の体積分率と反応に関与しないN₂の体積分率が径方向で実験結果とよく一致した。

Fig. 4.

Volume fraction of along the axial direction.

各Caseと実験のランス先端から距離に対するブローパイプ中心軸上のCO+CO₂の体積分率をFig.5に示す。Fig.6にz=0のxy平面におけるCO+CO₂の体積分率を示す。Fig.5からCase1，Case2ともに測定位置では実験結果とよく一致している。しかし，実験の測定点以外ではCase1とCase2でCO+CO₂の体積分率が異なる値を示している。Fig.5，Fig.6からEx-CPDモデルを用いたCase1の方が初期の揮発分放出速度が速いことが分かる。

Fig. 5.

Volume fraction at 300 mm from lance tip of CO+CO₂.

Fig. 6.

Volume fraction of CO+CO₂.

Ex-CPDモデルと1ステップモデルとの揮発分放出速度パラメータを比較するために，Ex-CPDモデルのシミュレーション結果から1ステップモデルにおけるA_v，E_v，Qファクターを見積もった。その結果をTable 5に示す。その方法はAppendix1に示した。

Table 5. Estimated devolatilization model parameter.

	Model	Av [1/s]	Ev [J/mol]	Qfac	Ref
Case1	Ex-CPD	4.639 × 106	5.17 × 105	1.17
Case2	1step	1.675 × 1017	2.35 × 105	1.2	37)

得られたA_v，E_vを用いて粒子温度に対する式(15)中のK_vをFig.7に示す。1ステップモデルを用いたCase2はEx-CPDモデルを用いたCase1と比較してK_v，つまり揮発分放出速度の粒子温度に対する感度が高い。一般的に，1ステップモデルのような単純なモデルのパラメータを汎用的に利用することは難しい。そのためモデル開発者は，解析対象の燃焼場を再現可能なパラメータをその都度決定する必要がある。本研究では，1ステップモデルを利用したCase2もEx-CPDモデルを用いたCase1と同程度に実験結果と一致した。しかし，Fig.7に示したように粒子温度の変化に対して，Case2のK_vは15桁程度変化するため，適用範囲を慎重に見極めた上でモデルパラメータを使用しなければ，すぐに過大評価または過小評価する危険性がある。さらに，炭種の変更や還元雰囲気や酸化雰囲気，CO₂分圧やH₂O分圧などの雰囲気条件の変更時に1ステップモデルのパラメータを事前に決定することは困難である。一方で，Ex-CPDは化学構造に基づいているため，揮発分放出パラメータは化学的に明確に決定することができる。そのため，Ex-CPDモデルは様々な燃焼場に対して汎用的に利用できるものと考えられる。今後，様々な条件下での実験結果とEx-CPDモデルを用いたシミュレーション結果を比較し，Ex-CPDモデルの高炉微粉炭吹込み条件下への適用性を検証していく。

Fig. 7.

Devolatilization parameter (K_v) of each cases.

5. 考察

Ex-CPDモデルの大きな特徴として，揮発分放出速度およびQファクターが粒子の昇温速度に依存して決定される。そこで，Case1の個々の粒子の揮発分放出挙動に着目し検討を行った。

粒子温度を1000 Kとしたときの粒径毎のK_vと粒子の昇温速度の関係をFig.8に示した。ここで粒子の昇温速度は，粒子温度が600 K，800 K，1000 K，1200 Kに到達する時間から最小二乗法を用いて求めた。Fig.8からK_vと粒子の昇温速度には正の相関がある。また，径の小さな粒子のK_vは粒子毎のばらつきが大きく，これは粒子毎の昇温速度のばらつきが大きいためだと考えられる。以下，40 μm以下の径の粒子を小径粒子とする。

Fig. 8.

Effect of the coal particle heating ratio on devolatilization parameter (K_v).

Fig.9にランス出口でのy軸に垂直な方向な面の流速ベクトルを示す。ランスの中心部から吹き込まれるキャリアガスのN₂はランス出口面に対して垂直な流れになるのに対し，ランス壁面の流れはコアンダ効果により大きな運動量を持つ周囲の熱風の流れに引き込まれる。z=0におけるxy断面でのN₂の体積分率と微粉炭粒子の位置をFig.10に示す。ランス先端から100 mmの地点でのyz平面における微粉炭粒子の分布をFig.11に示す。ランス中央部から放出された小径粒子はキャリアガスのN₂に囲まれながら昇温するため，昇温が緩やかになる。一方でランス壁面付近から放出された小径粒子はキャリアガスのN₂と一緒に熱風に引き込まれることで即座に昇温される。Fig.11に示すように小径粒子はz方向にばらつきがあり，このばらつきが小径粒子の昇温速度のばらつきの原因であると考えられる。一方でFig.10に示すように大きな径の粒子は，粒子自身の持つ慣性力が大きいためランス出口に対して垂直な方向へ運動し，キャリアガスのN₂と分離され熱風に昇温される。Fig.11に示すように大きな径の粒子ほど粒子位置にばらつきがない。そのため昇温速度が粒子毎に大きく変化せず揮発分の放出速度も差異が小さい。

Fig. 9.

Gas velocity vectors around the lance tip.

Fig. 10.

N₂ volume fraction distribution and particle distribution with particle diameters.

Fig. 11.

Particle distribution at 100 mm from lance tip.

Ex-CPDモデルを用いたミュレーションにより本研究で対象にした微粉炭燃焼場では，粒径毎に，小径粒子では粒子毎に燃焼挙動が異なることが示唆された。現在，微粉炭吹込みランスを2重管にして外管から酸素吹込み微粉炭粒子の周囲の酸素濃度を高める吹込み方法⁴²⁾，または外管に酸素ではなく易燃性の気体還元剤を吹き込むことで微粉炭を急速加熱する方法¹³⁾が検討・実用化されている。これらの複雑な微粉炭燃焼場では，条件により粒子の熱履歴が大きく異なると考えられる。また，微粉炭粒子周囲の酸素濃度や周囲ガスの温度が粒子毎に大きく異なり粒子個々の熱履歴も異なることが予想される。このような様々な微粉炭燃焼場を熱履歴の影響を加味しない1ステップモデルで汎用的な揮発分放出パラメータを決定することは困難である。しかしながら，粒子の個々の熱履歴に対して可変的に揮発分放出速度を決定することができるEx-CPDモデルならば，熱履歴に依存しない微粉炭特有のパラメータを利用するだけで粒子個々の熱履歴が大きく異なる燃焼場であっても汎用的に利用できる可能性がある。今後，上記の様々な条件下で実験結果とEx-CPDモデルの比較を実施し，Ex-CPDモデルの高炉への微粉炭吹込みへの適用性を検証していく。

6. 結言

本研究では，LESとEx-CPDモデルの高炉への微粉炭吹込みへの適用性を検証するために，まず最も簡易的な微粉炭吹込み条件でLESとEx-CPDモデルを組み合わせた数値シミュレーションとコークス充填層型小型燃焼炉を用いた燃焼実験結果を比較した。その結果，以下の知見を得た。

(1)LESおよびEx-CPDモデルを用いた数値シミュレーションは，微粉炭の炭素分のガス化率を示すCO₂とCOの体積分率の和，および反応に関与しないN₂の体積分率ついて実験結果を再現できることを確認した。

(2)Ex-CPDモデルを用いたシミュレーションでは，微粉炭の揮発分放出速度が粒径に加え，粒径が40 μm以下の微粉炭粒子では粒子毎に異なることが示された。

(3)Ex-CPDモデルは様々な高炉への微粉炭吹込み条件に対して，パラメータの変更なく汎用的に利用できる可能性がある。

Symbols

A_p：粒子表面温度[K]

A_v：頻度因子[1/s]

A_k：頻度因子

A₀：頻度因子

a：熱拡散率[m²/s]

C：炭素分[−]

C_d：粒子の抵抗係数[−]

c_p,p：粒子の比熱[J/(kg·K)]

D：化学種の拡散係数[m²/s]

D_p：粒子径[m]

E：活性化エネルギ[J/mol]

E_k：活性化エネルギ[J/mol]

E_v：活性化エネルギ[J/mol]

F：燃料の濃度[mol/m³]

g：重力加速度[m/s²]

h：エンタルピ[J/kg]

I：輻射強度[J/(m²·s)]

K1：モデル定数[−]

k_s：SGS成分の運動エネルギ[m²/s²]

K_v：揮発分放出度速度係数[J/mol]

L_v：蒸発潜熱[J/kg]

M_Ci：シミュレーション上の固定炭素量[kg]

M_C0：初期の固定炭素量[kg]

m_p：粒子質量[kg]

O：酸化剤の濃度[mol/m³]

P_Ai：ガス化剤の分圧[Pa]

Q_p：粒子のソース項[J/s]

Q_fac：Qファクター[−]

q_p”：粒子表面における熱流量[J/(m²·s)]

R：理想気体のガス定数[J/(mol·K)]

Re：レイノルズ数[−]

Re_p：粒子レイノルズ数[−]

S：生成項

T_g：気相温度[K]

T_p：粒子温度[K]

t：時間[s]

u：流速[m/s]

u_p：粒子の速度[m/s]

V：粒子から放出された揮発分質量[kg]

V*：粒子中の全揮発分量[kg]

x：空間座標

x_p：粒子の位置座標

Y：化学種の質量分率[−]

Greeks

γ：チャーの反応率[−]

δ：クロネッカーのデルタ

ε_p：放射率[−]

η：累積質量分率[−]

Θ_R：輻射温度[K]

μ：粘性係数[Pa·s]

ρ：密度[kg/m³]

σ：Stefan-Boltzman定数[J/(m²·s·k⁴)]

τ：応力テンソル[Pa]

φ：物理量

Ψ_i：細孔構造係数[−]

$\dot{\omega }$ ：反応速度[mol/(s·m³)]

演算子

$\overline{●}$ 空間フィルタ

$\widetilde{●}$ Favreフィルタ

Appendix1

本解析では揮発分の放出とチャーの表面反応は同時に生じる。そこで，揮発分の放出量のみを検討するためにCase1においてチャーの表面反応式(17)，(18)，(19)の反応を除いたシミュレーションを実施し，80個の微粉炭粒子の時間変化をモニタリングした。Fig.A1に時間と粒子温度の逆数の関係の例を示す。今回サンプリングした粒子の時間と粒子温度の逆数の関係は下に凸の2次関数または単調減少の1次関数で近似した。まず2次関数で近似を試みた後，式(A1)においてα₂<0の場合は，単調減少の1次関数(式(A2))で近似した。結果として，40個の粒子を式(A1)で40個の粒子を式(A2)で近似した。なお，式(A1)，式(A2)の係数は最小二乗法により決定した。

  

$$\frac{1}{T_p}={\alpha }_0+{\alpha }_{1}t+{\alpha }_2{t}^2ただし{\alpha }_0<0,{\alpha }_1<0,{\alpha }_2>0$$(A1)

  

$$\frac{1}{T_p}={\alpha }_0+{\alpha }_{1}tただし{\alpha }_0<0,{\alpha }_1<0$$(A2)

Fig. A1.

Time history of reciprocal particle temperature.

近似曲線と解析の結果の80個の微粉炭粒子の相関係数の平均値は0.996である。

粒子温度の逆数と時間の関係を2次関数で近似した場合，揮発分放出速度式(14)および式(15)は式(A3)となり，解析解は式(A4)または式(A5)となる。

  

$$\frac{1}{\left(V^{*}-V\right)}dV=A_v\text{exp}\left\{\frac{E_v}{R}\left(\frac{{\alpha }_1{}^2}{4{\alpha }_2}-{\alpha }_0\right)\right\}\text{exp}\left\{-\frac{E_v{\alpha }_2}{R}{\left(t+\frac{{\alpha }_1}{2{\alpha }_2}\right)}^2\right\}dt$$(A3)

  

$$-\text{log}\left(V^{*}-V\right)=A_v\text{exp}\left\{\frac{E_v}{R}\left(\frac{{\alpha }_1{}^2}{4{\alpha }_2}-{\alpha }_0\right)\right\}\sqrt{\frac{\pi R}{4E_v{\alpha }_2}}erf\left\{\sqrt{\frac{E_v{\alpha }_2}{R}}\left(t+\frac{{\alpha }_1}{2{\alpha }_2}\right)\right\}+C_0$$(A4)

  

$$V^{*}-V=C_1\text{exp}\left[-A_v\text{exp}\left\{\frac{E_v}{R}\left(\frac{{\alpha }_1{}^2}{4{\alpha }_2}-{\alpha }_0\right)\right\}\sqrt{\frac{\pi R}{4E_v{\alpha }_2}}erf\left\{\sqrt{\frac{E_v{\alpha }_2}{R}}\left(t+\frac{{\alpha }_1}{2{\alpha }_2}\right)\right\}\right]$$(A5)

ただし，C₁=exp(−C₀)は積分定数である。また，erfは誤差関数である。粒子温度の逆数と時間の関係を1次関数で近似した場合，解析解は式(A6)または式(A7)となる。

  

$$\text{log}\left(V^{*}-V\right)=\frac{A_{v}R}{E_v{\alpha }_1}\text{exp}\left(-\frac{E_v\left({\alpha }_0+{\alpha }_{1}t\right)}{R}\right)+C_0$$(A6)

  

$$V^{*}-V=C_1\text{exp}\left[\frac{A_{v}R}{E_v{\alpha }_1}\text{exp}\left(-\frac{E_v\left({\alpha }_0+{\alpha }_{1}t\right)}{R}\right)\right]$$(A7)

ただし，C₁=exp(−C₀)は積分定数である。A_v，E_vは次のようにして求めた。まず積分定数C₁を決定する。揮発分の放出が開始される時間ステップをt₀，揮発分の放出が終了する時間ステップをt_endとし，t₀≦t≦t_endまでの間に時間ステップがn個あるとする。各時間ステップにおける式(A5)の両辺の総和をとればC₁は式(A8)となる。また式(A7)の総和をとれば式(A9)となる。

  

$$C_1=\frac{n{V}^{*}-{\displaystyle \sum _{k=1}^n{V}_k}}{\text{exp}\left[-A_v\text{exp}\left\{\frac{E_v}{R}\left(\frac{{\alpha }_1{}^2}{4{\alpha }_2}-{\alpha }_0\right)\right\}\sqrt{\frac{\pi R}{4E_v{\alpha }_2}}\right]{\displaystyle \sum _{k=1}^n\text{exp}\left[erf\left\{\sqrt{\frac{E_v{\alpha }_2}{R}}\left(t_k+\frac{{\alpha }_1}{2{\alpha }_2}\right)\right\}\right]}}$$(A8)

  

$$C_1=\frac{n{V}^{*}-{\displaystyle \sum _{k=1}^n{V}_k}}{\text{exp}\left(-\frac{A_{v}R}{E_v{\alpha }_1}\right){\displaystyle \sum _{k=1}^n\text{exp}\left\{\text{exp}\left(-\frac{E_v\left({\alpha }_0+{\alpha }_1{t}_k\right)}{R}\right)\right\}}}$$(A9)

次にA_vを求める。t_k≦t≦t_k+1の間ではA_vkが一定と仮定し，式(A4)の時間ステップ間の差分を取れば式(A10)を導ける。また，式(A6)の時間ステップ間の差分を取れば式(A11)を導ける。

  

$$A_{vk}=\frac{\text{log}\left\{\raisebox{1ex}{$\left(V*-V_k\right)$}\!\left/ \!\raisebox{-1ex}{$\left(V*-V_{k+1}\right)$}\right.\right\}}{A_v\text{exp}\left\{\frac{E_v}{R}\left(\frac{{\alpha }_1{}^2}{4{\alpha }_2}-{\alpha }_0\right)\right\}\sqrt{\frac{\pi R}{4E_v{\alpha }_2}}\left[erf\left\{\sqrt{\frac{E_v{\alpha }_2}{R}}\left(t_{k+1}+\frac{{\alpha }_1}{2{\alpha }_2}\right)\right\}-erf\left\{\sqrt{\frac{E_v{\alpha }_2}{R}}\left(t_k+\frac{{\alpha }_1}{2{\alpha }_2}\right)\right\}\right]}$$(A10)

  

$$A_{vk}=\frac{\text{log}\left\{\raisebox{1ex}{$\left(V*-V_k\right)$}\!\left/ \!\raisebox{-1ex}{$\left(V*-V_{k+1}\right)$}\right.\right\}}{\frac{A_{v}R}{E_v{\alpha }_1}\left\{\text{exp}\left(-\frac{E_v\left({\alpha }_0+{\alpha }_1{t}_{k+1}\right)}{R}\right)-\text{exp}\left(-\frac{E_v\left({\alpha }_0+{\alpha }_1{t}_k\right)}{R}\right)\right\}}$$(A11)

式(A10)，式(A11)で求められたn-1個のA_vkの平均をA_vとする。

  

$$A_v=\frac{{\displaystyle \sum _{k=1}^{n-1}{A}_{v_k}}}{n-1}$$(A12)

最後にE_vは時刻t₀における揮発分放出量V’₀=0となるように求めた。上記で求めた揮発分放出量とCase1での揮発分放出量の時間変化の例をFig.A2に示す。なおCase1での揮発分放出量の時間変化と上記の方法で求めた揮発分放出量の時間変化の相関係数は0.97である。Qファクターは式(16)の通りFig.A2に示す全揮発分放出量V*と工業分析値の揮発分量Q_proの比から求めている。

Fig. A2.

Time history of released volatile matter.

文献

• 1)  T. Ariyama, M. Sato, R. Murai, K. Miyagawa, K. Nozawa and T. Kamijo: Tetsu-to-Hagané, 81(1995), 1114.
• 2)  M. Sato, R. Murai, T. Ariyama, A. Maki, A. Shimomura and K. Mori: Tetsu-to-Hagané, 84(1998), 37.
• 3)  M. Deguchi, K. Shibata and M. Fukuoka: Tetsu-to-Hagané, 76(1990), 691.
• 4)  Y. Iwanaga: Tetsu-to-Hagané, 77(1991), 71.
• 5)  K. Yamaguchi, T. Uno, T. Yamamoto, H. Ueno, N. Kon-no and S. Matsuzaki: Tetsu-to-Hagané, 82(1996), 641.
• 6)  M. Sato, R. Murai, T. Ariyama, A. Maki, A. Shimomura and K. Mori: Tetsu-to-Hagané, 85(1999), 717.
• 7)  K. Kunitomo, T. Orimoto, T. Nishimura, M. Naito and J. Yagi: Tetsu-to-Hagané, 90(2004), 190.
• 8)  J. He, M. Kuwabara and I. Muchi: Tetsu-to-Hagané, 72(1986), 1847.
• 9)  H. Aoki, H. Nogami, H. Tsuge, T. Miura and T. Furukawa: ISIJ Int., 33(1993), 646.
• 10)  K. Takeda and F.C. Lockwood: ISIJ Int., 37(1997), 432.
• 11)  R. Murai, M. Sato and T. Ariyama: ISIJ Int., 44(2004), 2168.
• 12)  Y. Shen, B. Guo, A. Yu, D. Maldonado, P. Austin and P. Zulli: ISIJ Int., 48(2008), 777.
• 13)  A. Murao, Y. Kashihara, K. Takahashi, N. Oyama, H. Matsuno and M. Sato: Tetsu-to-Hagané, 101(2015), 653.
• 14)  S. Badzioch and P.G.W. Hawksley: Ind. Eng. Chem. Process Des. Dev., 9(1970) 521.
• 15)  H. Kobayashi, J.B. Howard and A.F. Sarofim: Proc. Combust. Inst., 16(1977), 411.
• 16)  H.Watanabe, K.Tanno, Y.Baba, R.Kurose and S.Komori: Turbulence, Heat and Mass Transfer 6, Begell House, Danbury, CT, (2009), 1027.
• 17)  R. Kurose, H. Watanabe and H. Makino: KONA Powder Part. J., 27(2009), 144.
• 18)  M. Muto, H. Watanabe, R. Kurose, S. Komori, S. Balusamy and S. Hochgreb: Fuel, 142(2015), 152.
• 19)  S. Ahn, K. Tanno and H. Watanabe: Appl. Therm. Eng., 124(2017), 1194.
• 20)  T. Fletcher, A. Kerstein, R. Pugmire and D. Grant: Energy Fuels, 4(1990), 54.
• 21)  S. Niksa: Combust. Flame, 100(1995), 384.
• 22)  O. Stein, G. Olenik, A. Konenburg, F. Cavallo Marincola, B.M. Franchetti, A.M. Kempf, M. Ghiani, M. Vascellari and C. Hasse: Flow Turbul. Combust., 90(2012), 859.
• 23)  S. Umemoto, S. Kajitani, K. Miura, H. Watanabe and M. Kawase: Fuel Proc. Technol., 159(2017), 256.
• 24)  木円重雄，柳瀬眞一郎:乱流力学，朝倉書店，東京，(1999), 93.
• 25)  Y. Matsui, Y. Yamaguchi, M. Sawayama, S. Kitano, N. Nagai and T. Imai: Tetsu-to-Hagané, 92(2006), 919.
• 26)  P.Moin, K.Squires, W.Cabot and S.Lee: Phys. Fluids A, 3(1991), No.11, 2746.
• 27)  M.Germano, U.Piomelli, P.Moin and W.H.Cabot: Phys. Fluids A, 3(1991), No.7, 1760.
• 28)  D.K. Lilly: Phys. Fluids A, 4(1992), 633.
• 29)  A.A. Amdsen and F.H. Harlow: J. Comput. Phys, 6(1970), 322.
• 30)  C.T. Crowe, M.T. Sharma and D.E. Stock: J. Fluids Eng., 99(1977), 325.
• 31)  W.A. Fiveland: J. Thermophys. Heat Transf., 2(1988), 309.
• 32)  K. Lei, N. Taniguchi and T. Kobayashi: JSME Int. J., 45B(2002), 164.
• 33)  W. Zhang, K. Tainaka, S. Ahn, H. Watanabe and T. Kitagawa: Adv. Powder Technol., 29(2018), 2322.
• 34)  S. Bhatia and D. Perlmutter: AIChE J., 26(1980), 379.
• 35)  P.E.Desjardin and S.H.Frankel: Phys. Fluids., 10(1998), No.9, 2298.
• 36)  W.P. Jones and R.P. Lindstedt: Combust. Flame, 73(1988), 233.
• 37)  C.K. Westbrook and F.L. Dryer: Combust. Sci. Technol., 27(1981), 31.
• 38)  H. Watanabe, K. Tanno, H. Umetsu and S. Umemoto: Fuel, 142(2015), 250.
• 39)  Chemkin-Pro, Release 15113, Reaction design, Canonsburg, PA, (2013).
• 40)  G.P.Smith, D.M.Golden, M.Frenklach, N.W.Moriarty, B.Eiteneer, M.Goldenberg, C.T.Bowman, R.K.Hanson, S.Song, W.C.Gardiner, Jr., V.V.Lissianski, Z.Qin: GRI-Mech 3.0, GRI-Mech 3.0, http://www.me.berkeley.edu/gri_mech/, (accessed 2016-10-18).
• 41)  T. Maffei, A. Frassoldati, A. Cuoci, E. Ranzi and T. Faravelli: Proc. Combust. Inst., 34(2013), 2401.
• 42)  Z. Zhou, Q. Xue, C.LI, G.Wang, X.She and J.Wang: Appl. Therm. Eng., 123(2017), 1096.