Thermal Science & Engineering Vol.15 No.4 (2007)

LES による微粉炭燃焼解析
山本 研二†

岡崎 輝幸†

室田 知也†

谷口 正行†

Large Eddy Simulation of Pulverized Coal Combustion
Kenji YAMAMOTO†, Teruyuki OKAZAKI†, Tomoya MUROTA† and Masayuki TANIGUCHI†

Abstract
A numerical simulation program for pulverized coal combustion based on LES (Large eddy simulation) has beendeveloped. It was validated by comparing with experimental data of a pulverized coal jet flame ignited by a hightemperature (1420 K) flue gas jet. The difference between predicted and measured coal burnout on the centerline of thejets was less than 5 percent, which corresponds to experimental error. The predicted ignition position was in good
agreement with the experimental data. Clustering of coal particles was observed significantly when the Stokes numberof the particles is nearly equal to unity.
Key Words: Coal combustion, LES, Turbulent jet, Pulverized coal

記

号

Cc
CD
CoalH L
Cs
dp
E
f
g
I
Ic
m
n
Nu
Pr
q
Ra
Re
s
†

: スマゴリンスキー定数

[-]

St
S
T
u

: 粒子径

[m]

U

: 粒子流速[m/s]

x, y , z
X

: 座標

[m]

: 粒子の座標[m]: 粒子の比熱: 抗力係数: 熱分解の潜熱[J/kgK][-][J/kg]: 吸収効率因子: 体積分率[m3/m3]: 重力[m/s2]: ふく射強度[W/m2]: クラスタ指数

[1/m]

: 質量

[kg]

[-]

: ガス流速

[m/s]

: 生成項
: 温度

: 面積と体積の比率yuz: 粒子の数: ヌッセルト数: ストークス数

[K]

[1/m]

: 混合層の幅[m]: 軸方向流速差[m/s]: 微少量（=0.1 mm）[m]: 一般化した保存量[-]Γ: 拡散係数: プラントル数[-]: 累積重量分率

: 伝熱量

[W]

: 粘性係数[Pa s]: 吸収係数[1/m]: ガス密度[kg/m3]: 粒子密度[kg/m3]: レイノルズ数[-]: ふく射における長さ[m]p: ステファンボルツマン定数[%][W/m2K4]

受付日: 2007 年 3 月 14 日, 第 44 回日本伝熱シンポジウムより受付，担当エディター: 花村 克悟
株式会社日立製作所 電力・電機開発研究所 (〒312-8507 茨城県ひたちなか市堀口 832 番地の２)- 233 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)t: 粒子の緩和時間[s]: 乱流の特性時間

[s]

t

xi

(

ui )

xi

(Γ

xi

) S

(1)

添え字
b

: 黒体の値c: 対流[4]．運動量と圧力の定義点が異なるスタッガード格i: i=1,2,3 = x, y, z子を用いた．運動量の保存式において対流項の離散p: 粒子の値化には 2 次精度中心差分法を用い，その他の保存式

r

: ふく射

の対流項には，ハイブリッド法[4]を用いた．時間積

rel

: 相対値

分には，完全陰解法を用いた．時間積分の精度低下

t

: 乱流

を防ぐために，時間刻みを小さくした（クーラン数:

vol: ボラタイルの値

約 0.1）．乱流シュミット数と乱流プラントル数を 1.0

char: チャーの値と仮定した．

この方程式を有限体積法で離散化して計算する

2.2
１

緒

言

流動モデル

LES 乱流モデルを用いて運動方程式を解いた．サブグリッドスケール応力の評価に Smagorinsky モデ

微粉炭燃焼ボイラの設計にあたっては，(1)燃料のル[5]を用いた．Smagorinsky 定数 Cs を 0.1 に設定し燃焼が十分に進むこと，(2)NOx や CO 等の環境汚染

た．圧力と流速の補正に，SIMPLE 法[4]を用いた．

物質の排出が少ないことが要求される．定格負荷から部分負荷にわたってこの要求を満足し，かつ経済

2.3

粒子モデル

