特許 7302621 有機材料輸送特性評価装置及び有機材料輸送特性評価方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2023-06-26

(45)【発行日】2023-07-04

(54)【発明の名称】有機材料輸送特性評価装置及び有機材料輸送特性評価方法

(51)【国際特許分類】

   G16C 10/00 20190101AFI20230627BHJP

   G01N 19/00 20060101ALI20230627BHJP

   G01N 25/18 20060101ALI20230627BHJP

   G01N 11/00 20060101ALN20230627BHJP

【ＦＩ】

G16C10/00

G01N19/00 A

G01N25/18 L

G01N11/00 Z

【請求項の数】 5

(21)【出願番号】P 2021063594

(22)【出願日】2021-04-02

(65)【公開番号】P2022158584

(43)【公開日】2022-10-17

【審査請求日】2022-03-23

(73)【特許権者】

【識別番号】000003609

【氏名又は名称】株式会社豊田中央研究所

(74)【代理人】

【識別番号】110001210

【氏名又は名称】弁理士法人ＹＫＩ国際特許事務所

(72)【発明者】

【氏名】名児耶 彰洋

(72)【発明者】

【氏名】大庭 伸子

(72)【発明者】

【氏名】馬場 健

(72)【発明者】

【氏名】梶田 晴司

【審査官】塩田 徳彦

(56)【参考文献】

【文献】特開２０１９－１４８４４６（ＪＰ，Ａ）

【文献】特開２０２０－０８７４５６（ＪＰ，Ａ）

【文献】特開２０１４－１０６５５４（ＪＰ，Ａ）

【文献】米国特許出願公開第２００２／００６１５４０（ＵＳ，Ａ１）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ１６Ｃ １０／００－９９／００

Ｇ０１Ｎ １９／００

Ｇ０１Ｎ ２５／１８

Ｇ０１Ｎ １１／００

(57)【特許請求の範囲】

【請求項1】

有機材料の輸送特性を評価する有機材料輸送特性評価装置であって、
コンピュータを、
前記有機材料における流体の輸送特性Ｋが

【数1】

で表され、
輸送特性Ｋに対応する物理量Ｑ（ｔ）の時間微分であるフラックス［Ｑ（ｔ）ドット］において前記有機材料の分子間のフラックスの相関を除いて輸送特性Ｋを算出する手段として機能させることを特徴とする有機材料輸送特性評価装置。

【請求項2】

請求項１に記載の有機材料輸送特性評価装置であって、
前記有機材料の分子間の熱流相関を除いてアモルファス構造を有する直鎖高分子鎖材料の熱伝導度を計算することを特徴とする有機材料輸送特性評価装置。

【請求項3】

請求項１に記載の有機材料輸送特性評価装置であって、
前記有機材料の分子間の応力相関を除いて有機材料の粘性を計算することを特徴とする有機材料輸送特性評価装置。

【請求項4】

請求項１に記載の有機材料輸送特性評価装置であって、
前記有機材料の分子間の速度相関を除いて有機材料の拡散係数を計算することを特徴とする有機材料輸送特性評価装置。

【請求項5】

有機材料の輸送特性を評価する有機材料輸送特性評価方法であって、
コンピュータに、
前記有機材料における流体の輸送特性Ｋが

【数2】

で表され、
輸送特性Ｋに対応する物理量Ｑ（ｔ）の時間微分であるフラックス［Ｑ（ｔ）ドット］において前記有機材料の分子間のフラックスの相関を除いて輸送特性Ｋを算出させることを特徴とする有機材料輸送特性評価方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、有機材料輸送特性評価装置及び有機材料輸送特性評価方法に関する。

【背景技術】

【0002】

有機材料について熱伝導度、粘度及び拡散係数等の時間的な物理量の変化を示す輸送特性について、平衡分子動力学（ＥＭＤ）計算によるＧｒｅｅｎ－Ｋｕｂｏ（ＧＫ）公式を用いる方法や非平衡分子動力学（ＮＥＭＤ）計算による直接計算法が用いて評価が行われている。

【0003】

Ｇｒｅｅｎ－Ｋｕｂｏ（ＧＫ）公式を用いた平衡分子動力学（ＥＭＤ）計算を高速に行う手法が提案されている。例えば、必要な精度を得るための計算量を誤差解析によって決定する方法が提案されている（非特許文献１）。また、Ｇｒｅｅｎ－Ｋｕｂｏ（ＧＫ）公式を用いた分子動力学（ＥＭＤ）計算において熱流の自己相関関数（ＨＦＡＣＦ）をランダムウォークでモデル化し、調和解析を行って誤差要因となる周波数成分を分離することで、計算誤差の標準偏差を求める方法が提案されている（非特許文献２）。また、赤池情報量基準を用いたモデル精度評価と調和解析による方法が提案されている（非特許文献３）。また、Ｇｒｅｅｎ－Ｋｕｂｏ（ＧＫ）公式において、空間分割して計算することで不均一系を解析する方法が提案されている（特許文献１）。また、敵対的生成ネットワークを用いた分子動力学（ＥＭＤ）計算における負荷低減が提案されている（特許文献２）。

