JP5469696B2 - 高分子材料のエネルギーロスの計算方法 - Google Patents

Description

本発明は、実際の試験を行うことなく高分子材料のエネルギーロスを安定的に精度良く計算しうる高分子材料のエネルギーロスの計算方法に関する。

例えば、プラスチック又はゴムのような高分子材料において、そのエネルギーロスは製品の様々な特性に影響を及ぼす重要な物理量である。例えば、ゴム製品であるタイヤでは、エネルギーロスは燃費性能やグリップ性能に密接に係わっており、適切な制御が必要と考えられている。

従来、例えば、目的とするゴム材料を得るために、ゴム材料の配合設計、試作及びエネルギーロスの測定試験といった工程を含む開発サイクルが繰り返される。しかしながら、従来の開発サイクルでは、ゴムのエネルギーロスについては、実際に材料を試作して引張試験等を行う必要があり、効率的ではなかった。

そこで、近年では、コンピュータを用いた高分子材料のシミュレーションを行い、その結果からエネルギーロスを計算することが提案されている（例えば、下記特許文献１乃至２参照）。これらの方法では、概ね、次のようなステップａ乃至ｅが時系列的に行われる。
ａ）解析対象となる高分子材料のシミュレーション用の分子構造モデルの設定
ｂ）分子構造モデルにポテンシャルを定義
ｃ）分子構造モデルをセルにランダムに配置した高分子材料モデルを用いて構造緩和計算
ｄ）構造緩和計算後、高分子材料モデルに周期的な変形を与えて応力と歪とを計算
ｅ）計算結果からヒステリシスループを求め、その面積からエネルギーロスを計算

特開２００７−１０７９６８号公報 特開２００９−２７１７１９号公報

しかしながら、従来の方法では、上記ステップｄの計算途中において、計算が不安定となって異常終了してしまうことがあった。また、ステップｅで計算されたエネルギーロスの値が、実験値等に比べて著しく異なる場合があった。

発明者らは、鋭意研究の結果、ポテンシャルの固有振動数に基づいた分子構造モデルの振動周期と、高分子材料モデルの変形の周期とが実質的に一致した場合、計算が不安定となって異常終了すること等を突き止めた。

本発明は、以上のような実情に鑑み、案出なされたもので、ポテンシャルの固有振動数に基づく分子構造モデルの振動周期と、高分子材料モデルの変形の周期とを比較し、両者が実質的に一致していないときに変形計算を行うことを基本として、高分子材料のエネルギーロスを精度良く計算しうる計算方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、コンピュータを用いて高分子材料のエネルギーロスを計算する方法であって、前記高分子材料の分子構造を、原子又はその集合体を表す有限個の粒子モデルと、この粒子モデル間の相対位置を特定する結合鎖とを含んでモデル化したシミュレーション用の分子構造モデルを設定するステップ、前記分子構造モデルの各粒子モデル間にポテンシャルを定義するステップ、前記分子構造モデルをセルに配置した高分子材料モデルを設定し、予め定めた条件と前記ポテンシャルに基づいて構造緩和計算を行うステップ、前記構造緩和計算の後、前記高分子材料モデルに周期的な変形を与えて応力と歪とを計算する変形計算のステップ、及び前記変形計算の結果に基づいて、エネルギーロスを計算するステップを含み、前記高分子材料モデルの前記変形の周期と、前記ポテンシャルの固有振動数に基づく振動周期とを比較し、両者が実質的に一致していないと判断されたときに前記変形計算を行うことを特徴とする。

また請求項２記載の発明は、前記構造緩和計算は、前記高分子材料モデルの密度が、実際の高分子材料の密度よりも小さく定義されて行われる請求項１記載の高分子材料のエネルギーロスの計算方法である。

また請求項３記載の発明は、前記構造緩和計算は、圧力及び温度が一定の下で少なくとも１０（ps）以上行われる請求項１又は２に記載の高分子材料のエネルギーロスの計算方法である。

