JP6311835B2

JP6311835B2 - 反応機構生成方法及び反応機構生成装置

Info

Publication number: JP6311835B2
Application number: JP2017500640A
Authority: JP
Inventors: 宏俊平井
Original assignee: Toyota Central R&D Labs Inc
Current assignee: Toyota Central R&D Labs Inc
Priority date: 2015-02-19
Filing date: 2016-02-10
Publication date: 2018-04-18
Anticipated expiration: 2036-02-10
Also published as: JPWO2016133002A1; WO2016133002A1; US20180032704A1

Description

本発明は、反応機構生成方法及び反応機構生成装置の技術に関する。

製造工程において、生産性などを改良するためにプロセスシュミレーションが検討されることが多くなってきている。

例えば、化学反応を利用した製造過程においては、その反応がどのようにして起こり、どのようにすれば収率や反応速度を改良できるかが検討される。このようなプロセスシュミレーションを行うためには、反応メカニズムの詳細、反応速度の算出などが必要となる。

例えば、特許文献１及び非特許文献１には、分子動力学計算を用いて、原子間の結合情報を算出し、その原子間の結合情報の変化から、反応メカニズムを推定する方法が提案されている。

特許第４７１３８５６号公報

ＪｏｕｒｎａｌｏｆＭｏｌｅｃｕｌａｒＧｒａｐｈｉｃｓａｎｄＭｏｄｅｌｌｉｎｇ，５３（２０１４），１３−２２．

ところで、分子数の多い複雑系の反応の場合、例えばＡ＋Ｂ→Ｄ＋ＥとＣ→Ｆ＋Ｇ等の複数の素反応が同時に起こる場合があるが、特許文献１及び非特許文献１の方法では、このような素反応の解析が考慮されていない。したがって、特許文献１及び非特許文献１の方法では、分子数の多い複雑系の反応を解析する場合、見かけ上の多分子反応、例えばＡ＋Ｂ＋Ｃ→Ｄ＋Ｅ＋Ｆ＋Ｇを抽出してしまい、反応を精度よく解析することができない場合がある。

そこで、本発明の目的は、分子数の多い複雑系の反応を解析する場合でも、反応を精度よく解析することができる反応機構生成方法及び反応機構生成装置を提供することである。

本発明の反応機構生成方法は、反応系内の各分子を構成する原子について、時間ステップ毎に分子動力学計算を行うステップと、前記時間ステップ前後で前記反応系内に化学反応が起こった場合、前記化学反応に寄与した反応分子及び生成分子を特定するステップと、前記反応分子と前記生成分子との原子の種類の関連性に基づいて、前記関連性のある反応分子及び生成分子から構成された素反応を構築するステップと、前記構築した素反応の反応速度定数を算出するステップと、を含む。

また、前記反応機構生成方法において、前記反応分子及び前記生成分子を特定するステップでは、前記時間ステップ前後の分子動力学計算から得られる原子の位置情報から、前記化学反応に寄与した前記反応分子及び前記生成分子を特定することが好ましい。

また、前記反応機構生成方法において、前記素反応の反応速度定数を算出するステップでは、前記反応分子と同じ分子種の分子の個数の時間履歴を算出し、前記反応分子と同じ分子種の分子の個数の時間履歴に基づいて、前記素反応の反応速度定数を算出することが好ましい。

また、前記反応機構生成方法において、前記素反応を構築するステップにおいて構築された複数の素反応を含むベースモデルを作成するステップと、前記分子動力学計算中に起こった素反応の回数に基づいて、前記ベースモデルの前記複数の素反応を絞り込み、前記ベースモデルより素反応の数が少ない簡略化モデルを作成するステップと、を含むことが好ましい。

また、前記反応機構生成方法において、前記簡略化モデルの精度を評価するステップを含むことが好ましい。

また、前記反応機構生成方法において、前記素反応を構築するステップにおいて構築された複数の素反応を含むベースモデルを作成するステップと、前記分子動力学計算中に起こった素反応の回数、及び前記分子動力学計算において設定された複数の設定温度のうち、同一の素反応が存在する設定温度の個数に基づいて、前記ベースモデルの前記複数の素反応を絞り込み、前記ベースモデルより素反応の数が少ない簡略化モデルを作成するステップと、を含むことが好ましい。

また、本発明の反応機構生成装置は、反応系内の各分子を構成する原子について、時間ステップ毎に分子動力学計算を行う力学計算演算部と、前記時間ステップ前後で前記反応系内に化学反応が起こった場合、前記化学反応に寄与した反応分子及び生成分子を特定する分子特定部と、前記反応分子と前記生成分子との原子の種類の関連性に基づいて、前記関連性のある反応分子及び生成分子から構成された素反応を構築する素反応構築部と、前記構築した素反応の反応速度定数を算出する反応速度定数算出部と、を含む。

また、前記反応機構生成装置において、前記分子特定部は、前記時間ステップ前後の分子動力学計算から得られる原子の位置情報から、前記化学反応に寄与した反応分子及び生成分子を特定することが好ましい。

また、前記反応機構生成装置において、前記反応速度定数算出部は、前記反応分子と同じ分子種の分子の個数の時間履歴を算出し、前記反応分子と同じ分子種の分子の個数の時間履歴に基づいて、前記素反応の反応速度を算出することが好ましい。

また、前記反応機構生成装置において、前記素反応構築部において構築された複数の素反応を含むベースモデルを作成するベースモデル作成部と、前記分子動力学計算中に起こった素反応の回数に基づいて、前記ベースモデルの前記複数の素反応を絞り込み、前記ベースモデルより素反応の数が少ない簡略化モデルを作成する簡略化モデル作成部と、を含むことが好ましい。

また、前記反応機構生成装置において、前記簡略化モデルの精度を評価するモデル精度評価部を含むことが好ましい。

また、前記反応機構生成装置において、前記素反応構築部において構築された複数の素反応を含むベースモデルを作成するベースモデル作成部と、前記分子動力学計算中に起こった素反応の回数、及び前記分子動力学計算において設定された複数の設定温度のうち、同一の素反応が存在する設定温度の個数に基づいて、前記ベースモデルの前記複数の素反応を絞り込み、前記ベースモデルより素反応の数が少ない簡略化モデルを作成する簡略化モデル作成部と、を含むことが好ましい。

本発明によれば、分子数の多い複雑系の反応を解析する場合でも、反応を精度よく解析することができる。

