すす粒子モーメント計算へのガウシアン過程回帰の導入
研究代表者 九州大学大学院工学研究院機械工学部門 助教 安藤 詩音
1．緒言
多環芳香族炭化水素（PAH，Polycyclic Aromatic Hydrocarbon）やすすは気候変動を引き起こし得るのみならず，人体の呼吸器官に対しても悪影響をおよぼし得るため，排出の低減が強く求められている．これらは主に自動車用エンジンや産業炉内における炭化水素
燃料の不完全燃焼によって生じる．したがって，すす粒子の生成メカニズムを明らかにすることは重要であり，これまで多くの研究が進められてきた[1–3]．
すす粒子は PAH からの核生成の後，凝集や表面成長を繰り返すことで成長する．ナノ粒子の成長モデルの理論的枠組みは Smoluchowski[4–6]によって初めて構築され，粒子の凝集過程が Smoluchowski 方程式と呼ばれる支配方程式によって記述された．その後エア
ロゾルの研究領域を中心に洗練され，エアロゾルの速度分布や粒径分布が統計的に記述された[7]．同時に，工学分野でも応用されるようになり，Williams は Williams 方程式を構築し，衝撃噴流からの衝突や空力抵抗により成長する噴霧火炎中の液滴の挙動を記述した[8]．
Smoluchowski 方程式は無限個の微分方程式によって粒子集団の時間発展を表す．ゆえにこの方程式を厳密に解くことは不可能であり，いくつかの解析手法が提案されている．
モンテカルロ法は最も精度の高い計算手法の一つであり，stochastic particles と呼ばれる
ランダムにサンプリングされた粒子の集団の成長を追跡することで，粒子の粒径分布関数
を取得する方法である[9]．しかし，精度の高い粒径分布関数を取得するためには，多数の粒子を追跡する必要があり，計算コストは高い．ゆえに，予混合層流火炎などに適用が限定されている[10]．一方で，モーメント法は粒径分布関数を求める代わりに，数次の粒径分布関数のモーメントのみを算出する．そのために計算負荷が小さく，すす粒子の粒径分
布関数を評価する手法として最も広く用いられている[11]．例えば，Appel らは C2 炭化水
素燃料の層流予混合火炎に適用し，バーナー軸上のすす体積分率を高い精度で再現できることを示した[12]．また，Wang はエチレンの対流拡散火炎中でのすす生成のシミュレー
ションを実施した[13]．彼らは実験との比較から，すす表面成長におけるアセチレンの効
果の大きさを示した．近年では，3 次元でのディーゼルエンジン内のシミュレーション[14]や多成分燃料の燃焼シミュレーション[11]においても適用され始めている．
モーメント法の重要な特徴は凝集や表面成長に伴うソース項が輸送モーメントの関数になっておらず，分数次モーメントや負モーメントに依存することである．言い換えれば，これがクロージャー問題を生むことである．このソース項を閉じる手段によってモーメン
ト法は MOMIC（Moments with Interpolative Closure）[15,16]と QMOM
（the quadrature
method of moments）[17,18]の 2 つに分類される．MOMIC は Frenklach[16]によって提－ 247 －

