JP7059687B2

JP7059687B2 - 輸送係数の解析プログラムおよび輸送係数の解析方法

Info

Publication number: JP7059687B2
Application number: JP2018031932A
Authority: JP
Inventors: 穂高森; 龍介泉; 伸幸松林
Original assignee: Denso Corp; Osaka University NUC
Current assignee: Denso Corp; Osaka University NUC
Priority date: 2018-02-26
Filing date: 2018-02-26
Publication date: 2022-04-26
Anticipated expiration: 2038-02-26
Also published as: JP2019148446A

Description

本発明は、輸送係数の解析プログラムおよび輸送係数の解析方法に関するものである。

従来、緩和弾性率、粘度、電気伝導度、熱伝導度、拡散係数等の輸送係数の解析については、測定による方法の他、シミュレーションによる方法が提案されている。例えば非特許文献１では、イオン間の距離に基づいて電気伝導度を解析する方法が提案されている。

K.-M. Tu, R. Ishizuka, N. Matubayasi, "Spatial-decomposition analysis of electrical conductivity in concentrated electrolyte solution", The Journal of Chemical Physics 141, 044126 (2014)

接着面等の界面、特にナノオーダーの凹凸が形成された界面の特性が強調されるモデルでは、輸送係数の空間不均一性が系全体の特性に大きく影響すると予想され、このようなモデルの接着、電気伝導、熱伝導等の現象の解析は、産業上重要である。また、接着面、フィラー周りの電気特性、熱特性などの解析は、不具合の解明や新たな接着処理、フィラー処理等の実現、新規材料の創出等に役立つと考えられる。

しかしながら、非特許文献１に記載の方法では、このような不均一系の輸送係数を解析することができない。

本発明は上記点に鑑みて、不均一系の輸送係数を解析することが可能な輸送係数の解析プログラムおよび輸送係数の解析方法を提供することを目的とする。

上記目的を達成するため、請求項１に記載の発明では、輸送係数の解析プログラムであって、解析対象の系を複数の領域に分割する分割部（Ｓ２）と、系に含まれる原子の運動を計算する運動計算部（Ｓ３）と、運動計算部の計算結果に基づいて、複数の領域それぞれについて輸送係数に対応する物理量を算出する物理量算出部（Ｓ４）と、物理量の時間相関関数を算出する時間相関関数算出部（Ｓ５）と、時間相関関数に基づいて、複数の領域それぞれについて輸送係数を算出する輸送係数算出部（Ｓ６）と、を備え、複数の領域の数をＮ _ｓｕｂとし、時間をｔとし、自然数をｊ、iとし、複数の領域のうちｊ番目の領域の時間における輸送係数をＧ _ｊ（ｔ）とし、ｊ番目の領域の時間における応力テンソルのａｂ成分をσ _ｊ ^ａｂ（ｔ）とし、ｉ番目の領域の時間における応力テンソルのａｂ成分をσ _i ^ａｂ（ｔ）とし、輸送係数算出部は、下記式を用いてＧ _ｊ（ｔ）を算出し、算出したＧ _ｊ（ｔ）に基づいて、輸送係数を算出する。

このように、系を複数の領域に分割し、各領域について輸送係数に対応する物理量を算出することにより、不均一系の輸送係数を解析することができる。

また、請求項５に記載の発明では、輸送係数の解析方法であって、解析対象の系を複数の領域に分割すること（Ｓ２）と、系に含まれる原子の運動を計算すること（Ｓ３）と、原子の運動を計算することの結果に基づいて、複数の領域それぞれについて輸送係数に対応する物理量を算出すること（Ｓ４）と、物理量の時間相関関数を算出すること（Ｓ５）と、時間相関関数に基づいて、複数の領域それぞれについて輸送係数を算出すること（Ｓ６）と、を備え、複数の領域の数をＮ _ｓｕｂとし、時間をｔとし、自然数をｊ、iとし、複数の領域のうちｊ番目の領域の時間における輸送係数をＧ _ｊ（ｔ）とし、ｊ番目の領域の時間における応力テンソルのａｂ成分をσ _ｊ ^ａｂ（ｔ）とし、ｉ番目の領域の時間における応力テンソルのａｂ成分をσ _i ^ａｂ（ｔ）とし、下記式を用いてＧ _ｊ（ｔ）を算出し、算出したＧ _ｊ（ｔ）に基づいて、輸送係数を算出する。