本実施形態に係る反応機構生成装置の構成の一例を示すブロック図である。本実施形態の反応機構生成方法の一例を示すフロー図である。各原子間の結合有無を行列情報として表した図である。各原子間の結合有無を示す行列情報をブロック対角行列情報に変換した図である。特定した反応分子と生成分子の原子の関連性を説明する図である。反応分子と生成分子の原子の関連性の有無を行列情報として表した図である。特定の素反応が起こった時間に対する分子種Ａ及びＢの数の変化を示した図である。反応速度と温度との関係を示す図及びアレニウスの式を示す図である。本実施形態に係る反応機構生成装置の構成の一例を示すブロック図である。本実施形態の反応機構生成装置の動作の一例を示すフロー図である。水素分子と酸素分子をシミュレーションセル中に配置した状態を示す図である。実施例１において特定した分子数の時間変化を示す図である。実施例１で構築された素反応の一覧を示す図である。Ｈ_２＋Ｏ→ＯＨ＋Ｈの素反応における反応速度定数の計算結果を示す図である。３０００Ｋにおける素反応、反応速度定数、及び分子動力学計算中に起こった素反応の回数を含むベースモデルの一部を示す図である。素反応、パラメータ（Ａ、ｎ、Ｅａ）、素反応が起こった設定温度の個数を含む簡略化モデルの一部を示す図である。３_４簡略化モデルを用いた３０００Ｋにおけるモル分率の時間変化及び分子動力学計算により得られるモル分率の時間変化を示す図である。各簡略化モデル中の素反応の数及び化学種の数を示す図である。３_２簡略化モデル〜２０_５簡略化モデルの着火遅れ時間の計算結果を示す図である。１０_５簡略化モデルを用いた３０００Ｋにおけるモル分率の時間変化を示す図である。従来知られている３つの簡略化法を用いて得られた簡略化モデルにおける素反応数と化学種数の結果を示す図である。

以下、本発明の実施形態の一例について、図面に基づいて説明する。

図１は、本実施形態に係る反応機構生成装置の構成の一例を示すブロック図である。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろな形で実現できることは、本明細書に触れた当業者には理解されるところである。

図１に示すように、反応機構生成装置１００は、入力装置１０２および出力装置１０４と接続される。入力装置１０２は、反応機構生成装置１００で実行される処理に関係するユーザの入力を受けるためのキーボード、マウスなどであってもよい。入力装置１０２は、インターネットなどのネットワークやＣＤ、ＤＶＤなどの記録媒体から入力を受けるよう構成されていてもよい。出力装置１０４は、ディスプレイなどの表示機器やプリンタなどの印刷機器であってもよい。

反応機構生成装置１００は、反応条件設定部１０６と、初期条件設定部１０８と、分子動力学演算部１１０と、分子特定部１２０と、素反応構築部１２２と、反応速度定数算出部１１２と、データ保持部１１４と、表示制御部１１６とを備える。

反応条件設定部１０６は、入力装置１０２を介してユーザから入力された入力情報に基づいて、反応系内の初期分子種、各初期分子種の個数、反応系の温度、圧力、体積等の反応条件を設定する。

初期条件設定部１０８は、入力装置１０２を介してユーザから入力された入力情報に基づいて、分子を構成する原子の位置情報、原子の速度等の分子動力学計算に必要な初期条件を設定する。初期条件設定部１０８は、分子の個数が複数個設定されている場合には、例えば、分子の構造を壊さない範囲で、各原子の位置情報及び速度をランダムに設定する。

分子動力学演算部１１０は、設定された反応条件及び初期条件の下、ニュートンの運動方程式の時間積分により逐次的に計算する分子動力学計算を行い、時間ステップ毎に各原子の位置情報及び速度等を算出する。

分子特定部１２０は、分子動力学計算の時間ステップ前後で、反応系内に化学反応が起こった場合に、時間ステップ前後の分子動力学計算から得られる原子の位置情報等から化学反応に寄与した反応分子及び生成分子を特定する。反応系内に化学反応が起こったか否かは、後述するように、原子間の結合距離、結合次数、分子の数の変化等により判定する。

素反応構築部１２２は、特定した反応分子及び生成分子の原子の関連性に基づいて、原子の関連性のある反応分子及び生成分子から構成される素反応を構築する。反応分子及び生成分子の原子の関連性は、後述するように、反応分子を構成する原子と生成分子を構成する原子とを比較し、反応分子を構成する原子が生成分子に含まれているか否か等により決定する。

反応速度定数算出部１１２は、反応速度式を用いて、特定した素反応の反応速度定数を算出する。

データ保持部１１４は、素反応構築部１２２により構築された素反応及び反応速度定数算出部１１２により算出された反応速度定数をデータとして保持する。表示制御部１１６は、データ保持部１１４に保持された素反応及び反応速度定数を出力装置１０４に表示させる。

以下に、本実施形態に係る反応機構生成装置１００の動作を図２に示すフロー図に従って説明する。

ステップＳ１０では、反応条件設定部１０６により、反応を解析しようとする反応系の諸条件が設定される。ここで、反応系の諸条件としては、初期分子種、各初期分子種の個数、反応系の温度、圧力、体積等である。初期分子種としては、反応系の原料となる反応分子等であり、例えば、複数の原子から構成される分子、単原子分子、イオン、ラジカル等が挙げられる。各初期分子種の個数は特に制限されるものではないが、様々な反応経路の出現を容易にするため、複数個設定することが望ましい。

ステップＳ１２では、初期条件設定部１０８により、分子動力学計算で用いられる初期条件が設定される。初期条件は、反応系内の分子を構成する原子の位置及び速度等である。前述したように、分子の個数を複数個設定した場合には、分子の構造を壊さない範囲で、各原子の位置及び速度をランダムに設定する。また、複数の分子について分子動力学計算を行う場合には、各分子を適当な間隔をおいて配置させ、ランダムな方向に配向させる。

ステップＳ１４では、分子動力学演算部１１０により、反応系内の分子を構成する原子について、時間ステップ毎に分子動力学計算が行われる。分子動力学計算は、反応系の諸条件及び分子動力学計算で用いられる初期条件の下、下式（１）で表されるニュートンの運動方程式に従って計算する手法である。分子動力学演算部１１０では、時間ステップ（Δｔ）を定め、Δｔ後の原子の位置情報を、ニュートンの運動方程式の時間積分により逐次的に計算している。

ｍ_ｉ：ｉ番目の原子の質量
Ｒ_ｉ：ｉ番目の原子の位置（Ｘ_ｉ、Ｙ_ｉ、Ｚ_ｉ）
Ｆ_ｉ：ｉ番目の原子に働く力（Ｆ_Ｘｉ、Ｆ_Ｙｉ、Ｆ_Ｚｉ）

原子に働く力（Ｆ）は、例えば、第一原理計算等の非経験的手法、分子力場法等の経験的手法等により求められる。原子に働く力の計算精度は、分子動力学計算による原子の位置及び速度の計算結果に大きく影響するため、反応系で起こりうる化学反応を精度よく記述できる手法が好ましい。

このような分子動力学計算により、各原子の位置（Ｒ）が、予め設定した時間ステップ毎に求められる。以下、各原子の位置を位置情報と称する。また、時間ステップ毎の原子の位置情報から、時間ステップにおける原子の移動距離が分かるため、この移動距離から原子の速度が求められる。あるいは、速度に関する時間差分方程式を計算することで原子の速度が求められる。

ステップＳ１６では、分子特定部１２０により、分子動力学計算の時間ステップ前後（ｔ→ｔ＋Δｔ）で、反応系内に化学反応が起こったか否かが判断される。化学反応が起こったか否かは、例えば、以下の方法により判定される。（Ａ）各原子間の結合距離あるいは結合次数の大きさから、各原子間の結合の有無を評価し、各原子間の結合の有無が、分子動力学計算の時間ステップの前後で異なる場合に化学反応が起こったと判定する方法、（Ｂ）各原子の位置情報から各分子の分子種を特定し、何れかの分子種の分子の個数が分子動力学計算の時間ステップの前後で異なる場合に化学反応が起こったと判定する方法等である。以下、（Ａ）及び（Ｂ）の手法について具体的に説明する。