唱された手法でラグランジュ補間を 2 回用いることによって方程式を閉じる（詳細は後述）．MOMIC は計算効率性のみならず数値ロバスト性に長けており，3D-CFD との組み合わせにも適している[19]．一方で，補間手法により，計算精度が大きく影響を受けることが指摘されてきた[15]．初期の MOMIC モデルでは負の無限大オーダーのモーメントを用いて，負の分数次モーメントを補間したが[15]，後のモデルでは負の分数次モーメントを 0,1,2 次モーメントで補間する手法が提案されている[20,21]．これらの手法の精度がWick らによって検証されている[19]．しかし，補間スキームによる精度検証については十分実施されていない．ラグランジュ補間は考慮するモーメントのオーダーを多くするとル
ンゲ現象と呼ばれる数値不安定性を示すことも指摘されている[22]．これに対処するため，
Roy らはラグランジュ補間の代わりに Stineman 補間を用いて分数次モーメントを算出した[22]．その結果，ルンゲ現象を回避できる可能性はあるものの，精度としては大きな改善が見られなかったと結論付けている．
そこで，本研究ではラグランジュ補間の代わりにガウシアン過程（GPR）を用いること
で分数次モーメントおよび凝集項を算出する．ガウス過程はノンパラメトリックな回帰手法であり，非線形モデルに対しても適用可能である．したがって非線形的な関数の出力を推測する手法として，燃焼研究でも応用され始めている[23]．その特徴として，目標変数の期待値に加えて，信頼度が算出される．ゆえに，輸送モーメント補間時における誤差を評価することが可能となる．本研究では炭化水素燃料の 0 次元熱分解計算に MOMIC とガウス過程を組み合わせた新モデル(MOMIC-GP)を適用し，凝集項計算の誤差評価を実施した．計算結果は Alexiou らの衝撃波管実験と比較した[24]．
2．数値計算2-1．MOMIC
Smoluchowski 方程式からの MOMIC の導出過程については Frenklah[16]に詳細に解説されている．ここでは概略を述べる．
𝑟𝑟次モーメントは式(1)のように定義される．
�𝑀𝑀� � � 𝑚𝑚�� 𝑁𝑁�

(1)

���

核生成，凝集および表面成長を考慮したすす粒子のモーメント方程式は式(2)のように書け
る．
𝑑𝑑𝑑𝑑� 𝑑𝑑𝑑𝑑  𝑅𝑅�  𝐺𝐺�  𝑊𝑊�

(2)

ここで，𝐺𝐺� は凝集項であり，粒子クラス𝑖𝑖=1 から∞までの Smoluchowski から数学的に厳密に式(3)のように変形される．

�

���

�

�

��� ���

𝐺𝐺�   � ��� �� � 𝑚𝑚�� 𝑚𝑚���� 𝛽𝛽�� 𝑁𝑁� 𝑁𝑁� �(3)ただし𝛽𝛽�� は衝突係数である．粒子が希薄に分布し自由行程が十分長いとみなせるならば式－ 248 －

(4)のように書ける．
�𝐺𝐺�  �������𝑚𝑚�  𝑚𝑚� �

���

��







�



𝑚𝑚�� � 𝑚𝑚�� � �𝑚𝑚�� �  𝑚𝑚�� � � 𝑁𝑁� 𝑁𝑁�

(4)

式(4)は式(3)に代入して以下のように書ける．ただし，式(4)は(mi+mj)1/2 項を含むため，これ以上展開することができない．これに対処するため，MOMIC では補間を用いる．ここで gird function を以下のように定義する．
�𝑓𝑓�  �

���

�

�

���

�









�

�𝑚𝑚�  𝑚𝑚� � 𝑚𝑚�� � 𝑚𝑚�� � �𝑚𝑚�� �  𝑚𝑚�� � � 𝑁𝑁� 𝑁𝑁�

(5)gird function に対して，𝑙𝑙が非負整数の時，�𝑖𝑖  𝑗𝑗���項は展開できる．例えば，𝑙𝑙  であれば式(6)となり，分数次オーダーのモーメントの和積で表すことができる．
�𝑓𝑓�  �

���

�

�













𝑚𝑚�� � 𝑚𝑚�� �  𝑚𝑚�� � 𝑚𝑚�� �  𝑚𝑚�� � 𝑚𝑚�� � 𝑁𝑁� 𝑁𝑁�  𝜇𝜇�� 𝜇𝜇��  𝜇𝜇�� 𝜇𝜇��

���

(6)

