Title 引力を考慮したDEM法による多孔体構造作成シミュレー ション(粉体

Title
Author(s)
Citation
Issue Date
URL
引力を考慮したDEM法による多孔体構造作成シミュレー
ション(粉体物理の現状と展望,2006年度後期基礎物理学
研究所研究会)
林, 秀光
物性研究 (2007), 88(2): 297-300
2007-05-20
http://hdl.handle.net/2433/110795
Right
Type
Textversion
Departmental Bulletin Paper
publisher
Kyoto University
「粉体物理の現状と展望j
引力を考慮した DE
恥f法による多孔体構造作成シミュレーション
株式会社豊田中央研究所反応制御研究室
林 秀 光1
1 はじめに
ディーゼ、ル車から排出される粒子状物質 (PM)を除去するためのフィルタとして， D
P
F
(
D
i
e
s
e
l
P
a
r
t
i
c
u
l
a
t
eF
i
l
t
e
r
)が考案され，現在急速に普及しつつある. DPFの構造は複数提案されている
a
l
l
t
h
r
o
u
g
h型である.この型の DPFでは，エンジン側から流入した排気ガ
が，現在の主流は W
スは，多孔体で形成される流路壁を通過して初めて大気側へ流出できる仕組みになっている.排気
ガス中に浮遊している PMは，排気ガスが流路壁を通過する過程で多孔体内部に捕集される.こ
の様な WaUt
h
r
o
u
g
h型 DPFの動作原理から考えて，捕集効率や圧力損失など DPFの主要特性
凶
は，多孔体の構造に強く依存すると予想される.
iC粉末を焼結して作成される SiC-DPFを対象とし，紛体のパッキングにより流
本研究では， S
路壁の多孔体を作成した.粉のパッキング問題の解法としては，従来は，粉の幾何学的形状のみを
考慮してランダムにパッキングする事が行われていたが [
1
1，最近は， DEM法 (
D
i
s
t
i
n
c
tElement
Method)[
2
]により粉の運動を動力学的に追跡する事で現実に近いパッキングが可能になってきた.
ところで，粉の粒径が 100μm以下になると，パッキングの結果として作成される多孔体の P
o
r
o
s
i
t
y
が粒径の減少につれて増大する事が知られている.この現象は，粉の幾何学的形状のみを考慮し
た方法では再現できないが，粉と粉の聞に v
a
nd
e
rWaals力を考慮、した DEM法を用いると再現
できる事が示されており [
3
]，動力学的な多孔体作成方法の優位性を示している.そこで本研究で
and
e
rWaals力を考慮した DEM法を用いた紛体のパッキングにより多孔体を作成した.
は
， v
2 引力を考慮した D E
乱f法
DEM法は，弾性変形とエネルギー散逸を考慮して剛体の運動を追跡する計算手法である.その
基本方程式は，同IJ体に対する Newtonの運動方程式であり，並進運動と回転運動に関してそれぞ
[
2
1
.
dVi
L•生~
=T
.
d
t
m
.
_
-= .
l
1i
:
1
- d
t
明
.L
'
t
，
、‘，， r
i
咽
r
，，‘、
れ次の通りに書ける
ここで， miと Lは粒子 iの質量と慣性モーメント， V
iと ω4は粒子 iの重心の並進速度と重心周り
.
の角速度， F
iと Tiは粒子 iに働く力と重心周りのトルクである(図 1参照 )
粒子関で衝突を繰り返しながら重力により落下する粒子を考察する場合，粒子 iに働く力民と
トルク Tiは，次のように書ける.
F
i
=乞 (
F
i
j+時十時)+mig
E
-m
a
i
l
:e0905@mos
k
.t
y
t
l
a
b
s
.
c
o
.
j
p
.URL:h
t
旬:
jj
w
w
w
.
t
y
t
l
a
b
s
.
c
o
.
j
p
j
l
(
2
)
門
i
QU
“
っ
研究会報告
図1
:DEM法の方程式に含まれる各種ベクトル
i =乞 (
T
T
む+
T
L
)
(
3
)
J
ここで，F
i
j
. F，
む FL，mzgは，それぞれ，粒子 iが粒子 jから受ける接触カの法線成分，接触力
の接線成分， vand
e
rWaals力，粒子 iに働く重力である.また， Tむと Tらは，それぞれ，粒子 i
が粒子 jから受けるトルクと転がり摩擦である.
DEM法では，粒子 iが粒子 jから受ける接触力の法線成分
F
2を次式で表す.
F
2
=
[
;
E
J
E
d
/
3 rnEVR/
[
;
;(
い
(
4
)
一
ここで，Eは Young率 Y と P
o
i
s
s
o
n比 σを使って E = Y
/
(
lー σ2)と定義される量であり ，Rは
，
iRj/(九 +Rj)と書ける.また，ふは，
粒子 iと粒子 jの半径九と半径 Rjの調和平均であり長 =R
粒子 i と粒子 j の重なり長さを表しふ =max(~+ 均一 Irj - ril，
0) と定義されている .r
nは，法
線方向のダンピング係数である. Vijは粒子 iと粒子 jの重心の相対速度であり
る. llijは，粒子 jの重心から粒子 iの重心へ向かう単位ベクトノレで、あり，
Vij= Vi-Vj と書け
l
l
i
j= (ri- rj)/(
1ri- rj)
1
と定義される.
粒子 iが粒子 jから受ける接触力の接線成分 Fらは次式で与えられる.
F宇一一 μIF~lmin (
c
sと
，s
，
m
a
x
) V;j
一円叫
ん
ax
(
5
)
I
V
f
j
l
ここで， μ は摩擦係数，ごsは粒子 iと粒子 jが接触してから現時刻までの接触点の積算移動長さ，
最大積算移動長さと叩眠は
C
s，
m
砿
=
μ
[
(
2ー σ
)
/
2
/
(
1ー σ
)
]ふと定義される. V~j
は，粒子 i と粒
子 jの接触点の相対速度の接線成分であり，
・
+(
ω
4十 時 )xR
i
j
(
6
)
VL=vtj一(Vij l
l
i
j
)l
l
i
j
で与えられる.ここで，
~j = R
i
(
r
j-ri)/(~
R
i
j は，粒子 iの重心から粒子 jとの接触点へ向かうベクトノレで、あり，
+Rj)と書ける.接触点の積算移動長さとs の算出方法は，文献に詳しく書
かれていないが本研究では，運動方程式 (
1
)を積分する際の時間刻みを d
tとして，
値を粒子 iと粒子 jが接触を開始した時刻から積算して求めた.
- 298-
V
f
tの絶対
jd
「粉体物理の現状と展望j
vanderWaals力 FLは，次式で与えられる(文献 [
3
1およびその参考文献を参照).
64R?Rj(九十 ~+Rj )
I
4h+2
'
Rj
h+4~Rj)2 UZJ
F
Y
.
:= _Ha
り
ι
2十 2
6(
h
j
h
)
2
(
h
，2+2
Rih+2
R
(
7
)
ここで，Hα
，は
， Hamaker定数である. hは，粒子 iと粒子 jの表面聞の距離であり，カットオフ
距離 hm
i
nを考慮して，次式で算出される.
h= m
a
x
(
l
r
j-ril-(
払 +R
j
)，
hm
i
n)
粒子 iが粒子 jから受けるトルク
(
8
)
T
i
jは
，
Tむ =~j
xF
む
(
9
)
と書ける.また，転がり摩擦 Tらは，次式で導入した.
T:-ZJ ニーμr~ I
F
r
，
l，
Wι
1
Z
J1
I
叫
rl-~.
(
1
0
)
運動方程式 (
1
)は，力が速度に依存するため，数値的に積分する際には工夫が必要である.本研
P
r
e
d
i
c
t
o
rC
o
r
r
e
c
t
o
rV
e
l
o
c
i
t
yV
e
r
l
e
tMethod)を用いた.
究では，運動方程式の積分に PCVV法 (
PCVV法のアルゴリズムは以下の通りである 時刻 tでの粒子 iの位置，速度，角速度を r
j
O
)，v
j
O
)
?
J
)とし，これらを用いて計算した粒子 iに働く力とトルクを F ?とTjO)とする
まず，次式に
j
1
)，v
j
1
)
?
ω
j
1
)を求める.
より位置，速度，角速度の中間的な値 r
rj1)=rjO)+d
付 j hd
よ~r九
YZ jO)
(
1
1
)
t
F
j
O
)
vj1)=vjO)十 d
(
1
2
)
m
~l) =ω4(0)十 dtT4
(
0
)
I
(
1
3
)
次に，これらの中間的な値を使って求めた各粒子に働くカとトルクを F
j
I
)とT
j
1
)とする.さらに，
次式によりす)，
J
)を求める
d
t(
V
}
2
) = V
}
l
)+ _
F
)
O
)+F
}
l
)
)
η ¥
I
(
1
4
)
d2)=ω}1)+
(
T
)
O
) T}川
.笠
2 1 ¥+
- ZJ
(
1
5
)
乙 ，1
-~
j
1
)，v
f
)
?
ω
j
2
)とする.
最終的に，時期l
tt+dtで、の粒子 iの位置，速度，角速度を r
3 球のパッキングシミュレーション
本研究では，粉の形状は球とし，球の大きさはすべて問ーとして数値計算を実施した.数値計算
に用いた各種パラメータの値は，文献 [
3
1と同じ値を使用した.重力と垂直な方向の境界条件は周
期境界条件とし，その周期は粒子の直径 dの 1
0倍とした.重力に平行な方向の境界条件は，底面
-299-
研究会報告
は剛体壁(物性値は球と同一)とし，上面は自由境界条件とした.数値計算の初期条件は， P
o
r
o
s
i
t
y
が0
.
9
4となるように粒子をランダムに配置し，各粒子の速度と角速度はゼ.ロとした.
球の集団を自由落下させ，運動エネノレギーを失った後の球の集団は多孔体となる.多孔体の
P
o
r
o
s
i
t
yは，球の集団が含まれる領域の体積と球の外部となる空間の体積の比として容易に計算
o
r
o
s
i
t
yの値を他の計算や実験と比較
される.本研究のシミュレーション結果を用いて算出した P
(
a
)に示す.粉の粒径が大凡 100μm以下になるとパッキングの結果として得られた
した結果を図 2
o
r
o
s
i
t
y値が細密充填の場合の値 (
0
.
3
6
)より増加する事が知られている.これは， van
多孔体の P
(
a
)から，本
d
e
rWaals力の大きさが重力と比べて無視できなくなるためと考えられている.図 2
o
r
o
s
i
t
yの増大を
研究で開発された DEM法の計算機プログラムが，粉の粒径の減少にともなう P
正しく再現する事が分かる.
次に，パッキングシミュレーションの結果として得られた多孔体の細孔径分布を求める.ここで
は
， G
elbらの文献 [
4
Jに述べられている方法を用いた.この方法の要点は，次の通りである. (
1
)
細孔中のある点 (
a点とする)での細孔径を求めるとする. (
2
)細孔中の任意の点を中心とし a点を
3
)これらの球の半径の最大値を a点での細孔径
含み細孔壁と交わらない半径が最大の球を描く. (
(
b
)に示す.
とする.このようして求められた細孔径分布の計算結果を図 2
0
.
9
0
.
5
AMilewski
oF
e
ngandYu
A Y
ange
ta
l
.
・ 羽l
i
swork
0
.
8
>.0
.
7
e
‘
c
o
吉 0.
4
D
き0.3
。
4
0
.
6
τ
コ
0.
5
器 0.2
0.
4
‘
-0
~
.
1
0
仏
@
0
.
3
1
0
∞
o
1
α
)
(
)
10
100
PO開 s
i
z
e(μm)
P
a
r
t
i
c
l
ed
i
a
m
e
t
e
r(μm)
1000
図2
:球のパッキングシミュレーションの結果
(
a
)
P
o
r
o
s
i
t
yの球の直径依存性に関し，本研究の結果が， Yangらの計算結果 [
3
]及 び
M
i
l
e
w
s
k
iや臨時らの実験結果と比較されている. (
b
)粉の粒径が 20μm，50μm，100μm，
1000μmの場合に細孔径分布の計算結果が示されている.
参考文献
[
l
J 例えば; C
.Pan，M.H
i
l
p
e
r
t，andC
.
T
.Mi
1
l
e
r，P
h
y
s
.R
e
v
.E64(
2
0
0
1
)，0
6
6
7
0
2
.
[
2
]P
.
A
.C
u
n
d
a
l
l叩 dO
.
D
.
L
.S七r
a
c
k，
G
e
o
t
e
c
h
n
i
q
u
e29(
1
9
7
9
)，
4
7
.
.
P
.Zou，andA
.
B
.Yu，P
h
y
s
.R
e
v
.E62(
2
0
0
0
)，3ωO.
[
3
]R
.
Y
.Yang，R
[
4
]L
.
D
.GelbandK
.
E
.G山 b
i
n
s，
La
時 m
uir15(
1
9
9
9
)，
3
0
5
.
ハ
u
n
u
qJ