＜（Ａ）の方法＞
図３は、各原子間の結合有無を行列情報として表した図である。例えば、反応系内に配置された分子を構成する各原子に原子ＩＤ（１、２、３、４、５・・・）を付した反応系において、ある時間ｔにおける各原子の位置情報（Ｒ）の差を求める。各原子の位置情報の差の絶対値が、各原子間の結合距離となる。そして、図３の左側に示すように、各原子間の結合距離が閾値（Ｃｕｔｏｆｆ）を超えた場合には結合無しと評価し、各原子間の結合距離が閾値（Ｃｕｔｏｆｆ）より小さい場合には、結合有と評価する。そして、各原子間で結合有りと評価した場合を１とし、結合無しと評価した場合を０とし、図３の右側に示すような行列情報としてまとめる。また、時間ステップΔｔ後（時間ｔ＋Δｔ）における各原子の位置情報から、上記と同様にして結合有無を評価し、行列情報としてまとめる。そして、時間ｔにおける行列情報と時間ステップΔｔ後（時間ｔ＋Δｔ）の行列情報を比較し、各原子間の結合の有無が異なる場合に、化学反応が起こったと判断する。

また、第一原理計算手法や一部の分子力場法では、各原子間の結合次数を求めることができるが、各原子間の結合次数を用いる場合は、例えば、各原子間の結合次数が閾値を超える場合、結合あり「１」と評価し、閾値より小さい場合には、結合無し「０」と評価して行列情報を作成する。そして、前述と同様に、時間ｔにおける行列情報と時間ステップΔｔ後の行列情報を比較し、各原子間の結合の有無が異なる場合に、化学反応が起こったと判断する。

＜（Ｂ）の方法＞
後に示すように原子の位置情報から、反応系内の各分子の分子種を特定することができるが、これにより各分子種の分子の個数を各時間ステップにおいて求めることができる。時間ステップの前後を比較し、個数が変化している分子種があった場合に、化学反応が起こったと判断する。

分子動力学計算の時間ステップ前後において、化学反応が起こっていないと判断した場合にはステップＳ１８に進み、化学反応が起こったと判断した場合はステップＳ２０に進む。なお、化学反応が起こったか否かの判断は、分子動力学計算の時間ステップ毎に実行してもよいが、一般的に、化学反応が起こるのに要する時間よりも、分子動力学計算の時間ステップの方が短いため、複数の時間ステップ毎実行してもよい。

ステップＳ１８では、分子動力学計算の時間ステップ数が最大であるか否かを判断し、最大ステップ数未満であれば、再度、ステップＳ１４の分子動力学計算に戻り、時間ステップ数が最大ステップ数に達していれば終了する。なお、分子動力学計算が最大ステップ数であるか否かの判断は、分子特定部１２０により行ってもよいし、別途ステップ数判定部を設置して、当該ステップ数判定部により行ってもよい。また、最大ステップ数は、反応系内の分子の数等によって適宜設定されればよい。

ステップＳ２０では、分子特定部１２０により、分子動力学計算の時間ステップ前後に起こった化学反応に寄与した反応分子及び生成分子を特定する。分子を特定する方法としては、分子動力学計算により得られる原子の位置情報又は原子間の結合の有無を示す行列情報（結合情報）等から特定される。以下具体的に説明する。

ある時間ｔにおける原子の位置情報から、分子が特定され、同様に、時間ステップΔｔ後における原子の位置情報から、分子が特定される。そして、ある時間ｔで特定された分子のうち、時間ステップΔｔ後で特定された分子に含まれていない分子（例えば、原子ＩＤ１及び２の原子から構成される分子１、原子ＩＤ１５２及び１５３から構成される分子２、原子ＩＤ１２、１４０及び１４１から構成される分子３）が、化学反応に寄与した反応分子として特定される。また、時間ステップΔｔ後で特定された分子のうち、時間ｔで特定された分子に含まれていない分子（例えば、原子ＩＤ２から構成される分子４、原子ＩＤ１、１５２及び１５３から構成される分子５、原子ＩＤ１２から構成される分子６、原子ＩＤ１４０及び１４１から構成される分子７）が、化学反応に寄与した生成分子として特定される。

図４は、各原子間の結合有無を示す行列情報をブロック対角行列情報に変換した図である。各原子間の結合の有無を示す行列情報を用いて化学反応が起こったか否か判断した場合には、例えば、図４に示す行列情報Ａ１から、結合有とされた原子を並べ替えて、各分子に分割したブロック対角行列情報Ａ２を作成することが望ましい。このようなブロック対角行列情報を時間ｔにおける行列情報及び時間ステップΔｔ後における行列情報それぞれから作成し、時間ｔのブロック対角行列情報の分子のうち、時間ステップΔｔのブロック対角行列情報にない分子は、化学反応に寄与した反応分子として特定される。また、時間ステップΔｔ後のブロック対角行列情報の分子のうち、時間ｔのブロック対角行列情報の分子にない分子は、化学反応に寄与した生成分子として特定される。

このように特定された反応分子を左側に、生成分子を右側に記述することで反応式を構築することができる。しかし、分子数の多い複雑系の反応の場合、例えば、分子１＋分子２→分子４＋分子５と分子３→分子６＋分子７等の反応が同時に起こっているにも関わらず、ステップＳ２０から反応式を構築しようとすると、例えば、分子１＋分子２＋分子３→分子４＋分子５＋分子６＋分子７の複数の素反応が混合した反応式として抽出されてしまい、反応解析の精度が充分に保障されない場合がある。そのため、反応機構を精度よく解析するためには、真の素反応を構築する必要がある。

ステップＳ２２では、素反応構築部１２２により、特定した反応分子と生成分子の原子の関連性から素反応が構築される。図５は特定した反応分子と生成分子の原子の関連性を説明する図である。図５に示すように、反応分子として分子１、２、３が特定され、生成分子として分子４、５、６、７が特定された場合、各反応分子を構成する原子と、各生成分子を構成する原子とを比較し、反応分子に含まれる原子が生成分子の原子に含まれる場合、その反応分子と生成分子には関連性があると判断する。図５で示すように、例えば分子１の原子は分子４又は５の原子に含まれているが、分子６，７の原子には含まれていないため、分子１は分子４，５と関連性があるが、分子６，７とは関連性がないと判断することができる。分子２は、分子１と関連性がある分子５と関連性がある。これらの原子の関連性により、分子１，２，４，５から素反応が構成される。また、分子３は、分子６及び７と関連性があるが、分子４，５とは関連性がないため、分子３，６，７から素反応が構成されることになる。そして、反応分子を左側、生成分子を右側に記述することで、分子１＋分子２→分子４＋分子５と分子３→分子６＋分子７の素反応が構築される。抽出した素反応データは、反応速度定数算出部１１２に送られると共に、データ保持部１１４に格納される。