MOMIC ではまず，この分数次オーダーのモーメントを整数次モーメントのラグランジュ
補間によって評価する．本研究では Frenklah[20]を参照に，負分数次モーメントは最初の3 次のモーメントで，正分数次モーメントは非負整数次モーメント全てを使って求めた．
さらに，𝑙𝑙  の grid function については𝑙𝑙  ∙∙∙からラグランジュ二次補間によって求めた．以上に述べたように，MOMIC では凝集によるソース項を算出するために，分数次モーメントおよび grid function の補間を 2 回実施している．本研究では，この 2 回の分数次モーメントおよび grid function のそれぞれをラグランジュ補間とガウス過程の結果で比較し，ラグランジュ補間による誤差の評価を実施した．なお，表面成長項の計算においても，分数次モーメントが用いられるが[16]，MOMGP では表面成長項も GP による分数次モーメントを使って算出している．
2-2．GPGP についてはビショップら[25]に詳説されており，ここでは概略を述べる．
ガウス過程は入力𝑥𝑥から目標変数𝑦𝑦の分布を予測する確率的生成モデルの一つである．入力
𝑥𝑥と出力𝑦𝑦の関係が観測ノイズを用いて𝑦𝑦�  𝑓𝑓𝑥𝑥�   𝜖𝜖� のように書けるとする．ガウス過程では入力変数列𝑋𝑋  𝑥𝑥 𝑥𝑥∙∙∙ 𝑥𝑥𝑥𝑥に対する出力𝑦𝑦の値𝑦𝑦 𝑦𝑦∙∙∙ 𝑦𝑦𝑦𝑦の同時分布が多変量ガウ
ス分布に従う，すなわち𝑝𝑝𝐲𝐲𝑁𝑁 𝐾𝐾と仮定する．ただし，K はカーネル行列であり，その成分はカーネル関数式𝐾𝐾���  𝑘𝑘𝑥𝑥�  𝑥𝑥�� で与えられる．カーネル関数としてガウスカー

ネルを用いると𝑘𝑘𝑥𝑥�  𝑥𝑥��   𝜃𝜃� − 𝑥𝑥� − 𝑥𝑥�� � 𝜃𝜃� (7)と書ける．ガウス過程の学習は訓練データに含まれない入力𝑥𝑥��� での目標変数𝑦𝑦��� を予測することを目的とする．ガウス過程では𝑌𝑌に𝑦𝑦��� を追加した𝑦𝑦 ∗  𝑦𝑦 𝑦𝑦  𝑦𝑦��� もガウス分布に従うので－ 249 －

�𝑦𝑦

𝐲𝐲

���

𝑲𝑲
� � � � �
𝑘𝑘∗
�

(8)

𝑲𝑲∗
��
𝑘𝑘∗∗

ただし，𝑘𝑘∗  �𝑘𝑘𝑥𝑥���  𝑥𝑥�    𝑘𝑘𝑥𝑥���  𝑥𝑥� � かつ𝑘𝑘∗∗  𝑘𝑘𝑥𝑥���  𝑥𝑥��� である．ゆえに，𝒚𝒚が与えられたときの𝑦𝑦��� の条件付確率は以下のように求まる．

(9)𝑝𝑝𝑦𝑦��� 𝐲𝐲  ��𝒌𝒌𝑻𝑻∗ 𝑲𝑲�� 𝐲𝐲 𝑘𝑘∗∗ − 𝒌𝒌𝑻𝑻∗ 𝑲𝑲�� 𝒌𝒌∗ �

なお，カーネル関数式(7)にはハイパーパラメータ𝜃𝜃� および𝜃𝜃� が含まれているが，これは尤度が最大になるように最適化を実施した．
2-3．計算条件

本研究では炭化水素燃料の熱分解の 0 次元計算を実施した．Table 1 には計算条件を示
す．計算結果は Alexiou らの実験結果[24]と比較検証した．彼らは衝撃波管実験を実施しており，ゆえに本研究では定温定圧での計算を行った．初期温度は 1500 から 2500K まで変化させ，初期圧力は各条件で計測された圧力の平均値を使用した．計算時間は 1.0，1.5
および 2.0msec とした．
Table 1 Calculation conditions
Mixture (mol%)Pressure [kPa]