性を考慮した運転を実現するためには，火炉内の流ラグランジ座標系で重力と抗力を考慮して粒子

動・燃焼・伝熱挙動を解析する技術が必要となる．

の速度と座標を計算した[6]．この方法は，DPM とこれまでに我々は，k- 乱流モデルを使用した数値呼ばれる．

解析技術を開発した[1]．さらに性能の高い新型ボイラの設計には，乱流モデルの解析精度を高める必要

dU i , jがある．しかし，LES（ラージ・エディ・シミュレ

dtーション）法などの高精度な乱流モデルを用いた微粉炭燃焼解析はほとんど報告されていない[2]．本研

dX i , j

究の目的は，LES 乱流モデルを用いた微粉炭燃焼解dt

析プログラムを開発し，石炭を高温のガスで加熱し

3C D
(u i
4d p pU i , j )U rel , jgi

U i, j(2)(3)て燃焼させる実験を対象に燃焼率と着火位置の計算
精度を検証することである．

ここで， U i , j と X i , j は，j 番目の粒子の i 方向速度，
粒子位置である．また， C D , U rel , j は，抗力係数と２粒子とガスとの相対速度であり，以下の式で計算し数値計算方法た．
2.1 数値計算方法
本研究の解析で考慮する気相の支配方程式は，質CD量，運動量，エンタルピ，質量分率(N2, O2, CO2, CO,24(1 0.15 Re p

0.687

H2O, H2,CH4)の保存式である．また，粒子の運動は) / Re p

(4)

1/ 2

U rel , j

ラグランジ座標系で DPM (Discrete Particle Method)により計算した．気相と固相間の質量，運動量，エ(ui U i , j )

2

(5)

i

ンタルピの移動を考慮した．
気相での Favre 平均した保存量 の一般化した保存式は，次の式で表される[3,4]．

また，式(2), (3)の時間積分に 2 次精度ルンゲクッ
タ法を用いて，高精度化を図った．

- 234 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

2.4

伝熱モデル

Ra char

本研究では，粒子とガスのふく射，対流伝熱を計

f char

算した．

char

E charr2

ふく射伝熱については，散乱を無視したふく射輸

char

(9)

4 3
r
3

送方程式（式(6), (7)）を，DT 法（Discrete Transfer
Method）[7]により計算した．

(8)

対流伝熱については，粒子が球形と仮定してヌッ

dI
ds

Ib

RaI

RaI b

セルト数 Nu を計算した[10]．

(6)

10 3

Tg4

Re p

Nu

(7)

1 Re p1

2 0.52Re1p/ 3 Pr 1 / 3

(10)

103

1 つの面から出すレイの数を 36 本に設定した．計算

Nu

時間の短縮のため，10 タイムステップに 1 回ふく射

2 0.57 Re1p/ 2 Pr 1 / 3(11)の計算をした．ふく射の計算をしない場合は，エン
タルピのソース項を前回と同じ値に設定した．
吸収係数( Ra )は，ガス分子(CO2，H2O)や固体粒子(チャーや飛灰等)の濃度により決まると仮定した．

ふく射，対流伝熱，熱分解による潜熱を考慮して固体粒子の温度を式(12)で計算した．

初めにガスの吸収係数( Ra g )の計算方法について説明する．燃焼ガスは，N2, O2, CO2, H2O 等のガス組成dm p C c TPを含んでいる．この中で 2 種類以上の原子を含むも

dtのは，輻射エネルギーを吸収または放射するが，単原子分子はほとんど輻射に影響しない．また，ある特定の波長だけで輻射エネルギーの吸収や放射が生

2.5

qc

qr

dmvol
CoalH L
dt

(12)

燃焼モデル

じることも良く知られている．原理に忠実にガスのガスの燃焼は Eddy Dissipation Model[11]により評

輻射エネルギーの伝達を求めるためには，波長毎に価した．また，粒子の熱分解反応は，2 反応競合モ

エネルギーの移動を計算しなければならず，膨大なデル[12]により評価した．反応速度は，レーザによ計算を必要とする．簡便な計算方法として，波長依る微粉炭粒子の着火実験[13]により測定した．チャ