【先行技術文献】

【特許文献】

【0004】

【文献】特開２０１９－１４８４４６号公報

【文献】特開２０２０－８７４５６号公報

【非特許文献】

【0005】

【文献】R. E. Jones, J. Chem. Phys. 136, 154102 (2012)

【文献】L. Ercole, Scientific Reports, 7, 15835 (2017)

【文献】L. S. Oliveira, Phys. Rev. E 95, 023308 (2017)

【発明の概要】

【発明が解決しようとする課題】

【0006】

ところで、従来技術の全般において、調和解析による誤差の分離において閾値となるパラメータが必要となる。そのため、事前検討やパラメータの設定等のための前処理が必要であり、材料探索の自動化のためにより簡便な装置や方法が望まれている。

【0007】

非特許文献１に記載の技術では、目的の精度を得るための最小の計算量を求めることができるが、本質的な誤差の低減には繋がらない。また、調和解析や誤差の分離のための解析が必要である。また、非特許文献２に記載の技術では、調和解析による誤差の分離等の解析が必要であり、自動計算探索には不向きである。また、特許文献１に記載の技術では、有機材料の輸送特性の評価を高速化する処理についての記載はない。また、特許文献２に記載の技術では、機械学習を用いているために事前の学習が必要であり、多くの場合には計算速度が低く、様々な材料の自動探索には不向きである。

【課題を解決するための手段】

【0008】

本発明の１つの態様は、有機材料の輸送特性の評価する有機材料輸送特性評価装置であって、前記輸送特性に対する前記有機材料の分子間のフラックスの相関を除いて前記輸送特性を算出することを特徴とする有機材料輸送特性評価装置である。

【0009】

ここで、Ｇｒｅｅｎ－Ｋｕｂｏ公式を用いて前記輸送特性を算出する際に、前記有機材料の分子間のフラックスの相関を除いて前記輸送特性を算出することが好適である。

【0010】

また、Ｇｒｅｅｎ－Ｋｕｂｏ公式を用いて、分子間の熱流相関を除いてアモルファス構造を有する直鎖高分子鎖材料の熱伝導度を計算することが好適である。

【0011】

また、Ｇｒｅｅｎ－Ｋｕｂｏ公式を用いて、分子間の応力相関を除いて有機材料の粘性を計算することが好適である。

【0012】

また、Ｇｒｅｅｎ－Ｋｕｂｏ公式を用いて、分子間の速度相関を除いて有機材料の拡散係数を計算することが好適である。

【0013】

本発明の別の態様は、有機材料の輸送特性の評価する有機材料輸送特性評価方法であって、前記輸送特性に対する前記有機材料の分子間のフラックスの相関を除いて前記輸送特性を算出することを特徴とする有機材料輸送特性評価方法である。

【発明の効果】

【0014】

本発明によれば、有機材料の輸送特性をより高速に評価することが可能になる。

【図面の簡単な説明】

【0015】

【図1】本発明の実施の形態における有機材料輸送特性評価装置の構成を示す図である。

【図2】ポリエチレンに対する熱流の相関関数の評価結果を示す図である。

【図3】ポリエチレンに対する熱伝導度の評価結果を示す図である。

【図4】ポリスチレンに対する熱伝導度の評価結果を示す図である。

【図5】熱伝導度の評価結果を示す図である。

【図6】典型的な直鎖高分子材料の単量体の分子構造を示す図である。

【図7】単成分系における応力ゆらぎの時間相関の評価結果を示す図である。

【図8】２成分系における応力ゆらぎの時間相関の評価結果を示す図である。

【図9】混合系における混合比に応じた成分毎の拡散係数の評価結果を示す図である。

【発明を実施するための形態】

【0016】

本発明の実施の形態における有機材料輸送特性評価装置１００は、図１に示すように、処理部１０、記憶部１２、入力部１４、出力部１６及び通信部１８を含んで構成される。処理部１０は、ＣＰＵ等の演算処理を行う手段を含む。処理部１０は、記憶部１２に記憶されている有機材料輸送特性評価プログラムを実行することによって、本実施の形態における有機材料輸送特性評価方法における機能を実現する。記憶部１２は、半導体メモリやメモリカード等の記憶手段を含む。記憶部１２は、処理部１０とアクセス可能に接続され、有機材料輸送特性評価プログラム、その処理に必要な情報を記憶する。入力部１４は、情報を入力する手段を含む。入力部１４は、例えば、ユーザからの入力を受けるキーボード、タッチパネル、ボタン等を備える。出力部１６は、ユーザから入力情報を受け付けるためのユーザインターフェース画面（ＵＩ）等の有機材料輸送特性評価装置１００での処理結果を出力する手段を含む。出力部１６は、例えば、ユーザに対して画像を呈示するディスプレイを備える。通信部１８は、情報通信網１０２を介して、外部の情報処理装置との情報の通信を行うインターフェースを含んで構成される。通信部１８による通信は有線及び無線を問わない。