また請求項４記載の発明は、前記変形計算は、慣性半径がセルの大きさの０．５倍未満の分子構造モデルのみを対象として行われる請求項１乃至３のいずれかに記載の高分子材料のエネルギーロスの計算方法である。

また請求項５記載の発明は、前記高分子材料モデルの前記変形の周期Ｔｓと、前記ポテンシャルの固有振動数に基づく振動周期Ｔｐとの比Ｔｓ／Ｔｐが０．８〜１．２の範囲にないと判断されたときに前記変形計算を行う請求項１乃至４のいずれかに記載の高分子材料のエネルギーロスの計算方法である。

本発明では、高分子材料の分子構造を、原子又はその集合体を表す有限個の粒子モデルと、この粒子モデル間の相対位置を特定する結合鎖とを含んでモデル化したシミュレーション用の分子構造モデルを設定するステップ、前記分子構造モデルの各粒子モデル間にポテンシャルを定義するステップ、前記分子構造モデルをセルに配置した高分子材料モデルを設定し、予め定めた条件と前記ポテンシャルに基づいて、構造緩和計算を行う構造緩和ステップ、前記構造緩和計算の後、前記高分子材料モデルに周期的な変形を与えて応力と歪とを計算する変形計算のステップ、及び前記変形計算の結果に基づいて、エネルギーロスを計算するステップを含む。

発明者らの実験によれば、ポテンシャルの固有振動数に基づいた分子構造モデルの振動周期と、高分子材料モデルの変形の周期とが一致した場合、計算が不安定となって異常終了する他、計算結果に大きな精度の低下が表れることが判明した。この理由については、明確にはなっていないが、概ね、高分子材料モデルが、分子構造モデルのポテンシャルの固有振動数から計算される振動周期と、同じ周期で変形した場合、変形を繰り返すにつれて原子間距離の揺らぎの振幅が大きくなり、原子に異常に大きな力がかかることが原因ではないかと推測される。

そこで、本発明では、高分子材料モデルの変形計算を行うステップに、高分子材料モデルの変形の周期と、ポテンシャルの固有振動数に基づく振動周期とを比較し、両者が実質的に一致していないときに前記変形計算を行うこととした。このため、本発明において、高分子材料モデルの変形計算は、常に、高分子材料モデルの変形の周期と、ポテンシャルの固有振動数に基づく振動周期とが異なった状態で行われるため、計算の異常終了や計算精度の低下を確実に防止できる。

ポリブタジエンの構造式である。 本実施形態のシミュレーション方法の処理手順を示すフローチャートである。 本実施形態の分子構造モデルの部分平面図である。 ３重合度の分子構造モデルの斜視図である。 高分子材料モデルを説明する斜視図である。 構造緩和計算の時間と分子構造モデルの密度との関係を示すグラフである。 本実施形態の変形計算のフローチャートである。 （ａ）は変形計算を説明するセルの線図、（ｂ）は歪の変化を説明するグラフである。 変形計算によって得られた応力と歪との関係を示すグラフである。 （ａ）、（ｂ）は、計算のセルと分子構造モデルの慣性半径との関係を説明する線図である。

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の高分子材料のエネルギーロスを計算する方法は、コンピュータを用いて高分子材料のエネルギーロスが計算される。

前記高分子材料としては、少なくともゴム、樹脂又はエラストマーを含む。本実施形態では、高分子材料として、図１に示されるように、cis-１，４ポリブタジエン（以下、単にポリブタジエンという。）が用いられる。ポリブタジエンは、メチレン基（−ＣＨ_２−）とメチン基（−ＣＨ−）とからなるモノマー｛−［ＣＨ_２−ＣＨ＝ＣＨ−ＣＨ_２］−｝を重合して得られる。ただし、ポリブタジエン以外の他の高分子材料についても同様に実施できるのは言うまでもない。

図２には、本実施形態のシミュレーション方法の処理手順の一例が示されている。本実施形態では、先ず、高分子材料としての分子構造モデルをコンピュータに設定する工程が行われる（ステップＳ１）。