存性を考慮してガスの吸収係数を決める

ーの酸化・ガス化反応は，別途実験で測定したデー

WSGG(Weighted Sum of Gray Gasses)[8]モデルがあタを用いて開発したモデル[14]により評価した．

り，このモデルを使用した．本研究ではふく射に影
３

響するガスとして CO2 と H2O のみを考慮した．

解析対象

次にチャーの吸収係数の計算方法について説明す
る．微粉炭燃焼においては，チャーの平均粒径は 40

微粉炭の燃焼実験で測定した燃焼率とガス温度m 程度である．ガス温度を 1500 ℃と仮定すれば，[15]を用いて計算精度を検証する．図 1 に実験装置

輻射能の最大となる波長は 1.63 m となり，粒径パの形状を示す．
また，
図 2 に解析の境界条件を示す．

ラメータ(= d/ )は 70 以上になる．この条件では，中心に設置した直径 7 mm の 1 次ノズルから室温

幾何光学が適用できる[9]．チャーの吸収係数を式(8)（300 K）の空気と石炭を噴出する（図 2 ①）．流速で計算する． f char はチャーの体積分率である．

は 10 m/s である．石炭と空気の質量比率である C/Achar は体積を投影面積に変換する係数で，式(9)では 0.5 である．

表される．最後に， Echar は，吸収効率因子であり，
文献(8)の方法で計算した．

1 次ノズルの外側に設置した正方形ノズルから1420 K の高温ガスを噴出する（図 2 ②）．出口の流- 235 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

Table 1 Coal analysis.

速は 4.8 m/s である．この高温ガスにより，石炭を加熱し着火させる．

Proximate analysis (wt%, dry basis)
Volatile matter(VM)
31.1
Fixed carbon(FC)
54.0
Ash
14.9
Fuel ratio (FC/VM)
1.74

本実験装置では，ダクト出口にガスを吸引するブ
ロアが設置されており，一定の流量のガスを吸引す
る．このため，ノズル出口からフードまでの開放空間で室内の空気を吸引する(図 2 ③)．数値計算では，
ここに流速の境界条件を設定する必要がある．本研

Ultimate analysis (wt%, dry ash-free basis)
C
83.1
H
4.6
O
9.9
N
1.9
S
0.5

究では，ブロア出口の流量から，ノズル入口の流量を引いたものが，開放空間で吸引した空気量と考えた．この部分の吸引の流速が同じと仮定し，開放空間における流入流速を 0.16 m/s と決めた．
実験に使用した石炭は，燃料比 1.74 の瀝青炭である．石炭の性状を表 1 に示す．

石炭の粒径分布を図 3 に示す．初めに，実験で測定した粒径分布に近くなる Rosin-Rammler の近似式

Cycloneを決めた．この累積質量分率( )の近似式は以下の式

Pulverized coal + primary air

Exhaust hoodで表される．

Preheated surrounding gas

500mm

Open
area

BlowerHoneycomb

100

Secondary propane
Aluminabeads

Primary
propane量割合を計算した(表 2)．150

Air+N2Honeycomb
catalyst

粒子の代表径は，体面積平均により計算した(165

100mm

100

Honeycomb

90

100mm

300mm
500mm

③

7mm

x

Inlet nozzle(a)Whole②

Cumulative Fraction (%)m)．

zm 以下の代表粒子径40，62.5，113 m を代表粒径とした．150 m 以上のFig. 1 Experimental apparatus.

Velocity=0.16m/s③
Tg=300K
ρ=1.17kg/m3
N2=0.768kg/kg
O2=0.232kg/kg(6)は，範囲の最小と最大の中間とした．すなわち，15，

Fire brick

Outletd p 1.3)50

30-50 m，50-75 m，75-150 m，150 m 以上の質

Nozzle