【0017】

有機材料輸送特性評価装置１００は、有機材料輸送特性評価プログラムを実行可能な情報処理装置であれば様々なものを適用できる。例えば、有機材料輸送特性評価装置１００としては、据置型又は携帯型パーソナルコンピュータ（ＰＣ）、タブレット型コンピュータ等が使用できる。

【0018】

以下、本実施の形態における有機材料輸送特性評価方法について説明する。有機材料輸送特性評価方法は、有機材料輸送特性評価装置１００において以下の処理を実行する有機材料輸送特性評価プログラムを実行することによって実現される。

【0019】

流体における輸送特性Ｋはそれに対応する物理量Ｑ（ｔ）とその時間微分であるフラックス∂Ｑ（ｔ）／∂ｔ（式中ではＱ（ｔ）にドットを打って示す）によって以下のように表される。

【数1】

【0020】

［熱伝導度の評価］
有機材料の輸送特性として熱伝導度を評価対象とする場合、数式（２）で表されるフラックスを数式（１）に適用することで数式（３）のように熱伝導度κが得られる。ここで、ｋＢはボルツマン定数、ｅ_ｉはｉ番目の原子のエネルギー、ｖ_ｉはｉ番目の原子の速度及びＳ_ｉはｉ番目の原子の応力を示す。

【数2】

【数3】

【0021】

系の熱流Ｊ（ｔ）は原子毎の寄与Ｊｉの和として表される。この和を、同一分子に属する原子について和をとったＪｍｏｌで置き換え、異なる分子に属する原子間の相関を無視することで、数式（３）の熱伝導度κは数式（４）のように熱伝導度κの近似式として表すことができる。ここで、数式（４）に置いてダッシュを付したΣ記号は各分子についての和を取る演算子を意味する。

【数4】

【0022】

高分子の熱伝導度κにほとんど寄与しない分子間の熱流相関を無視することによって、そこに含まれるランダムな誤差が取り除かれる。これによって、有機材料輸送特性評価における統計誤差が低減され、複数の初期配置に対して平均を算出する必要がなくなり、評価における計算負荷を大幅に削減することができる。

【0023】

［実施例１］
Ｇｒｅｅｎ－Ｋｕｂｏ公式を用いて計算した熱伝導度に対する同一分子内及び異なる分子間の熱流の相関関数の寄与について検討した。

【0024】

アモルファス構造の直鎖高分子材料の熱伝導度κについて、Ｇｒｅｅｎ－Ｋｕｂｏ（ＧＫ）公式を用いた平衡分子動力学（ＥＭＤ）計算による有機材料輸送特性評価を行った。

【0025】

ポリマー１分子の長さを約１２０Åとし、原子数Ｎとして５万原子からなるアモルファス構造を作成して分子動力学（ＭＤ）計算を行った。ポリマーの初期構造の作成、及び、分子動力学（ＭＤ）計算にはプログラムＥＭＣ（“Enhanced Monte Carlo (EMC),” http://montecarlo.sourceforge.net）、及び、ＬＡＭＭＰＳ（S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” Journal of Computational Physics 117, 1 (1995).）を用いた。力場には水素原子を省略するＵｎｉｔｅｄ Ａｔｏｍ力場であるＯＰＬＳ－ＵＡを用いた。タイムステップには２．０ｆｓを用いた。構造の平衡化には、Ｎｏｓe－Ｈｏｏｖｅｒの熱浴（W. G. Hoover, “Canonical dynamics: Equilibrium phase-space distributions,” Physical Review A 31, 1695 (1985).）、及び、Ｐａｒｒｉｎｅｌｌｏ－Ｒａｈｍａｎの圧力浴（M. Parrinello and A. Rahman, “Polymorphic transitions in single crystals: A new molecular dynamics method,” Journal of Applied Physics 52, 7182 (1981).）による温度、圧力一定のＮＰＴアンサンブルを用いて、３００Ｋ、１ａｔｍで５００ｐｓ保持して平衡化を行った。熱伝導度κの評価には、Ｇｒｅｅｎ－Ｋｕｂｏ （ＧＫ）公式を用いた。このとき、体積とエネルギーを一定にしたＮＶＥアンサンブルを用いて３００Ｋにおいて５００ｐｓの熱流の計算を行った。ＨＦＡＣＦの期待値計算のための積分時間を４６０ｐｓ、ＧＫ公式における積分時間を４０ｐｓとした。

【0026】