図３及び図４には、分子構造モデル２の模式図が示されている。該分子構造モデル２は、前記ポリブタジエンを、コンピュータ上での分子動力学計算で取り扱えるように、設定される。従って、分子構造モデルの実態は、数値データである。また、図３及び図４に示されている分子構造モデル２は、解析対象のポリブタジエンの分子構造の高分子鎖の１単位を表している。

前記分子構造モデル２は、ポリブタジエンのモノマーであるブタジエンをモデル化したモノマーモデル３が、重合度ｎで連結されて構成される。なお、図４には、重合度が３の分子構造モデル２が示されている。

本実施形態のモノマーモデル３は、全ての原子（Ｃ及びＨ）がコンピュータのシミュレーションにおいて粒子モデル４として取り扱われる"Full Atom model"で構成される。ただし、例えば水素原子を陽として取り扱わない"United Atom Model"に従ってモデル化されたものでも良い。この場合、ＣＨ、ＣＨ２は、それぞれ水素原子の効果を織り込んだ原子の集合体として取り扱われる。

前記分子構造モデル２に含まれる粒子モデル４は、分子動力学計算に基づいたコンピュータシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル４には、質量、体積、直径、電荷及び／又は初期座標などのパラメータが定義される。これらの各パラメータは、上述のように、数値情報としてコンピュータに記憶される。

また、粒子モデル４、４間には、結合鎖５が定義される。結合鎖５は、粒子モデル４、４間を拘束するものである。結合鎖５には、例えば、平衡長とバネ定数とが定義されたバネとして取り扱われるもので、これらの情報もコンピュータに入力される。

個々のモノマーモデル３は、３次元構造を有し、各粒子モデル４、４間の結合長さである結合長、結合鎖５を介して連続する３つの粒子モデル４がなす角度である結合角、及び、結合鎖５を介して連続する４つの粒子モデル４において、隣り合う３つの粒子が作る二面角などがそれぞれ定義される。また、各モノマーモデル３は、慣例に従い、外力又は内力を受けることによって、上記結合長、結合角及び二面角が変化する。これにより、分子構造モデル２は、その三次元構造を変化させることができる。

以上のようなモデル化は、例えば（株）ＪＳＯＬ社製のＪ−ＯＣＴＡというソフトウエアを用いて処理することができる。

次に、本実施形態では、分子構造モデル２の各粒子モデル４、４間にポテンシャルが定義され、コンピュータに入力される（ステップＳ２）。

前記ポテンシャルは、粒子間の距離乃至角度の関数であって、２つの粒子モデル４、４間に作用する力を計算する際に用いられる。主なポテンシャルとしては、結合鎖５で連結された粒子モデル４、４間のボンドストレッチのポテンシャル、連続する３つの粒子モデル４で構成されるベンディングのポテンシャル、連続する３乃至４つの粒子モデル４で構成されるトーションのポテンシャル、及び、互いに連結されていない粒子モデル４、４間のファンデルワールス力のポテンシャルなどが定義される。これらのポテンシャルは、慣例に従い、距離又は結合角の関数として具体的に設定され、コンピュータに入力される。

次に、図５に示されるように、上記のように初期設定された分子構造モデル２を予め定めたセル（空間）Ｓにランダムに複数個配置した高分子材料モデルＰを設定し、予め定めた条件、前記ポテンシャル、及び分子動力学計算に基づいて、構造緩和計算が行われる（ステップＳ３）。このような構造緩和計算を行うことによって、セルＳ内の分子構造モデル２の人為的な配置による影響を無くして、その構造を安定化ないし準安定化させることができる。

前記セルＳは、解析対象の高分子材料の微小構造部分に相当する。本実施形態のセルＳは、微小な立方体として定義され、その中に３重合度の分子構造モデル２が複数本ランダムに初期配置されている。

前記構造緩和計算を行う条件等は、種々定めることができる。本実施形態では、特に好ましいものとして、高分子材料モデルＰの密度が、実際のポリブタジエンの密度よりも小さく設定されて構造緩和計算が行われる。前記密度は、一つのセルＳの体積と、該セル内に配置された分子構造モデル２の全質量との比で表される。

