微粒子吸着フィルター内の内部構造のモデル化と燃焼解析による再生システムの検討
研究代表者

1.

名古屋大学大学院工学研究科助教授

山本和弘

はじめにディーゼル車はエネルギー効率が良く，大きな動力を得られることから，ヨーロッパを中心に今後も需要が多く見込まれる 1)．しかし，排気ガス中に含まれる有害成分が多く，特に粒子状物質(PM)は人体に悪影響を及ぼすため，対策が急務となっている 2)．近年，フィルター（Diesel Particulate Filter，DPF）を用いて，排ガス中の PM を除去することが提案されている．多孔体構造を持つセラミックのフィルターに排気ガスを通し，PM を内部に吸着させて除去する装置である．フィルターは一定期間使用すると内部に PM が目詰まりし，再生する過程が必要である．PM を効率的に捕集・除去するためには，実験に加えて数値解析によりデータを蓄積することが望ましい．しかし，フィルターの内部構造は非常に複雑であるため，N-S 方程式など各種の保存式を解く従来の手法では，内部の現象を模擬することは容易ではない．そこで我々は格子ボルツマン法（LBM）に着目し，多孔体内の流れを数値的に検討してきた．これまでに，流路内に障害物を配置して多孔体内部を簡略化するモデルを提案し，流れを計算した 3)．また，実際の多孔体の内部構造を CT で求め，すすの燃焼を解析している 4）．
LBM では，粒子の分布関数を用いて格子の状態を記述し，分布関数に対する発展方程式をもとに流れを解析する 5)．この方法では，ミクロな視点から現象をモデル化しているため境界条件の記述が容易であり，複雑な境界面を持つ多孔体内の流れに対して特に有効である．これまでに行われた粒子の堆積計算の例として，Filippova らの研究が挙げられる 6)．流れを LBM で,また粒子をラグランジュ的に解析し,流れにより運ばれた径の均一な粒子が固体壁に吸着する様子を模擬した．しかしこの方法では，流れと粒子の運動を別々に扱う必要があり，また粒子の微細構造まで考慮しなければならない．DPF の内部構造の最適化やシステム設計を行うためには，多孔体内部にすすが堆積する目詰まり現象を考慮し，どの程度の圧力損失があればフィルターを再生する必要があるか，など事前に現象を模擬し，フィルターの特性を予測することが重要となる．
そこで本研究では，格子ボルツマン法によりすす粒子を含む流れを計算した．ここでは，粒子濃度の分布関数をもとにすすの堆積を計算した 7)．これにより，すすが堆積することによる流れの変化を模擬することができる．

2.

解析方法

2.1 格子ボルツマン法（LBM）
ここでは，LBM の計算方法について簡単に説明する．LBM では格子ガス法（LGA）と同様，粒子の並進と衝突の過程から流れを解析する．粒子の分布関数は，１時間ステップで１格子だけ分布関数が移動する並進過程と，分布関数が衝突により平衡分布へ緩和する衝突過程により記述される．今回は３次元の流れを計算するため，d3q15 モデル 8)（図１）を用いた．
空間は立方格子で離散化され，粒子の速度ベクトルはそれぞれ以下のようになる．

[e 1e2e31= c  0 0e4− 100010e5e6e7e8e9e 10e 11e 12e 13e 140− 1000100− 1111− 1− 1− 111− 1− 1− 111− 11− 11− 11− 1− 1e 15 ]− 11100 0 

ここで，c =δx /δt であり，δx とδt は格子間隔と時間間隔である．多孔体内では流速が小さいため，密度変化を考慮しない非圧縮のモデルを用いた．この場合，圧力 p に対する分布関数をもとに流れを計算する．圧力 p の分布関数に対する発展方程式と平衡分布はそれぞれ以下のようになる．

1pα ( x + eα δ t , t + δ t ) − pα ( x , t ) = − [ pα ( x , t ) − pαeq ( x , t )]τ ( e ⋅ u) 9 ( eα ⋅ u)2 3 u 2  pαeq = wα  p + p0 3 α 2 +−2 c42 c 2  c( 2)(3)

ここで pαeq は平衡分布関数であり，ｗα＝1/9(α＝1～6)，1/72(α＝7～14)，および 2/9(α＝15) である．
圧力 p と速度 u は以下の式から求められる．圧力 p と速度ベクトル u は次の式から求められる．

p = ∑ pα(4)αu = ∑ eα pα / p0(5)αまた，p0 =ρ0 cs2 であり，ρ0 は非圧縮性を仮定した際に一定とした密度，cs（=√1/3）は音速である．動粘性係数νは緩和時間τを用いて以下のように表される．