Ar 99 %/ C6H5CH3 1 %258

Ar 99 %/ C6H5CH3 0.7 %, n-C7H16 0.3 %273

Ar 99 %/ C6H5CH3 0.7 %, n-C7H16 0.3 %276

Ar 99 %/ C6H5CH3 0.7 %, n-C7H16 0.3 %

288

3．結果
上述したようにカーネル関数（式(7)）はハイパーパラメータを 2 つ（𝜃𝜃� ，𝜃𝜃� ）含んでおり，これらは尤度が最大化するように最適化する必要がある．しかしながら，熱分解シミュレーションは非定常な現象であり，モーメントは時間の経過とともに変化する．そのため，最適化を時間刻みごとに実行すると，CPU 時間が増大してしまう．そのため，まず粒子成長計算を含まない，すなわちモーメント法を組み込まない計算を実施して各化学種の挙動を調査した．Fig. 1 にはトルエンおよび核生成に関わる PAH であるベンゾ[ghi]ペリ
レンおよびピレンのモル分率を示す．さらに，Fig. 1 にはピレンの変化率も示した．初期
条件はアルゴン 99%＋トルエン 1%（モル分率），初期温度 2000K である．Fig. 1 に示されるように，計算開始直後にトルエンは熱分解を開始し，同時にピレンの生成が始まる．
本研究で使用した反応モデルでは，C1-C4 分子からのベンゼンまでの分子成長に加えて，コロネンまでの PAH の生成および成長を含んでいる[26]．したがって，ピレンの消費が始
まると同時にベンゾ[ghi]ペリレンの生成が始まる．コロネンも生成され始めるが，その増加は緩やかである．しばらくすると，0.2ms 頃にベンゾ[ghi]ペリレンおよびピレンのモル－ 250 －

分率も減少し始めるが，コロネンは 2ms ま
でモル分率が増大し続けた．したがって，
全化学種のモル分率が 2ms 以内に一定になることはなかった．しかし，例えば，コロネンのモル分率の変化率は 0.5ms 以降，
10-4[1/sec]以下であり，それ以降も減少し続ける．上述のようなモル分率の変化は，他の初期条件においても観察された．ゆえに，0.5ms 以降のモーメント値もほぼ一定の値を示すものと予想される．そこで，本研究では 0.5ms におけるモーメントにガウス過程を適用した際の最適なハイパー
パラメータを用いて，MOMIC-GP による
すす生成計算を実施した．
次に分数次モーメントの評価におよぼす補間および回帰手法の影響を調査する．
Fig. 2 には Fig. 1 と同様の初期条件（アルFig. 1 Mole fractions of

ゴン 99%+トルエン 1%，初期温度 2000K）(a) C6H5CH3, (b) pyrene,

における分数次モーメントの値と，それらとラグランジュ補間および GP によって求めた分数次モーメントの値を示す．Fig. 2には凝集項および表面項の計算に必要な(c) BGHIPER, and (d) coronene
without considering soot formation（In(d), time derivative of mole fraction of

分数次モーメント（M-4/6,M-3/6,~M31/6）をプ

coronene was also plotted）