前記セルＳに多数本の分子構造モデル２を初期配置する場合、もし、分子構造モデル２同士が重なっていると、分子動力学計算において、異常な力が生じコンピュータでの計算が不安定になるおそれがある。しかし、本実施形態のように、分子構造モデル２の密度を、実際のポリブタジエンの密度よりも小さい条件を与えて構造緩和計算を開始することにより、分子構造モデル２同士が重なる確率を減じ、計算を安定させることができる。本実施形態では、初期配置において、分子構造モデルの密度が０．０１ｇ／cm³で設定されている。

また、前記構造緩和計算は、圧力一定、温度一定、及び周期境界条件の下、少なくとも１０ps（ピコ秒）以上行われることが望ましい。図６には、上記条件にて構造緩和計算を行った結果、分子構造モデル２の密度と時間との関係が示されている。構造緩和計算を少なくとも１０ps以上行うことにより、高分子材料モデルの密度が実際の値に収束、安定する。従って、次のステップで行われる高分子材料モデルＰの変形計算が安定する。

次に、前記構造緩和計算を終え安定した高分子材料モデルＰに、周期的な変形を与えて応力と歪とを計算する変形計算が行われる（ステップＳ５）。

変形計算ステップの処理手順の一例は、図７に示されている。本実施形態では、先ず、コンピュータが、分子構造モデル２に定義されたポテンシャルの固有振動数に基づく振動周期Ｔｐを計算する（ステップＳ５１）。本実施形態のように、ポテンシャルが複数定義されている場合、このステップＳ５１では、それぞれのポテンシャルに基づいて、それぞれの振動周期Ｔｐが計算される。

例えば、いま、Ｃ−Ｃ結合間のボンドスットレッチのポテンシャルＵが、下式（１）で定義されているとする。
Ｕ＝０．５＊ｋ＊（ｘ−ｘ0）^２ …（１）
ここで、"ｋ"は、バネ定数、"ｘ"はＣ−Ｃ結合間のボンド長、"ｘ0"は結合鎖の平衡長である。

この場合、前記ボンドスットレッチに関する固有振動数ｆは、下式（２）で計算することができる。
ｆ＝１／（２π）＊√（ｋ／Ｍ） …（２）
ここで、"π"は円周率、"Ｍ"はＣ原子の質量である。

従って、ボンドスットレッチのポテンシャルの固有振動数ｆに基づく振動周期Ｔｐは、下式（３）で計算される。
Ｔｐ＝１／ｆ …（３）

次に、コンピュータは、これから行おうとする高分子材料モデルＰの変形計算の周期Ｔｓを設定する（ステップＳ５２）。本実施形態の変形計算では、図８（ａ）に示されるように、高分子材料モデルＰに、初期の状態から１軸の引張変形を与えた後、それと逆方向に同一の歪で圧縮変形を与え、再び、初期の状態へと戻される工程を１周期とする引張、圧縮変形をシミュレーションする。図８（ｂ）には、このような高分子材料モデルＰの周期的な変形を表すものとして、歪と時間との関係が示されており、本実施形態では、歪が正弦波で変化するが、三角波のように変化するものでも良い。また、周期Ｔｓの他、歪の大きさやポアソン比なども設定される。

次に、コンピュータは、ポテンシャルの固有振動数に基づく振動周期Ｔｐと、高分子材料モデルＰに定義された変形の周期Ｔｓとを比較し、両者が実質的に一致しているか否かを判断する（ステップＳ５３）。

そして、ポテンシャルの固有振動数に基づく振動周期Ｔｐと、高分子材料モデルＰに定義された変形の周期Ｔｓとが実質的に一致していると判断された場合（ステップＳ５３でＹ）、高分子材料モデルＰの変形計算の周期Ｔｓが再度、設定される（ステップＳ５２）。