図６は、反応分子と生成分子の原子の関連性の有無を行列情報として表した図である。図６に示すように、反応分子と生成分子の原子の関連性において、関連性がある場合を１、関連性が無い場合を０として、関連性のある分子をまとめてブロック対角化した行列情報を作成してもよい。図６の各ブロックＣ，Ｄが真の素反応に対応している。そして、図６に示すブロック毎に、反応分子を左側、生成分子を右側に記述することで、分子１＋分子２→分子４＋分子５と分子３→分子６＋分子７の素反応を構築することができる。

ステップＳ２４では、反応速度定数算出部１１２により、構築した素反応の反応速度定数が算出される。ある分子種Ａの反応速度式は、各素反応の反応速度定数（ｋ）を用いて、下式（２）の上段として表される。さらに、この式は、式（２）の下段のように、同時に進行する様々な素反応による寄与の和として表される。

また、各素反応の寄与について、反応速度式を差分形式で表すと、下式（３）の一次反応は、下式（４）で表され、下式（５）の二次反応は、下式（６）で表される。

ここで、Ｎは反応系内の各分子種の分子の数であり、Ｖは反応系の体積であり、δｔは特定の素反応が１回起こるのに要した時間である。そして、特定の素反応が１回起きた場合の寄与を考えるので、分子種Ａが反応分子である場合には、式（４）のＮ_Ａ（ｔ＋δｔ）−Ｎ_Ａ（ｔ）は−１となり、式（４）の下段の式に帰着する。また、二次反応の場合で、反応分子の分子種ＡとＢが同一の分子である場合、Ｎ_Ａ（ｔ＋δｔ）−Ｎ_Ａ（ｔ）は−２となるが、その場合には、反応式の左辺に１／２の因子が定義上現れるため相殺されて、式（６）の下段の式に帰着する。すなわち、他の素反応の寄与がない場合には、式（２）から直接反応速度定数ｋを求めることができる。

しかし、他の素反応の寄与がある場合には、δｔの間にも他の素反応によって、各分子の数は変化し得るため、式（２）から、反応速度定数ｋを求めるより、以下で示す式（７）、（８）から反応速度定数ｋを求める方が、精度の高い計算結果が得られる。

図７は、特定の素反応が起こった時間に対する分子種Ａ及びＢの数の変化を示した図である。式７及び式８は、式４及び式６の分母（Ｎ_Ａδｔ、Ｎ_ＡＮ_Ｂδｔ）を、図７に示すδｔの時間区間（ｔ_{ｓｔａｒｔ}からｔ_ｅｎｄまで）における各反応分子の数の時間積分値、すなわち反応分子の数の時間履歴に置き換えたものである。一般的に各素反応は分子動力学計算中に複数回起こりうるため、反応速度定数の計算を素反応が起こる度に実行し、複数の得られた値に対して統計平均を取ることで、反応速度定数の計算精度をさらに向上させることができる。また、素反応が起こる度に、δｔの時間区間で反応分子の個数の時間積分値を計算し、反応速度定数を計算してもよいが、反応分子の個数の時間積分値を全時間区間で求めておき、素反応が起こった回数で割ることで、素反応１回当たりの反応速度定数の平均値を求めてもよい。算出した反応速度定数データをデータ保持部１１４に格納する。

反応速度定数の計算後、ステップＳ１８に戻り、分子動力学計算の時間ステップ数が最大であるか否かを判断し、最大ステップ数未満であれば、再度、ステップＳ１４の分子動力学計算に戻り、時間ステップ数が最大ステップ数に達していれば終了する。

本実施形態では、分子動力学計算中に逐次素反応を抽出する場合を例示したが、これに制限されるものではなく、分子動力学計算を設定した最大ステップ数まで実行し、各原子の位置情報を格納しておき、その後の処理として、各原子の位置情報を利用して素反応の抽出を行っても良い。

一般に反応速度定数は温度に依存するため、熱浴を用いた温度一定の分子動力学計算を行うことで、定めた温度における反応速度定数を求めることができる。図８は反応速度と温度との関係を示す図及びアレニウスの式を示す図である。図８に示すように、複数の温度での反応速度定数を求め、求めた反応速度定数をアレニウスの式にフィッティングさせることで、反応速度パラメータを算出することが可能となる。

以上のように、本実施形態では、反応系が複数の反応分子を含む場合でも、反応分子と生成分子との原子の関連性から素反応を構築するステップを含むため、見かけ上の多体反応を抽出してしまうリスクを減らして、真の素反応を構築することができるため、反応を精度良く解析することができる。

反応速度定数の計算方法としては、遷移状態理論等の他の計算手法を用いてもよい。しかし、遷移状態探索等の計算コストの大きい演算処理を必要とせず、容易に反応速度定数を計算することができる等の点から、前述したように、反応分子の数の時間履歴を用いて、分子動力学計算の結果から直接反応速度定数を求めることが好ましい。

このようにして得られた素反応及び反応速度定数は、反応プロセスのシミュレーションに用いることができる。また、流体力学計算と連成させることで、流れや拡散の影響を取り入れることも可能となる。また、実験データや経験則を用いることなく、反応機構（一連の素反応及び反応速度定数の集合）を構築することができるため、反応機構が未知の系に対してもシミュレーションの実行が可能となる。

また、このようにして得られた素反応及び反応速度定数は、化学過程と物理過程をモデル化した三次元数値流体力学シミュレーション等に利用することも可能である。

しかし、三次元数値流体力学シミュレーションを用いて、例えばエンジン燃焼シミュレーションを実行しようとすると、燃料の化学反応のモデル化がネックとなる。例えば、炭化水素燃料の燃焼のような複雑な反応をモデル化すると、膨大な数の化学種と素反応を含んだものとなるため、化学反応を考慮した三次元数値流体力学シミュレーションを実行する際には、モデルに含まれる素反応の数に比例して解くべき方程式の数が増加する問題がある。したがって、三次元数値流体力学シミュレーションの計算時間を短縮する等のために、反応機構モデルの簡略化が求められている。

そこで、反応機構モデルの簡略化方法について、以下説明する。

図９は、本実施形態に係る反応機構生成装置の構成の一例を示すブロック図である。図９に示す反応機構生成装置１０１は、反応機構モデル形成部１２４、ベースモデル作成部１２６、簡略化モデル作成部１２８、モデル精度評価部１３０、データ保持部１１４、表示制御部１１６を有している。ここに示す各ブロックは、ハードウエア的には、コンピュータのＣＰＵ（central processing unit）をはじめとする素子や機械装置で実現でき、ソフトウエア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウエア、ソフトウエアの組合せによっていろいろな形で実現できることは、本明細書に触れた当業者には理解されるところである。

図９に示す反応機構モデル形成部１２４は、図示しないが、図１に示す反応条件設定部１０６と、初期条件設定部１０８と、分子動力学演算部１１０と、分子特定部１２０と、素反応構築部１２２と、反応速度定数算出部１１２と、を備えており、複数の設定温度での素反応の構築及び反応速度定数の算出を実行する。素反応の構築及び反応速度定数の算出については前述した通りであるので、その説明を省略する。

図９に示すベースモデル作成部１２６は、各設定温度での各素反応、各素反応の反応速度定数、及び各素反応について分子動力学計算中に起こった素反応の回数を含むベースモデルを作成する。

