シミュレーション

広島国際学院大学研究報告，第 3
2巻 (
1
9
9
9
)，21-31
2
1
凹凸火炎面近傍における多成分分子拡散過程の
シミュレーション
(分子動力学法の適用の試み)
黄
樹偉
(乎成 1
1年 1
0月1
4日受理)
A
p
p
l
i
c
a
t
i
o
noft
h
eM Dr
n
e
t
h
o
dt
ot
h
eS
i
r
n
u
a
t
i
o
noft
h
e
MulticornponentGasD
i
f
f
u
s
i
o
n
nearaWrinkledF
l
a
r
n
e
l
e
t
ShuweiHUANG
(ReceivedOctober14，1999)
Ana
t
t
e
m
p
twas mede t
oa
p
p
l
yt
h
em
o
l
e
c
u
l
a
rdynamics methodt
o
t
h
es
i
m
u
l
a
t
i
o
no
ft
h
emulticomponentg
a
sd
i
f
f
u
s
i
o
nphenomenonn
e
a
ra
o
l
e
c
u
l
e
s
w
r
i
n
k
l
e
df
l
a
m
e
l
e
to
fs
t
o
i
c
h
i
o
m
e
t
r
i
candl
e
a
np
r
e
m
i
x
t
u
r
e
s
. Allm
i
nt
h
er
e
c
t
a
n
g
u
l
a
rc
a
l
c
u
l
a
t
i
o
nc
e
l
lweret
a
k
e
nt
ober
i
g
i
db
a
l
l
sandhave
J
ones p
o
t
e
n
t
i
a
l between any two m
o
l
e
c
u
l
e
s
. Among t
h
e
t
h
e Lennardfthemi
nt
h
ed
i
r
e
c
t
i
o
no
ff
l
a
m
e
s
i
xs
i
d
earoundt
h
ec
a
l
c
u
l
a
t
i
o
nc
e
l
l，twoo
p
r
o
p
a
g
a
t
i
o
nweret
a
k
e
nt
ob
em
i
r
r
o
rr
e
f
l
e
c
t
i
o
nb
o
u
n
d
a
r
i
e
sandt
h
eo
t
h
a
r
f
o
u
rt
obep
e
r
i
o
d
i
cb
o
u
n
d
a
r
i
e
s
. Thedynamicse
q
u
a
t
i
o
nwass
o
l
v
e
dby
u
s
i
n
gt
h
eV
e
r
l
e
tmethod.
Keywords
i
m
u
l
a
t
i
o
n，W
r
i
n
k
l
e
dF
l
a
m
e
l
e
t，Gas
M
o
l
e
c
u
l
a
r Dynamics，S
D
i
f
f
u
s
i
o
n
.
1.まえがき
エンジン内での予混合乱流火炎は，乱れにより変形された層流火炎片 (
F
l
a
m
e
l
e
t
) の集合とし
て考えることができる?層流火炎片での局所燃焼速度，変形および火炎伝播には，乱流特性のほ
流火炎の燃焼機構を解明する
か，混合気の組成ガスの拡散特性も影響を及ぼしている。従って，舌L
には，乱流により変形された凹凸火炎片での熱拡散と反応物質拡散過程を解明する必要がある O し
2
2
樹偉
黄
かしながら，火炎凹凸部の最小尺度が O.lmm以下であることから，実測は極めて困難で、ある O
一方，近年極く小さい空間内での物質構造の研究では分子シミュレーシヨンが有力な手法として
用いられ，その有効性が明らかにされている?従って，本研究では火炎凹凸部での分子拡散過程
への，この手法の適用を試みた。分子シミュレーションにおいては，現在モンテカルロ法と分子動
去の 2
つの手法があるが，モンテカルロ法では計算結果がサンプルの数や乱数列の性質に依存
力学 i
する 12)ので，本研究では分子動力学法を用いた。
2
. 分子動力学法の導入
2.1 計算の対象
ほとんどの分子動力学法によるシミュレーションは単純物質を対象としているが，本研究では，
系の構成粒子が多成分であるため，計算の対象を混合物質系まで拡張した。作成したプログラムは
プロパン，エタン，メタン，水素，酸素，窒素など多成分ガス系に適用でき，必要に応じ構成成分
の追加ができる。混合気中の燃料ガスはプロパン，エタン，及びメタンのどれか 1つを選択し，必
要に応じてさらに水素添加量(添加の体積割合で指定される)を指定することもできる O 図 1は計
算条件の設定画面を示す。図 1から，プログラムでは，未燃混合気の組成，温度，圧力，セル中の
全分子数を指定するだけで済む。既燃ガスの組成，温度は，定圧断熱完全燃焼として算出される O
図中の STEP係数は後述するが，保存間隔はこの値のステップ数ごとに計算結果を保存する設定
値である。また，計算セルの形状は長方形である O
随議選総録亘書一一一一一一一一一一一一一ー一一一ー一一ー一一一一一一一一一一一一 _
I
G
l
x
l
新規計算
」斗
データを保存する場所を指定して下さ~\
~!\\ì. 1
!日間DATA~C3H84附
炭化水素
'
"CH4
水素割合同
当愛比
1
¥
.
0
圧 力 匂 加 ) 1¥.0
温 度 以)1273
分子数(岱) ド蹴
STEP係数
105
保存間続
15
r C2H6
r C3H8
ヱ
。j 」
ユ
。j
間 三J
J
∞
ュJJ
1
ょJ
1
ょJ
2
ょj2
0
ょj抑
2 ~J
斗 4叩
∞1 ~..J
』百 5
1
ユ
j J
~1O
セjし
形i
丈守享さを 1とする)
高さ
長さ
i
司
1~
J
斗
1~
J
ょ]32
図 1 新規計算時の計算条件設定画面
16
凹凸火炎面近傍における多成分分子拡散過程のシミュレーション
2
3
2.2 粒子構造モデル
気相の分子拡散過程を対象としたため，分子同士の距離は原子同士のそれより遥かに大きく，分
子の内部構造および原子聞の相互作用による拡散過程への影響は無視でき，分子を変形しない剛体
粒子とした。
2.3 粒子同士の相E作用モデル
ここで粒子は分子で，分子同士の相互作用については，現在，完全剛体衝突とするハードコア型，
斥力場だけを持つソフトコア型，及び引力と斥力場を有する Lennard-Jonesポテンシャル型 (L
Jm
!
. 図 2参照)などのモデルが用いられており t3i 本研究では，モデル定数の実験データが比較
的豊富な L-J型モデルを用いる。粒子の受ける合力は力切断距離以内の粒子からそれぞれの作用
の総和として考慮する。即ち， η にある
j 番目の分子から
η にある
1 番目のものにはたらくポテ
Y
i
j)と力 F
i
}は，それぞれ次式で与えられる。
ンシャル φ(
φ(Yij)=
(
1
)
3
( b ¥ U ) @ 2 0比
パ
U
2
市 山 山 崎
P
o
t
e
n
t
i
a
l
一
一 Force/e
o
i
i
/
ア
日
)
恒
吉一 1
・
.
o
f
/
3
!.I
J
a
.
2
¥I
3
1
.5
2
2
.
5
C
e
n
t
e
rS
e
p
a
r
a
t
i
o
n (σ)
0
.
5
3
図 2 LenardJ
onesポテンシャルと相互作用力
F
i
jニ
-
v
'i
φ(
Y
i
j)
(
2
)
F
z
=
F
t
f
ここにみ二=入
入 Y
i
j
二
(
3
)
I
;
u
lで
V
'iは子z方向の勾配である o Fiは t
番目の分子が受ける合力であ
i
jとη
σ について，分子対
るO また， ε
1 と 1が同種物質の場合，
L-J方程式中の ε (エネルギ一
長さ次元を持つモデル定数で，コア直径と称され，通常，分子直径
次元を持つモデル定数)と σ (
とされる)をそのまま用い，異種分子対では，次のように与える 141
ムZ
(
4
)
σrt(σi+σJ
(
5
)
ε!J
2
4
黄
樹偉
2.4 運動方程式の解き方
式(
2
)と(
3
)より，時間 tでの全ての粒子の位置から i番目の粒子が受ける合力がわかるので ，d
t
後の粒子の位置は次式より求まる O
F
;
dZ;，
d
tZ m;
(
6
)
式 中 間 は i番目の分子の質量である O 式(
6
)の差分化による計算は，窒素分子のパラメータを基準
e
r
l
e
t法¥3そ用いて行った。
値に無次元化し， V
2.5 境界条件
プログラムの有効性を確認する場合，境界条件は凹凸火炎面なしで基本セルの 3方向ともに周期
境界条件またはともに鏡面対称境界とした。実際のシミュレーションでは図 3に示す境界条件とし
た。図中，左右の境界だけが鏡面対称境界で，その他の境界は周期境界である O
既燃側
未燃側
対
利γーーーーー一ーーーーーー
火炎伝ぱ‘
鏡面対祢
z
i
/
凹凸火炎面
図 3 境界条件の設定
2.6 初期条件
系全体が均一混合気の場合，各組成成分の全ガス分子の平均温度が等しく，各成分の分子速度が
正規分布となるように，粒子の質量に応じて初期速度の大きさを与える O また，初期位置と運動方
向は統計的に均ーかっ等方になるように与える O 最後に，系の質量中心の初期運動量がゼロになる
ように各粒子の初速度を再調整した O
系が既燃側と未燃側に分かれる場合，系の圧力が同一で，両側の体積が等しい条件でそれぞれの
温度とガス組成に応じて各種の粒子の数を決定する O
3 プログラムの設計と主要パラメータの決定
分子動力学法によるシミュレーションの特徴は計算量が大きいことと，データ量が多いことであ
る。系の粒子数を N とすれば，計算量は既ねがに比例する。大型コンピュータを用いても，数千
個の粒子を初期条件から実用的な時間まで計算するには長い CPU時間を要する O 特に本研究では
PCを用いるため，より長い計算時間を必要とする。これらのことを考慮して，次のようなプログ
ラム仕様を決定した。
3.1 プログラムの主要仕様
(
1
) 断続計算ができること
O
このことから，初期設定条件，中間計算結果をファイルに保存するほ
かさらに，計算が中断された場合，中断時の状態をそのままファイルに保存する O 継続計算時に
2
5
凹凸火炎面近傍における多成分分子拡散過程のシミュレーション
はこれをファイルから読み出す。
(
2
) 計算とデータ処理は別々のプログラムで行うこと 。これにより，一度にメモリにロードするコー
ド量が小さく，より大きいメモリ空間をデータの配置に残すことができる 。
(
3) プログラムが操作しやすいこと 。分子動力学法で扱う物理量は通常の物理量とはかなり異なる
ため，入力ミスを回避するため ，ほとんどの必要な物理量はデータベ ースとしてプログラムに組
み込んでおり ，図 1に示すように，簡単な入力で計算を行うことができる 。
3.2 プログラムの構造と入出力フ ァイル
以上の仕様により，プログラムの構造を図 4に示すように設計した。 図 4から ，各モジュ ール聞
のデータ交換はファイルを介して行っていることがわかる 。
保存ファイル
主 要 モ デュ ー ノ
レ
バ計算条件
じ
初期設定条件
入力部
ト ップ メ ニ ュ ー
中断時の状態
L
A
ST
.S
A
V
動力学計算部
粒子位置
XYZ.D
A
T
分子運動
表示部
I
X" I
粒子速度
V
X
Y
Z.D
A
T
: 司エネ
ル
"
， '- 統 計 値
S
U
M.
D
A
T
図4
l
CONFIG.INI
II
I
I
プログラムの構造と使用するファイル
3.3 主要パラメータの決定
(
1
) 切断距離 r
c
r
cは 2粒子問の相互作用が実効的になくなると見なせる粒子問の最短距離であり，これを大きく
すればより正確に計:
算できるが，計算量が増加する 。図 2からわかるように，九は分子のコア直径
σの 2-4倍程度である 。本研究では混合ガスを主要対象とするので， これを窒素分子のコア直径
の1
0倍と大きめにとったが， 273K， O
.lMPaでは平均して約 1
.5
個の分子がこの半径九の球の中に
入るので，計算量の増加は小 さい。
(
2) 時間ステップ L
1tと STEP係数
運動方程式(
6)
を解く際，L
1t
が大きすぎると解が不安定となり，小さすぎると変位量子+.
1;の計
算時には桁落ちの問題が起 こ りやすい。 同じ L
1t
でも ，系の温度，圧力及び粒子の質量により解の
安定性 (
断熱系では系の運動エネルギー とポテンシャルエネルギーの和であるハミ jレトニアンが保
存されれば解が安定である )が変わる 。従って，計算条件毎に L
1t
の適正範囲が存在する 。 これを
調べるため，基本セル中を対向運動する 2粒 子 (
水素分子)の場合の，系のハミルトニアンの時間
2
6
黄
変化を計算して ，L
1tの適正範囲を得た。
樹偉
T を水素分子が与えられた温度における平均速度で自分
01-5MPa，温度が 2000K以下の範囲内
の半径分ほど移動するのに要する時間とすれば，圧力が 0.
1t= (
0.001-0.
2
)r時の解が安定であることが確認されたので，L
1t=(
0
.
0
5ー0
.
1)rとした。
で，L
0.05-0.
1)を，本研究では STEP係数 (
図 1参照) と称する 。 また，実際のシミュ
ここの係数 (
レーション時には，系のハミルトニアンを計算しているので，解の安定性を随時に!照合できる 。
4 試行計算結果
本研究では与えられた初期状態粒子分布から，計算終了 (
いつでも継続計算できるが) までの各
粒子の運動軌跡を記録することにより，分子運動を再現することができる 。図 5と 6は，それぞれ
0
0
0ステ ップ目の分子位置を示し，初期設定条件は図 7に示す。 また，図 8に示してい
初期状態と 3
•-
図 5 初期状態の 1例 (
計算条件は図 7に示す)
・
・
田
園
・
・
・
・
・
・
•・
図 6 計算の途中結果の 1例 (
3
0
0
0ステ ップ目，計算条件は図 7に示す)
2
7
凹凸火炎面近傍における多成分分子拡散過程のシミュレーション
計算条件
i
'
i
i
2
炭化水ヨ"ti8(ー〉
来燃倒分子構感
水素割合{石一 一 (%)
α嗣 8 f
a
7
一
一
H20 ~
ゼ ル 幅 下 在 日-C
m
l
当量比 「 一 ← }
H2
f
a
7
一
一
C02 同 一 一
セ)~高さ戸8E-08
01
[石「ー
制
制空
同59
H2
r
o
-
切断臨ド荏
02
f
3
5
一
E
寺問S'
e
p
N2
同τ一
E
力 『 一 一 白加}
未 燃 温 度 防γー は }
全分子畝
既餓倒分子情戚
1
4
0
0
0一{⑩
τ
「←}
セル形状医
8
r
o
:
セ J~.IJê
Cm
l
fT9妄o8Cm
l
a Cml
1
2
_
5
8E
1
3 ¥
S
e
C
J
S
t
e
p
鮒 同 一- (-l
図 7 図 5と 6の計算条件
100
33280
80
33260
60
33240
40
(
'
"33220
z
20
'
" 33200
。
533m
-20
均四
P
o
t
e
n
t
i
a
lEnergy
O
c
w
33160
-40
33140
u0
33120
-80
1
∞o
2
A生
nU
n
u
n
u
数
子
L
vu
変
・
目
勾
，
問
時
の
、
ネ
ギ
レ
エ
の
系
図8
∞o
関口
。申 位
C
0
田島
33100
(NZUFO 旬 ) 同 国LωC凶
33300
-100
2500
∞o
3
d
t= 2.
5
8E1
3
s
ec
)
るように，系の全エネルギーはほとんど変化せず，解が安定していることがわかった。 図 5と 6か
0と1
1に示
ら抽出した混合気中の水素分子と酸素分子並びに局所当量比の分布を，それぞれ図 9， 1
す。
5
. 計算結果の検討
a)
からわかるように，水素分子の数が少ないため，初期状態では，未燃側の空間に等確率的
図 9(
に配置しでも ， あまり均等になっていないことがわかった。 これを改善するために，全分子数を大
b)
からは， 3
0
0
0ステップ目には水素分子の分布は初期状態
幅に増大させる必要がある 。 また，図 9(
より幾分均等になったことがわかった。 しかし，火炎面に沿って，未燃側に突き出た火炎(図中の
下方で，以降，凸状火炎と称、し，逆の場合は凹状火炎と称する ) と凹状火炎とでは，水素拡散の差
は未だ不明確であるため，今後分子総数を増やすとともにステップ数に比例する動力学時間を長く
する必要がある 。
0では ，酸素分子の密度分布は水素分子のそれに比べ幾分均等になっていることを見受けたも
図1
のの， まだ不十分と考えている 。火炎面近傍では，酸素分子密度の勾配は動力学時間につれて小さ
2
8
黄
樹偉
くなったことがわかった。 凸状火炎面の近傍では，酸素分子が既燃側に多く拡散したことがわかっ
た。 このことは，そこでの局所当量比が変化していることの可能性を示唆している 。
図1
1によると，局所当量比の変化はまだ十分明確なものではないが，凸状火炎下方の既燃側近傍
では，幾分増大したことがわかった。 このことから，凹状火炎に比べ，反応物質の不足成分 (
希薄
混合気では燃料分子)の凸状火炎への拡散がより速いことを意味する 。
Di
s
t
r
i
b此 i
ono
fH2NumberDensi
tyOn
i
ti
a
ls
t
a
t
e，d
i
m
e
n
s
i
o
n
l
e
s
s
l
。
2
3
2 345 6 7 8 9
(
a) 初期状態
、
D
is廿 i
bu
t
i
o
no
fH2NumberDens
i
ty(
3000stepsa
f
t
e
r
，d
i
m
e
rsi
onl
ess
l
234 5 6 789
(
b
)3
0
0
0ステップ目の時
図 9 水素分子の無次元数密度の分布 (
図 9-11の軸の数字はサブセルの座標)
凹凸火炎面近傍における多成分分子拡散過程のシミュレーション
2
9
D
i
s
t
ri
b
u
t
i
o
no
f02NumberDensi
t
y(
l
ni
t
i
a
lS
t
a
t
e.Di
mens
i
o
n
l
e
s
s
l
2~
a56
7
r
1
81
9 2~1
7 8 9 1
01
11
'
;
-1
'
31
41
'
51
'
;
;
'
2
2 ;~4~'57s'2T;'8 29
(
a) 初期状態
∞ostepsafter.dimensionlessl
Di
s
t
r
i
b
u
t
i
o
no
f02NumberD
e
n
s
i
t
y(
3
23456 7 8 9
(
b
)3
0
0
0ステップ自の時
図1
0 酸素分子の無次元数密、度の分布
d
o5
黄
30
樹偉
、
D
i
s
t
ri
b
u
t
i
o
no
fE
q
u
i
v
a
l
e
rtRat
i
o(
In
i
t
i
a
lS
t
a
t
e)
.1
.2
-1
.4
・
ロ1ー1.2
0.
8
1
ロ0.
6-0.
8
ロ0.
4-0.
6
・
0.
20
.
4
ロ0-0.
2
。
2
3
5
2 3 456 7 8 9
(
a
) 初期状態
∞oSteps杭 er)
Di
st
r
i
b
u
t
i
o
no
fEqui
v
a
l
e
n
tR
a
t
i
o(
3
.1
.
2
ー
1
.4
・
ロ1
1.
2
0.
8
1
00
.
6-0.
8
ロ0.4-0.
6
0.
20
.
4
・
2
3
5
2 3 4 5 6 7 8 9 1
01
1
(
b
) 3000ステップ目の時
図1
1 局所当量比の分布
6
. あとがき
本研究で作成したプログラムは，凹凸火炎面近傍での分子拡散過程をシミュレ ー トすることは一
通りにできたが，以下のいくつかの点を改良 ・再検討する必要がある 。
(
1
) 実際の火炎面近傍での分子拡散過程には化学反応を伴うので，これをより真実に再現するには，
反応性分子が未燃側から既燃側へ移動したとき，既燃物質の分子になることを考慮する必要があ
る。 この場合，系の粒子数が変化しうるアルゴリズムに変える必要がある 。
凹凸火炎面近傍における多成分分子拡散過程のシミュレーション
3
1
(
2
) 境界条件について，ほとんどの研究では，周期境界か自由境界を用いているが，火炎伝播の場
合，伝播方向ではこれらの条件は適用できなし、。本研究では鏡面対称としたが，計算セルの火炎
伝播方向での長さを十分長くする必要があった。今後より実際に近い境界条件を検討したい。
(
3
)
凹凸火炎面近傍での分子拡散過程をシミュレーションする目的は，凸状火炎面と凹状火炎面と
では，拡散過程がどのように違うかを調べるためである。しかしながら，現在では PCを用いる
ため，実際の凹凸火炎の尺度より遥かに小さい尺度のセルしか計算できない。この尺度の差によ
り計算結果にどんな影響を与えているかは未だ未調査である。今後， WS環境でのシミュレーシ
ヨンを行う予定である O
(
4
)
計算結果の抽出について，本報では，温度分布の抽出は行っていないが，今後これを追加する
予定である。
7年度の黄研究室の卒業生諸君に謝意を表する O
最後に，本研究に協力した 9
参考文献
1.城戸・黄，燃焼研究， N
o
.
9
2， 1
9
9
3，p
l
2
0
2
. 上回顕，コンピュータシミュレーシヨン，朝倉書!苫， 1
9
9
0
3
. 大津映二，片岡洋右，分子動力学法とモンテカルロ法，講談社サイエンテイフイク， 1
9
9
4
.
C
.，P
r
a
u
s
n
i
t
z，1
.
M
. andSherwood，T
.
K
.，TheP
r
o
p
e
r
t
i
e
so
fGasesandL
i
q
u
i
d
s，
4
. R
e
i
d，R
McGraw-Hill，1997