他方、ポテンシャルの固有振動数に基づく振動周期Ｔｐと、高分子材料モデルＰに定義された変形の周期Ｔｓとが実質的に一致していないと判断された場合（ステップＳ５３でＮ）、高分子材料モデルＰについて、周期Ｔｓで変形計算が行われる（ステップＳ５４）。なお、ステップＳ５３では、ポテンシャルの固有振動数に基づく振動周期Ｔｐが複数ある場合、これらの全てについて、実質的に一致するか否かが判断される。

発明者らの実験によれば、分子構造モデル２のポテンシャルの固有振動数に基づいた振動周期Ｔｐと、高分子材料モデルＰの変形計算の周期Ｔｓとが実質的に一致した場合、計算が不安定となり計算が異常終了する他、計算結果に大きな精度低下が表れることが判明した。この理由については、さらなる研究が必要であるが、高分子材料モデルＰが、分子構造モデル２のポテンシャルの固有振動数から計算される振動周期Ｔｐと、同じ周期Ｔｓで変形した場合、変形を繰り返すにつれて原子間距離の揺らぎの振幅が大きくなり、原子に異常に大きな力が作用することが原因と推測されている。なお、高分子材料モデルの変形計算を繰り返さず、１周期のみ行うようにすれば、上記不具合は解消されるようにも見える。しかしながら、分子動力学を用いて周期変形を行う場合、変形開始後しばらくの間は、変形による影響が応力に現れないか、又は現れにくい状態になる。これは、計算系が変形に馴染む過渡期の様な状態である。従って、変形の周期Tpが固有振動数Tsに一致しているか否かに係わらず、変形計算開始後に変形の影響が見えない時間がある。このため、変形を１周期だけに限って計算することは現実的ではない。

そこで、本発明では、高分子材料モデルＰの変形計算を行うに先立ち、高分子材料モデルＰの変形の周期Ｔｓと、分子構造モデル２のポテンシャルの固有振動数に基づく振動周期Ｔｐとを比較し、両者が実質的に一致していないことを条件として、高分子材料モデルＰの変形計算が行われる。このため、本発明で行われる高分子材料モデルＰの変形計算では、常に、高分子材料モデルの前記変形の周期Ｔｓと、ポテンシャルの固有振動数に基づく振動周期Ｔｐとが本質的に異なったものとなり、コンピュータでの計算の異常終了や計算精度の低下などの不具合が確実に防止される。

なお、本実施形態において、前記ステップ５３の判断、即ち、各周期Ｔｓ及びＴｐが実質的に一致するか否かの判断は、記高分子材料モデルＰの変形の周期Ｔｓと、ポテンシャルの固有振動数に基づく振動周期Ｔｐとの比Ｔｓ／Ｔｐが０．８〜１．２の範囲にある場合には、実質的に一致すると判断し、この範囲にない場合には、実質的に一致しないと判断している。ただし、このような態様に限定されるものではない。

図９には、このような変形計算により得られた高分子材料モデルＰの応力−歪曲線の１ループが示されている。本発明では、コンピュータによる計算の異常終了なしに、このような応力−歪曲線を得ることができ、かつ、その精度も高いものとなる。そして、高分子材料モデルＰのエネルギーロスを求めるには、応力−歪曲線が描くヒステリシスループの面積を計算すれば良い。

前記変形計算の好ましい態様して、前記変形計算は、慣性半径がセルＳの大きさの０．５倍未満の分子構造モデルのみを対象として行われるのが良い。

図１０（ａ）には、慣性半径ＲｇがセルＳの大きさＬ（この例では一つのセルＳの一辺の長さ）の０．５倍以上の場合が示されている。分子構造モデルを用いた変形計算では、セルＳが前後左右上下に連続するいわゆる周期境界条件が与えられる。周期境界条件では、小さな系である一つのセルＳの一方の境界が、そのセルの反対側の境界と繋がっているものとして取り扱われる。この条件の下では、分子構造モデル２がセルＳの一方の境界Ｓａから外側にはみ出した部分は、セルＳの他方の境界ＳｂからセルＳ内に進入する。この際、慣性半径ＲｇがセルＳの大きさＬの０．５倍以上の分子構造モデル２では、同一の分子構造モデルの粒子モデル同士が重なることがあり、ポテンシャルによる異常な力が働く場合がある。一方、図１０（ｂ）に示されるように、慣性半径ＲｇがセルＳの大きさＬの大きさの０．５倍未満の分子構造モデル２の場合、上述のような干渉がなく、精度の良い計算が行える。