SiC honeycomb(

Rosin-Rammler の式を使って，粒径を 30 m 以下，

Heater

Cooling
water

1 exp

②

Velocity＝4.8m/s ①
Tg=1420K
Velocity＝10.0m/s
N2=0.760kg/kg
Tg=300K
O2=0.101kg/kg
N2=0.768kg/kg
CO2 =0.093kg/kg O2=0.232kg/kg
H2 O=0.045kg/kg C/A＝0.5

Measurement
Rosin-Rammler

80
70
60
50
40
30
20
10
0

Ⅰ
1

Ⅱ ⅢⅣ Ⅴ

10
100
Particle diameter (μm)(b)Inlet nozzle

Fig. 3 Particle diameter distribution.

Fig. 2 Boundary conditions.

- 236 -C 2007 The Heat Transfer Society of Japan○

1000

Thermal Science & Engineering Vol.15 No.4 (2007)

Table 2 Coal particle diameter.
Particle name
Particle diameter(μm)
Mass fraction (wt%)

I
15
37.5

II
40
27.6

III
62.5
16

Ⅳ
113
16.7

を図 6 に示す．点 A(入り口)の温度は 1420 K で，点B (z=100 mm)の温度は 1390 K であった．高温ガスの

Ⅴ
165
2.2

温度が低下するのは，ふく射による放熱の影響である．
z=150 から 200 mm までの間で，1 次空気噴流近傍

本研究で使用したメッシュを図 4 に示す．メッシに外側の高温ガスよりも温度の高い場所（C, D, E）

ュ数は約 100 万である．特に 1 次ノズル付近のガスが形成された．これは，熱分解により発生した可燃と粒子の混合が重要であるので，ここにメッシュをガスに着火して温度が上昇したことを示す．この位多く配置し，最小メッシュ間隔を 1 mm に設定した．

置を着火位置 IP(Ignition Position)と定義する．計算のまた，噴流の噴出方向を z 方向と定義する．

IP は約 140 mm であった．その下流では，さらにガス温度の高い場所がある（F, G, H）．温度の上昇する

300mm

過程を調べるために，計算結果のアニメーションをmm500作成した．これを分析すると，C, D, E などの火種が生じた後，この火種が下流に流されるにつれて温度が上昇し，最終的に連続的な火炎へと変化する．こ400

れらの計算で確認される現象は実験[15]と一致した．
実験の着火位置を高速度カメラ画像の輝度の高くなった場所と定義すると，着火位置は 130 mm±20％300であった．計算と実験の着火位置は，実験のばらつ
きの範囲内で一致した．
200

mm
500
100

400

Z=0

Preheatedsurrounding gas

Pulverized coal +primary air

300
Fig. 4 Numerical mesh.
４200

解析結果と考察
計算によるガス温度の瞬時値を図 5 に示す．石炭100と空気を供給する 1 次空気の温度は 300 K であり，
石炭はノズルの出口近くではすぐに着火しない．z=0から 150 mm までに，1 次空気の温度は，外側の高温ガスと混合することで上昇する．しかし，この場
所では混合するだけなので，中心軸上のガス温度は，

300

(K)

Z=0
1600

外側の高温ガスよりも低くなる．
次に表示のカラーレンジを変えた計算のガス温

Fig. 5 Calculated gas temperature.

度と，実験で測定した燃焼火炎の高速度カメラ画像- 237 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

1600 (K)

1300

mm
500
O2 Mass fraction

H

G

400
Hi-speed video
camera image
300

400

H

F
D

C

B

IP

A

300

G

F

E

mm
500

CO Mass fraction

E

200

200

C

D
100 IP

100

Z=0
Calculation

0.0

ExperimentFig. 6 Flame structure.

(kg/kg)Z=0(kg/kg)

0.232

0.0

0.02

Fig. 7 Mass fraction of oxygen and carbon monoxide.

次に，数値計算のガス濃度を調べた（図 7）．O2 濃さらに図 8 は，粒子径により粒子の分散に差があ度は，1 次空気のノズル近傍で高い．その後，周囲ることも示している．例えば，粒子径 165の高温ガスとの混合のために低下する．着火位置よ合，粒子は分散し，局所的に粒子が集合している場

りも下流では，高温ガスの酸素濃度（0.1 kg/kg）よ所はない．しかし，40，63，112.5 m の粒子の場合りも低い場所が生じる．熱分解により発生した可燃

は，A, B, C, D のような粒子の集合している場所がガスが，酸素と反応したためである．

ある．本研究では，粒子の集合をクラスタリングと

CO 濃度の高い場所は，1 次空気の噴流の外側に点在している．CO の高い場所は，熱分解反応の生じ

m の場呼ぶ．クラスタリングを定量的に評価するために，
以下のクラスタ指数（ I C ）を考えた．

ている場所である．z=0～150 mm の範囲で CO 濃度が高い位置は，半径 5～10 mm であった．ノズル出

IC

口の半径が 3.5 mm であるので 1 次空気噴流の外側に拡散した石炭粒子が加熱され，熱分解反応により

1 n nn2 k 1 j 1

1
1/ 2

3

X i, j

X i,k

(7)

2

i 1

CO が発生したと考えられる．また，1 次ノズル出口から z=150 mm までは，CO 濃度の最大値は変化しな
い．これは，CO の生成速度と消滅速度がつりあっていることを示す．しかし，着火の開始した点（図粒子が 1 点に集まっていると I C の値が大きくなる．

7 C, D, E）から CO 濃度が高くなった．ガスの着火は 0 割防止の微少量であり，本研究では 0.1 mm とによりガス温度が上昇し，熱分解とチャーの酸化反した．クラスタ評価式の計算結果を図 9 に示す．図

応が加速されたためと考えられる．

8 の見た目と同様，粒径 63 m の粒子で最もクラスで追跡した粒子の位置を示す．径の小さい粒子は加タリングが発生していることが分かった．
クラスタリングは，非定常的に粒子，ガスの動き

熱されやすいので熱分解率は高い．逆に径の大きなを計算する LES を用いたために計算できた．時間平

粒子は加熱に時間がかかるので，熱分解率が低い．

均化した RANS（Reynolds Averaged Navier-Stokes）

これらの傾向は，実験結果[15]と一致する．

方程式を解いた場合には粒子のクラスタリング（図図 8 は粒子の熱分解率の瞬時値を示す．点は計算- 238 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)8），局所的に発生する火種（図 6）を予測できない．