図２は、ポリエチレン（ＰＥ）に対して熱流（ＨＦ）の相関関数の計算結果を示す。図２（ａ）は、同一分子内の熱流の自己相関関数（ＡＣＦ）を示す。図２（ｂ）は、異なる分子間の熱流の相互相関関数（ＣＣＦ）を示す。なお、図２（ａ）及び図２（ｂ）において、１０^２ｆｓ～１０^４ｆｓの期間における計算結果を拡大して内挿図として示した。

【0027】

これらの和を時間に対して積分することで熱伝導度κが得られる。図２（ａ）から、評価開始時から１０^３ｆｓ以下の期間において、分子内の熱流の自己相関関数（ＡＣＦ）が比較的大きな値をとり、熱伝導度に大きく寄与することがわかった。一方、図２（ｂ）から、異なる分子間における熱流の相互相関関数（ＣＣＦ）は、評価開始時から１０^３ｆｓ以下の期間において分子内の熱流の自己相関関数（ＡＣＦ）よりも値が十分に小さく、１０^３ｆｓ以上の期間においてランダムな誤差を含むことがわかった。

【0028】

図２（ａ）及び図２（ｂ）から、熱流の相関関数を分子内の自己相関関数（ＡＣＦ）と分子間の相互相関関数（ＣＣＦ）に分けることで、熱伝導度に主に寄与する成分と誤差成分とを分離できる。本実施例では、熱流の相関関数を分子内の自己相関関数（ＡＣＦ）のみから熱伝導度を求める。これによって、分子間の相互相関関数（ＣＣＦ）の誤差の影響を避けることができる。一方、従来法では、複数の初期構造について計算を行い、平均をとる必要がある。

【0029】

図３は、ポリエチレン（ＰＥ）に対してＧｒｅｅｎ－Ｋｕｂｏ（ＧＫ）公式を用いて熱伝導度を計算した結果を示す。図３において、横軸はＧｒｅｅｎ－Ｋｕｂｏ（ＧＫ）公式の計算に用いる時間積分のカットオフ時間を示す。図３（ａ）は、ポリエチレン（ＰＥ）の２０個の初期配置について熱伝導度を計算した結果の平均値を示す。図３（ａ）において、太実線は熱伝導度に対する同一分子内の自己相関関数（ＡＣＦ）の寄与を示し、破線は熱伝導度に対する異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与を示す。また、図３（ａ）において、細実線は、同一分子内の自己相関関数（ＡＣＦ）の寄与と異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与との和を示し、これが従来法において計算された熱伝導度を示すことになる。

【0030】

図３（ｂ）は、ポリエチレン（ＰＥ）の２０個の初期配置の各々について熱伝導度の同一分子内の自己相関関数（ＡＣＦ）の寄与を計算した結果を示す。図３（ｃ）は、ポリエチレン（ＰＥ）の２０個の初期配置の各々について熱伝導度の異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与を計算した結果を示す。

【0031】

図３（ａ）において、ポリエチレン（ＰＥ）の熱伝導度（細線）の計算結果は２×１０^４ｆｓ以上の期間において増加した。これは、熱伝導度に対する異なる分子間の熱流の相互相関関数（ＣＣＦ）のランダムな誤差に起因する統計誤差に起因する。図３（ａ）の太実線で示すように、異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与を無視して、同一分子内の自己相関関数（ＡＣＦ）の寄与のみを用いた場合、１０^４ｆｓ以上の期間において熱伝導度は一定値に収束した。

【0032】

図４は、ポリスチレン（ＰＳ）に対してＧｒｅｅｎ－Ｋｕｂｏ（ＧＫ）公式を用いて熱伝導度を計算した結果を示す。図４において、横軸はＧｒｅｅｎ－Ｋｕｂｏ（ＧＫ）公式の計算に用いる時間積分のカットオフ時間を示す。図４（ａ）は、ポリスチレン（ＰＳ）の２０個の初期配置について熱伝導度を計算した結果の平均値を示す。図４（ａ）において、太実線は熱伝導度に対する同一分子内の自己相関関数（ＡＣＦ）の寄与を示し、破線は熱伝導度に対する異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与を示す。また、図４（ａ）において、細実線は、同一分子内の自己相関関数（ＡＣＦ）の寄与と異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与との和を示し、これが従来法において計算された熱伝導度を示すことになる。

【0033】

図４（ｂ）は、ポリスチレン（ＰＳ）の２０個の初期配置の各々について熱伝導度の同一分子内の自己相関関数（ＡＣＦ）の寄与を計算した結果を示す。図４（ｃ）は、ポリスチレン（ＰＳ）の２０個の初期配置の各々について熱伝導度の異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与を計算した結果を示す。

【0034】

図４（ａ）に示すように、ポリスチレン（ＰＳ）の熱伝導度の計算結果では、１０^４ｆｓ以上の期間において、太実線で示した熱伝導度に対する同一分子内の自己相関関数（ＡＣＦ）の寄与、及び、細実線で示した同一分子内の自己相関関数（ＡＣＦ）の寄与と異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与との和のいずれも一定値に収束した。すなわち、ポリスチレン（ＰＳ）の熱伝導度の計算では、異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与は十分小さく、無視できる程度であった。

