大分工業高等専門学校紀要

第 42 号(平成 17 年 11 月)

数値シミュレーションによる廃熱ボイラの気流混合最適化
菊川裕規 1 ，中尾里志 2 ，吉田尚平 2 ，亀谷敏博 31 機械工学科 ， 2 機械工学 科平成16年度卒業生， 3 日鉱金属株式会社

The exhaust-heat recovery boiler of a copper refining plant has the role of which oxidizes the sulfide
contained gas by blowing air from the exhaust gas discharged in self-melting furnace. When observing
an internal phenomenon about the high temperature and large scale phenomenon, the real scale
experiment is difficult. Then, in order to solve the problem, the numerical simulation was used. In thisresearch, the numerical simulation of the flow blowing with the exhaust gas in a boiler is carried out.
The main purpose of this research is to search the optimal condition for the air entry rate making an
optimizing mixed state by visualizing those gases. As a result of calculating by a numerical simulation,
it turned out that the circulation style is deeply related to the air current mixture in a boiler. If the big
circulation style in a boiler is obtained, it will be thought that air current mixture of an exhaust gas andair is promoted.
Key Words: Computational Fluid Dynamics , Heat Recovery Boiler , Navier-Stokes equation

本研究ではボイラ内における排ガスと吹き込 み空

１．はじめに気の流れを数値シミュレーションし，それらの気流
銅 精 錬工 場の廃 熱 回収 ボイラ は ，次 の２つ の 目的

混合を可視化することで，理想的な混合状態を作りから使用される．

だすための流入速度の最適条件を求める事を目的と（１）壁面裏に設置 された 配管内 を流れ る水を輻射する．

熱で加熱し，水蒸気でタービンを動かす．

解析対象である廃熱回収ボイラの概要を Fig.1 に（２）廃熱回収ボイ ラ手前 の自溶 炉から 排出さ れる示す．この廃熱ボイラは長さが約 36m,高さが約 18m

硫化物を酸化させフィルタで回収して除去する．

という大きさである.

銅精錬過程の自溶炉内では

２．解析手法

CuFeS 2 + SiO2 + O2 →廃熱回収ボイラのように，実物が大きく内部 が高Cu 2 S ⋅ FeS + 2 FeO ⋅ SiO2 + SO2 + 反応熱

温となるような物について内部の現象を観察す る場の反応が起こってお り，左 辺第１ 項は銅 精鉱，右辺第１項はマット（銅品位 65%の銅），右辺第２項はス
ラグ（酸化鉄と珪酸 の化合 物）を 表す． 廃熱回収ボ
イラでは，自溶炉か ら排出 される この排 ガスから熱を回収すると共に， マット （硫化 物）を ，空気を吹き込むことで酸化させる役目がある．

合，実物での実験による現象の観察は難しい． 取り扱いを容易にするため，寸法が原型と相似な模 型を用いた実験がしばしば利用される．しかし，模 型実験では模型作製の手間がかかる他，相似則の緩 和によりいくつかの現象に誤差が生じる可能性があ る．
また，高温状態における熱の浮力の再現や混合 状態の可視化を模型実験で行うことは難しい．
そこで，模型実験の問題点を解決するため， コン
ピュータを使い実物と同寸法の３次元モデルを 作成
し，境界条件や保存則式を考慮させた上で数値 シミ
ュレーションを行う．これにより，模型作製の コス
ト削減を図り，モデル形状の変更や解析条件を 変化させた時のデータを得ることが容易となり，時 間短縮も可能である．
数値シミュレーションを行うために用いるの が，計 算 流 体 力 学

Fig.1 Schematics of Heat Recovery Boiler

CFD(Computational

Fluid

Dynamics)である．CFD では，流体の運動を表現するいくつかの重要な方程式として，ナビエ・ストー−29−ク ス の 方 程 式 (運 動 量 保 存 式 )や 連 続 の 式 (質 量 保 存

各面に作用する力を３方向成分について整理すると，

式)等の支配方程式を離散化し，コンピュータを使っ

x 方向についてはて高速に演算して流 れ場や 温度場 等を数 値的に解析する手法である．今回は，ADAPTIVE(RESEARCH

社の流体解析ソフト CFD2000 を用いた．

ｙ，ｚ方向についても同様に

３．支配方程式及び離散化法
３-１∂σ x ∂τ yx ∂τ zx++)dxdydz
∂x
∂y
∂z