ストークス数が大きいときには粒子はガスの動きに

すなわち，着火位置の予測には，LES による非定常

追従しないが，ストークス数が小さい場合は，ガス場の計算が有効である．

の動きに追従する．また，ストークス数が 1 に近い場合，粒子濃度が高くなることが報告されている[16]．本研究では以下の式でストークス数を定義しInitial particle diameter
15μm

40μm

63μm

112.5μm

た．

165μm
mm
500

dp

p

(8)

18

D
C

2

400

St

(9)
t

300

B
A

y
uz

t

200

(10)

ここで， y は混合層の幅， u z は軸方向流速の速度差である．z=100 mm における値（ y =20 mm，100

u z =2.5 m/s）から，乱流の特性時間をτt=8 ms とした．また，粒子の密度を 1400 kg/m3，周囲のガスは，
Devolatilization ratio

0

100%

Z=0

500 ℃の空気を仮定した．ストークス数が 1 に近い
場合に，クラスタ指数が高く，従来の研究と同じ結論を得た．

Fig. 8 Devolatalization ratio.
15
Clustering index (1/m)

Clustering index (1/m)

15
14
13
12

14
13
12
11

11

10
10

0

50
100
150
Initial particle diameter（μm)

200Fig. 9 Dependence of clustering index on particle
diameter for.

0

1

2

4
Stokes number(-)

68Fig. 10 Dependence of clustering index on Stokes
number.
図 11 に，中心軸上の燃焼率の検証結果を示す．実

さらに，乱流中の粒子の挙動が，ストークス数 St

験の粒子サンプリングに使用したプローブは，内径（粒子の緩和時間( )と乱流の特性時間( t )の比）でが 10 mm であるので，計算でも直径 10 mm の領域整理できるか調べた(図 10)．乱流の特性時間は，速の平均値を使用した．計算結果は，実験の誤差（±

度勾配の逆数であり，ガスの混合時間に関係する．

5％）内で実験値と一致した．

- 239 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

参考文献

40