【0035】

以上のように、異なる分子間の熱流の相互相関関数（ＣＣＦ）を無視することによって、そこに含まれるランダムな誤差が取り除かれ、複数の初期配置に対して平均を算出することなく、熱伝導度を高い精度で評価することができることが分かった。また、計算負荷を大幅に削減することができることが分かった。

【0036】

図５は、熱伝導度について典型的な直鎖高分子材料の単量体を繋ぎ合わせたときの分子の長さ依存性を示す。図５（ａ）は、比較例として従来法によって計算された熱伝導度の分子の長さ依存性を示す。図５（ｂ）は、実施例において異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与を無視して計算された熱伝導度の分子の長さ依存性を示す。図５において、横軸は計算に用いた分子の長さ（Å）を示す。また、図５において、丸印は２０個の分子の初期配置に対して計算した熱伝導度の平均値を示す。また、図５において、誤差バーは計算結果の標準偏差を示す。

【0037】

なお、図６は、計算に用いた典型的な直鎖高分子材料の単量体の分子構造を示す。図６において、星印は単量体の結合位置を表す。

【0038】

図５（ａ）に示すように、従来法では２０個の初期配置に対する個々の計算値の標準偏差（誤差棒）が大きくなった。したがって、複数の初期配置に対して計算を行い、平均値を算出することで熱伝導度の傾向が得られた（丸印）。ただし、平均値が滑らかな曲線には乗らなかったことから、図５（ａ）に示した結果には誤差が残っていると推察される。

【0039】

一方、図５（ｂ）に示すように、実施例の有機材料輸送特性評価方法を適用した２０個の初期配置に対する個々の計算値ではほとんどばらつきがなかった。したがって、１つの初期配置に対する個々の計算値でも十分な精度で熱伝導度を評価することができた。実施例では、異なる分子間の熱流の相互相関関数（ＣＣＦ）の寄与を無視したため、熱伝導度はわずかに過小評価されたものの、各直鎖高分子材料の熱伝導度の序列は従来法の結果とよく一致した。したがって、実施例の有機材料輸送特性評価方法においても計算自動探索には十分な精度であった。また、実施例では、標準偏差は従来法の１／１０程度に低減され、すなわち統計誤差が１／１００程度に低減された。

【0040】

［時間－応力ゆらぎ相関関数］
粘性ηは、Ｇｒｅｅｎ－Ｋｕｂｏ公式を用いて評価することができる。ある時点における系の応力の非対角成分ｐ^ａｂ（ａｂ＝｛ｘｙ，ｘｚ，ｙｘ｝）は数式（５）で表すことができる。ここで、Ｖは系の体積、Ｎは系の全原子数、ｍ_ｉはｉ番目の原子の質量、ｖ_ｉはｉ番目の原子の速度、ｒ_ｉｊはｉ番目の原子とｊ番目の原子の間の距離、ｆ_ｉｊはｉ番目の原子とｊ番目の原子の間の相互作用である。

【数5】

【0041】

この応力に対する時間相関関数を用いて数式（６）を定義する。ここで、ｋＢはボルツマン定数、Ｔは温度である。また、＜・・・＞は初期時間ｔ０を変えてサンプリングしたアンサンブル平均操作を意味する。

【数6】

【0042】

φ（ｔ）に対する時間積分は数式（７）と表すことができ、数式（８）にて粘性ηを求めることができる。

【数7】

【数8】

【0043】

数値計算上は数式（８）での時間ｔ→∞を扱うことはできないが、時間ｔを長く取ることによってΦ（ｔ）が有限値に収束することから、一定時間に対する計算で打ち切って粘性ηを評価することができる。

【0044】

しかしながら、収束が遅く、数値計算上の結果のばらつきも大きくなるため、複数のサンプリングによる統計的平均処理やモデル式へのフィッティングを行うことで精度を高める方法が従来から採用されている。

【0045】

本実施の形態における有機材料輸送特性評価装置１００では、せん断弾性定数Ｇ_∞が数式（７）を用いて数式（９）と対応付けられることを利用する。

【数9】

【0046】

数式（９）で表される量の統計誤差や計算負荷が大幅に小さいことを利用して、数式（１０）のアレニウス式を用いて粘性ηを算出する。ここで、Ｔ_ｂは系の沸点、α，α’はパラメータ、η_０は沸点における粘性である。

【数10】

【0047】

数式（５）は、原子毎の応力を示す数式（１１）を用いて数式（１２）と書き直すことができる。ここで、Ｎｍｏｌは系に含まれる分子数、ｐバー（ｐ上にバーを付した値）は分子毎の応力である。分子毎の応力ｐバーは、ある分子ｋを構成する原子の応力の和を意味する。

【数11】

【数12】

【0048】