ロットしており，GP については期待値と共に 95%信頼区間も同時に示している．ここで，計算初期からラグランジュと GP を使った場合では 0.5ms ですでにモーメントの値が異なっている．ここでは，手法の比較のために，分数次モーメントは計算初期からラグランジュ補間を用いた，すなわち MOMIC を
0.5ms まで使用して求めた値としている．Fig. 2 に示されるように，オーダーの増大と共にモーメントの値はほぼ単調に増大した．ラグランジュ補間と GP で求められる分数次モーメントを比較すると，両者はよく一致しており，Fig. 2 ではそれほど差が見られなかった．さらに詳細に両者を比較するために，それらの差を算出した．Fig. 3 には GP で求められたモーメントをベースとしてラグランジュ補間と GP の差を示している．95%信頼区
間には収まらないものの，Fig. 2 と同様に，オーダー0 から 5 までにおいて差は小さく，
せいぜい 20%以内に収まった．一方で，外挿区間，すなわち，オーダー-4/6 から-1/6 および 31/6 においては両者の差は大きくなり，25%付近まで増大した．外挿区間における差の原因はラグランジュ補間の直線性にあると考えられる．Roy らはルンゲ現象を回避するた－ 251 －

めのラグランジュ補間に変わる新たな補間手法として Stineman 補間の適用を提案した．その精度を評価するために，彼らはMOMIC より正確な Discrete sectionalmodel で分数次モーメントを求め，それぞれの結果と比較した[22]．その結果，両者にはそれほど差が無いものの，両者とも
DSM よりも負オーダーのモーメントが小さくなったことを発見した．彼らはこの原因を Stineman 補間やラグランジュ補間の直線性にあると考察しており，実際の負オーダーにおけるモーメントの非線形性をFig. 2 Reduced whole moments

捉えることができなかったものと結論付(black) and fractional moments

けた．本研究では，GP はこの非線形性を

calculated by the Lagrange

より正確に捉えることができており，ゆえinterpolation (red) and GP (blue)に負オーダーにおいて偏差が負になったものと考えられる．
以下では上記のような GP によって分数次モーメントの非線形性を捉えることができたことがどのようにすす生成計算の精度に影響するかについて調査した．Fig.
4 にはラグランジュおよび GP を適用した
0 次から 5 次までの整数次モーメントを示す．初期条件は Fig. 1 と同様の条件であり，2msec における値を示している．ただし Fig. 2, Fig. 3 と異なり，GP は計算開始
時刻から適用されていることに留意して頂きたい．上述のように，分数次モーメン
トは凝集，PAH 凝固，表面反応の 3 段階で使用される．ここでは，どの段階で分数

Fig. 3 Difference between thefractional moments calculated by the

次モーメントを使用したことが結果に影

Lagrange interpolation and GP

響したのかを明らかにするため，凝集項の算出に GP を適用した場合を GPa，凝集項・PAH 凝集項に GP を適用した場合を GPb，
凝集項・PAH 凝集・表面反応全てに適用した場合を GPc，ラグランジュ補間を使用した場
合を Lag と名付けた．Fig. 4(b)に示されるように，1 次モーメントは全ての項に GP を適用しても影響は無かった．式(1)に示されるように，1 次モーメントは単位体積当たりに含－ 252 －

Fig. 4 Whole moments calculated by Lagrange interpolation and GP. GPa, GPb
and GPc correspond to the case where GP is applied to the aggregation,
aggregation and condensation, aggregation and condensation and surface reaction
term, respectively（Bars show the 95 % credible interval）
まれる炭素の個数を示している．ゆえに，凝集や PAH 凝集・表面反応が変化しても炭素の数密度は影響を受けないため，1 次モーメントは変化しなかったと考えられる．実験との検証に使われる特性の一つはすす転化率であり，燃料中の炭素がすす粒子に転化した割合である．これは 1 次モーメントから計算される．補間手法が 1 次モーメントにおよぼす影響は無かったため，すす転化率に対しても影響は無かった．ここで，本モデルの精度を検証するため，すす転化率を実験値(Alexiou)と比較する（Fig. 5）．Alexiou の論文はすす転化率を測定した正確な時刻が記述されていないため，1.0msec から 0.5msec ごとに－ 253 －