Coal burnout(%)

30

20
Calc.
Exp.

10

0
0

50

100

150 200 250 300
Axial position (mm)

350

400

Fig. 11 Validation of coal burnout.
５

結論

LES（ラージ・エディ・シミュレーション）法を用いた微粉炭燃焼解析プログラムを開発し，石炭を高温のガスで加熱して燃焼させる実験を対象に燃焼率と着火位置について計算精度を検証した．
1) ガス温度から判断した計算の着火位置と，高速度
カメラ画像の着火位置は実験のばらつきの範囲内（±20％）で一致した．
2)計算した中心軸上の燃焼率は，実験値と実験の誤差範囲（±5％）内で一致した．
3) ストークス数（=粒子の緩和時間/乱流の時定数）が 1 に近い場合に，クラスタ（粒子の濃度むら）を確認した．
これは，従来の研究結果と一致する．
また，クラスタを評価するクラスタ評価式を提案
し，クラスタ指数が大きいほどクラスタが発生することを確認した．
4) CO 濃度の最大値は，着火するまでは変化しない
が，着火後に熱分解反応が加速されたことで高くなった．CO 濃度が最大となる位置は中心軸から
半径 5～10 mm であり，1 次空気ノズルの半径 3.5mm よりも外側であった．外側に拡散した石炭粒子が加熱され，熱分解反応により CO が発生した．
5) 径の大きな粒子ほど，出口における熱分解率が低いことは，実験結果と一致した．

[1] 山本 研二他，”微粉炭燃焼ボイラの燃焼解析技術”, 日立評論, 82-2 (2000), 27-30.
[2] R. Kurose and H. Makino, “Large-eddy simulationof a solid-fuel jet flame” , Combustion and Flame,
135 (2003), 1-16.
[3] 日本機械学会編, 燃焼の数値解析 (2001) , 丸
善.
[4] S. V. Patanker 原著, コンピュータによる熱移動と流れの数値解析 (1985) , 森北出版.
[5] S.B. Pope, Turbulent Flows (2000), 587-604,
Cambridge university press.
[6] E.E. Michaelides, “On the drag coefficient and thecorrect integration of the equation of motion ofparticles in gases”, Journal of Fluids Engineering,
110 (1988), 339-341.
[7] F.C. Lockwood et al., “A New Radiation Solution
Method for Incorparation in General Combustion
Prediction Procedures”, Eighteenth Symposium(Int.) on Combustion (1981), 1405 – 1414.
[8] T.F. Smith et al., “Evaluation of coefficients for theweighted sum of gray gases model”, ASME Journalof Heat Transfer, 104 (1982), 602-608.
[9] C. Kim and N. Lior, “Easily computable good
approximations for spectral radiative properties ofthe particulate-gas components and mixture inpulverized coal combustors”, Fuel, 74 (1995),
1891-1902.
[10] 一色 尚次，北山 直方，伝熱工学 (1984), 85,
森北出版.
[11] B.F. Magnussen et al., “On Mathematical
Modelling of Turbulent Combustion with Special
Emphasis on Soot Formation and Combustion”,
Sixteenth Symposium (Int.) on Combustion (1976),
719-729.
[12] H. Kobayashi et al., “Coal Debolatilization at High
Temperatures”, Sixteenth Symposium (Int.) on
Combustion (1976), 411-425.
[13] 岡崎 洋文他, “レーザ加熱による石炭粒子の熱分解過程のモデル化”, 日本機械学会論文集
62-598, B (1996), 374-379.
[14] 山本 研二他, ” 微粉炭燃焼場におけるチャーの 表 面 反 応 モ デ ル ”, 日 本 機 械 学 会 論 文 集
65-637, B (1999), 260-266.
[15] M. Taniguchi et al., ”Pyrolysis and Ignition
Characteristics of Pulverized Coal Particle”,
Journal of Energy Resources Tec., 123 (2001),
32-38.
[16] 雷 康斌他，”固気混相乱流における固体粒子の密度分布”, ながれ, 23 (2004) , 385-393.

- 240 -C 2007 The Heat Transfer Society of Japan○