数式（９）のせん断弾性定数Ｇ∞を求める際、数式（６）の相関関係部分は数式（１３）として表すことができる。なお、数式（１３）の近似式は、単一分子種の純溶媒のときに等号が成立する。

【数13】

【0049】

このように、全系の応力に関する時間相関を分子毎の応力の時間相関に変換することによって、初期時刻ｔ_０の違いだけでなく、分子数での統計的な平均をとることが可能となる。したがって、従来の方法に対して数値計算上のばらつきを抑制することができる。

【0050】

［粘度及び拡散係数の評価］
数式（１３）に基づいて、多成分系における粘度及び拡散係数を評価する方法について検討する。以下では、説明を簡単にするために成分Ａ及び成分Ｂからなる液体の２成分系について検討する。

【0051】

成分Ａ及び成分Ｂのモル比がｒ_Ａ：ｒ_Ｂ（ここでｒ_Ａ＋ｒ_Ｂ＝１）であるとすると、数式（１３）は数式（１４）及び数式（１５）で表される。なお、以下、初期時刻を示すｔ_０は省略する。また、数式（１６）の様に表記を省略した。

【数14】

【数15】

【数16】

【0052】

せん断弾性定数Ｇ_∞は数式（１７）及び数式（１８）で表される。数式（１８）に示すように、せん断弾性定数Ｇ_∞は、成分ｘのモル分率とせん断弾性定数Ｇ_∞，xの線形和で表すことができる。

【数17】

【数18】

【0053】

ここで、粘性予測のための最も簡潔なアレニウス式（数式（１９））を用いると、粘性ηは数式（２０）及び数式（２１）で表される。

【数19】

【数20】

【数21】

【0054】

これは、全系の粘性ηが、系を構成する成分毎の比率と粘性を用いて評価できることを示している。さらに、経験的な混合溶媒粘性の予測式（例えば、Viswanath, D. S. et al., “Viscosity of Liquids: Theory, Estimation, Experiment, and Data” (Springer, 2007).の式(5.28)）と一致する。

【0055】

粘性ηと拡散係数Ｄとの関係は、数式（２２）のＳｔｏｋｅｓ－Ｅｉｎｓｔｅｉｎの関係として表すことができる。ここで、Ｒは拡散粒子径、λは理論的には６（Ｓｔｉｃｋ条件）又は４（Ｓｌｉｐ条件）である。λの値は、拡散粒子と溶媒との関係で決定され、実験的には定数にならない場合もある。

【数22】

【0056】

数式（２２）で示されるように、拡散係数Ｄと粘性ηとは反比例関係にある。拡散係数Ｄは、さらに数式（２３）で表される。ここで、ｃは比例定数である。また、数式（２１）を用いた。

【数23】

【0057】

数式（２３）は、全系の拡散係数Ｄが系を構成する成分ｘの粘性ηｘと成分比率を用いて表されることを意味する。そこで、Ｍ成分系に一般化すると、ｘ成分の拡散係数Ｄｘは数式（２４）で表され、全系の拡散係数Ｄは数式（２５）で表される。

【数24】

【数25】

【0058】

［実施例２］
本実施例では、単成分液体にて、分子応力相関の和表式（数式（１３））が全応力相関（数式（６））と変わらないか確認した。

【0059】

分子動力学（ＭＤ）計算は、ＬＡＭＭＰＳプログラム（version 5 Jun 2019）（Plimpton, S., “Fast Parallel Algorithms for Short-Range Molecular Dynamics”, J. Comput. Phys. 117, 1-19 (1995), http://lammps.sandia.gov.）を用いて実施した。分子力場パラメータは、ＧＡＦＦ（Wang, J. et al., “Development and testing of a general amber force field”, J. Comput. Chem. 25, 1157-1174 (2004). doi:10.1002/jcc.20035.）を適用し、電荷はＲＥＳＰ法（Bayly, C. I. et al., “A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model”, J. Phys. Chem. 97, 10269-10280）でＡｎｔｅｃｈａｍｂｅｒプログラム（Case, D. A. et al., “AMBER 2018”, University of California, San Francisco, 2018.）を用いてアサインした。なおＲＥＳＰ電荷の算出には、Ｇａｕｓｓｉａｎ０９プログラム（Frisch, M. J. et al., “Gaussian 09, Revision E.01”, Gaussian Inc. Wallingford CT.）によるＨＦ／６－３１Ｇ／／Ｂ３ＬＹＰ／６－３１＋＋Ｇ＊＊レベルの結果を利用した。

【0060】

初期構造には、Ｐａｃｋｍｏｌプログラム（Martinez, L. et al., “PACKMOL: A package for building initial configurations for molecular dynamics simulations”, J. Comput. Chem. 30, 2157-2164 (2009). doi:10.1002/jcc.21224.）を用いて、構成分子群をランダム配置で充填した立方セルを用いた。セルの大きさは、約１ｇ／ｃｍ^３となるように設定した。セル中の分子数は、モル比に合わせて６００個含まれるようにした。