図９に示す簡略化モデル作成部１２８は、ベースモデル作成部１２６により作成されたベースモデルに存在する素反応の絞り込みを行い、ベースモデルより素反応の数が少ない簡略化モデルを作成する。簡略化モデルの作成については、以下で詳述する。

図９に示すモデル精度評価部１３０は、簡略化モデル作成部１２８により作成された簡略化モデルの精度を評価する。モデル精度の評価は、例えば、簡略化モデルの化学特性とベースモデルの化学特性とを比較して、その差が許容範囲内にあるかを確認することにより行われる。そして、許容範囲内にある場合には、その簡略化モデルを簡略化した反応機構モデルとして、データ保持部１１４に送信する。

以下に、本実施形態の反応機構生成装置１０１の動作の一例を図１０に示すフロー図に従って説明する。

ステップＳ５０では、ベースモデル作成部１２６により、例えば、設定温度Ｔ_１〜Ｔ_Ｎ（例えば、３０００Ｋ、３２５０Ｋ・・・）毎に、各素反応（素反応Ａ、素反応Ｂ、素反応Ｃ・・・）、各素反応の反応速度定数（ｋ１、ｋ２、ｋ３・・・）、分子動力学計算中に起こった各素反応の回数（ａ１、ａ２、ａ３・・・）を含む素反応リストＬ_１〜Ｌ_Ｎとして表されるベースモデルが作成される。

各素反応、各素反応の反応速度定数は、反応機構モデル形成部１２４により構築された素反応及び計算された反応速度定数である。また、分子動力学計算中に起こった素反応の回数は、反応機構モデル形成部１２４により構築された素反応のうち同一の素反応の数に相当する。例えば、設定温度Ｔ１（例えば３０００Ｋ）において、反応機構モデル形成部１２４により構築された素反応Ａが１０個存在する場合、設定温度Ｔ１において分子動力学計算中に起こった素反応Ａの回数は１０となる。同一素反応の個数のカウントは反応機構モデル形成部１２４（例えば素反応構築部１２２）により行われても良いし、ベースモデル作成部１２６により行われても良いが、いずれにしろカウントされた数が、各設定温度において分子動力学計算中に起こった各素反応の回数となる（以下、単に素反応の回数と称する場合がある）。

なお、ベースモデルは、必ずしも分子動力学計算中に起こった素反応の回数を含まなくても良い。この場合、ベースモデルには、同一素反応を１つの素反応としてまとめることなく、全てベースモデルに含ませる必要がある。例えば、設定温度Ｔ１（例えば３０００Ｋ）において、反応機構モデル形成部１２４により構築された素反応Ａが１０個存在する場合、ベースモデルの素反応リストＬ１には、１０個の素反応Ａが存在することになる。

次に、ステップＳ５２では、簡略化モデル作成部１２８により、上記作成されたベースモデルの中で、各設定温度における素反応の回数が所定回数（ａ）以上である素反応を抽出し（所定回数（ａ）未満の素反応を削除し）、各設定温度における素反応の回数がａ以上の素反応を含む暫定簡略化モデルＭ１が作成される。例えば、設定温度（Ｔ_１）の素反応リストＬ_１において素反応Ａの回数（ａ１）が１０、素反応Ｂの回数（ａ２）が４、素反応Ｃの回数（ａ３）が２で、所定値（ａ）を３に設定した場合、素反応Ｃが削除され、素反応Ａ及びＢを含む素反応リストＬ_１’が作成される。その他のリストＬ_２〜Ｌ_Ｎも同様に行う。このようにして作成された暫定簡略化モデルＭ１は、例えば設定温度（Ｔ_１）〜（Ｔ_Ｎ）毎に、素反応の回数が所定回数以上（例えば３回以上）起こった各素反応、及び各素反応の反応速度定数を含む素反応リストＬ_１’〜Ｌ_Ｎ’として表されるモデルとなる。なお、ベースモデルに素反応の回数が含まれていない場合、ベースモデルに含まれる同一の素反応の数をカウントし、その数が所定回数以上である素反応を抽出し、暫定簡略化モデルＭ１が作成される。

次に、ステップＳ５４では、簡略化モデル作成部１２８により、暫定簡略化モデルＭ１の中から、同一素反応が存在する設定温度の個数が所定個数（ｂ）以上である素反応を抽出した暫定簡略化モデルＭ２が作成される。例えば、素反応Ａが、５つの設定温度に存在（５つの素反応リストに存在）し、素反応Ｂが３つの設定温度に存在（３つの素反応リストに存在）し、素反応Ｃが４つの設定温度に存在（４つの素反応リストに存在）し、所定個数（ｂ）を４に設定した場合、素反応Ｂが削除され、素反応Ａ及び素反応Ｃを含む暫定簡略化モデルＭ２が作成される。すなわち、図１０に示す暫定簡略化モデルＭ２は、所定回数以上の反応が起こった素反応であり、且つ所定個数以上の設定温度に存在する素反応（すなわち、所定個数以上の設定温度で反応が起こった素反応）を含むモデルとなる。なお、暫定簡略化モデルＭ２は、例えば設定温度（Ｔ_１）〜（Ｔ_Ｎ）毎の素反応リストとして表されるモデルであってもよい。本実施形態では、暫定簡略化モデルＭ２を簡略化モデルとして、モデル精度評価部１３０に送信する。

このように、本実施形態では、分子動力学計算中に起こった素反応の回数及び同一の素反応が存在する設定温度の個数に基づいて、ベースモデル中の素反応を絞り込み、ベースモデルより素反応の数を減らした簡略化モデルを作成している。しかし、ベースモデル中の素反応の絞り込みはこれに制限されるものではない。例えば、ベースモデル中の素反応毎に、素反応の回数を合計し、その合計回数が所定回数以上である素反応を抽出したモデルを簡略化モデルとしてもよい。したがって、簡略化モデル作成部１２８では、少なくとも分子動力学計算中に起こった素反応の回数に基づいて、ベースモデル中の素反応を絞り込んだ簡略化モデルが作成されればよい。ベースモデル中の素反応の絞り込みは、簡略化モデルの精度の点等から、各設定温度において分子動力学計算中に起こった素反応の回数及び同一の素反応が存在する設定温度の個数に基づいてベースモデル中の素反応を絞り込んだ簡略化モデルが作成されることが好ましい。

簡略化モデル作成部１２８により作成された簡略化モデルは、各素反応における反応速度パラメータ（アレニウスパラメータ）Ａ，ｎ，Ｅａを含むことが好ましい。反応速度パラメータは、ベースモデル又は暫定簡略化モデルＭ１等を用いて、図８に示すアレニウスの式のフィッティングを行うことにより求められる。反応速度パラメータを算出するために、反応速度パラメータ算出手段を別途設けても良いし、簡略化モデル作成部１２８又はベースモデル作成部１２６に、反応速度パラメータを算出する機能を付加してもよい。また、反応機構モデル形成部１２４において、予め反応速度パラメータを算出している場合には、ベースモデルに、予め算出した反応速度パラメータを含ませても良い。