なお、各構成要素等に付された括弧付きの参照符号は、その構成要素等と後述する実施形態に記載の具体的な構成要素等との対応関係の一例を示すものである。

第１実施形態が適用される不均一系の概略構成を示す図である。図１のII部分の拡大図である。輸送係数の解析処理のフローチャートである。輸送係数の解析結果を示す図である。

以下、本発明の実施形態について図に基づいて説明する。なお、以下の各実施形態相互において、互いに同一もしくは均等である部分には、同一符号を付して説明を行う。

（第１実施形態）
第１実施形態について説明する。ここでは、図１に示す系について、輸送係数の１つである緩和弾性率を算出する場合について説明する。

図１に示す系は、金属層１の表面に接着剤２が塗布され、接着剤２に対して金属層１とは反対側の領域が真空とされたものである。金属層１は、ここではアルミニウムで構成されており、接着剤２は、ここでは炭素、酸素、水素で構成された樹脂とされている。

図２に示すように、金属層１と接着剤２との界面付近は空間的に不均一な構成とされており、このような構成の系については、従来の方法では輸送係数の解析が困難である。なお、接着剤２のうち金属層１との界面から離れた部分は均一な構成とされているが、この部分の輸送係数にも、金属層１との界面付近の特性が影響すると考えられる。

本実施形態では、図３に示すステップＳ１～Ｓ６を順に実行することで、図１に示す系の緩和弾性率を算出し、界面付近の特性、および、界面から離れた部分の特性の解析を可能とする。なお、図３に示す処理は、例えばＬＡＭＭＰＳ、Ｐｙｔｈｏｎ等を用いた解析プログラムによって実行することができるが、他のソフトウェアおよびプログラミング言語を用いて解析を行ってもよい。

ステップＳ１では、解析対象の系のデータを取得する。このデータには、図２に示すような金属層１、接着剤２を構成する複数の原子の種類および初期位置の他、原子の初速度等が含まれている。なお、ステップＳ１では、このようなデータを解析プログラムの外部から取得してもよいし、解析プログラム内で設定してもよい。

ステップＳ２では、系を複数の領域に分割する。例えば、図１の紙面左右方向、奥行き方向、上下方向をｘ方向、ｙ方向、ｚ方向として、系をｘ方向およびｙ方向に２０分割し、ｚ方向に１００分割する。

ステップＳ３では、ＭＤ（分子動力学）法等に基づいて系に含まれる複数の原子の運動を計算し、計算の開始から終了までの複数の時刻における各原子の位置および速度を取得する。

ステップＳ４では、各時刻における各原子の位置および速度に基づいて、各時刻における各領域の輸送係数に対応する物理量を算出する。本実施形態では、各領域の応力テンソルを算出し、緩和弾性率に対応する物理量としてせん断応力を用いる。

ステップＳ５では、せん断応力のＴＣＦ（時間相関関数）を算出する。そして、ステップＳ６では、ステップＳ５で求めたＴＣＦに基づいて各領域の緩和弾性率を算出する。ステップＳ５およびステップＳ６では、具体的には、次のようにして緩和弾性率を求める。

すなわち、複数の領域の数をＮ_ｓｕｂとし、時間をｔとし、複数の領域のうちｊ番目の領域の時間ｔにおける緩和弾性率をＧ_ｊ（ｔ）とする。また、ｊ番目の領域の時間ｔにおける応力テンソルのａｂ成分をσ_ｊ ^ａｂ（ｔ）とし、ｋ番目の領域の時間ｔにおける応力テンソルのａｂ成分をσ_ｋ ^ａｂ（ｔ）とし、数式１を用いてＧ_ｊ（ｔ）を求める。応力テンソルのａｂ成分は、例えば応力テンソルのｘｙ成分、ｙｚ成分、ｚｘ成分等である。なお、ＴＣＦの算出にはmulti－tau相関法などのスムージング手法を用いてもよい。

ステップＳ５では、ステップＳ４で算出された各時刻の応力テンソルについて順次計算を行ってもよいが、処理効率の向上のために、各時刻の応力テンソルについて並列処理を行うことが望ましい。

例えば、ステップＳ４において、応力テンソルを１ｆｓ毎に算出し、算出結果のデータを１００ｐｓ毎に出力する。そして、ステップＳ５にて、１００ｐｓ毎に分割された複数のデータを取得し、複数のデータそれぞれについて並列にＴＣＦを算出する。なお、ＴＣＦを算出する際には、遅れ時間に応じて、計算対象としているデータに、このデータに続く次の１００ｐｓのデータを追加して用いる。そして、１００ｐｓ毎のＴＣＦの計算結果を平均化する。