ν=2τ − 1 δ 2x6 δt(6 )すす濃度は，これまでの燃焼解析 9)におけるスカラー量と同様に，すす濃度に対する分布関数をもとに決定した．その発展方程式を示す．

FC ,α ( x + cα δ t , t + δ t ) − FC ,α ( x , t ) = −1[ FC ,α ( x , t ) − FCeq,α ( x , t )]τs(e ⋅ u) 9 (eα ⋅ u)2 3 u 2 −FCeq,α = wαYC 1 + 3 α 2 +2 c42 c2 c(8)(7)

これらの分布関数を積算することによりすす濃度（質量分率）が求められる．

YC = ∑ FC ,α(9)α

またすすの濃度場を計算するのに必要な緩和時間τc は，以下の式によりすすの物質拡散係数（Dc）に対応する値を用いる．

2τC − 1 δ 2x6 δtDC =(10)ただし正確な値は不明なため，本計算では式(6)の動粘性係数と同じ値を用いた．また今回は粒子状物質としてすすのみを考慮した．実際の粒子の吸着現象では，吸着表面の特性，重力，さえぎり効果やブラウン拡散が複雑に作用するため，どの程度のすすが吸着するかを簡単には予測できないが 10)，本解析では障害物の壁面に衝突したすすは，その表面にすべて吸着するものとした．またすすの堆積の影響は，物体周辺に吸着したすす濃度を計算し，各格子点のすすの濃度が１となった時点で空隙から障害物に変化させた．これにより，すすが堆積して流路が狭まる様子を模擬できる．
次に燃焼場の計算方法について説明する．実際の燃焼場ではさまざまな反応が存在し，かつ，温度上昇による急激な密度変化もあるため，燃焼計算には高次精度の離散化や安定なアルゴリズムが必要なことが多い．
また噴霧燃焼などのように異なる相が存在することもあるため，新しい計算手法である LBM により燃焼場を模擬することは，現象をよりミクロな視点からモデル化したり，計算の適用範囲を広げる意味でも有意義であると思われる．本解析では，既に我々が提案した燃焼場による LBM の解析方法を用いる 9)．そこでは，従来の流れ場に対する分布関数に加えて，以下のような温度・濃度に対する分布関数を与えることにより，温度場と濃度場を決定した．

Fs ,α ( x + cα δ t , t + δ t ) − Fs ,α ( x , t ) = −1[ Fs ,α ( x , t ) − Fseq,α ( x , t )] + wα Qs ,τs(e ⋅ u) 9 (eα ⋅ u)2 3 u 2 −Fseq,α = wα s 1 + 3 α 2 +422 c2 ccs = T, Yi ( i = C, O2 , CO2 )これらの分布関数を積算することにより温度と化学種の濃度が求められる．

T = ∑ FT ,α(13)Yi = ∑ FYi,α(14)ααs = T , Yi(12)(11)それぞれの緩和時間は，熱拡散率と各化学種の拡散係数により決定される．反応に関連する項 Qs は，計算している系と対象とする実際の場との相似性から与えられる．例えば LBM における化学種 i の質量分率に対する反応項 QYi は，実際の反応速度ωi を用いて以下のように与えられる．

QYi = ( ωi ) LBM = ωi {ρ oU 0 / L} {ρ 0U 0 / L}LBM−1(15)

ここで U0 と L は，それぞれの座標系における代表速度と代表長さであり，LBM の添え字をつけた場合は，LBM の座標系におけるそれぞれの値を示す．また化学種の反応速度は，Lee らが提案した 1 段の反応式を用いた 11)．これにより酸素はすすと反応して消費され，温度は増加することになる．

147510315212111496813

図 1 Velocity directions in d3q15 model.

2.2 座標系と境界条件
図 2 に座標系を示す．透過性の多孔体内部を解析するため，系内に一辺が D の立方体（障害物とした格子数）を流路中央にランダムに配置した．これにより空隙率と濡れ面積を自由に変えることができるが，この簡略化したモデルでも多孔体の内部流れが模擬できることは文献 3 において確認している．図の左よりすすを含んだ空気が一様に流入し，障害物に衝突したすすが表面に吸着する．後で詳しく説明するが，すすの吸着が進み，格子点のすす濃度が１を超えた格子点は次の時間ステップから障害物に変化する．障害物を配置した長さを L とし，計算領域は流路長さ 5cm×幅 1cm×奥行き 1cm とした．格子数は，x 方向（流れ方向）が 201，y 方向が 41，z 方向が 41 であり， L は 2.5cm とした．また境界条件としては，物体の表面で流速がゼロとなる Non-Slip の境界条件（Bounce-back 則 12））を用いて多孔体壁を模擬する．また出口の圧力を一定(大気圧)とし，流路の上下左右４つの側面を対称境界とした．なお，流入速度(Uin)は 20cm/s とした．
次に，すすの境界条件であるが，入口では流入するすすの濃度 Yc,in の平衡値，出口では十分に発達しているとしてすす濃度の x 方向勾配を０とした．また，障害物に衝突したすすの分布関数は，壁面に衝突したすすが全て吸着するとして 0 に置き換えた．このようにして固体壁にすすが吸着していく様子を模擬し，すすが堆積した領域が広がることで固体壁が成長することを考慮した．
また，燃焼計算では，物体の表面にすすが均一（質量分率は 0.1）に吸着しているものとし，壁面の温度を 1200K と高い値に設定することによりすすを燃焼させた．燃焼反応によりすすの濃度がゼロとなった格子点は，次のタイムステップから流体が進入できる格子点とみなすことにした．