次に、ステップＳ５６では、モデル精度評価部１３０により、簡略化モデルの精度が評価される。精度の評価は、例えば、簡略化モデルの反応特性と、ベースモデルの反応特性とを比較し、その差が許容範囲内であるか否かにより評価される。そして、反応特性の差が許容範囲内である場合には、簡略化モデルは簡略化した反応機構モデルとして、データ保持部１１４に送信され、反応特性の差が許容範囲を超える場合には、上記設定した所定値（分子動力学計算中に起こった素反応の回数及び同一素反応が存在する設定温度の個数）を再設定して、再度簡略化モデルを作成する。所定値の再設定は、所定値再設定部等を設けても良いし、入力装置１０２から作業者等が直接入力してもよい。反応特性の差が許容範囲を超える場合、簡略化モデルの作成において、反応特性に影響を与える重要な素反応が削除されていると考えられるため、再設定される所定値は、素反応の数が前の簡略化モデルの素反応の数より多くなるように設定される必要がある。

モデル精度評価部１３０により求められる反応特性は、例えば、化学動力学計算により求められるものであれば特に制限されるものではなく、例えば、モル分率の時間履歴や着火遅れ時間等が挙げられる。反応特性の計算には、例えば、比熱、エントロピー、エンタルピー等の熱力学データ、反応速度パラメータ（アレニウスパラメータ）等が必要である。熱力学データは、モデル精度評価部１３０に予め記憶させておくことが望ましいが、反応特性の計算の際に、入力装置１０２から作業者により入力されてもよい。

モル分率の時間履歴を反応特性として用いる場合、ベースモデルにおけるモル分率の時間履歴は、図７に示すように予め作成された分子数の時間履歴をモル分率に変換して、それを援用してもよいし、簡略化モデルと同様に、化学動力学計算により求めても良い。着火遅れ時間は、ベースモデル、簡略化モデルを用いて、それぞれ化学動力学計算により求められる。

また、反応特性の差が許容範囲内である場合であっても、例えば、許容範囲内である簡略化モデルをベースに、更なる素反応の簡略化を行ってもよい。以下、その一例について説明する。

まず、ベースモデル作成部１２６に上記簡略化モデルが送信される。この際、素反応を絞り込むために設定される所定値（分子動力学計算中に起こった素反応の回数及び同一の素反応が存在する設定温度の個数）が再設定される。以下、分子動力学計算中に起こった素反応の回数を制限する所定値をａ、同一の素反応が存在する設定温度の個数を制限する所定値をｂとする。

再設定される所定値は、簡略化モデル作成部１２８により新たに作成される簡略化モデルの素反応の数（及び化学種の数）が、ベースモデル作成部１２６に送信された簡略化モデルの素反応の数（及び化学種の数）より少なくなるように設定される。原則、所定値ａ，ｂの値が大きいほど、素反応の数（及び化学種の数）は減少するため、例えば、最初に設定された所定値ａが３、所定値ｂが４の場合、再設定の際には、所定値ａを３より大きい値（例えば１０）、又は所定値ｂを４より大きい値（例えば５）、或いはそれらの両方を行う。そして、簡略化モデル作成部１２８により、再設定された所定値を満たさない素反応が削除され、新たな簡略化モデルが作成される。なお、所定値ａ，ｂはそれぞれ１つに限らず、複数設定して、新たな簡略化モデルを複数作成してもよい。

次に、簡略化モデル作成部１２８により新たに作成された簡略化モデルの精度が、モデル精度評価部１３０により評価される。精度評価は、ベースモデル作成部１２６により作成された最初のベースモデルの反応特性と新たに作成された簡略化モデルとを比較し、その差が予め設定した許容範囲を満たすか否かを判断することが望ましいが、１つ前の簡略化モデルと新たに作成された簡略化モデルとを比較して、その差が予め設定した許容範囲を満たすか否かを判断しても良い。そして、許容範囲を満たす場合には、新たに作成された簡略化モデルを簡略化した反応機構モデルとして、データ保持部１１４に送信する。

また、新たな簡略化モデルを複数作成した場合には、モデル精度評価部１３０において、許容範囲を満たす簡略化モデル全てをデータ保持部１１４に送信してもよいが、例えば、許容範囲を満たす簡略化モデルのうち、素反応の数（及び化学種の数）が最も少ない簡略化モデルを、簡略化した反応機構モデルとして、データ保持部１１４に送信してもよい。

本実施形態によれば、モデル精度を維持しつつ、素反応の数（及び化学種の数）を減らすことが可能となるため、当該反応機構モデルを用いた三次元数値流体シミュレーションの実行が容易となる。例えば、簡略化した反応機構モデルを用いることで、エンジン燃焼シミュレーションに掛かる時間を短縮することが可能となる。

以下に、実施例を挙げて本発明をより具体的に説明するが、本発明は以下の実施例に制限されるものではない。

（実施例１）
図１１は、水素分子と酸素分子をシミュレーションセル中に配置した状態を示す図である。実施例１では、１辺が２５オングストロームの正方体としたシミュレーションセル中に、水素分子６６分子、酸素分子３３分子をランダムな位置にセットした。次に、反応速度定数を計算したい温度となるように、各原子にランダムな初期速度を与えた。初期原子位置、初期原子速度のもと、上式（１）で表されるニュートンの運動方程式に従って、各分子を構成する原子の位置情報を逐次計算した。運動方程式は時間に関する差分法によって計算した。運動方程式の解法アルゴリズムには様々なものがあるが、実施例１では速度ベルレ法を用いた。計算の時間ステップは０．１フェムト秒とした。また、反応系の温度を一定に保つためにＮｏｓｅ−Ｈｏｏｖｅｒサーモスタットを用いた。また、各原子に働く力を計算する方法として第一原理法、半経験的方法、分子力場法等があるが、実施例１では反応分子力場法ＲｅａｘＦＦを使用した。

実施例１では、反応分子力場法の計算過程で得られる結合次数を用いて、各原子間の結合の有無を判定した。結合次数の閾値として、Ｈ−Ｈを０．５５、Ｏ−Ｏを０．６５、Ｏ−Ｈを０．４とし、結合次数が閾値以上の場合に結合ありと判定した。結合有りを「１」、結合無しを「０」として行列情報を作成した。行列情報の更新は１０００時間ステップ毎に行った。この行列情報を図４に示すようにブロック対角化して、原子集団を分子に分割した。このようにして、１０００時間ステップ毎に反応系に含まれる分子を特定した。図１２に実施例１で特定した分子数の時間変化をまとめた。図１２では、Ｈ_２、Ｏ_２、ＯＨ、Ｈ_２Ｏのみを示したが、解析結果には、Ｈ_２Ｏ_２、ＨＯ_２等の分子も特定された。そして、１０００時間ステップ毎の計算結果において、原子の結合情報の変化前の分子を反応分子として左辺に、変化後の分子を生成分子として右辺に記述して反応式を構築した。この際、例えば、Ｈ＋Ｏ_２＋ＯＨ→ＯＨ＋Ｏ＋Ｈ＋Ｏ_２のような複数の素反応が混合した反応式が抽出されることを抑制するために、反応分子と生成分子の原子の関連性に基づいて、図５及び６に示すような素反応の解析を行った。図１３に、実施例１で構築された素反応の一覧を示す。

