(57)【要約】【課題】  シミュレーション精度を向上させる。【解決手段】  コンピュータを用いて、フィラーが配合されたゴム材料の伸長変形を解析するための方法である。このシミュレーション方法は、散逸粒子動力学法に基づいて、フィラーを複数の第１粒子１１で表現したフィラーモデル６を入力する工程、散逸粒子動力学法に基づいて、ゴム材料の分子鎖を複数の第２粒子１２で表現したポリマーモデル７を入力する工程、第１粒子１１と第２粒子１２との間に、Lenard-JonesポテンシャルＵLJ（ｒij）を定義する工程、コンピュータが、仮想空間に配置されたフィラーモデル６及びポリマーモデル７を用いて分子動力学計算を行うことにより、構造緩和されたゴム材料モデルを生成する工程、及び、ゴム材料モデルの伸長変形を計算する工程を含む。【選択図】図７
2020085692.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを用いて、フィラーが配合されたゴム材料の伸長変形を解析するための方法であって、
散逸粒子動力学法に基づいて、前記フィラーを複数の第１粒子で表現したフィラーモデルを、前記コンピュータに入力する工程、
前記散逸粒子動力学法に基づいて、前記ゴム材料の分子鎖を複数の第２粒子で表現したポリマーモデルを、前記コンピュータに入力する工程、
前記第１粒子と前記第２粒子との間に、Lenard-Jonesポテンシャルを定義する工程、
前記コンピュータが、前記コンピュータの仮想空間に配置された前記フィラーモデル及び前記ポリマーモデルを用いて分子動力学計算を行うことにより、構造緩和されたゴム材料モデルを生成する工程、及び
前記コンピュータが、前記ゴム材料モデルの伸長変形を計算する工程を含むゴム材料のシミュレーション方法。
【請求項２】
隣接する前記第１粒子間に、結合角ポテンシャルを定義する工程をさらに含む、請求項１記載のゴム材料のシミュレーション方法。
【請求項３】
前記結合角ポテンシャルは、下記式（１）で定義される、請求項２記載のゴム材料のシミュレーション方法。
【数１】
2020085692.tif 000009
ここで、
Ｕ（θ）：結合角ポテンシャル
θ：隣接する３つの第１粒子間がなす結合角の外角
ｋ：フィラーのバネ定数
θ0：平衡状態の隣接する３つの第１粒子間がなす結合角の外角
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、フィラーが配合されたゴム材料の伸長変形を解析するためのシミュレーション方法に関する。
【背景技術】
【０００２】
下記特許文献１は、散逸粒子動力学法（Dissipative Particles Dynamics Method：ＤＰＤ法）に基づいて、フィラーが配合された高分子材料を解析するためのシミュレーション方法を提案している。下記特許文献１のシミュレーション方法では、複数の粒子で構成された３次元構造のフィラーモデルと、複数の粒子で構成された直鎖状のポリマーモデルとが用いられている。ポリマーモデルの粒子と、フィラーモデルの粒子との間には、ソフトコアポテンシャルが定義されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１３−１０８９５１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ソフトコアポテンシャルは、隣接するフィラーモデルの粒子と、ポリマーモデルの粒子との間の斥力相互作用のみが定義されている。このため、上記特許文献１のシミュレーション方法では、ポリマーがフィラーの表面に吸着される現象を計算することができない。従って、シミュレーション精度の向上には、さらなる改善の余地があった。
【０００５】
本発明は、以上のような実状に鑑み案出されたもので、シミュレーション精度を向上しうるゴム材料のシミュレーション方法を提供することを主たる目的としている。
【課題を解決するための手段】
【０００６】
本発明は、コンピュータを用いて、フィラーが配合されたゴム材料の伸長変形を解析するための方法であって、散逸粒子動力学法に基づいて、前記フィラーを複数の第１粒子で表現したフィラーモデルを、前記コンピュータに入力する工程、前記散逸粒子動力学法に基づいて、前記ゴム材料の分子鎖を複数の第２粒子で表現したポリマーモデルを、前記コンピュータに入力する工程、前記第１粒子と前記第２粒子との間に、Lenard-Jonesポテンシャルを定義する工程、前記コンピュータが、前記コンピュータの仮想空間に配置された前記フィラーモデル及び前記ポリマーモデルを用いて分子動力学計算を行うことにより、構造緩和されたゴム材料モデルを生成する工程、及び、前記コンピュータが、前記ゴム材料モデルの伸長変形を計算する工程を含むことを特徴とする。
【０００７】
本発明に係る前記ゴム材料のシミュレーション方法において、隣接する前記第１粒子間に、結合角ポテンシャルを定義する工程をさらに含んでもよい。
【０００８】
本発明に係る前記ゴム材料のシミュレーション方法において、前記結合角ポテンシャルは、下記式（１）で定義されてもよい。
【数１】
2020085692.tif 000003
ここで、
Ｕ（θ）：結合角ポテンシャル
θ：隣接する３つの第１粒子間がなす結合角の外角
ｋ：フィラーのバネ定数
θ0：平衡状態の隣接する３つの第１粒子間がなす結合角の外角
【発明の効果】
【０００９】
本発明のゴム材料のシミュレーション方法は、散逸粒子動力学法に基づいて、フィラーを複数の第１粒子で表現したフィラーモデルを入力する工程、前記散逸粒子動力学法に基づいて、前記ゴム材料の分子鎖を複数の第２粒子で表現したポリマーモデルを入力する工程、前記第１粒子と前記第２粒子との間に、Lenard-Jonesポテンシャルを定義する工程、コンピュータの仮想空間に配置された前記フィラーモデル及び前記ポリマーモデルを用いて分子動力学計算を行うことにより、構造緩和されたゴム材料モデルを生成する工程、及び、ゴム材料モデルの伸長変形を計算する工程を含んでいる。
【００１０】
本発明のシミュレーション方法では、前記散逸粒子動力学法に基づいて、前記フィラーモデル及び前記ポリマーモデルが表現されるため、全原子分子動力学法に基づくシミュレーション方法に比べて、大きな時空間を扱うことができる。
【００１１】
本発明のシミュレーション方法では、前記Lenard-Jonesポテンシャルによって、前記第１粒子と前記第２粒子との間の斥力相互作用だけでなく、引力相互作用を定義することができる。これにより、本発明のシミュレーション方法では、フィラーモデルとポリマーモデルとの吸着を表現することができるため、ゴム材料モデルの伸長変形を高い精度で計算することができる。シミュレーションすることができる。したがって、本発明のシミュレーション方法は、シミュレーション精度を向上しうる。
【図面の簡単な説明】
【００１２】
【図１】ゴム材料のシミュレーション方法を実行するためのコンピュータの一例を示す斜視図である。
【図２】ゴム材料の分子鎖の一例を示す構造式である。
【図３】ゴム材料のシミュレーション方法の処理手順の一例を示すフローチャートである。
【図４】フィラーモデル、及び、ポリマーモデルが配置された仮想空間の一例を示す概念図である。
【図５】フィラーモデルの一例を示す図である。
【図６】ポリマーモデルの一例を示す図である。
【図７】フィラーモデル、及び、ポリマーモデルの相互作用ポテンシャルの一例を説明する概念図である。
【図８】本発明の他の実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。
【図９】本発明の他の実施形態のフィラーモデルの一例を示す図である。
【図１０】結合角ポテンシャルＵ（θ）と結合角の外角θとの関係の一例を示すグラフである。
【図１１】実施例及び比較例Ａ、Ｂの応力と歪との関係を示すグラフである。
【図１２】（ａ）は、伸長変形前のゴム材料モデルを示す斜視図、（ｂ）は、伸長変形後の実施例のゴム材料モデルを示す斜視図、（ｃ）は、伸長変形後の比較例Ｂのゴム材料モデルを示す斜視図である。
【発明を実施するための形態】
【００１３】
以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態のゴム材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある。）では、コンピュータを用いて、フィラーが配合されたゴム材料の伸長変形が解析される。
【００１４】
図１は、ゴム材料のシミュレーション方法を実行するためのコンピュータの一例を示す斜視図である。コンピュータ１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含んでいる。この本体１ａには、例えば、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリ、磁気ディスクなどの記憶装置、及び、ディスクドライブ装置１ａ１、１ａ２が設けられている。また、記憶装置には、本実施形態のシミュレーション方法を実行するためのソフトウェア等が予め記憶されている。
【００１５】
フィラーとしては、例えば、シリカ、カーボンブラック、又は、アルミナ等が採用される。ゴム材料の分子鎖としては、例えば、天然ゴム、合成ゴム、又は、樹脂等が採用される。
【００１６】
図２は、ゴム材料の分子鎖２の一例を示す構造式である。本実施形態の分子鎖２としては、図２に示したcis-１，４ポリイソプレン（以下、単に「ポリイソプレン」ということがある。）が例示される。ポリイソプレンを構成する分子鎖２は、メチン基等（例えば、−ＣＨ＝、＞Ｃ＝）、メチレン基（−ＣＨ２−）、及び、メチル基（−ＣＨ３）によって構成されるイソプレンのモノマー（イソプレン分子）３が、重合度ｎで連結されて構成されている。なお、分子鎖には、ポリイソプレン以外のものが用いられてもよい。
【００１７】
図３は、ゴム材料のシミュレーション方法の処理手順の一例を示すフローチャートである。本実施形態のシミュレーション方法では、先ず、コンピュータの仮想空間（以下、単に「仮想空間」ということがある。）が、コンピュータ１に定義される（工程Ｓ１）。図４は、フィラーモデル６、及び、ポリマーモデル７が配置された仮想空間４の一例を示す概念図である。
【００１８】
仮想空間（セル）４は、ゴム材料の一部に対応させたものである。仮想空間４は、少なくとも互いに向き合う一対の面５、５、本実施形態では、互いに向き合う三対の面５、５を有しており、直方体又は立方体（本例では、立方体）として定義されている。
【００１９】
各面５、５には、周期境界条件が定義されている。このような仮想空間４が用いられることにより、後述の粗視化分子動力学計算において、例えば、後述のフィラーモデル６及びポリマーモデル７について、一方側の面５ａから出て行ったフィラーモデル６及びポリマーモデル７の一部が、他方側の面５ｂから入ってくるように計算することができる。従って、仮想空間４は、一方側の面５ａと、他方側の面５ｂとが連続している（繋がっている）ものとして取り扱うことができる。
【００２０】
仮想空間４の一辺の各長さＬ１は、適宜設定することができる。本実施形態の長さＬ１は、後述のポリマーモデル７の拡がりを示す量である慣性半径（図示省略）の３倍以上が望ましい。これにより、後述の粗視化分子動力学計算では、周期境界条件による自己のイメージとの衝突の発生を防いで、ポリマーモデル７の空間的拡がりを適切に計算することができる。また、仮想空間４の大きさは、例えば１気圧で安定な体積に設定される。これにより、仮想空間４は、解析対象のゴム材料の少なくとも一部の体積を定義することができる。仮想空間４は、コンピュータ１に記憶される。
【００２１】
次に、本実施形態のシミュレーション方法では、散逸粒子動力学法に基づいて、フィラーを表現したフィラーモデル６が、コンピュータ１に入力される（工程Ｓ２）。図５は、フィラーモデル６の一例を示す図である。
【００２２】
フィラーモデル６は、フィラー（図示省略）を複数の第１粒子１１で表現している。本実施形態のフィラーモデル６は、第１粒子１１、１１間を拘束する第１ボンドポテンシャルＰ１が定義されることにより、三次元構造を有している。
【００２３】
第１粒子１１は、仮想空間４（図４に示す）の内部で凝集するように配置されている。第１粒子１１の配置は、例えば、実際のゴム材料を撮像した電子線透過画像のフィラーの一次粒子の位置に基づいて行うことができる。これにより、フィラーモデル６は、実際のフィラーの形状を精度よく表現することができる。
【００２４】
第１粒子１１は、後述の分子動力学計算において、運動方程式の質点として取り扱われる。即ち、第１粒子１１には、質量、直径、電荷又は初期座標などのパラメータが定義される。
【００２５】
第１ボンドポテンシャルＰ１は、下記式（２）のポテンシャルＵHarmonicとして定義される。
【００２６】
【数２】
2020085692.tif 000004
ここで、
ｋ：粒子間のばね定数
ｒij：粒子間の距離
ｒ0：平衡長  
なお、距離ｒij及び平衡長ｒ0は、粒子の中心間の距離として定義される。
【００２７】
上記式（２）で定義される第１ボンドポテンシャルＰ１は、Harmonic型である。Harmonic型とは、いわゆる線形バネが定義され、平衡長ｒ0からのずれの大きさに比例した復元力が働くポテンシャルである。このような第１ボンドポテンシャルＰ１は、距離ｒijの上限値が、平衡長ｒ0に限定されない（即ち、伸びきり長が設定されない）。各変数は、論文１（ Robert D. Groot and Patrick B. Warren 著、「Dissipaticle dynamics: Bridging the gap between atomistic and mesoscopic simulation」、J. Chem Phys. vol.107, No.11, 15 September 1997、p4423-4435）に基づいて、適宜設定することができる。フィラーモデル６は、コンピュータ１に記憶される
【００２８】
次に、本実施形態のシミュレーション方法では、散逸粒子動力学法に基づいて、分子鎖２を表現したポリマーモデル７（図４に示す）が、コンピュータ１に入力される（工程Ｓ３）。図６は、ポリマーモデル７の一例を示す図である。
【００２９】
ポリマーモデル７は、分子鎖２（図２に示す）を複数の第２粒子１２で表現している。本実施形態のポリマーモデル７は、複数個の第２粒子１２と、第２粒子１２、１２間を接続する第２結合鎖１４とから構成された三次元構造を有している。
【００３０】
第２粒子１２は、ポリマーのモノマー又はモノマーの一部分をなす構造単位を置換したものである。本実施形態の第２粒子１２は、任意の個数分（例えば、数十個分）のモノマー３を構造単位として、１個の第２粒子１２に置換される。これにより、ポリマーモデル７には、複数（例えば、１０〜５００個）の第２粒子１２が設定される。
【００３１】
第２粒子１２は、後述の分子動力学計算において、運動方程式の質点として取り扱われる。即ち、第２粒子１２には、例えば、質量、直径、電荷又は初期座標などのパラメータが定義される。
【００３２】
第２結合鎖１４は、各第２粒子１２、１２間を拘束するものである。第２結合鎖１４は、第２ボンドポテンシャルＰ２で定義される。第２ボンドポテンシャルＰ２は、下記式（３）のポテンシャルＦSijで定義される。
【００３３】
【数３】
2020085692.tif 000005
ここで、
Ｃ：ばね定数
ｒij：粒子間の距離
ｒs：平衡結合長
ｎij：ｒij／｜ｒij｜（単位ベクトル）
なお、距離ｒij及び平衡結合長ｒsは、粒子の中心間の距離として定義される。
【００３４】
上記式（３）で定義される第２ボンドポテンシャルＰ２は、ＤＰＤ型のボンドポテンシャルである。ＤＰＤ型のボンドポテンシャルは、結合長がガウス分布に従う伸縮ポテンシャルであり、ＤＰＤシミュレーションで一般的に使用される。各変数は、例えば、論文２（Satoru Yamamoto 著、「Mesoscopic Simulation for the Drying Process of Polymer Films」、日本レオロジー学会誌、vol.32, No.5, p295-301）に基づいて、適宜設定することができる。
【００３５】
仮想空間４内に配置されるポリマーモデル７の個数については、仮想空間４の大きさ等に基づいて、適宜設定することができる。ポリマーモデル７の個数の一例としては、５０〜１０００本程度である。ポリマーモデル７は、コンピュータ１に記憶される。
【００３６】
次に、本実施形態のシミュレーション方法では、隣接するフィラーモデル６の第１粒子１１、１１間、及び、隣接するポリマーモデル７の第２粒子１２、１２間に、相互作用ポテンシャルＰ３、Ｐ４が定義される（工程Ｓ４）。図７は、フィラーモデル６、及び、ポリマーモデル７の相互作用ポテンシャルの一例を説明する概念図である。
【００３７】
隣接するフィラーモデル６の第１粒子１１、１１間の相互作用ポテンシャルＰ３、及び、隣接するポリマーモデル７の第２粒子１２、１２間の相互作用ポテンシャルＰ４は、散逸粒子動力学法に基づいて、下記式（４）のポテンシャルＦcijによって定義される。
【００３８】
【数４】
2020085692.tif 000006
ここで、
ａij：ポテンシャルの強度
ｒij：粒子間の距離
ｒc：カットオフ距離
ｎij：ｒij／｜ｒij｜（単位ベクトル）
なお、距離ｒij及びカットオフ距離ｒcは、各粒子の中心間の距離として定義される。
【００３９】
上記式（４）で定義される相互作用ポテンシャルＰ３、Ｐ４は、粒子間の距離ｒijが、予め定められたカットオフ距離ｒc以下になったときに、粒子間に斥力を作用させる。さらに、粒子間の距離ｒijが予め定めたカットオフ距離ｒcに等しい場合、相互作用ポテンシャルＰ３、Ｐ４はゼロであり、粒子間に斥力は生じない。従って、相互作用ポテンシャルＰ３、Ｐ４は、隣接するフィラーモデル６の第１粒子１１、１１間、及び、隣接するポリマーモデル７の第２粒子１２、１２間の斥力相互作用のみが定義される。各変数は、例えば、論文２に基づいて適宜設定することができる。相互作用ポテンシャルＰ３、Ｐ４は、コンピュータ１に記憶される。
【００４０】
次に、本実施形態のシミュレーション方法では、第１粒子１１と第２粒子１２との間に、Lenard-JonesポテンシャルＵLJ（ｒij）が定義される（工程Ｓ５）。Lenard-Jonesポテンシャルは、下記式（５）で定義される。
【００４１】
【数５】
2020085692.tif 000007
ここで、
ｒij：粒子間の距離
ｒc：カットオフ距離
ε：粒子間に定義されるＬＪポテンシャルの強度
σ：粒子の直径に相当
なお、距離ｒij及びカットオフ距離ｒcは、各粒子の中心間の距離として定義される。
【００４２】
Lenard-JonesポテンシャルＵLJ（ｒij）は、第１粒子１１と第２粒子１２との間の距離ｒijが２1/6σから０に近づくほど、その値が無限大に大きくなる。このとき、Lenard-JonesポテンシャルＵLJ（ｒij）は、第１粒子１１と第２粒子１２との間に斥力相互作用を生じさせることができる。
【００４３】
Lenard-JonesポテンシャルＵLJ（ｒij）は、第１粒子１１と第２粒子１２との間の距離ｒijが２1/6σであるときに最小となり、その値が−εとなる。即ち、Lenard-JonesポテンシャルＵLJ（ｒij）が最小となるのは、第１粒子１１と第２粒子１２とが当接したときである。
【００４４】
Lenard-JonesポテンシャルＵLJ（ｒij）は、第１粒子１１と第２粒子１２との間の距離ｒijが２1/6σよりも大きくなるほど、その値が大きくなって０に近づく。このとき、Lenard-JonesポテンシャルＵLJ（ｒij）は、第１粒子１１と第２粒子１２との間に引力相互作用を生じさせることができる。
【００４５】
このように、Lenard-JonesポテンシャルＵLJ（ｒij）は、第１粒子１１と第２粒子１２との間の斥力相互作用だけでなく、引力相互作用を定義することができる。Lenard-JonesポテンシャルＵLJ（ｒij）の各変数は、例えば、論文３（ Kurt Kremer & Gary S. Grest 著、「Dynamics of entangled linear polymer melts: A molecular-dynamics simulation」、J. Chem Phys. vol.92, No.8, 15 April 1990、p5057-5086）に基づいて、適宜設定することができる。Lenard-JonesポテンシャルＵLJ（ｒij）は、コンピュータ１に記憶される。
【００４６】
次に、本実施形態のシミュレーション方法では、コンピュータ１が、構造緩和されたゴム材料モデルを生成する（工程Ｓ６）。工程Ｓ６では、図４に示した仮想空間４に配置されたフィラーモデル６及びポリマーモデル７を用いて、散逸粒子動力学法に基づく分子動力学計算が行われることにより、構造緩和されたゴム材料モデル２０が生成される。
【００４７】
本実施形態の分子動力学計算では、例えば、仮想空間４について所定の時間、フィラーモデル６、及び、ポリマーモデル７が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻でのフィラーモデル６、及び、ポリマーモデル７の動きが、シミュレーションの単位時間毎に追跡される。構造緩和の計算は、例えば（株）ＪＳＯＬ社製のソフトマテリアル総合シミュレーター（Ｊ−ＯＣＴＡ）に含まれるＣＯＧＮＡＣ、又は、ＶＳＯＰを用いて処理することができる。これにより、工程Ｓ６では、構造緩和されたゴム材料モデル２０を生成することができる。散逸粒子動力学法に基づくゴム材料モデル２０は、コンピュータ１に記憶される。
【００４８】
本実施形態のシミュレーション方法では、上記式（５）で定義されるLenard-JonesポテンシャルＵLJ（ｒij）により、第１粒子１１と第２粒子１２との間の斥力相互作用、及び、引力相互作用を定義することができる。これにより、本実施形態のシミュレーション方法は、散逸粒子動力学法に基づくシミュレーションにおいて、フィラーモデル６とポリマーモデル７との吸着を表現することができる。
【００４９】
このような吸着（物理吸着）は、フィラー（図示省略）と分子鎖２（図２に示す）との間に働くファンデルワールス相互作用及び静電相互作用によって、分子鎖２がフィラーの表面に吸着される現象であり、吸着の強さが、ゴム材料の強度に影響すると考えられている。従って、本実施形態のシミュレーションでは、フィラーモデル６とポリマーモデル７との吸着を表現できるため、シミュレーション精度を向上しうる。
【００５０】
次に、本実施形態のシミュレーション方法では、コンピュータ１が、ゴム材料モデル２０の伸長変形を計算する（工程Ｓ７）。工程Ｓ７では、図４に示したゴム材料モデル２０を予め定められた方向に引っ張る単軸引張試験が計算される。ゴム材料モデル２０を引っ張る方向については、適宜選択することができ、本例では、ｚ軸方向に引っ張っている。
【００５１】
ゴム材料モデル２０の変形計算は、ｚ軸方向において、ゴム材料モデル２０の一端（図４に示した仮想空間４の一方側の面５ａ）、及び、ゴム材料モデル２０の他端（図５に示した仮想空間４の他方側の面５ｂ）が互いに離間するように、ゴム材料モデル２０の伸長が計算される。ゴム材料モデル２０の伸長計算は、分子動力学計算に基づいて計算される。
【００５２】
工程Ｓ７では、ゴム材料モデル２０の伸長計算により、フィラーモデル６、及び、ポリマーモデル７の熱運動が計算される。このような熱運動は、ゴム材料モデル２０に与えられた歪み、ボンドポテンシャルＰ１、Ｐ２（図５及び図６に示す）、相互作用ポテンシャルＰ３、Ｐ４（図７に示す）、Lenard-JonesポテンシャルＵLJ（ｒij）、及び、運動方程式に基づいて計算される。
【００５３】
変形条件については、適宜設定することができる。変形条件の一例としては、ゴム材料モデル２０に与えられる歪みが０．１〜０．３程度であり、また、１歪み当たりの変形速度Ｖ１（図４に示す）が１０００〜１００００τ程度である。
【００５４】
工程Ｓ７では、ゴム材料モデル２０の変形開始から変形終了までの間、分子動力学計算の単位ステップ（ＭＤステップ）毎に、ゴム材料モデル２０の変形に伴う物理量が計算される。物理量としては、適宜採用することができる。本実施形態の物理量は、ゴム材料モデル２０の応力が含まれる。工程Ｓ７では、ゴム材料モデル２０に与えられる歪毎に、ゴム材料モデル２０の伸長方向の応力成分が計算される。これにより、本実施形態のシミュレーション方法では、ゴム材料モデル２０の応力と歪との関係を示す応力−歪曲線を取得することができる。この応力−歪曲線に基づいて、ゴム材料の強度を評価することができる。
【００５５】
本実施形態のシミュレーション方法では、散逸粒子動力学法に基づいて、フィラーモデル６及びポリマーモデル７が表現されるため、全原子分子動力学法に基づくシミュレーション方法に比べて、大きな時空間を扱うことができる。また、本実施形態のシミュレーション方法では、上記式（５）で定義されるLenard-JonesポテンシャルＵLJ（ｒij）により、散逸粒子動力学法に基づくシミュレーションにおいて、フィラーモデル６とポリマーモデル７との吸着を表現することができる。従って、本実施形態のシミュレーション方法では、ゴム材料モデル２０の強度、及び、物理量の計算を精度よく行うことができ、シミュレーション精度を向上させることができる。
【００５６】
次に、本実施形態のシミュレーション方法では、コンピュータ１が、伸長変形したゴム材料モデル２０の物理量が良好か否かを判断する（工程Ｓ８）。良好か否かを判断するための物理量の基準については、ゴム材料に求められる強度等に応じて適宜設定することができる。
【００５７】
工程Ｓ８において、ゴム材料モデル２０の物理量が良好と判断された場合（工程Ｓ８において、「Ｙ」）、ゴム材料モデル２０に設定された諸条件に基づいて、ゴム材料が製造される（工程Ｓ９）。他方、工程Ｓ８において、物理量が良好なゴム材料モデル２０が存在しないと判断された場合（工程Ｓ８において、「Ｎ」）、ゴム材料モデル２０に設定された諸条件が変更され（工程Ｓ１０）、工程Ｓ１〜工程Ｓ８が再度実施される。これにより、本実施形態のシミュレーション方法は、フィラーとポリマーとの相互作用に基づいて、強度に優れるゴム材料を開発することができる。
【００５８】
図５に示したように、本実施形態のフィラーモデル６の第１ボンドポテンシャルＰ１は、上記式（２）のポテンシャルＵHarmonicで定義されるため、第１粒子１１、１１間の拘束力が弱い。このため、図４に示したゴム材料モデル２０の伸長変形を計算する工程Ｓ７において、第１粒子１１、１１間で引き伸ばされ、フィラーモデル６が大きく変形（膨張）する場合がある。このように膨張したフィラーモデル６は、周期境界条件による自己のイメージとの衝突によって大きく変形し、フィラーモデル６の形状を維持できなくなる場合がある。このようなフィラーモデル６の膨張及び変形を防ぐために、フィラーモデル６には、隣接する第１粒子１１、１１間の拘束力が高められるのが望ましい。図８は、本発明の他の実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。図９は、本発明の他の実施形態のフィラーモデル６の一例を示す図である。なお、この実施形態において、これまでの実施形態と同一の構成については、同一の符号を付し、説明を省略することがある。
【００５９】
この実施形態のシミュレーション方法では、隣接する第１粒子１１、１１間に、結合角ポテンシャルＵ（θ）が定義される（工程Ｓ１１）。結合角ポテンシャルＵ（θ）は、隣接する第１粒子１１、１１間を拘束するためのものである。本実施形態の結合角ポテンシャルＵ（θ）は、例えば、下記式（１）で定義される。
【００６０】
【数１】
2020085692.tif 000008
ここで、
Ｕ（θ）：結合角ポテンシャル
θ：隣接する３つの第１粒子間がなす結合角の外角
ｋ：フィラーのバネ定数
θ0：平衡状態の隣接する３つの第１粒子間がなす結合角の外角
【００６１】
平衡状態の結合角の外角θ0については、フィラーモデル６の構造に応じて適宜設定することができる。図１０は、結合角ポテンシャルＵ（θ）と結合角の外角θとの関係の一例を示すグラフである。図１０の結合角ポテンシャルＵ（θ）は、平衡状態の結合角の外角θ0が６０°に設定されている。
【００６２】
結合角ポテンシャルＵ（θ）は、結合角の外角θが６０°又は−６０°のときに０となる。一方、結合角ポテンシャルＵ（θ）は、結合角の外角θが６０°又は−６０°から乖離するほど大きくなる。また、結合角ポテンシャルＵ（θ）は、ばね定数ｋが大きいほど、隣接する第１粒子１１、１１間を強固に拘束することができる。
【００６３】
このように、結合角ポテンシャルＵ（θ）は、フィラーモデル６の隣接する３つの第１粒子１１間がなす結合角の外角θが、平衡状態の結合角の外角θ0とは異なる場合に、結合角の外角θを平衡状態の結合角の外角θ0に戻そうと作用を定義できる。これにより、フィラーモデル６の隣接する３つの第１粒子１１間がなす結合角の外角θを、平衡状態の結合角の外角θ0に維持することができる。
【００６４】
このような結合角ポテンシャルＵ（０）により、隣接する第１粒子１１、１１を、平衡状態を維持するように強固に拘束できるため、図４に示したゴム材料モデル２０の伸長変形を計算する工程Ｓ７において、フィラーモデル６の膨張及び変形するのを防ぐことができる。従って、この実施形態のシミュレーション方法は、フィラーモデル６とポリマーモデル７との吸着を効果的に表現することができるため、シミュレーション精度を向上させることができる。
【００６５】
工程Ｓ１１では、フィラーモデル６の構造に応じて、平衡状態の結合角の外角θ0が異なる複数種類の結合角ポテンシャルＵ（θ）が定義されてもよい。本実施形態では、平衡状態の結合角の外角θ0が０°に設定された結合角ポテンシャルＵ（θ）と、平衡状態の結合角の外角θ0が６０°に設定された結合角ポテンシャルＵ（θ）とが定義されている。これにより、フィラーモデル６の変形を効果的に防ぐことができる。結合角ポテンシャルＵ（θ）は、コンピュータ１に記憶される。なお、結合角ポテンシャルＵ（θ）は、隣接する第１粒子１１、１１間を拘束できるものであれば、上記式（１）に限定されない。例えば、フィラー形状を固定化するポテンシャルが定義されることにより、フィラー粒子が連結して房状の凝集体になったフィラーモデルや、表面が非結晶のアモルファス構造を有するフィラーモデルを定義することができる。
【００６６】
以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。
【実施例】
【００６７】
構造緩和されたゴム材料モデルが生成され、ゴム材料モデルの伸長変形が計算された（実施例、比較例Ａ、比較例Ｂ及び比較例Ｃ）。実施例、比較例Ａ、及び、比較例Ｂでは、散逸粒子動力学法に基づいて、フィラーを複数の第１粒子で表現したフィラーモデル、及び、分子鎖を複数の第２粒子で表現したポリマーモデルが入力された。一方、比較例Ｃでは、Kremer-Grestの粗視化モデルに基づいて、フィラーモデル及びポリマーモデルが入力された。
【００６８】
実施例及び比較例Ｃの第１粒子と第２粒子との間には、Lenard-Jonesポテンシャルが定義された。さらに、実施例及び比較例Ｃでは、隣接する第１粒子間に、結合角ポテンシャルが定義された。
【００６９】
比較例Ａ、Ｂの第１粒子と第２粒子との間には、従来の散逸粒子動力学法に基づいて、ソフトコアポテンシャルが定義された。さらに、比較例Ａでは、隣接する第１粒子間に、結合角ポテンシャルが定義された。一方、比較例Ｂでは、Ｂでは、隣接する第１粒子間に、結合角ポテンシャルが定義されなかった。
共通仕様、及び、各ポテンシャルのパラメータについては、次のとおりである。
仮想空間：
一辺の長さＬ１：１７．２４σ
体積：５１２１σ3
フィラーモデル：
フィラーモデルの個数：４個
１個のフィラーモデルを構成する第１粒子の個数：７９１個
フィラーモデルの第１粒子の密度：３（１／σ3）
ポリマーモデル：
ポリマーモデルの個数：１２０個
１個のポリマーモデルを構成する第２粒子：１００個
架橋密度：０．０４（１／σ3）
構造緩和計算：
時間刻み：０．０１
ステップ数：１,０００,０００
変形計算：
歪み：０．０１τ
合計ステップ：１００,０００
合計時間：１,０００τ
ひずみ速度：０．００１（１／τ）
伸長速度：０．０１７２（σ／τ）
ポアソン比：０．５
フィラーモデルの第１粒子間のボンドポテンシャル（ＵHarmonic）：
平衡長ｒ0：０．７８σ
ばね定数ｋ：１００（ε／σ2）
ポリマーモデルの第２粒子間のボンドポテンシャル（Ｆsij）：
平衡結合長ｒs：０．８６σ
ばね定数Ｃ：４
相互作用ポテンシャル（ソフトコアポテンシャル）Ｐ３、Ｐ４：
カットオフ距離ｒc：１．０σ
ポテンシャルの強度ａij：２５
Lenard-Jonesポテンシャル（実施例、比較例Ｃ）：
カットオフ距離：２．５σ
ポテンシャルの強度ε：１．０σ
第１粒子及び第２粒子の直径に相当する値σ：１．０σ
第１粒子間の結合角ポテンシャル（実施例、比較例Ｂ）：
フィラーのばね定数ｋ：２００ε
平衡状態の結合角の外角θ0：０°、６０°
【００７０】
図１１は、実施例及び比較例Ａ、Ｂの応力と歪との関係を示すグラフである。実施例では、歪０．７付近から応力が起ち上がる挙動を確認できた。このような挙動は、ゴム材料の伸長変形時の挙動に近似する。一方、比較例Ａ、Ｂでは、歪０．７付近から応力が起ち上がる挙動を確認できなかった。従って、実施例では、Lenard-Jonesポテンシャルにより、フィラーモデルとポリマーモデルとの吸着を表現することができた。さらに、実施例では、結合角ポテンシャルにより、フィラーの硬さを表現することができ、シミュレーション精度を向上させることができた。また、実施例では、その計算時間を、比較例Ｃの計算時間の１／５０に短縮することができた。したがって、実施例は、計算コストを大幅に低減することができた。
【００７１】
図１２（ａ）は、伸長変形前のゴム材料モデルを示す斜視図である。図１２（ｂ）は、実施例の伸長変形後のゴム材料モデルを示す斜視図である。図１２（ｃ）は、比較例Ｂの伸長変形後のゴム材料モデルを示す斜視図である。図１２（ａ）〜（ｃ）では、ポリマーモデルを省略して示している。
【００７２】
図１２（ｂ）に示した実施例の伸長変形後のフィラーモデルは、第１粒子間の結合角ポテンシャルにより、図１２（ａ）に示した伸長変形前のフィラーの形状を維持できた。一方、図１２（ｃ）に示した比較例Ｂの伸長変形後のフィラーモデルは、フィラーモデルの形状が大きく変形し、図１２（ａ）に示した伸長変形前のフィラーの形状を維持できなかった。従って、実施例のゴム材料モデルは、シミュレーション精度を向上させることができた。
【符号の説明】
【００７３】
６      フィラーモデル
７      ポリマーモデル
１１      第１粒子
１２      第２粒子
【図１】
2020085692.tif 000010
【図２】
2020085692.tif 000011
【図３】
2020085692.tif 000012
【図４】
2020085692.tif 000013
【図５】
2020085692.tif 000014
【図６】
2020085692.tif 000015
【図７】
2020085692.tif 000016
【図８】
2020085692.tif 000017
【図９】
2020085692.tif 000018
【図１０】
2020085692.tif 000019
【図１１】
2020085692.tif 000020
【図１２】
2020085692.tif 000021