ObstacleLDExitInlet

図 2 Coordinate and obstacle regions.

3.

計算結果および考察

3.1 すすの堆積による流れの変化
濃度 0.1 ですすを流入させた場合の t = 0.13s のすす濃度と速度ベクトルを図 3 に示す．障害物の大きさ Dは 3 とし，吸着前の空隙率は 0.95 である．流路中央の x-y 断面における速度ベクトルとすすの２次元分布も同時に示した．また，異なる時刻 t = 0.19s の結果を図 4 に示した．これによると，堆積が進むにつれ流れが変化している様子が確認される．特に，流路前方ですすの堆積が顕著となっている．通常，DPF 内部で起こる吸着現象を直接可視化することができないため，粒子の吸着場所を実験により特定できないが，高速度顕微鏡カメラを用いて，DPF を輪切りにした部分を撮影した例がある 13)．それによると，フィルター表面壁面付近のくぼみにすすが吸着する様子が報告されている．この実験結果は，すすの堆積領域が流路前方となった本解析結果によく対応していた．
次に，t = 0.19s における圧力分布を，すすが堆積する前の結果と比較したものを図 5 に示す．すすが堆積している領域においては非常に大きな圧力降下があり，すすが堆積していない下流の領域では圧力降下が小さいことがわかる．また各時刻において，表面に吸着するすすの平均堆積量を求めた．流入させるすすの濃度を 0.1，0.05，0.01 と変化させて比較した結果を図 6 に示す．計算回数が進むにつれ，堆積量が急激に増加している．これは，すすが堆積する面積も拡大していることが原因と思われる．
また，今回の流れ場を理論と比較した．多孔体内の流れにおいて，流路の流動抵抗は理論的に検討されており，経験式として以下のような Ergun の式 14) がよく知られている．Ergun の式に使われている f は流動抵
抗，Re はレイノルズ数，Dp は等価直径である．本研究と同様の計算条件ですすがない場合の流れ解析を行った結果，この Ergun の式と計算結果がよく一致することを確認している．

f = 150 / Re + 1 . 75(16 )堆積が進んでいる過程で，どのようにレイノルズ数および流動抵抗が変化しているのかを調べた．一例として，図６と図７で示した時刻が 0.13s と 0.19s の結果を，Ergun の式と比較して図 7 に示す．ただし式 16 で必要な圧力勾配としては，図 5 で示したようにすすが堆積している領域内の急激な圧力勾配(Case1)と物体を配置した領域内における平均の圧力勾配 (Case2)の２条件で比較を行った．その結果，Case1 では Ergun の式より高い流動抵抗となった．逆に Case2 では流動抵抗が Ergun の式より低い値を示した．
Ergun の式は，定常を仮定して導出された経験式であり，すす粒子が堆積して時々刻々流れが変化しているような場では，この式は適応できないと言える．したがって，すす粒子が吸着している流れ場では，フィルター内部の平均的な圧力損失（圧力降下）をもとに Ergun の式から流動抵抗を単純に予測するだけでは不十分であり，本解析のような非定常計算が必要であることがわかった．

図3

Distribution of soot with velocity vector at t = 0.13 s.

図4Distribution of soot with velocity vector at t = 0.19 s.

p/pout1.8
after deposition
before deposition

1.6
1.4
1.2
1.0
0.8
00.5

1.0

1.5
x/L

2.0

図 5 Distributions of pressure with soot deposition.

Mean accumulated soot0.2
Yc=0.1

Yc=0.050.1

0

Yc=0.01

0.05

0.10

0.15

0.20
t, s0.25

図 6 Mean accumulated soot with different inlet condition.

f103102

Ergun eq.
Case1, t=0.25sCase1, t=0.37s
Case2, t=0.25s
Case2, t=0.37s

101

100 110102Re103

図 7 Variations of friction factor with Reynolds number compared with Ergun equation.