2.0msec のすす転化率を算出した．Fig. 5に示されるように，実験値も計算も1900K から 2000K 程度で最大値を示している．これはベル特性と呼ばれ[27]，本モデルが定性的にはすす生成挙動を再現できていることを示している．しかしながら，計算値は 1800K 以下では実験値を上回り，1800K 以上では下回ってしまっており，定量的には再現できたとは言えない．これは核生成反応や表面反応速度定数の不正確性に原因があると考えられるが，これらは今後調査する予定である．

Fig. 5 Soot yields obtained by MOMIC
and experimental measurement [24]

1 次モーメント以外の整数次モーメン
トについては補間手法により差が見られた．Fig. 4 に示されるように，0 次モーメントはGP を適用した方がラグランジュを適用した場合よりも小さくなり，2 から 5 次モーメン
トは大きくなった．加えて，Lag と GPa の差は GPa と GPb および GPb と GPc の差に比べて大きかった．このことから，補間手法の影響は PAH 凝集項や表面反応項に対してよりも凝集項に対して大きく影響をおよぼすことが示唆された．Lag に比べて 0 次モーメ
ントが小さくなった原因は 0 次の凝集項（式(4)）の補間に必要なグリッドファンクション
が 19/6 次以下の分数次モーメントの和積で構成されるためである．Fig. 3 に示されるように，GP は 0 次から 4 次程度では Lag とのモーメント差が小さいものの，負のオーダーでは大きく上回る．したがって，GP で算出される凝集項は Lag よりも大きくなる．凝集の活性はすす粒子の個数を減少させるため，0 次モーメントは Lag に比べて小さくなったと考えられる．同様にして，0 次から 4 次の分数次モーメントで構成されるグリッドファン
クションを含む 2 次 3 次の凝集項はそれほどモーメント値に影響せず，かつ，95%信頼区
間も小さかった（Fig. 4(c)(d)）．それに対して，31/6 次モーメントを含む 4 次 5 次は Lagに比べて GPa はかなり大きくなった（Fig. 4(e)(f)）．さらに，PAH 凝集によってすす粒子の個数は変化しないため，0 次の PAH 凝集項は 0 であるにもかかわらず，0 次モーメントでは GPa よりも GPb の方が小さくなった．これは 1 次の PAH 凝集項の計算に 0 次，2/6
次，4/6 次モーメントを含むためで，これらが Lag よりも大きい（Fig. 3）ために PAH 凝
集項は GPa よりも GPb の方が大きく評価される．このことはすなわち，PAH がより多く
PAH 凝集に消費され，核生成に使用される PAH の量が減少することを示唆しており，したがって GPa よりも GPb の 1 次モーメントは小さくなったものと考えられる．2 次以上の PAH 凝集項についても，同様に核生成に消費される PAH が減少したことが原因で GPa
よりも GPb のモーメントが小さくなったものと考えられる．一方で，GPb と GPc に大きな差は見られなかった．表面反応項でも分数次モーメントは使用されるものの，0 次から－ 254 －

5 次までの内挿区間の分数次モーメントであり，負や 31/6 次などの外挿区間のモーメントは使用されないことが原因と考えられる．加えて，表面反応項自体が凝集項や PAH 凝集
項に比べて 3 オーダーほど小さかったことも原因と考えられる．
すす転化率以外に実験値との検証に使われる指標としてはすす粒子径が挙げられる．すす粒子径の計算精度におよぼす補間および回帰手法の影響を調査した．モーメント法においては 2/6 次モーメントから計算される．Fig. 6 には Table 1 に示した条件下でのすす粒子径の算出結果を示す．Fig. 6 においてはどの条件下でもラグランジュ補間で求めた 2/6
次モーメントから算出した．単調ではないものの，2000K 付近で最大値を示す傾向がどの条件でも観察された．また，トルエンの初期モル分率が小さくなるのに伴い，粒子径が小さくなったことから，燃料の違いによるすす粒子形成の違いを定性的に本モデルは再現できたといえる．補間による影響については，どの条件でも GPb,GPc>GPa>Lag となった．
これは上述のように，GP は凝集項をラグランジュよりも大きく評価するためであると考えられる．以上のことから，補間手法はすす転化率の算出には影響しないものの，すす粒子径の評価には影響をおよぼすことが明らかとなった．