【0061】

分子動力学（ＭＤ）計算では、まずエネルギー最小化後、時間刻み０．２５ｆｓにて初期構造緩和（ＮＶＥとＮＶＴアンサンブル）を実施し、その後時間刻み１．０ｆｓ、１ａｔｍ下でのＮＰＴアンサンブルにて体積も緩和して平衡構造を作成した。得られた平衡構造から、本計算を時間刻み１．０ｆｓのＮＶＴアンサンブルにて１ｎｓ（１，０００，０００ｓｔｅｐ）行い、解析に用いた。温度は室温とした。長距離相互作用の計算にはＰＰＰＭ法（精度パラメータは１０^?４）を用い、相互作用のカットオフ値は１０Åとした。

【0062】

原子・分子ごとの応力算出には、ＬＡＭＭＰＳプログラムのｃｏｍｐｕｔｅ ｓｔｒｅｓｓ／ａｔｏｍコマンドと、ｃｏｍｐｕｔｅ ｒｅｄｕｃｅ／ｃｈｕｎｋコマンドを利用した。前者は原子ごとの応力を計算する機能、後者は設定したまとまり（ここでは分子単位）で合算して出力する機能である。

【0063】

図７は、ポリカーボネイト（ＰＣ）に対する応力に対する時間相関をプロットした結果を示す。図７（ａ）は短時間領域における計算結果を示し、図７（ｂ）は長時間領域における計算結果を示す。図７（ａ）及び図７（ｂ）において、濃実線は全応力相関を示し、薄実線は分子毎の応力相関（図では計算結果から－０．０５した値）を示す。また、図７（ｃ）は、図７（ｂ）における全応力相関と分子毎の応力相関の差を示す。

【0064】

図７（ａ）及び図７（ｂ）に示すように、短時間領域及び長時間領域のいずれにおいても全応力相関と分子応力相関和の振動挙動は全体的によく一致した。しかしながら、図７（ｃ）に示したように、全応力相関と分子応力相関和にはわずかに差が見られた。時間ｔ＝０のときの全応力相関と分子毎の応力相関の差異は－０．０００９［ｍＰａｓ／ｆｓ］であった。相関時間が進むにつれ、徐々に誤差が大きくなり、時間ｔ＝３０付近で最大となり、その後小さくなったが０近傍で一定の幅をもって振動が続いた。当該誤差は、実施例において用いた近似によるもの、すなわち分子間の相互相関を無視していることによるものと推察される。しかしながら、本実施例における有機材料輸送特性評価方法では、ごく短い相関時間（１０ｆｓ以下）の情報が重要であり、当該時間では誤差も相対的に小さいことから、拡散係数の評価に使用できると考えられる。

【0065】

［実施例３］
本実施例では、２成分溶液系における応力について評価した。図８は、ポリカーボネイト（ＰＣ）と１，２－ジメトキシエタン（ＤＭＥ）を５０％：５０％で混合した系に対する応力ゆらぎの時間相関を示す。図８（ａ）は短時間領域における計算結果を示し、図８（ｂ）は長時間領域における計算結果を示す。図８（ａ）及び図８（ｂ）において、濃実線は全応力相関を示し、薄実線は分子毎の応力相関（図では計算結果から－０．０５した値）を示す。さらに、図８（ａ）では、薄破線はポリカーボネイト（ＰＣ）のみの成分和、薄鎖線は１，２－ジメトキシエタン（ＤＭＥ）のみの成分和を示す。また、図８（ｃ）は、図８（ｂ）における全応力相関と分子毎の応力相関の差を示す。

【0066】

図７に示した単成分系での議論と同様に、２成分系でも分子応力相関の総和は、短時間領域及び長時間領域ともに全応力相関によく一致した。図８（ｃ）に示した両者の差では、時間ｔ＝０のときの全応力相関と分子毎の応力相関の差異は－０．００１２［ｍＰａｓ／ｆｓ］であった。両者の差は、大きさ及び振動幅において少しの差異は見られたものの、おおむね単成分系と同様の挙動を示した。

【0067】

構成成分毎の計算結果では、ポリカーボネイト（ＰＣ）と１，２－ジメトキシエタン（ＤＭＥ）の成分に差異があることが分かった（なお、時刻ｔ＝０での一致は偶然である）。ポリカーボネイト（ＰＣ）の成分の挙動は、モル比が５０％であるために単成分での応力ゆらぎにて振幅がちょうど半分になったことに対応する。すなわち、溶媒分子の種類は違ったとしても、個々の分子の寄与を足し合わせれば全系の応力相関を示すだけではなく、溶媒種類毎に足し合わせて成分寄与として評価できる。

【0068】