なお、Ｇ_ｊ（ｔ）を用いて、数式２により、時間ｔにおける系全体の緩和弾性率Ｇ（ｔ）を算出することができる。ここで、Ｖは系全体の体積、ｋはボルツマン定数である。数式１、数式２の導出方法については後述する。

このようにして、系の輸送係数を求めることができる。なお、ステップＳ２は分割部に対応し、ステップＳ３は運動計算部に対応し、ステップＳ４は物理量算出部に対応し、ステップＳ５は時間相関関数算出部に対応し、ステップＳ６は輸送係数算出部に対応する。

図３に示す処理によって、例えば図４に示すような結果が得られる。なお、図４において、黒丸は、金属層１と接着剤２との界面付近の輸送係数の算出結果であり、実線は、この算出結果の近似曲線である。また、黒四角は、接着剤２のうち金属層１との界面から離れた部分の輸送係数の算出結果であり、破線は、この算出結果の近似曲線である。

図４に示す結果では、金属層１と接着剤２との界面付近が、接着剤２のうち金属層１との界面から離れた部分に比べて、Ｇ_ｊ（ｔ）の減衰が遅く、弾性的になっている。金属層１の表面処理の方法や、接着剤２を構成する樹脂の種類等を様々に変化させて図３に示す処理を実行し、得られた結果を比較することで、用途に適した緩和弾性率が得られる材料等の条件を導き出すことができる。

数式１、数式２の導出方法について説明する。均一な系の緩和弾性率は、数式３で表され、せん断応力のＴＣＦに対応することが知られている。

また、原子にポテンシャルＵが働く場合、マクロな応力は、各原子からの寄与の和として、数式４のように表されることが知られている。ｖ_ｉ１は、系に含まれる複数の原子のうちｉ_１番目の原子の速度であり、ｒ_ｉ１ｉ２はｉ_１番目の原子とｉ_２番目の原子との間の距離である。

ここで、系全体の応力のうちｉ_１番目の原子からの寄与は、local atomic stress elementsと呼ばれる以下の量ｓ_ｉ１ ^ａｂで表される。

ｓ_ｉ１ ^ａｂを用いると、系全体の応力は数式６で表される。なお、ｓ_ｉ１ ^ａｂとしてはＬＡＭＭＰＳ等において具体的に定義されたものを用いればよい。

系を複数の領域に分割し、各領域についてｓ_ｉ１ ^ａｂの和を求めるようにすると、数式６は数式７のように変形される。Ｎ_ａはｊ番目の領域に含まれる原子の個数、ｓ_ｉ１，ｊ ^ａｂはｊ番目の領域に含まれる複数の原子のうちｉ_１番目の原子のlocal atomic stress elementsである。

ここで、ｊ番目の領域の応力をσ_ｊ ^ａｂとすると、σ_ｊ ^ａｂは、数式８のようにｓ_ｉ１、ｊ ^ａｂの和で表されるので、系全体の応力σ^ａｂは、σ_ｊ ^ａｂを用いて数式９のように表される。

数式９を数式３に代入して変形すると、緩和弾性率Ｇ（ｔ）は、数式１０のように表される。

このように、系全体の緩和弾性率Ｇ（ｔ）は、領域ごとの相関関数の和で表される。そこで、数式１０のうち相関関数で表された部分を各領域の緩和弾性率Ｇ_ｊ（ｔ）とし、数式１のように定義すれば、系全体の緩和弾性率Ｇ（ｔ）は、数式２に示すようになる。

以上説明したように、本実施形態では、系を複数の領域に分割し、各領域について輸送係数に対応する物理量を算出することで、各領域の輸送係数を算出することができる。これにより、不均一系について、界面付近の特性、および、界面から離れた部分の特性を解析することが可能となる。

（他の実施形態）
なお、本発明は上記した実施形態に限定されるものではなく、特許請求の範囲に記載した範囲内において適宜変更が可能である。

例えば、上記第１実施形態では緩和弾性率の解析方法について説明したが、粘度、電気伝導度、熱伝導度、拡散係数等の他の輸送係数についても、第１実施形態と同様に、ステップＳ１～Ｓ６の処理により解析することができる。