－ 255 －

Fig. 6 Soot particle diameter calculated by Lagrange interpolation, GPa, GPb
and GPc. (a) Ar 99 %/ C6H5CH3 1 %, (b) Ar 99 %/ C6H5CH3 0.7 %, n-C7H16 0.3 %,(c) Ar 99 %/ C6H5CH3 0.5 %, n-C7H16 0.5 %, and (d) Ar 99 %/ C6H5CH3 0.3 %, nC7H16 0.7 %4．結言
すす生成計算におよぼす補間手法の効果について明らかにするため，ラグランジュ補間の代わりに GP を適用したモーメント法（MOMIC-GP）を開発した．GP により分数次モーメントを計算し，算出した分数次モーメントは凝集項・PAH 凝集項および表面反応項の計算に利用された．MOMIC-GP で 0 次元のトルエン・ヘプタン混合気の熱分解計算を実施し，ラグランジュ補間を適用した MOMIC と比較した．初期圧力および温度はそれぞれ
258-288kPa，1500-2500K である．さらに，計算結果を衝撃波管実験の結果と比較検証した．化学種の濃度変化が 0.5ms でほぼ定常的になったため，ハイパーパラメータは 0.5ms
時に最適化し，同じ値のハイパーパラメータを全計算時間中に利用した．0-5 次の分数次モーメントはラグランジュと GP とでは大きな差は見られなかったものの，負の分数次モーメントおよび 31/5 次モーメントでは 25%程度の差が見られた．これは GP が分数次モーメントの非線形性を再現できたことを示唆している．また，負の分数次モーメントは GP－ 256 －

の方が大きく評価するため，0 次の凝集項を大きく評価した．このことから 0 次モーメン
トは小さく評価し，一方ですす粒子径は大きく評価することになった．また，凝集項が 0となる 1 次モーメントについては MOMIC-GP と MOMIC ではほとんど差が見られなかった．
謝辞
本研究は公益財団法人 JFE21 世紀財団の技術研究助成により実施されたことをここに記し，深謝いたします．
参考文献
[1]

J.B. Howard, Carbon addition and oxidation reactions in heterogeneous
combustion and soot formation, Symp. Combust. 23 (1991) 1107–1127.

[2]

H. Richter, J.B. Howard, Formation of polycyclic aromatic hydrocarbons andtheir growth to soot-a review of chemical reaction pathways, 2000.

[3]

A. D’Anna, Combustion-formed nanoparticles, Proc. Combust. Inst. 32 I (2009)593–613.

[4]

M. von. Smoluchowski, Drei vortrage uber diffusion, brownsche bewegung und
koagulation von kolloidteilchen, Zeitschrift Fur Phys. 17 (1916) 557–585.

[5]

M. V. Smoluchowski, An experiment on mathematical theorization of coagulation
kinetics of the colloidal solutions, Zeitschrift Fur Phys. Chemie. 92 (1917) 129–
168.

[6]

S. Chandrasekhar, Stochastic problems in physics and astronomy, Rev. Mod.
Phys. 15 (1943) 1.

[7]

G. M. Hidy, J. R. Brock, Topics in Current Aerosol Research, 1st ed., 1971.

[8]

F.A. Williams, Spray combustion and atomization, Phys. Fluids. 1 (1958) 541–
545.

[9]

D.E. Rosner, R. McGraw, P. Tandon, Multivariate population balances via
moment and Monte Carlo simulation methods: An important sol reactionengineering bivariate example and “mixed” moments for the estimation ofdeposition, scavenging, and optical properties for populations of nonsphericalsuspended particles, Ind. Eng. Chem. Res. 42 (2003) 2699–2711.

