生命科学関連特許情報
| タイトル： | 公開特許公報(A)_粘弾性体・フィラー複合体の動的挙動のシミュレーション方法、シミュレーション装置およびシミュレーションプログラム |
| 出願番号： | 2012123688 |
| 年次： | 2013 |
| IPC分類： | G06F 17/50,G06F 19/00,G01N 19/00 |
特許情報キャッシュ
大熊 孝広 JP 2013250673 公開特許公報(A) 20131212 2012123688 20120530 粘弾性体・フィラー複合体の動的挙動のシミュレーション方法、シミュレーション装置およびシミュレーションプログラム 株式会社ブリヂストン 000005278 杉村 憲司 100147485 冨田 和幸 100119530 山口 雄輔 100164448 大熊 孝広 G06F 17/50 20060101AFI20131115BHJP G06F 19/00 20110101ALI20131115BHJP G01N 19/00 20060101ALI20131115BHJP JPG06F17/50 612AG06F17/50 612GG06F19/00 110G01N19/00 A 12 2 ＯＬ 20 5B046 5B046AA04 5B046JA04 5B046JA09 本発明は、シミュレーション方法、シミュレーション装置およびシミュレーションプログラムに関し、より詳細には、粘弾性体中にフィラーが分散した材料（以下、「粘弾性体・フィラー複合体」と称する）の動的挙動をシミュレーションするシミュレーション方法、シミュレーション装置およびシミュレーションプログラムに関するものである。 従来、粘弾性体、例えばゴムにカーボンブラックやシリカ等のフィラーを配合すると補強効果があることが知られており、ゴムにフィラーを配合した材料（以下、「ゴム・フィラー複合体」と称する）が各種タイヤやコンベヤベルト等のゴム製品に適用されている。このゴム・フィラー複合体を構成するフィラー部分は、複数のフィラーが連なった複雑なネットワーク構造を形成しており、ゴム・フィラー複合体において、ある程度偏ったフィラー分布を示すことが知られている。 このような複雑なフィラー構造を有するゴム・フィラー複合体の高い性能を実現するためには、変形が与えられた際の、ゴム・フィラー複合体の動的挙動を把握することが肝要である。そこで、これまで、こうしたゴム・フィラー複合体の動的挙動をシミュレーションするための様々な手法が提案されている。 例えば、特許文献１には、有限要素法（ＦＥＭ）を用いて、ゴム製品の粘弾性応答性能を予測する方法について記載されている。また、特許文献２には、粗視化分子動力学法を用いて、ゴムとカーボンブラックとを少なくとも含むゴム材料を評価する技術について記載されている。 しかし、特許文献１に記載された方法では、固体表面の境界に対して、ノンスリップ条件等の適切な境界条件を満足させるような計算が必要となり、ゴム・フィラー複合体のように境界が多い材料においては、計算コストが膨大になる問題がある。また、ＦＥＭ等においては、解析対象を小さな要素に分割するためのメッシュを作成するが、このメッシュをフィラーの形状に合わせて設定する必要があり、フィラーの位置が時々刻々変化する場合には、そのメッシュの再作成の計算コストが膨大となることも問題である。 また、特許文献２に記載された方法では、ＦＥＭにおいて必要となる境界条件を考慮する必要はないが、マクロな物性を評価するためには、膨大な計算コストがかかり、ゴム・フィラー複合体のような大きな系のマクロな物性を評価することは現実的には不可能である。 一方、非特許文献１には、液体の溶媒中にコロイド粒子が分散したコロイド分散系において、コロイド粒子を剛体ではなく粘性の高い液体ドメインと見なして、コロイド分散系の動的挙動を計算する、いわゆる流体粒子ダイナミクス（Ｆｌｕｉｄ Ｐａｒｔｉｃｌｅ Ｄｙｎａｍｉｃｓ，ＦＰＤ）法が提案されている。このＦＰＤ法においては、コロイド粒子の高粘性ドメインと周囲の粘性液体の物性とが滑らかに接続されるため、境界条件をあらわに解くことなく、流体効果を考慮したコロイド分散系の動的挙動を高効率にシミュレーションすることが可能である。 さらに、特許文献３には、非特許文献１に記載された方法を発展させて、荷電粒子が液体の溶媒中に分散した荷電コロイド分散系の動的挙動をシミュレーションする方法について記載されており、電気泳動ディスプレイにおける電極間の帯電顔料粒子のコロイド分散液において生じるコロイド電気泳動現象のシミュレーションが行われている。 このように、ＦＰＤ法により、コロイド分散系の動的挙動を高効率にシミュレーションできていることから、上記ＦＰＤ法をゴム・フィラー複合体の動的挙動のシミュレーションに適用して、高性能なゴム・フィラー複合体の設計に役立たせることが期待される。特開２００５−４９３３３号公報特開２００６−６４６５８号公報特開２００５−３０８７８２号公報Ｈ．Ｔａｎａｋａ ａｎｄ Ｔ．Ａｒａｋｉ，Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．８５（６），１３３８（２０００） しかしながら、非特許文献１に記載されたＦＰＤ法においては、溶媒として粘性の液体を用いることが想定されているため、ＦＰＤ法の、ゴム・フィラー複合体のような、溶媒が粘弾性体の場合に対する適用は一切検討されて来なかった。そのため、どのような条件の下でシミュレーションを行うことができるのか不明となっていた。 そこで、本発明の目的は、ゴムをはじめとする粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動を簡便にシミュレーションできるシミュレーション方法、シミュレーション装置およびシミュレーションプログラムを提供することにある。 発明者は、上記課題を解決する方途について鋭意検討した。上記ＦＰＤ法において、溶媒が従う流体の運動方程式であるナビエ・ストークス方程式は、液体の溶媒および固体であるコロイド粒子の粘度をパラメータとして含み、コロイド粒子の粘度を液体である溶媒の粘度よりも高く設定することにより、粒子内に流れが入り込まず、コロイド粒子の境界条件が加味される。しかしながら、上記ＦＰＤ法においては、液体の溶媒中に剛体粒子であるフィラーを充填させても、フィラーの弾性的性質を取り込むことができず、フィラーによるゴム弾性の向上効果を評価できないことが判明した。そこで、この問題を解決するためには、溶媒である粘弾性体の弾性的な特性を加味した式を用いてシミュレーションを行う必要があるとの認識にたどり着き、本発明を完成させるに至った。 すなわち、本発明の要旨構成は以下の通りである。（１）粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動をシミュレーションする方法であって、ａ）前記粘弾性体および前記フィラーの物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力ステップと、ｂ）入力された前記計算条件の下で前記粘弾性体が従う流体の運動方程式を解く演算ステップと、ｃ）前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出ステップと、ｄ）算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新ステップと、ｅ）前記ステップｂ）〜ｄ）の繰り返し処理により求めた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップと、を含み、前記粘弾性体が従う流体の運動方程式は、前記粘弾性体の弾性的性質を加味したものであることを特徴とするシミュレーション方法。（２）前記粘弾性体が従う流体の運動方程式は以下の式（Ａ）で表される、前記（１）に記載のシミュレーション方法。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆｐは外場、Πは連続体の応力テンソル、Πｊ（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Πｊ（ｆ）は連続体の応力テンソルのフィラー部分の成分、Ｎｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、τｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、ηｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、ηｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。（３）前記式（Ａ）における粘度ηｊ（ｆ）およびηｊ（ｂ）、緩和時間τｊ（ｆ）およびτｊ（ｂ）は下記の式（Ｂ）を満足する、前記（２）に記載のシミュレーション方法。（４）前記データ出力ステップは、前記求められたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに含む、前記（１）〜（３）のいずれか一項に記載のシミュレーション方法。（５）粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動をシミュレーションする装置であって、前記粘弾性体および前記フィラーの物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力部と、入力された前記計算条件の下で前記粘弾性体が従う流体の運動方程式を解く演算器と、前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出器と、算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新器と、前記演算器、前記移動量算出器および前記配置更新器により求められた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力器とを有する解析部とを備え、前記ゴムが従う流体の運動方程式は、前記ゴムの弾性的性質を加味したものであることを特徴とするシミュレーション装置。（６）前記粘弾性体が従う流体の運動方程式は以下の式（Ａ）で表される、前記（５）に記載のシミュレーション装置。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆｐは外場、Πは連続体の応力テンソル、Πｊ（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Πｊ（ｆ）は連続体の応力テンソルのフィラー部分の成分、Ｎｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、τｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、ηｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、ηｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。（７）前記式（Ａ）における粘度ηｊ（ｆ）およびηｊ（ｂ）、緩和時間τｊ（ｆ）およびτｊ（ｂ）は下記の式（Ｂ）を満足する、前記（６）に記載のシミュレーション装置。（８）前記解析部は、前記得られたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析器をさらに有する、前記（５）〜（７）のいずれか一項に記載のシミュレーション装置。（９）粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動を求める装置として構成するコンピュータに、ａ）前記粘弾性体および前記フィラーの物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力ステップと、ｂ）入力された前記計算条件の下で前記粘弾性体が従う流体の運動方程式を解く演算ステップと、ｃ）前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出ステップと、ｄ）算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新ステップと、ｅ）前記ステップｂ）〜ｄ）の繰り返し処理により求めた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップとを実行させ、前記流体の運動方程式は、前記粘弾性体の弾性的性質を加味したものであることを特徴とするシミュレーションプログラム。（１０）前記粘弾性体が従う流体の運動方程式は以下の式（Ａ）で表される、前記（９）に記載のシミュレーションプログラム。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆｐは外場、Πは連続体の応力テンソル、Πｊ（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Πｊ（ｆ）は連続体の応力テンソルのフィラー部分の成分、Ｎｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、τｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、ηｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、ηｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。（１１）前記式（１）における粘度ηｊ（ｆ）およびηｊ（ｂ）、緩和時間τｊ（ｆ）およびτｊ（ｂ）は下記の式（Ｂ）を満足する、前記（１０）に記載のシミュレーションプログラム。（１２）前記データ出力ステップは、前記求められたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに含む、前記（９）〜（１１）のいずれか一項に記載のシミュレーションプログラム。 本発明によれば、粘弾性体・フィラー複合体の動的挙動を簡便にシミュレーションできるシミュレーション方法、シミュレーション装置およびシミュレーションプログラムを提供することができる。非特許文献１に記載のＦＰＤ法により求めたゴム・フィラー複合体の物性を示す図である。本発明のシミュレーション方法の一実施形態を示すフローチャートである。（ａ）はＲ＝（ηｆ／τｆ）／（ηｂ／τｂ）とフィラー効果との関係を示す図であり、（ｂ）はτｆ／τｂとフィラー効果との関係を示す図である。本発明のシミュレーション装置の一実施形態を示す図である。（ａ）はフィラーの充填率とＧｕｔｈ−Ｇｏｌｄ因子との関係を示す図であり、（ｂ）は周波数と振幅および位相シフトとの関係を示す図である。本発明および非特許文献１に記載の方法により得られたゴム・フィラー複合体の弾性率を示す図である。（シミュレーション方法） 以下、図面を参照して、本発明の実施形態について説明する。 本発明の粘弾性体・フィラー複合体の動的挙動をシミュレーションする方法は、ａ）粘弾性体およびフィラーの物性値、並びに粘弾性体中におけるフィラーの初期配置を含む、動的挙動を求めるための計算条件を入力する入力ステップと、ｂ）該入力された計算条件の下で粘弾性体が従う流体の運動方程式を解く演算ステップと、ｃ）流体の運動方程式の解に基づいてフィラーの移動量を算出する移動量算出ステップと、ｄ）算出されたフィラーの移動量を用いてフィラーの配置を更新する配置更新ステップと、ｅ）ステップｂ）〜ｄ）の繰り返し処理により求めた粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップとを含む。ここで、粘弾性体が従う流体の運動方程式は、粘弾性体の弾性的性質を加味したものであることが肝要である。なお、本発明において、「フィラー」とは、カーボンブラックやシリカ等、粘弾性体の補強効果を図るために粘弾性体に配合される充填剤であり、粒子状のものを意味している。ここで、フィラーの形状は球形に限定されない。また、「粘弾性体・フィラー複合体」とは、粘弾性体中にフィラーが分散された材料を意味している。以下の説明においては、粘弾性体がゴムの場合を例に説明するが、それに限定されない。 上述のように、非特許文献１に記載されたＦＰＤ法は、特許文献１に記載されたＦＥＭのような境界条件をあらわに解く必要なしに、粘性液体溶媒の流体効果を考慮したコロイド分散系の動的挙動を高効率にシミュレーションすることを可能にする。このＦＰＤ法において、溶媒である粘性液体の運動は、下記のナビエ・ストークス方程式で表される。ここで、Ｄ／Ｄｔはラグランジュ微分（＝∂／∂ｔ＋ｖ・∇）、Ｆｐは力場、ρは質量密度、ｖは速度、ｐは圧力、ηは粘度、ζは熱遥動力である。 ＦＰＤ法を、ゴム・フィラー複合体の動的挙動のシミュレーションに適用する場合、上記式（１）において、粘度ηの値を、溶媒として粘性液体を用いる場合よりも大きな値に設定することにより、ゴムのフィラーによる補強性を近似的に取り込むことができると考えられる。そこで、式（１）を用いて、ゴム・フィラー複合体の動的挙動をシミュレーションしてその粘弾性特性を調べた。具体的には、ゴムの粘度に対して熱エネルギーは非常に小さいと考えられるためゼロとし、２次元のモデルに剛体フィラーの充填率（面積分率）を０％〜７％まで変化させて、周期せん断条件の下、シミュレーションを行った。無次元化は、密度ρ=１、粘度η＝０．５、フィラー半径Ｒ＝３．８となる条件で行った。得られた結果を図１に示す。ここで、図１（ａ）はフィラー充填率と貯蔵弾性率Ｇ’との関係を、図１（ｂ）はフィラー充填率と損失弾性率Ｇ’’との関係をそれぞれ示している。ここで、貯蔵弾性率Ｇ’および損失弾性率Ｇ’’は、後述するように、動的粘弾性を特徴付ける量であり、貯蔵弾性率Ｇ’および損失弾性率Ｇ’’は、粘弾性体の弾性的性質、粘性的性質をそれぞれ特徴付ける量である。これらの図から明らかなように、フィラー充填率の増加とともに損失弾性率Ｇ’’は増加する一方、貯蔵弾性率Ｇ’は、フィラー充填率にかかわらず０であり、ＦＰＤ法では粘弾性体の弾性的性質を取り込めていないことが分かる。すなわち、ＦＰＤ法によりシミュレーションを行うと、ゴムにフィラーを充填しても弾性率を求めることができず、適切な結果が得られないことが判明した。 このように、発明者は、従来のＦＰＤ法では、粘弾性体の動的挙動のシミュレーションを正しく行えないことを受けて、溶媒としてのゴムの弾性的な性質を取り込むことができる流体の運動方程式を用いてシミュレーションを行う必要があるとの認識にたどり着いたのである。そこで、発明者は、試行錯誤の末、以下に示す、本発明のゴム・フィラー複合体の動的挙動のシミュレーション方法に用いる流体の運動方程式を導出するに至った。 すなわち、まず、非圧縮性粘性流体の運動量保存則は、以下の式（２）で与えられる。ここで、Πは連続体の応力テンソルであり、と表され、また、ｖは、連続の条件から以下の式（４）の関係を満たす。 本発明においては、粘弾性体の弾性的性質を加味するために、鋭意検討した結果、上記式（２）における応力テンソルΠの構成方程式を、フィラー部分とゴム部分とに分けて考え、Ｍａｘｗｅｌｌモデルに基づいて下記の式（５）および（６）で書き換えることにより、粘弾性体の弾性的性質を加味できることを見出した。 ここで、Πは連続体の応力テンソル、Ｄは連続体の歪み速度テンソル、Ｎｂはゴム部分の粘弾性Ｍａｘｗｅｌｌ要素数、Ｎｆはフィラー部分の粘弾性Ｍａｘｗｅｌｌ要素数、ηはＭａｘｗｅｌｌ要素における粘性要素の粘度、τはＭａｘｗｅｌｌ要素の緩和時間であり、η／τはＭａｘｗｅｌｌ要素における弾性要素の弾性率を表す。また、∇は、上対流時間微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ Ｔｉｍｅ Ｄｅｒｉｖａｔｉｖｅ）であり、である。よって、である。また、Ｄは、である。上記式（５）および（６）は、Ｍａｘｗｅｌｌモデルに基づいたものであるが、原理的には、非線形性を含むＭａｘｗｅｌｌモデル以外のモデルに基づいた粘弾性構成則を導入することも可能である。また、φは、フィラーの位置で１を返しフィラー以外の位置で０を返す連続的な関数である。 上記式（２）における応力テンソルΠの構成方程式を、フィラー部分とゴム部分とに分けて記載するに当たり、その分け方として幾つかの方法が考えられる。発明者らは、上記（５）および（６）式を導出するに至る前に、様々な式について検討した。例えば、応力テンソルΠの構成方程式を、下記の式（１１）および（１２）のように表現することもできる。 しかし、上記式（１１）および（１２）を用いると、式（５）および（６）と比べて、移動するフィラーの界面近傍において、粘弾性緩和の特性に不均一が生じて結果が悪化する上に、計算の収束性も悪化してシミュレーションを行うことが困難になることが分かった。そこで、こうした問題を生じることのない式について鋭意検討した結果、上記式（５）および（６）のように、Ｄのみにφを考慮した式を導出するに至ったのである。これらの式の妥当性は、後述するフィラー効果を検証することにより示される。 ここで、φとして、フィラーが１つで球形である場合、例えば、あるいは、を用いることができる。ここで、ｒはフィラー中心からの距離、Ｒはフィラー半径、ａはフィラーとゴム部分を接続する界面の厚さ、δｘは離散格子の空間刻みである。 また、フィラーが球形でない場合には、上記式（１３）や（１４）を用いる代わりに、フィラーの形状を適切に表現した関数を用いるか、あるいはフィラーの位置においては１を返しフィラー以外の位置で０を返すように構成することもできる。このようにして、フィラーが球形でない場合にも対応することができる。 こうして、溶媒としてのゴム中でフィラーを粘弾性特性の大きく異なるドメインとして表現した流体の運動方程式（６）を導出することができた。この式（６）と、式（２）、（４）および（５）（特許請求の範囲における式（Ａ）に対応する）を連立させて解くことにより、粘弾性体・フィラーの動的挙動をシミュレーションすることができる。式（６）においては、前記応力Πの構成方程式の弾性成分に対する粘性成分の比が、フィラーが存在する位置でフィラーの形状に合わせて変化するため、フィラーとゴムとの境界条件をあらわに扱わずに計算コストを下げ、フィラー分散材料の動的挙動のシミュレーションが可能となるのである。また、粘弾性体であるゴムからの応力とフィラー間の相互作用力によって、フィラーの重心の位置が変化するため、フィラー同士の相対位置が変化した場合や、流れるような大きな変形を受けた場合についてシミュレーションできるようになる。以下、図２のフローチャートを参照して、本発明のシミュレーション方法の各ステップを説明する。ここでも、粘弾性体がゴムの場合を例に説明するが、これに限定されない。 まず、ステップＳ１において、ゴムおよびフィラーの物性値、並びにゴム中におけるフィラーの初期配置を含む、動的挙動を求めるための計算条件を入力する。ここで、ゴムおよびフィラーの物性値とは、具体的には、ゴムの粘度、弾性率、その周波数特性、フィラーのサイズ、個数、粘度、弾性率、その周波数特性である。 また、ゴム・フィラー複合体の動的挙動をシミュレーションするための計算条件とは、具体的には、フィラーの配置の更新回数や時間および空間のメッシュ、システムに与える周期せん断ひずみの大きさおよび周波数等である。フィラーの配置は、ランダムに重心の座標を指定することにより配置する。 次いで、ステップＳ２において、ステップＳ１にて入力された計算条件の下で、溶媒であるゴムが従う流体の運動方程式（６）を解く。具体的には、まず、モデルポテンシャルＵ（ｒ）を用いて、各フィラーに働く力Ｆｉを計算する。このモデルポテンシャルＵ（ｒ）としては、例えば、下記の式（１５）に示すレナード・ジョーンズ（ＬＪ）型のポテンシャルを用いることができる。ここで、εおよびａはフィラー間の相互作用の大きさを決定するパラメータである。 こうしたモデルポテンシャルＵ（ｒ）を用いることにより、各粒子に働く力Ｆｉは、Ｕ（ｒ）の勾配を計算することにより（Ｕ（ｒ）をフィラーの位置ベクトルｒｉで偏微分することにより）求めることができる。ここで、Ｎはフィラーの数であり、「力Ｆｉ」は、ベクトル量である。こうして求められた各フィラーに働く力Ｆｉに、i番目のフィラーのφi（ｒ）を掛け合わせて全てのフィラーについて足し合わせることにより、流体の運動方程式（６）における力場Ｆpが得られる。ここで、Ｖｉは、要素ｉの体積である。この力場Ｆ（ｒ）を式（６）に代入して式（６）を解くことにより、各フィラーの速度ｖｉが求められる。また、「速度ｖｉ」は、ベクトル量である。 続いて、ステップＳ３において、上記流体の運動方程式の解に基づいてフィラーの移動量を算出する。具体的には、ステップＳ２において運動方程式（６）を解いて算出された各フィラーの速度ｖｉを下記式（１８）に代入することにより、Δｔ後（例えば、０．０１秒後）の各フィラーの移動量Δｒｉを算出する。ここで、「移動量Δｒｉ」はベクトル量である。 その後、ステップＳ４において、得られたフィラーの移動量を用いてフィラーの配置を更新する。すなわち、ステップＳ３において得られた、Δｔ後の各フィラーの移動量Δｒｉを、各フィラーの座標に足し合わせることにより、各フィラーの配置を更新する。上記ステップＳ２〜Ｓ４を繰り返し行うことにより、ゴム・フィラー複合体の動的挙動をシミュレーションできる。 粘弾性を評価する際には、システムに指定の周期せん断の境界条件を設定して上記連立方程式（２）、（４）〜（６）を解く。境界条件としては、実際の速度ｖからせん断速度を引いた量に、例えば、Ｌｅｅｓ−Ｅｄｗａｒｄｓ法による周期境界条件を課す。フィラーが境界を跨ぐ際にも、その位置、速度をせん断速度に合わせて周期境界条件を課して計算を行う。 最後に、ステップＳ５にて、ステップＳ２〜Ｓ４の処理を所定の繰り返し回数だけ行ったかを判定し（例えば、Δｔ＝０．０１秒で１００秒のシミュレーションを行う場合には１００００回）、行っている場合には、得られたゴム・フィラー複合体の動的挙動に関するデータを出力して終了し、行っていない場合にはステップＳ２に戻る。ここで、ゴム・フィラー複合体の動的挙動に関するデータは、各繰り返し回におけるフィラーの配置や速度等のデータを含んでいる。 こうして、ステップＳ２〜Ｓ４の処理を繰り返し行うことにより、ゴム・フィラー複合体の動的挙動をシミュレーションすることができる。 ここで、式（６）におけるゴムの粘度ηおよび弾性率τとフィラー効果との関係について調べる。図３は、ゴム・フィラー複合体の粘度および弾性率とフィラー効果との関係を示しており、図３（ａ）は、Ｒ１＝（ηｆ／τｆ）／（ηｂ／τｂ）とフィラー効果との関係、図３（ｂ）は、Ｒ２=τｆ／τｂとフィラー効果との関係を示している。本発明において、「フィラー効果」とは、フィラーにより理論的に期待される、貯蔵弾性率および損失弾性率の増分に対する、貯蔵弾性率および損失弾性率の実際の増分（シミュレーションにより得られた増分）を意味しており（理論値＝１＋２．５ｇ、ｇは充填率）、Ｒ１はフィラー部、ゴム部におけるそれぞれのＭａｘｗｅｌｌ要素における弾性要素の弾性率の比を、ηｆ／τｆおよびηｂ／τｂは、それぞれフィラー部におけるＭａｘｗｅｌｌ要素における弾性要素の弾性率、ゴム部におけるＭａｘｗｅｌｌ要素における弾性要素の弾性率を意味している。また、上記結果は、ゴム・フィラー複合体に５％のひずみ、および周波数０．０１の振動が与えられた条件下で得られたものである。シミュレーションの数値は、密度ρ=１、粘度η＝０．５、フィラー半径Ｒ＝３．８となる条件で無次元化した。 これらの図から明らかなように、Ｒ１またはＲ２の値の増加とともに、フィラー効果は単調に増加する。すなわち、Ｒ１またはＲ２の値は大きいほどフィラー効果が大きいことを意味している。例えば、Ｒ１を５０以上、Ｒ２を１００以上とすることにより、９０％以上のフィラー効果が得られ、流体の運動方程式中において剛体としてのフィラーを表現できていると考えることができる。 以上の本発明の粘弾性体・フィラー複合体の動的挙動のシミュレーション方法において、データ出力ステップにて、求められたデータに基づいて粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに行うことにより、粘弾性体・フィラー複合体の粘度や弾性率等の物性値を求めることができる。粘弾性体としてゴムを例に説明すると、具体的には、ゴム・フィラー複合体の物性値は、動的粘弾性の貯蔵弾性率Ｇ’および損失弾性率Ｇ’’として求める。このＧ’およびＧ’’は以下のように求める。すなわち、大きさ（振幅）γ０、周波数ωである周期せん断ひずみγ０ｓｉｎ（ωt）を、システムのある方向（ｘ方向）にせん断流れを、その方向に直交する方向（ｙ方向）に勾配を付けて作用させて、その時の平均の応力σを次式（１９）から計算することにより求めることができる。 この応力がとなった時の、Ａ／γ０、Ｂ／γ０がそれぞれ、Ｇ’、Ｇ’’である。これらＧ’およびＧ’’は、完全弾性体ではＧ’≠０、Ｇ’’＝０であり、完全粘性体ではＧ’＝０、Ｇ’’≠０である。また、ｔａｎδ＝Ｇ’’／Ｇ’が小さいほど、弾性体に近い粘弾性材料を意味する。 こうして、従来のＦＰＤ法では実質的に不可能だった、粘弾性体・フィラー複合体のマクロな物性を求めて評価することができ、ゴム等の粘弾性材料の設計に有効利用することができる。（シミュレーション装置） 次に、本発明の粘弾性体・フィラー複合体の動的挙動のシミュレーション装置について説明する。図４は、本発明の粘弾性体・フィラー複合体の動的挙動のシミュレーション装置の一実施形態を示す図である。この装置１は、入力部１１と、解析部１２と、表示部１３とを備える。また、解析部１２は、演算器２１と、移動量算出器２２と、配置更新器２３と、データ出力器２４とを有する。以下、粘弾性体がゴムの場合を例に説明するが、これに限定されない。 入力部１１は、ゴムおよびフィラーの物性値、並びにゴム中におけるフィラーの初期配置を含む、動的挙動をシミュレーションするための計算条件を入力する。この入力部１１としては、例えばキーボードやマウス等を用いることができる。 解析部１２は、入力部１１により入力されたゴムおよびフィラーの物性値、並びに動的挙動を求めるための計算条件下でゴムが従う流体の運動方程式である式（６）を解いて、フィラーの移動量を算出することにより、ゴム・フィラー複合体の動的挙動をシミュレーションする。以下、解析部１２の各構成について説明する。 演算器２１は、入力部１１により入力されたフィラーの配置の下で、ゴムが従う流体の運動方程式を解く。ここで、ゴムの弾性的性質を加味された流体の運動方程式である式（６）は、記憶装置（図示せず）に予め格納しておき、流体の運動方程式を解く際に記憶装置から読み出されるように構成されている。 移動量算出器２２は、演算器２１により得られた流体の運動方程式の解に基づいて、フィラーの移動量を算出する。 配置更新器２３は、移動量算出器２２により算出されたフィラーの移動量を用いて前記フィラーの配置を更新する。 データ出力器２４は、演算器２１、移動量算出器２２および配置更新器２３により求められたゴム・フィラー複合体の動的挙動に関するデータを出力する。 上記演算器２１、移動量算出器２２、配置更新器２３およびデータ出力器２４からなる解析部１２は、パーソナルコンピュータ中のＣＰＵやメモリ、記憶装置等により構成することができる。 表示部１３は、データ出力器２４により出力されたゴム・フィラー複合体に関するデータを表示する。この表示部１３としては、例えばディスプレイを用いることができる。 また、ゴムが従う流体の運動方程式（６）において、図３に示した、Ｒ＝（ηｆ／τｆ）／（ηｂ／τｂ）とフィラー効果との関係、およびτｆ／τｂとフィラー効果との関係から、Ｒを５０以上、τｆ／τｂを１００以上とすることにより、９０％以上のフィラー効果が得られ、流体の運動方程式中において剛体としてのフィラーを表現できる。 こうして、本発明のシミュレーション装置により、粘弾性体・フィラー複合体の動的挙動を簡便にシミュレーションすることができる。 以上の本発明のシミュレーション装置において、解析部１２は、上記演算器２１、移動量算出器２２、配置更新器２３により求められた粘弾性体・フィラー複合体の動的挙動に関するデータに基づいて、粘弾性体・フィラー複合体の粘弾性特性を解析する粘弾性解析器をさらに有することができる。こうして、従来のシミュレーション装置では実質的に不可能だった、粘弾性体・フィラー複合体のマクロな物性を求めて評価することができ、ゴム等の粘弾性材料の設計に有効利用することができる。（シミュレーションプログラム） 上述した本発明のシミュレーション装置１として、コンピュータを好適に用いることができ、そのようなコンピュータは、図４に示した手段の各機能を実現する処理内容を記述したプログラムを、当該コンピュータの記憶部に格納しておき、当該コンピュータの中央演算処理装置（ＣＰＵ）によってこのプログラムを読み出して実行させることで実現することができる。こうして、本発明のシミュレーションプログラムをコンピュータに実行させることにより、粘弾性体・フィラー複合体の動的挙動を簡便にシミュレーションすることができる。＜シミュレーション方法の妥当性の検証＞ まず、本発明のシミュレーション方法の妥当性を検証した。ここで、粘弾性体としてゴムを想定した。液体中や弾性体中に剛体フィラーを希薄に分散させた場合、その流体効果として、全体の粘度や弾性率がフィラーを入れた数に応じて変化することがＥｉｎｓｔｅｉｎ則やＧｕｔｈ−Ｇｏｌｄ則として知られている。理論的には小さいフィラーを分散させた場合には、通常レイノルズ数が十分小さく、慣性の影響を無視したＳｔｏｋｅｓ近似で挙動を解析することが可能であり、希薄な極限として、フィラーを１つ考慮した際の結果と一致する。 粘弾性媒質中のフィラーの流体効果についても同様のことが成立すると考えられる。そこで、ゴムにフィラーを配合した際に、弾性率が上昇する効果の１つであるＧｕｔｈ−Ｇｏｌｄ則が成立していることを検証した。すなわち、上記法則が成立している場合には、フィラーの面積分率ｇに対して、となることが期待される。ここで、ωは周波数、Ｇ’およびＧ’’はゴム・フィラー複合体の粘弾性体の粘弾性係数であり、Ｇ’は貯蔵弾性率、Ｇ’’は損失弾性率である。また、Ｇ’ｂｕｌｋおよびＧ’’ｂｕｌｋはフィラーがない場合のゴムの粘弾性係数である。 図５（ａ）は、本発明のシミュレーション方法による、フィラーの充填率とＧｕｔｈ−Ｇｏｌｄ因子Ａとの関係を示す図である。ここで、この関係は、ゴム・フィラー複合体に対して１０％のひずみ、周波数０．０１の振動を与えることにより得られたものである。数値は、密度ρ＝１、粘度η＝０．５、フィラー半径Ｒ＝３．８となる条件で無次元化した。この図から明らかなように、フィラーの充填率の増大とともに、Ｇｕｔｈ−Ｇｏｌｄ因子が増大しており、ゴム・フィラー複合体の応力が増大していることが分かる。 また、図５（ｂ）は、本発明のシミュレーション方法による、ゴム・フィラー複合体に与える周波数とＧｕｔｈ−Ｇｏｌｄ因子および位相シフトとの関係を示している。振幅Ａと位相シフトΔは、フィラーがない場合の応力σ０（ｔ）に対して、フィラーがある場合の応力σ（ｔ）を、σ０（ｔ）＝Ａσ０（ｔ＋Δ／ω）でフィッティングして得られた量である。ここで、ゴム・フィラー複合体に与えられた歪みは１０％であり、フィラーの面積分率ｇは約８％、τｂは１００である。まず、Ｇｕｔｈ−Ｇｏｌｄ因子と周波数との関係を見ると、０．０２程度まではほとんど周波数に対する依存性は見られないが、０．０３を超えると急激に増大することが分かる。また、位相シフトは、周波数の増大とともに単調に増大することが分かる。 この周波数依存性は、フィラー部の慣性項が影響していると考えられる。ゴムなどの場合、慣性の影響は小さく、実験的に得られる周波数をシミュレーションするには、周波数をかなり小さくして長時間のシミュレーションを行う必要がある。レイノルズ数をおよそρＲ２ω／η≒３０ωと考えると、この条件下では、３０ω＜＜１のときに周波数依存性は十分小さく、ゴム・フィラー複合体の粘弾性を求めることができていることが分かる。＜ゴム・フィラー複合体の粘弾性特性の評価＞ 本発明のシミュレーション方法（発明例）および非特許文献１のＦＰＤ法（比較例）により、ゴム・フィラー複合体の粘弾性特性を調べた。具体的には、２次元のモデルに対してフィラーを適度に凝集させ、２０％の面積分率で配置、フィラー間にＬＪ相互作用をε＝０．５として計算した。周波数０．０１に対してひずみ３％、１０％、３０％での、粘度および弾性率を求めた。数値は、密度ρ=１、粘度η＝０．５、フィラー半径Ｒ＝３．８となる条件で無次元化した。また、Ｒ１＝５０、Ｒ２＝１００である。得られた結果を図６に示す。この図から明らかなように、貯蔵弾性率Ｇ’にはフィラー以外の寄与として粘弾性体からの寄与が含まれ、Ｇ’’もそれに合わせて変化しており、発明例においては、ゴム・フィラー複合体の粘弾性特性が正しく求められていることが分かる。１ シミュレーション装置１１ 入力部１２ 解析部１３ 表示部２１ 演算器２２ 移動量算出器２３ 配置更新器２４ データ出力器Ｓ１ 入力ステップＳ２ 演算ステップＳ３ 移動量算出ステップＳ４ 配置更新ステップＳ５ データ出力ステップ 粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動のシミュレーション方法であって、ａ）前記粘弾性体および前記フィラーの物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力ステップと、ｂ）入力された前記計算条件の下で前記粘弾性体が従う流体の運動方程式を解く演算ステップと、ｃ）前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出ステップと、ｄ）算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新ステップと、ｅ）前記ステップｂ）〜ｄ）の繰り返し処理により求めた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップと、を含み、前記粘弾性体が従う流体の運動方程式は、前記粘弾性体の弾性的性質を加味したものであることを特徴とするシミュレーション方法。 前記粘弾性体が従う流体の運動方程式は以下の式（Ａ）で表される、請求項１に記載のシミュレーション方法。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆｐは外場、Πは連続体の応力テンソル、Πｊ（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Πｊ（ｆ）は連続体の応力テンソルのフィラー部分の成分、Ｎｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、τｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、ηｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、ηｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。 前記式（Ａ）における粘度ηｊ（ｆ）およびηｊ（ｂ）、緩和時間τｊ（ｆ）およびτｊ（ｂ）は下記の式（Ｂ）を満足する、請求項２に記載のシミュレーション方法。 前記データ出力ステップは、前記求められたデータに基づいて前記ゴム・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに含む、請求項１〜３のいずれか一項に記載のシミュレーション方法。 粘弾性体中にフィラーが分散したゴム・フィラー複合体の動的挙動をシミュレーションする装置であって、 前記粘弾性体および前記フィラーの物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力部と、 入力された前記計算条件の下で前記粘弾性体が従う流体の運動方程式を解く演算器と、前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出器と、算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新器と、前記演算器、前記移動量算出器および前記配置更新器により求められた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力器とを有する解析部と、を備え、前記粘弾性体が従う流体の運動方程式は、前記粘弾性体の弾性的性質を加味したものであることを特徴とするシミュレーション装置。 前記粘弾性体が従う流体の運動方程式は以下の式（Ａ）で表される、請求項５に記載のシミュレーション装置。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆｐは外場、Πは連続体の応力テンソル、Πｊ（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Πｊ（ｆ）は連続体の応力テンソルのフィラー部分の成分、Ｎｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、τｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、ηｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、ηｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。 前記式（Ａ）における粘度ηｊ（ｆ）およびηｊ（ｂ）、緩和時間τｊ（ｆ）およびτｊ（ｂ）は下記の式（Ｂ）を満足する、請求項６に記載のシミュレーション装置。 前記解析部は、前記得られたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析器をさらに有する、請求項５〜７のいずれか一項に記載のシミュレーション装置。 粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動を求める装置として構成するコンピュータに、ａ）前記粘弾性体および前記フィラーの物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力ステップと、ｂ）入力された前記計算条件の下で前記粘弾性体が従う流体の運動方程式を解く演算ステップと、ｃ）前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出ステップと、ｄ）算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新ステップと、ｅ）前記ステップｂ）〜ｄ）の繰り返し処理により求めた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップと、を実行させ、前記流体の運動方程式は、前記粘弾性体の弾性的性質を加味したものであることを特徴とするシミュレーションプログラム。 前記粘弾性体が従う流体の運動方程式は以下の式（Ａ）で表される、請求項９に記載のシミュレーションプログラム。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆｐは外場、Πは連続体の応力テンソル、Πｊ（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Πｊ（ｆ）は連続体の応力テンソルのフィラー部分の成分、Ｎｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、τｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、ηｊ（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、ηｊ（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。 前記式（Ａ）における粘度ηｊ（ｆ）およびηｊ（ｂ））、緩和時間τｊ（ｆ）およびτｊ（ｂ）は下記の式（Ｂ）を満足する、請求項１０に記載のシミュレーションプログラム。 前記データ出力ステップは、前記求められたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに含む、請求項９〜１１のいずれか一項に記載のシミュレーションプログラム。 【課題】ゴム等の粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の簡便な動的挙動のシミュレーション方法、シミュレーション装置およびシミュレーションプログラムを提供する。【解決手段】ａ）粘弾性体およびフィラーの物性値、並びに粘弾性体中におけるフィラーの初期配置を含む、動的挙動をシミュレーションするための計算条件を入力する入力ステップと、ｂ）入力された計算条件の下で粘弾性体が従う流体の運動方程式を解く演算ステップと、ｃ）流体の運動方程式の解に基づいてフィラーの移動量を算出する移動量算出ステップと、ｄ）算出されたフィラーの移動量を用いてフィラーの配置を更新する配置更新ステップと、ｅ）ステップｂ）〜ｄ）の繰り返し処理により求めた粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップとを含み、粘弾性体が従う流体の運動方程式は、粘弾性体の弾性的性質を加味したものであることを特徴とする。【選択図】図２