(

ナビエ・ストークスの運動方程式

∂τ xy
∂x

+

∂σ y
∂y

+

∂τ zy
∂z

)dxdydz子の運動を詳しく調 べなけ ればな らず， 流体粒子の∂τ xz ∂τ yz ∂σ z++
)dxdydz∂y∂x∂z持つ物理量（３つの速度成分 u, v, w, 密度ρ，圧と表される．ここで単位質量に作用する体積力の x, y,

力 p 及び温度 ｔ ）を求めることにより流れの様子を

z,方向の成分 をそれぞれ Ｘ , Ｙ , Ｚ とし，ニュートン

調べる事ができる．

の第２法則(運動方程式)を流体粒子に適用すると

流 れ 場全 体の様 子 を詳 しく調 べ るた めには 流 体粒(

こ れ らの 物理量 は ，流 体の運 動 を支 配する ３つの法則（質量保存則， 運動量 保存則 ，エネ ルギー保存ρ∂τ yx ∂τ zx∂σDu)= ρX + ( x +
+
Dt
∂x
∂y
∂z

ρ

∂τ xy ∂σ y ∂τ zy
Dv= ρY + (
)
+
+
Dt
∂x
∂y
∂z

ρ

∂τ yz ∂σ z∂τDw
= ρZ + ( xz +
+
)
∂x
∂y
∂z
Dt

則）と状態方程式か ら求め られる ．ここ で取り扱う流体が音速より遅い場 合は 圧縮性の 影響 は省略され，流れは非圧縮性流れ として 取り扱 うこと ができる．
また温度一定の非圧 縮性流 体では 密度ρ は一定であり，粘性係数μも一 定とし てよい ため状 態方程式やエネルギー保存則と は無関 係に調 べる事 ができる．

と表される．せん断応力τについては

従って，非圧縮性・ 等エン トロピ ーの流 れは流体粒

dβ
∂u ∂v
= µγ xy = µ ( + )
dt
∂y ∂x
∂v ∂w
τ yz = τ zy = µ ( + )
∂z ∂y
∂w ∂u
τ zx = τ xz = µ ( + )∂x ∂zτ xy = τ yx = µ子の持つ３つの速度 成分と 圧力か ら決定 され，質量保存則・運動量保存則より求められる．
粘 性 流体 が流動 す る場 合，流 体 粒子 の伸縮 ・せん断変形などにより流 体中に 内部応 力が発 生する．ここで Fig.2 のように各辺の長さが dx,dy,dz,の微小六面体を考える．
こ の 微小 六面体 に 作用 する外 力 は流 体に均 等 に作の式が成立する．ここで dβ ，γはそれぞれ，せん断用する体積力であり ，内部 応力は 流体粒 子の単位面

歪・せん断歪速度を表す．また引張り応力σについ積あたりに作用する 表面力 と等し くなる ．また表面

ては

∂u 2
− µ∇ ⋅ v∂x 3∂v 2
σ y = − p + 2µ − µ∇ ⋅ v
∂y 3
∂w 2
σ z = − p + 2µ
− µ∇ ⋅ v∂z 3

力は各面に垂直な応 力と面 に平行 な２つ のせん断応

σ x = − p + 2µ

力に分けられる．微小六面体に働く表面力の x,y,z,方向成分は，

σ x , σ y , σ z , τ xy , τ yx , τ yz , τ zy , τ zx , τ xzの９つである．

の式が成立する．従って x, y, z 方向の運動方程式は∂u ∂u ∂u ∂u 1 ∂p ∂2u ∂2u ∂2u 1 ∂ ∂u ∂v ∂w+u +v +w = −+ν( + + ) + ν ( + + ) + f∂t ∂x ∂y ∂z ρ ∂x ∂x2 ∂y2 ∂z2 3 ∂x ∂x ∂y ∂z∂v ∂v ∂v ∂v 1 ∂p ∂2v ∂2v ∂2v 1 ∂ ∂u ∂v ∂w+u +v +w = −+ν( + + ) + ν ( + + ) + f∂t ∂x ∂y ∂z ρ ∂y ∂x2 ∂y2 ∂z2 3 ∂y ∂x ∂y ∂z∂w ∂w ∂w ∂w 1 ∂p ∂2w ∂2w ∂2w 1 ∂ ∂u ∂v ∂w+u +v +w = −+ν( + + ) + ν ( + + ) + f∂t ∂x ∂y ∂z
ρ ∂z ∂x2 ∂y2 ∂z2 3 ∂z ∂x ∂y ∂zとなる．これらを一つの式にまとめると

1
⎛ ∂u
⎞
⎛
⎞
+ u・∇u ⎟ = −∇p + µ ⎜ ∇ 2 u + ∇(∇・u )⎟ + f
3
⎝ ∂t
⎠
⎝
⎠

ρ⎜
Fig.2 Inside Power Distributionの形で表される．

−30−

ρ・u・p・μはそれぞ れ， ある位置 と時 刻におけるこの代数方程式に一般的かつ効率的な数値解法を適

流体の密度・速度・ 圧力・ 粘性率 である ．方程式は用して解を求める．

運動量保存則に基づ いてお り，左 辺は運 動量の時間今回解析を行ったボイラを直交座標系で簡略化し的・空間的変化を表 し，右 辺は圧 力およ び粘性によて作成したモデルを Fig.3 に示す．Fig.3 における

る運動量輸送を表し ている ．流体 を非圧 縮性流体とInlet は流入口,Outlet は流出口である.

見なすと， ∇ ⋅ u = 0

より

∂u
1
µ
+ (u ⋅ ∇)u = − ∇p + ∇ 2 u + fρρ∂t

４．解析及び結果
流入口 1 からは排ガスの主成分である二酸化硫黄を流入し，流入口 2 からは空気を流入した．廃熱回

これを微分形式で表すと，

∂u i
∂u
∂ ui
1 ∂p
+uj i = −
+ν
+ fi
∂t
∂x j
ρ ∂xi
∂x j

収ボイラには，排ガスに含まれる硫化物を酸化 させ2る目的があるので，ボイラ内でいかに効率良く ２つの気流を混合するかが問題である．そこで，２ つのとなる．これは粘性 流体の 運動に 対して 一般的に成流入口の流速を変化させボイラ内のガス濃度分 布を立する方程式で，こ の方程 式はナ ビエ・ ストークス

求める事で，最適な流入条件を推測した．

方程式である．
左 辺 第１ 項は速 度 の時 間変化 を 表す 非定常 項，第２項は各位置での流 速分布 を表す 移流項 である．この移流項が非線形性 を持っ ている ．右辺 第１項は圧力の変化を表す圧力 項，第 ２項は 運動及 び圧力エネ
ルギーが拡散される 変化率 を表し ，数値 計算において分布を平滑化させ 安定さ せる効 果を持 つ拡散項，
第３項は重力等の外 力項で ある． 移流と 粘性の強さの比率はレイノルズ 数と呼 ばれる 無次元 数であり，レイノルズ数がある 閾値を 越える と微小 なかく乱が移流項の非線形性に より拡 大して いくこ とで流れ場は非定常な乱流となる．
流 体 のシ ミュレ ー ショ ンは， こ の方 程式を 数値的に解く（離散化する）ことで流体の挙動を予測する．
３-２

Fig.3 Simplified Model for analysis

エネルギー方程式

固 定 され た微小体 積内 の流体 に 熱力 学第一 法則を適用すると，次式が得られる．
∂ ( ρH ) ∂ ( ρu i H )∂ κ ∂H
∂p
∂p
+
=
(
)+
+ ui
+ Φ + Q + S H ,P∂t∂xi∂xi Cp ∂xi∂t∂xi

ここに H は局所的な密度で平均化された静的な流体のエンタルピーを表し，κ は熱伝導率，Cp は質量で 平 均 化 し た定 圧比 熱 ， S H , p は エ ン タ ル ピ ー生 成・

吸収項を表す． Φ はストークスの分子散逸関数で

Φ = µ[(∂u i ∂u j∂u
2 ∂u k+)−δ ij ] i∂x j ∂xi
3 ∂x k∂x j

Fig.4 Velocity Vector for Result A (XZ-plane)(Inlet 1:4.4m/s, Inlet 2:7.5m/s)で表される．
今回解析に用いた流体解析ソフト CFD2000 では，
ナビエ・ストークス 方程式 ，及び 連続の 式，エ ネル
ギー方程式等の支配 方程式 を，有 限体積 法を用 いて離散化する事で解析を行う．
３-３

離散化手法

有限体積法とは，解析領域を多数のセルに分割し，個々のセルにおいて 支配方 程式を 積分し て全領 域において保存則を満た す離散 化式を 導く手 法であ る．

−31−

Fig.5 Velocity Vector (XZ-Plane)Fig.6 Distribution of Air Density for Result A(XZ-plane)

Fig.8 Velocity Vector for Result B (XZ-plane)(Inlet 1:4.4m/s, Inlet 2:15.0m/s)

Fig.7 Distribution of Air Density for Result A(XY-plane)
４-１

Result A

ま ず ， 作 成 し た モ デ ル の 精 度 を 検 証 す る た め ，銅精錬工場技術開発セン ター が行った 解析 と比較した．
各流入口の流入条件を同じとし，流入口 1 の流速を4.4m/s，流入口 2 の流速を 7.5m/s とした時の解析結

Fig.9 Distribution of Air Density for Result B(XZ-plane)

果を Result A とする．得られた流速ベクトルを Fig.4に示す．銅精錬工場 技術開 発セン ターの 解析結果を
Fig.5 に示す.
その結果，技術開発センターの解析結果と同様に，ボイラ上部に循環流 が現れ ，各断 面での 流速分布もFig.10 Distribution of Air Density for Result B(XY-plane)

同様の結果が得られ た．従 って， 今回作 成した簡略モデルの解析精度は 十分で あると 考えら れる．得られた空気濃度分布を Fig.6 及び Fig.7 に示す．
Fig.6 及び Fig.7 における数字は,ボイラ内の混合
ガス成分中の空気の割合を示す.Fig.6 からボイラ上部に 0.7（混合ガス成分中の空気の割合が 70%）という領域が広がっている 事が 見受けら れる .更に,全体的に 0.6 という領域が広い事から,この条件では気流混合が行われていないと言える.また,Fig.7 から 0.2〜0.6 という領域が複雑に存在している事からも,こFig.11 Velocity Vector for Result C (XZ-plane)(Inlet 1:4.4m/s, Inlet 2:22.5m/s)の条件では気流混合 が行わ れてい ないと 言う事ができる.
４-２

Result B

次に，流入口 1 の流速は 4.4m/s のまま変化させず，流入口 2 の流速を 2 倍の 15.0m/s に変化させた時の解析結果を Result B とする．得られた流速ベクトル
を Fig.8 に示す．Fig.8 より，循環流の位置が Fig.4
よりも下の方へ移動 してい る事が 見受け られる．空気濃度分布を Fig.9 及び Fig.10 に示す．Fig.9 かFig.12 Distribution of Air Density for Result C(XZ-plane)

ら,0.5 という領域がボイラ内の大部分を占めている事が見受けられる.また,Fig.10 から 0.2〜0.6 という領域が入り組んで存在している事が分かる.
以上の結果より,Result A よりも Result B の方が気流混合に関しては良い条件であると言える.

−32−

Fig.13 Distribution of Air Density for Result C(XY-plane)

４-３

Result C

次に,流入口 1 の流速は 4.4m/s のまま変化させず，流入口 2 の流速を 3 倍の 22.5m/s に変化させた時の解析結果を Result C とする．得られた流速ベクトル
を Fig.11 に示す．Fig.11 より，循環流の位置が Fig.8
よりも下の方へ移動 してい る事が 見受け られる ．空気濃度分布を Fig.12 及び Fig.13 に示す．Fig.12 から,0.5 という領域がボイラ内の大部分を占めている事が分かる.Fig.13 から 0.3〜0.7 という領域が存在

Fig.14 Velocity Vector for Result D (XZ-plane)(Inlet 1:4.4m/s, Inlet 2:30.0m/s)

しているが大部分は 0.6 という領域である事が分かる.また,Fig.12 及び Fig13 から,Fig.6 及び Fig.7 に比べて全体的に濃度が 均一 になって いる 事が分かる．
以上の結果から,Result C の方が Result A よりも気流混合に関しては良い条件であると言える.
４-４

Result D

次に,流入口 1 の流速は 4.4m/s のまま変化させず，流入口 2 の流速を 4 倍の 30.0m/s に変化させた時の解析結果を Result D とする．得られた流速ベクトル

Fig.15 Distribution of Air Density for Result D(XZ-plane)

を Fig.13 に示す．Fig.14 より，循環流の位置が Fig.11
よりも上の方へ移動 してい る事が 見受け られる ．空気濃度分布を Fig.15 及び Fig.16 に示す．Fig.15 から,0.5 という領域がボイラ内の大部分を占めている事が分かる.Fig.16 から,0.3〜0.6 という領域が存在

Fig.16 Distribution of Air Density for Result D(XY-plane)

しているが,0.5 という領域が大部分を占めている事が分かる.また,Fig.6 及び Fig.7 に比べて全体的に濃度が均一になっている事が分かる．以上の結果から,Result D の方が Result A よりも気流混合に関しては良い条件であると言える.
４-５

Result E

次に,流入口 1 の流速を 2 倍の 8.8m/s に変化させ，流入口 2 の流速を 7.5m/s のまま変化させなかった時の解析結果を Result E とする．得られた流速ベクト

Fig.17 Velocity Vector for Result E (XZ-plane)(Inlet 1:8.8m/s, Inlet 2:7.5m/s)

ルを Fig.17 に示す． Result A〜D では循環流がボイラ内に 1 つ存在していたが,Result E では Fig.17
よりボイラ上部に 1 つ,ボイラ内下部に 1 つ見受けられる．空気濃度分布を Fig.18 及び Fig.19 に示す．
Fig.18 から,ボイラ内の大部分を 0.4 という領域が占めている事が分かる.Fig.19 から,0.1〜0.3 という領域 が 存 在 し て い る 事 が 分 か る . ま た ,Fig.18 及 び
Fig.19 から,Fig.5 及び Fig.6 に比べて全体的に濃度が均一になっていない事が分かる．以上の結果より,Result E は Result A よりも気流混合に関しては

Fig.18 Distribution of Air Density for Result E(XZ-plane)

悪い条件であると言える.

Fig.19 Distribution of Air Density for Result E(XY-plane)

−33−

５.結論
今回の解析結果では,流入口 2 の流速を上げるほど，
循環流の位置が下の方へ移動する傾向が見られ た．
また，混合ガス濃度分布に関しては，Result A では
ボイラ内上部に空気濃度の高い混合ガスの領域 があり，全体的にも若干空気濃度が高いのに対して ，流入口 2 の流速を増加させた場合の Result C では，ボFig.20 Velocity Vector for Result F (XZ-plane)(Inlet 1:13.2m/s, Inlet 2:7.5m/s)

イラ内上部の空気濃度の高い混合ガス領域が無 くなり，全体的にも排ガスと空気の混合が均一の領 域が広がっている事から，ボイラ内の気流混合が促 進されたと考えられる．
逆に,流入口 1 の流速を上げると循環流が無くなり,ボイラ内の気流が大きく変わった.また,混合ガス濃度分布に関しても,気流混合が行われておらず,全体
的に排ガス濃度が高くなっている.
以上の結果から，循環流がボイラ内の気流混合に深く関係しており，ボイラ内に広い循環流が得 られ

Fig.21 Distribution of Air Density for Result F(XZ-plane)

れば，排ガスと空気の気流混合が促進されると 考えられる．また,今回の解析では，流入口 2 の流速を流
入口 1 の流速の３倍に設定した場合が最も望ましい
気流混合が得られた．

Fig.22 Distribution of Air Density for Result F(XY-plane)
４-６

謝辞
本研究は日鉱金属株式会社佐賀関精錬所との 共同研究として行った．ここに記して謝意を表する．

Result F

次に,流入口 1 の流速を 3 倍の 13.2m/s に変化させ，流入口 2 の流速を 7.5m/s のまま変化させなかった時

参考文献の解析結果を Result F とする．得られた流速ベクト(1)棚橋隆彦,はじめての CFD-移流拡散方程式-,

ルを Fig.20 に示す．Fig.20 より，Result A〜E で見
受 け ら れ た 循 環 流 が 無 く な り ,ボ イ ラ 内 の 流 れ が 大

コロナ社,(1996 年)(2)内田茂男,流体力学序説 基礎理論から応用まで,

きく変わっている事が分かる．空気濃度分布をFig.21 及び Fig.22 に示す．Fig.21 から,ボイラ内上

森北出版株式会社,(1990 年)(3)大槻義彦,div,grad,rot,…,共立出版株式会社,部に 0.7,0.8 という領域が存在しており,ボイラ内の大 部 分 は 0.3 と い う 領 域 が 占 め て い る 事 が 分 か(1993 年)
(4)ランダウ＝リフシッツ,流体力学 1,東京図書株式

る.Fig.22 から,0.1〜0.2 という領域がボイラ内の大会社,（1970 年）

部 分 を 占 め て い る 事 が 分 か る . ま た ,Fig.21 及 び
Fig.22 より,Fig.5 及び Fig.6 に比べて全体的に濃度が均一になっていない事が分かる．この結果か
ら,Result F は Result A よりも気流混合に関しては悪い条件であると言える.

−34−（2005.9.30 受付）