[10]

M. Balthasar, M. Frenklach, Monte-Carlo simulation of soot particle coagulation
and aggregation : the effect of a realistic size distribution, Proc. Combust. Inst.
30 (2005) 1467–1475.

[11]

Y. Wang, A. Raj, S.H. Chung, Soot modeling of counterflow diffusion flames of
ethylene-based binary mixture fuels, Combust. Flame. 162 (2015) 586–596.

－ 257 －

[12]J. Appel, H. Bockhorn, M. Frenklach, c, Combust. Flame. 121 (2000) 122–136.

[13]

H. Wang, D.X. Du, C.J. Sung, C.K. Law, EXPERIMENTS AND NUMERICALSIMULATION ON SOOT FORMATION IN OPPOSED-JET ETHYLENEDIFFUSION FLAMES, 1996.

[14]

B. Zhong, S. Dang, Y. Song, J. Gong, 3-D simulation of soot formation in a directinjection diesel engine based on a comprehensive chemical mechanism andmethod of moments, 7830 (2012).

[15]

M. Frenklach, S. Harris, J, Aerosol Dynamics Modeling Using the Method ofMoments, J. Colloid Interface Sci. 118 (1987) 252–261.

[16]

M. Frenklach, Method of moments with interpolative closure, Chem. Eng. Sci. 57(2002) 2229–2239.

[17]

R. McGraw, Description of aerosol dynamics by the quadrature method of
moments, Aerosol Sci. Technol. 27 (1997) 255–265.

[18]

D.L. Marchisio, R.O. Fox, Solution of population balance equations using thedirect quadrature method of moments, J. Aerosol Sci. 36 (2005) 43–73.

[19]

A. Wick, M. Frenklach, H. Pitsch, Systematic assessment of the Method ofMoments with Interpolative Closure and guidelines for its application to sootparticle dynamics in laminar and turbulent flames, Combust. Flame. 214 (2020)
450–463.

[20]

M. Frenklach, H. Wang., Detailed mechanism and modeling of soot particle
formation, in: Soot Form. Combust. Springer, Springer, Berlin, Heidelberg, 1994:
pp. 165–192.

[21]

A. Kazakov, M. Frenklach, Dynamic Modeling of Soot Particle Coagulation andAggregation: Implementation With the Method of Moments and Application to
High-Pressure Laminar Premixed Flames, Combust. Flame. 114 (1998) 484–501.

[22]

S.P. Roy, P.G. Arias, V.R. Lecoustre, D.C. Haworth, H.G. Im, A. Trouvé,
Development of high fidelity soot aerosol dynamics models using method of
moments with interpolative closure, Aerosol Sci. Technol. 48 (2014) 379–391.

[23]

S. Nakaya, Y. Yamaki, F. Ren, M. Tsue, Effects of the discharge frequency on the
dielectric barrier discharge ignition behaviors for lean methane–air mixtures atvarious pressure values, Combust. Flame. 234 (2021) 111650.

[24]

A. Alexiou, A. Williams, Soot formation in shock-tube pyrolysis of toluene-nheptane and toluene-iso-octane mixtures, Fuel. 74 (1995) 153–158.

[25]

Cristopher M. Bishop, Pattern Recognition and Machine Learning, Springer,
Berlin, Heidelberg, 2006.

[26]

Y. Wang, A. Raj, S.H. Chung, A PAH growth mechanism and synergistic effect－ 258 －

on PAH formation in counterflow diffusion flames, Combust. Flame. 160 (2013)
1667–1676.
[27]

H. Yui, Y. Namai, K. Iwata, O. Imamura, J. Hashimoto, K. Akihama, すす粒子生成モデルにおける表面反応および核形成反応がすす生成特性に及ぼす影響, 自動車技術会論文集. 48 (2017) 1207–1212.

－ 259 －