すなわち、系全体の物理量を数式６のように各原子からの寄与の和で表すことができれば、寄与の和の部分を数式７のように領域ごとの和に分割して、系全体の物理量を数式９のように各領域の物理量の和を用いて表すことができる。また、輸送係数は、数式３と同様に、対応する物理量のＴＣＦを用いて表されるので、数式９に相当する式を、輸送係数を表す式のＴＣＦの部分に代入して変形すれば、各領域の物理量のＴＣＦの和を用いて輸送係数を表すことができる。そして、各領域についてのＴＣＦの部分を各領域の輸送係数とすれば、第１実施形態と同様に不均一系の輸送係数を解析することができる。

Ｓ２分割部
Ｓ３運動計算部
Ｓ４物理量算出部
Ｓ５時間相関関数算出部
Ｓ６輸送係数算出部

Claims

輸送係数の解析プログラムであって、
解析対象の系を複数の領域に分割する分割部（Ｓ２）と、
前記系に含まれる原子の運動を計算する運動計算部（Ｓ３）と、
前記運動計算部の計算結果に基づいて、前記複数の領域それぞれについて前記輸送係数に対応する物理量を算出する物理量算出部（Ｓ４）と、
前記物理量の時間相関関数を算出する時間相関関数算出部（Ｓ５）と、
前記時間相関関数に基づいて、前記複数の領域それぞれについて前記輸送係数を算出する輸送係数算出部（Ｓ６）と、を備え、
前記複数の領域の数をＮ _ｓｕｂとし、
時間をｔとし、
自然数をｊ、iとし、
前記複数の領域のうちｊ番目の前記領域の時間における前記輸送係数をＧ _ｊ（ｔ）とし、
ｊ番目の前記領域の時間における応力テンソルのａｂ成分をσ _ｊ ^ａｂ（ｔ）とし、
ｉ番目の前記領域の時間における応力テンソルのａｂ成分をσ _i ^ａｂ（ｔ）とし、
前記輸送係数算出部は、

の式を用いて前記Ｇ _ｊ（ｔ）を算出し、算出した前記Ｇ _ｊ（ｔ）に基づいて、前記輸送係数を算出する輸送係数の解析プログラム。
時間における前記系の全体の前記輸送係数をＧ（ｔ）とし、
前記系の全体の体積をＶとし、
ボルツマン定数をｋとし、
前記輸送係数算出部は、前記Ｇ _ｊ（ｔ）を

の式に代入することにより、前記輸送係数を算出する請求項１に記載の輸送係数の解析プログラム。
前記輸送係数は、緩和弾性率である請求項１または２に記載の輸送係数の解析プログラム。
前記時間相関関数算出部は、並列処理を用いて前記時間相関関数を算出する請求項１ないし３のいずれか１つに記載の輸送係数の解析プログラム。
輸送係数の解析方法であって、
解析対象の系を複数の領域に分割すること（Ｓ２）と、
前記系に含まれる原子の運動を計算すること（Ｓ３）と、
前記原子の運動を計算することの結果に基づいて、前記複数の領域それぞれについて前記輸送係数に対応する物理量を算出すること（Ｓ４）と、
前記物理量の時間相関関数を算出すること（Ｓ５）と、
前記時間相関関数に基づいて、前記複数の領域それぞれについて前記輸送係数を算出すること（Ｓ６）と、を備え、
前記複数の領域の数をＮ _ｓｕｂとし、
時間をｔとし、
自然数をｊ、iとし、
前記複数の領域のうちｊ番目の前記領域の時間における前記輸送係数をＧ _ｊ（ｔ）とし、
ｊ番目の前記領域の時間における応力テンソルのａｂ成分をσ _ｊ ^ａｂ（ｔ）とし、
ｉ番目の前記領域の時間における応力テンソルのａｂ成分をσ _i ^ａｂ（ｔ）とし、

の式を用いて前記Ｇ _ｊ（ｔ）を算出し、算出した前記Ｇ _ｊ（ｔ）に基づいて、前記輸送係数を算出する輸送係数の解析方法。
時間における前記系の全体の前記輸送係数をＧ（ｔ）とし、
前記系の全体の体積をＶとし、
ボルツマン定数をｋとし、
前記Ｇ _ｊ（ｔ）を

の式に代入することにより、前記輸送係数を算出する請求項５に記載の輸送係数の解析方法。
前記輸送係数は、緩和弾性率である請求項５または６に記載の輸送係数の解析方法。
前記時間相関関数を算出することでは、並列処理を用いて前記時間相関関数を算出する請求項５ないし７のいずれか１つに記載の輸送係数の解析方法。