3.2

すすの燃焼解析

一例として，すす燃焼時の温度分布を図 8 に示す．図では物体の大きさ D が３（実際の長さは 1.5mm），空隙率εが 0.9 の場合の結果である．図の右側から空気が流入してくることにより，酸素が物体表面にあるすす近傍に輸送される．これにより壁面に付着しているすすは酸素と壁面上で反応し，壁面近くの温度は高く，酸素濃度は低くなる．図では３つの温度領域を示しており，入り口近くの温度が 600K の領域，その後温度が 1100K となる領域，壁面近くの燃焼が活発な 1210K の領域を示す．これによると，物体表面にあるすすが燃焼することで温度が上昇していることがよくわかる．
また，温度分布，酸素・すすの質量分率分布，反応速度分布を詳しく比較したものをそれぞれ図 9 に示す．
壁面を高温にして 0.63 秒後の xy 断面の結果である．時間がたつと燃焼反応によりすすの濃度は次第に減少していく．計算初期においてすすは物体表面に均一に分布させていたが，図 9(c)でわかるように，この時点で壁面に存在していたかなりのすすが，反応により消費されてなくなっていることがわかる．

図 8 Temperature field, ε=0.9, D=3.

DPF の再生過程において，局所的に急激な反応が起こると非常に温度の高い領域が存在することになり，DPF の内部構造が破壊されてしまう恐れがある．酸素濃度を変えることで，すすの反応速度を制御することができるため，本研究でも流入する空気中の酸素濃度を変えてすすの燃焼反応を計算した．すべての物体の表面における平均の反応速度を計算し，その時間変化を調べた．その結果を図 10 に示す．酸素濃度は，空気中の質量分率 0.233 と 0.1 で行った．これによると，全体的に酸素濃度が低い場合には反応速度が小さくなるものの，反応が長い間持続していることがわかる．これは酸素濃度が高いもののほうが早い段階ですすの酸化が完了することに対応していると考えられるので妥当な結果である．

図 9 Distributions of (a) Temperature, (b) Oxygen concentration,(c) coot concentration, and (d) reaction rate field; ε=0.9, D=3.

図 104.

Mean reaction rate of soot as a function of time with oxygen mass fraction of 0.1 and 0.233 in flow.

まとめ
格子ボルツマン法により，多孔体内部にすす粒子が堆積していく様子を３次元で解析し，DPF の目詰まり現象を模擬した．その結果，すすが堆積するにつれて流れは大きく変化し，圧力降下もより大きくなった．
すすが堆積していくと，吸着する面積も同時に増加するため，すすの堆積量も急激に増大した．特に，すすの堆積は上流部分で顕著であり，流入するすすの濃度が大きいほどすすの堆積量も増加した．また，粒子が目詰まりして流れが時間的に変化しているため，流路内の平均圧力分布から求めた流動抵抗は Ergun の式とは異なることがわかった．また燃焼計算では，すすの温度を高くすることにより，すすが周囲にある酸素と反応し，表面から徐々に燃焼してなくなっていく様子が計算できた．以上により，フィルターの再生過程が模擬できた．

参考文献
1) 浅海靖男，日本燃焼学会誌第 44 巻 127 号，25-32 (2002)．
2) 自動車 NOx・PM 法 (2002)．
3) 山本和弘，日本流体力学会誌ながれ 23, 295-302, (2004) ．
4) K. Yamamoto, N. Takada, and M. Misawa, Proc. Combust. Inst. 30, 1509-1515 (2005)5) S. Chen and G. D. Doolen, Annual Rev.Fluid Mech, 30, 329-364 (1998).
6) O. Filippova and D. Hanel, J. Comput. Physics B, Nos.1&2, pp. 197-200 (2003)．
7) B. Chopard, A. Masselot, and A. Dupuis, Computer Physics Communications 129, 167-176 (2000).
8) Y. H. Qian, D. d’Humie`res and P. Lallemand, Europhys. Lett. 17, 479–484 (1992).
9) K. Yamamoto, X. He, and G. D. Doolen, J. Statist. Phys., 7(1/2), 367-383 (2002).
10) 林秀光，R&D Review of Toyota CRDL Vol.38 No.1, 17-25 (2003)．
11) K. B. Lee, M. W. Thring and J. M. Beer, Combustion and Flame 6, 137-145 (1962).
12) Q. Zou and X. He, Phys. Fluids 9 (6), 1591-1598 (1997).
13) K. Hanamura and T. Tanaka, Proc. JSAE Annual Congress, 20055487, 23-26 (2005).
14) R.B. Bird, W.E. Stewart and E.N. Lightfoot, Transport Phenomena, Wiley, New York, (1960).