次に、得られた素反応の反応速度定数を上式（７）、（８）を用いて計算した。実施例１では、２５００Ｋ、２７５０Ｋ、３０００Ｋ、３２５０Ｋ、３５００Ｋ、４０００Ｋの温度におけるＨ_２＋Ｏ→ＯＨ＋Ｈの素反応に対して反応速度定数を計算した。その結果を図１４に示す。図１４に示す反応速度定数は、１０回の独立した分子動力学計算を行い、統計平均を取った後の値である。なお、参考のため、実験データ等を元に構築された反応機構による値も示した。図１４に示すように、本実施例１で得られた反応速度定数の結果は、実験データ等を元に構築された反応機構による値を良く再現していると言える。

（実施例２）
実施例２では、１辺が２５オングストロームの正方体としたシミュレーションセル中に、メタン分子５０分子、酸素分子１００分子をランダムな位置にセットし、３０００Ｋ、３２５０Ｋ、３５００Ｋ、３７５０Ｋ、４０００Ｋの設定温度で、実施例１と同様に分子動力学計算（１０回）を行い、素反応、分子動力学計算中に起こった素反応の回数、反応速度定数を求めた。図１５に、３０００Ｋにおける素反応、反応速度定数、及び分子動力学計算中に起こった素反応の回数を含むベースモデルの一部を示す。

所定値ａを３に設定し、各設定温度における素反応、素反応の回数、反応速度定数を含むベースモデルから、素反応の回数が３回未満の素反応を削除し、素反応の回数が３回以上の素反応を抽出した暫定簡略化モデルを作成した。なお、素反応の回数の少ない素反応の反応速度定数は統計誤差が大きく値の信頼性が低い。上記のように、３回以上起こった素反応に絞り込むことで、反応速度定数の誤差の大きい素反応を削除することができる。

次に、素反応の回数が３回以上の素反応に絞り込んだ暫定簡略化モデルを用いて、図８に示すアレニウス式のフィッティングを行い、反応速度パラメータ（Ａ、ｎ、Ｅａ）を求めた。また、暫定簡略化モデルの各素反応それぞれにおいて、同一の素反応が存在する設定温度の個数を算出した。

図１６に、素反応、パラメータ（Ａ、ｎ、Ｅａ）、同一の素反応が存在する設定温度の個数を含む簡略化モデルの一部を示す。所定値ｂを４に設定し、図１６に示す簡略化モデルから、同一の素反応が存在する設定温度の個数が４個未満の素反応を削除し、４個以上の素反応を抽出した３_４簡略化モデルを作成した。以下、素反応の回数をａ以上、且つ同一素反応が存在する設定温度の個数をｂ以上として、素反応を絞り込んだ簡略化モデルをａ_ｂ簡略化モデルと称する。

図１７に、３_４簡略化モデルを用いた３０００Ｋにおけるモル分率の時間変化及び分子動力学計算により得られるモル分率の時間変化を示す。対象とする化学種をメタン、酸素、ホルムアルデヒド、水、一酸化炭素、二酸化炭素とした。３_４簡略化モデルを用いたモル分率の時間変化は、化学動力学計算により求めた。また、分子動力学計算により得られるモル分率の時間変化は、実施例２で行った分子動力学計算の計算において、各分子数の時間変化を求め、その結果をモル分率に換算することにより求めた。

図１７に示すように、３_４簡略化モデルを用いたモル分率の時間変化は、分子動力学計算により得られるモル分率の時間変化をよく再現した結果となった。具体的には、全時間において、３_４簡略化モデルを用いたモル分率は、分子動力学計算により得られるモル分率に対して５．７％以下に収まっており、高い精度の簡略化モデルが得られたと判断できる。

３_４簡略化モデルは、素反応数２５３、化学種数６４であった。一方、ベースモデルは、素反応数３７０、化学種数９２であった。すなわち、３_４簡略化モデルは、ベースモデルと比較して、モデル精度を維持しつつ、素反応数及び化学種数が簡略化された反応モデルであると言える。

次に、３_４簡略化モデルをベースとして、さらに、素反応数の簡略化を実行した。

まず、素反応の回数の所定値を３〜４０、素反応が存在する設定温度の数を２〜５の範囲で設定し、前述の素反応の絞り込みと同様にして、簡略化モデル（３_２簡略化モデル〜４０_５簡略化モデル）を作成した。

図１８に、各簡略化モデル中の素反応の数及び化学種の数をまとめた。なお、図１８における折れ線の一番左の素反応の数、化学種の数は、最初に設定したベースモデルに対する素反応の数、化学種の数であり、一番左から２番目以降の素反応の数、化学種の数が、簡略化モデル（３_２簡略化モデル〜４０_５簡略化モデル）に対応している。

次に、素反応を絞り込んだ各簡略化モデルに対して、定容断熱条件下での着火遅れ時間（温度１０００Ｋ、１５００Ｋ、２０００Ｋ、２５００Ｋ、３０００Ｋ）を化学動力学計算により求めた。各初期温度から１０００Ｋ上昇した時点を着火時点とみなし、着火遅れ時間を計算した。

図１９に、３_２簡略化モデル〜２０_５簡略化モデルの着火遅れ時間の計算結果を示す。図１９に示す着火遅れ時間は、ベースとなる３_４簡略化モデルの着火遅れ時間を基準とし、その着火遅れ時間に対するその他の簡略化モデルの着火遅れ時間の比率（％）で表している。なお、図１８から明らかなように、３_２簡略化モデル、３_３簡略化モデル、３_５簡略化モデルは、素反応の数、化学種の数はいずれも、３_４簡略化モデルより多いので、本来であれば、着火遅れ時間の評価は行わなくてよい。

３_４簡略化モデルの着火遅れ時間に対する３_４簡略化モデル以外の簡略化モデルの着火遅れ時間の比率の許容範囲を２０％以下に設定し、簡略化モデルの精度を評価したところ、５_４簡略化モデル、５_５簡略化モデル、１０_３簡略化モデル、１０_４簡略化モデル、１０_５簡略化モデルが許容範囲内であった。なお、３_３簡略化モデル、３_５簡略化モデルも許容範囲内であるが、前述したように素反応の数、化学種の数はいずれも、３_４簡略化モデルより多いので、評価対象外である。

５_４簡略化モデル、５_５簡略化モデル、１０_３簡略化モデル、１０_４簡略化モデル、１０_５簡略化モデルを簡略化した反応機構モデルとして特定してもよいが、これらの中で、素反応の数及び化学種の数が最も少ない１０_５簡略化モデルを反応モデルとして特定することが望ましい。なお、実施例２では、３_４簡略化モデルの着火遅れ時間と比較しているが、当然最初に設定したベースモデルの着火遅れ時間と比較してもよい。

図２０に、１０_５簡略化モデルを用いた３０００Ｋにおけるモル分率の時間変化を示す。また、比較のため、３_４簡略化モデルを用いたモル分率の時間変化及び分子動力学計算により得られるモル分率の時間変化を示す。対象とする化学種をメタン、酸素、ホルムアルデヒド、水、一酸化炭素、二酸化炭素とした。

図２０に示すように、１０_５簡略化モデルを用いたモル分率の時間変化は、分子動力学計算により得られるモル分率の時間変化をよく再現した結果となった。具体的には、全時間において、１０_５簡略化モデルを用いたモル分率は、分子動力学計算により得られるモル分率に対して１４．８％以下に収まっていた。