［実施例４］
本実施例では、多成分系での拡散係数Ｄの評価を検証した。多成分系にて成分毎の拡散係数Ｄは、ＮＭＲ測定法による実験が可能であり、報告例が多く存在する。本実施の形態において検証対象とした系は、ポリカーボネイト（ＰＣ）と１，２－ジメトキシエタン（ＤＭＥ）の２成分系及びポリカーボネイト（ＰＣ）とジエチルカーボネート（ＤＥＣ）の２成分系とした。

【0069】

成分毎の拡散係数Ｄを計算する際、従来法では分子間の応力相関も考慮した全応力相関を用い、実施例では分子間の応力相関は無視して分子毎の応力相関を用いた。従来法では、ＭＤ解析として一般的なＥｉｎｓｔｅｉｎ法（平均二乗変位の傾き）を用いた。なお、計算条件は、解析に用いたサンプリングステップ数を実施例より大幅に長く（２０倍）した以外、実施例と同等とした。

【0070】

図９（ａ）は、ポリカーボネイト（ＰＣ）と１，２－ジメトキシエタン（ＤＭＥ）の２成分系に対する成分毎の拡散係数Ｄを示す。図９（ｂ）は、ポリカーボネイト（ＰＣ）とジエチルカーボネート（ＤＥＣ）の２成分系に対する成分毎の拡散係数を示す。図９（ａ）及び図９（ｂ）では、ポリカーボネイト（ＰＣ）に対する１，２－ジメトキシエタン（ＤＭＥ）又はジエチルカーボネート（ＤＥＣ）の構成比率に対する成分毎の拡散係数Ｄを示す。黒四角印と太実線は実験値（早水紀久子, “3．リチウム電池用電解液の構造と輸送現象のNMR による研究”, Electrochemistry 81, 995-1000 (2013). doi:10.5796/electrochemistry.81.995.よるＮＭＲ測定値）、白三角印と鎖線は従来法による計算結果、白丸印と破線は本実施例による計算結果を示す。

【0071】

純溶媒では、拡散係数Ｄはポリカーボネイト（ＰＣ）＜１，２－ジメトキシエタン（ＤＭＥ）、ポリカーボネイト（ＰＣ）＜ジエチルカーボネート（ＤＥＣ）であることが知られている。ただし、図９に示したように、混合溶媒になると単体では低かったポリカーボネイト（ＰＣ）の拡散係数Ｄは高くなり、逆に単体では高かった１，２－ジメトキシエタン（ＤＭＥ）やジエチルカーボネート（ＤＥＣ）の拡散係数Ｄは低くなり、濃度とも相関することが報告されている。なお、溶媒の種類によって傾向は異なり、ポリカーボネイト（ＰＣ）－１，２－ジメトキシエタン（ＤＭＥ）系の場合は各成分の拡散係数Ｄに差が見られるが、ポリカーボネイト（ＰＣ）－ジエチルカーボネート（ＤＥＣ）系の場合は各成分の拡散係数が近い値を示す。

【0072】

図９（ａ）及び図９（ｂ）に示したように、従来法によって計算された各成分の拡散係数Ｄは、溶媒の濃度に応じた変化傾向はポリカーボネイト（ＰＣ）－１，２－ジメトキシエタン（ＤＭＥ）系、ポリカーボネイト（ＰＣ）－ジエチルカーボネート（ＤＥＣ）系ともに再現された。しかしながら、従来法によって計算された各成分の拡散係数Ｄは、実験値に対して定量的に大幅に過小評価された。

【0073】

一方、本実施例では、ポリカーボネイト（ＰＣ）を０％として１，２－ジメトキシエタン（ＤＭＥ）又はジエチルカーボネート（ＤＥＣ）の純溶媒とした場合、それぞれの成分の拡散係数Ｄは実験値と一致した。これに対して、ポリカーボネイト（ＰＣ）の拡散係数Ｄについては評価の精度が若干低下した。したがって、ポリカーボネイト（ＰＣ）－１，２－ジメトキシエタン（ＤＭＥ）系では、濃度変化に対して１，２－ジメトキシエタン（ＤＭＥ）の成分に対する拡散係数Ｄはよい一致をみせたが、ポリカーボネイト（ＰＣ）の成分に対する拡散係数Ｄは概ね過小評価となった。ポリカーボネイト（ＰＣ）－ジエチルカーボネート（ＤＥＣ）系では、前述のように２成分の溶媒の拡散係数Ｄが一致する傾向があり、ポリカーボネイト（ＰＣ）の成分及びジエチルカーボネート（ＤＥＣ）の成分ともに実験値に対して僅かに低い傾向を示した。しかしながら、従来法に比べて実験値に対する定量性は大幅に改善した。さらに、実施例では従来法に比べて計算コストは１／２０程度であるにもかかわらず、計算値のばらつきは１／５～１／１０の確度で算出できた。

【符号の説明】

【0074】

１０ 処理部、１２ 記憶部、１４ 入力部、１６ 出力部、１８ 通信部、１００ 有機材料輸送特性評価装置。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】