本発明の効果を確認するために、上述の実施形態に準じて、コンピュータを使用したエネルギーロスの計算が行われた。結合長のボンドストレッチのポテンシャルのばね定数ｋ＝９１８［ε／σ²］、粒子モデルの質量Ｍ＝１．０［ｍ］、固有振動数ｆ＝４．８２［１／τ］、振動周期は０．２０７［τ］となる（ε、σ、ｍ及びτは、シミュレーションにおける単位エネルギー、単位長さ、単位質量及び単位時間を表している。）。

高分子材料モデルの変形の周期Ｔｓを０．２０７［τ］とし、ポテンシャルの固有振動数に基づく振動周期と一致させて計算を行ったところ（比較例）、コンピュータの計算が、途中で異常終了した。

一方、本発明に従い実施例では、上記条件では、高分子材料モデルの変形計算は行われることがない。また、高分子材料モデルの変形の周期Ｔｓを、ポテンシャルの固有振動数に基づく振動周期Ｔｐの１０倍、即ち２．０７［τ］として計算を行ったところ（実施例）、図９に示したヒステリシスループで表されるような応力と歪との関係を得ることができた。

また、高分子材料モデルの変形の周期Ｔｓを種々異ならせて計算の完了の有無を確認したところ、周期の比Ｔｓ／Ｔｐが０．８〜１．２の範囲において、計算の異常終了することがあった。

２ 分子構造モデル
３ モノマーモデル
４ 粒子モデル
５ 結合鎖
Ｐ 高分子材料モデル
Ｓ セル

Claims (5)

1. コンピュータを用いて高分子材料のエネルギーロスを計算する方法であって、
   前記高分子材料の分子構造を、原子又はその集合体を表す有限個の粒子モデルと、この粒子モデル間の相対位置を特定する結合鎖とを含んでモデル化したシミュレーション用の分子構造モデルを設定するステップ、
   前記分子構造モデルの各粒子モデル間にポテンシャルを定義するステップ、
   前記分子構造モデルをセルに配置した高分子材料モデルを設定し、予め定めた条件と前記ポテンシャルに基づいて構造緩和計算を行うステップ、
   前記構造緩和計算の後、前記高分子材料モデルに周期的な変形を与えて応力と歪とを計算する変形計算のステップ、及び
   前記変形計算の結果に基づいて、エネルギーロスを計算するステップを含み、
   前記高分子材料モデルの前記変形の周期と、前記ポテンシャルの固有振動数に基づく振動周期とを比較し、両者が実質的に一致していないと判断されたときに前記変形計算を行うことを特徴とする高分子材料のエネルギーロスの計算方法。
2. 前記構造緩和計算は、前記高分子材料モデルの密度が、実際の高分子材料の密度よりも小さく定義されて行われる請求項１記載の高分子材料のエネルギーロスの計算方法。
3. 前記構造緩和計算は、圧力及び温度が一定の下で少なくとも１０（ps）以上行われる請求項１又は２に記載の高分子材料のエネルギーロスの計算方法。
4. 前記変形計算は、慣性半径がセルの大きさの０．５倍未満の分子構造モデルのみを対象として行われる請求項１乃至３のいずれかに記載の高分子材料のエネルギーロスの計算方法。
5. 前記高分子材料モデルの前記変形の周期Ｔｓと、前記ポテンシャルの固有振動数に基づく振動周期Ｔｐとの比Ｔｓ／Ｔｐが０．８〜１．２の範囲にないと判断されたときに前記変形計算を行う請求項１乃至４のいずれかに記載の高分子材料のエネルギーロスの計算方法。

Images