１０_５簡略化モデルは、素反応数１０９、化学種数３７であった。すなわち、１０_５簡略化モデルは、素反応数３７０、化学種数９２のベースモデル、素反応数２５３、化学種数６４の３_４簡略化モデルと比較して、モデル精度を維持しつつ、大幅に素反応数及び化学種数が簡略化された反応モデルであると言える。

（参考例）
３_４簡略化モデルをベースとして、従来知られているＤＲＧ（ＤｉｒｅｃｔｅｄＲｅｌａｔｉｏｎＧｒａｐｈ）法、ＤＲＧＥＰ（ＤＲＧｗｉｔｈＥｒｒｏｒＰｒｏｐａｇａｔｉｏｎ）法、ＤＲＧＥＰＳＡ（ＤＲＧＥＰａｎｄＳｅｎｓｉｔｉｖｉｔｙＡｎａｌｙｓｉｓ）法の３つの簡略化法を実行した。この３つの簡略化法においては、実施例２と同様に、得られる簡略化モデルの定容断熱条件下での着火遅れ時間が、３_４簡略化モデルの着火遅れ時間に対して２０％以下となるように、モデルの簡略化を実行した。また、１０_５簡略化モデルをベースとして、上記３つの簡略化法を実行した。

図２１に、従来知られている３つの簡略化法を用いて得られた簡略化モデルにおける素反応の数と化学種の数の結果を示す。図２１に示すように、３_４簡略化モデルをベースとした場合、ＤＲＧ法では、素反応数１４１、化学種数２５、ＤＲＧＥＰ法では、素反応数１３１、化学種数２２、ＤＲＧＥＰＳＡ法では、素反応数１１２、化学種数２１となった。実施例２では、３_４簡略化モデルをベースとして得られた１０_５簡略化モデルは、素反応数１０９、化学種数３７であったので、実施例２の方が、従来の簡略化法と比較して、素反応の数の点で、より簡略化を行うことができていると言える。

また、１０_５簡略化モデルをベースとした場合、ＤＲＧ法では、素反応数９５、化学種数２６、ＤＲＧＥＰ法では、素反応数９１、化学種数２５、ＤＲＧＥＰＳＡ法では、素反応数８２、化学種数２４となった。この結果から、実施例２における簡略化と、従来の簡略化法を組み合わせることで、さらに素反応数を簡略化した反応モデルを作成することが可能となる。

１００，１０１反応機構生成装置、１０２入力装置、１０４出力装置、１０６反応条件設定部、１０８初期条件設定部、１１０分子動力学演算部、１１２反応速度定数算出部、１１４データ保持部、１１６表示制御部、１２０分子特定部、１２２素反応構築部、１２４反応機構モデル形成部、１２６ベースモデル作成部、１２８簡略化モデル作成部、１３０モデル精度評価部。

Claims

反応系内の各分子を構成する原子について、時間ステップ毎に分子動力学計算を行うステップと、
前記時間ステップ前後で前記反応系内に化学反応が起こった場合、前記化学反応に寄与した反応分子及び生成分子を特定するステップと、
前記反応分子と前記生成分子との原子の種類の関連性に基づいて、前記関連性のある反応分子及び生成分子から構成された素反応を構築するステップと、
前記構築した素反応の反応速度定数を算出するステップと、を含む反応機構生成方法。
前記反応分子及び前記生成分子を特定するステップでは、前記時間ステップ前後の分子動力学計算から得られる原子の位置情報から、前記化学反応に寄与した前記反応分子及び前記生成分子を特定することを特徴とする請求項１記載の反応機構生成方法。
前記素反応の反応速度定数を算出するステップでは、前記反応分子と同じ分子種の分子の個数の時間履歴を算出し、前記反応分子と同じ分子種の分子の個数の時間履歴に基づいて、前記素反応の反応速度定数を算出することを特徴とする請求項１記載の反応機構生成方法。
前記素反応を構築するステップにおいて構築された複数の素反応を含むベースモデルを作成するステップと、
前記分子動力学計算中に起こった素反応の回数に基づいて、前記ベースモデルの前記複数の素反応を絞り込み、前記ベースモデルより素反応の数が少ない簡略化モデルを作成するステップと、を含むことを特徴とする請求項１記載の反応機構生成方法。
前記簡略化モデルの精度を評価するステップを含むことを特徴とする請求項４記載の反応機構生成方法。
前記素反応を構築するステップにおいて構築された複数の素反応を含むベースモデルを作成するステップと、
前記分子動力学計算中に起こった素反応の回数、及び前記分子動力学計算において設定された複数の設定温度のうち、同一の素反応が存在する設定温度の個数に基づいて、前記ベースモデルの前記複数の素反応を絞り込み、前記ベースモデルより素反応の数が少ない簡略化モデルを作成するステップと、を含むことを特徴とする請求項１記載の反応機構生成方法。
前記簡略化モデルの精度を評価するステップを含むことを特徴とする請求項６記載の反応機構生成方法。
反応系内の各分子を構成する原子について、時間ステップ毎に分子動力学計算を行う力学計算演算部と、
前記時間ステップ前後で前記反応系内に化学反応が起こった場合、前記化学反応に寄与した反応分子及び生成分子を特定する分子特定部と、
前記反応分子と前記生成分子との原子の種類の関連性に基づいて、前記関連性のある反応分子及び生成分子から構成された素反応を構築する素反応構築部と、
前記構築した素反応の反応速度定数を算出する反応速度定数算出部と、を含むことを特徴とする反応機構生成装置。
前記分子特定部は、前記時間ステップ前後の分子動力学計算から得られる原子の位置情報から、前記化学反応に寄与した反応分子及び生成分子を特定することを特徴とする請求項８記載の反応機構生成装置。
前記反応速度定数算出部は、前記反応分子と同じ分子種の分子の個数の時間履歴を算出し、前記反応分子と同じ分子種の分子の個数の時間履歴に基づいて、前記素反応の反応速度を算出することを特徴とする請求項８記載の反応機構生成装置。
前記素反応構築部において構築された複数の素反応を含むベースモデルを作成するベースモデル作成部と、
前記分子動力学計算中に起こった素反応の回数に基づいて、前記ベースモデルの前記複数の素反応を絞り込み、前記ベースモデルより素反応の数が少ない簡略化モデルを作成する簡略化モデル作成部と、を含むことを特徴とする請求項８記載の反応機構生成装置。
前記簡略化モデルの精度を評価するモデル精度評価部を含むことを特徴とする請求項１１記載の反応機構生成装置。
前記素反応構築部において構築された複数の素反応を含むベースモデルを作成するベースモデル作成部と、
前記分子動力学計算中に起こった素反応の回数、及び前記分子動力学計算において設定された複数の設定温度のうち、同一の素反応が存在する設定温度の個数に基づいて、前記ベースモデルの前記複数の素反応を絞り込み、前記ベースモデルより素反応の数が少ない簡略化モデルを作成する簡略化モデル作成部と、を含むことを特徴とする請求項８記載の反応機構生成装置。
前記簡略化モデルの精度を評価するモデル精度評価部を含むことを特徴とする請求項１３記載の反応機構生成装置。