JP2014164555A - 粘弾性体・フィラー複合体の動的挙動のシミュレーション方法、シミュレーション装置およびシミュレーションプログラム - Google Patents

Abstract

【課題】ゴム等の粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の簡便な動的挙動のシミュレーション方法、シミュレーション装置およびシミュレーションプログラムを提供する。
【解決手段】粘弾性体・フィラー複合体のモデルに基づいてシミュレーションを行うための、粘弾性体、フィラー、および拘束層の物性値、並びにシミュレーションのための計算条件を入力する入力ステップと、入力された計算条件の下で流体の運動方程式を解く演算ステップと、流体の運動方程式の解に基づいてフィラーの移動量を算出する移動量算出ステップと、算出されたフィラーの移動量を用いてフィラーの配置を更新する配置更新ステップと、粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップとを含み、粘弾性体および拘束層が従う流体の運動方程式は、粘弾性体および拘束層の弾性的性質を加味する。
【選択図】図２

Description

本発明は、シミュレーション方法、シミュレーション装置およびシミュレーションプログラムに関し、より詳細には、粘弾性体中にフィラーが分散した材料（以下、「粘弾性体・フィラー複合体」と称する）の動的挙動を簡便にシミュレーションできるシミュレーション方法、シミュレーション装置およびシミュレーションプログラムに関するものである。

従来、粘弾性体、例えばゴムにカーボンブラックやシリカ等のフィラーを配合すると補強効果があることが知られており、ゴムにフィラーを配合した材料（以下、「ゴム・フィラー複合体」と称する）が各種タイヤやコンベヤベルト等のゴム製品に適用されている。このゴム・フィラー複合体を構成するフィラー部分は、複数のフィラーが連なった複雑なネットワーク構造を形成しており、ゴム・フィラー複合体において、ある程度偏ったフィラー分布を示すことが知られている。

このような複雑なフィラー構造を有するゴム・フィラー複合体の高い性能を実現するためには、変形が与えられた際の、ゴム・フィラー複合体の動的挙動を把握することが肝要である。そこでこれまで、こうしたゴム・フィラー複合体の動的挙動をシミュレーションするための様々な手法が提案されてきた。

例えば、特許文献１には、有限要素法（Ｆｉｎｉｔｅ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ，ＦＥＭ）を用いて、ゴム製品の粘弾性応答性能を予測する方法について記載されている。
また、特許文献２には、粗視化分子動力学法を用いて、ゴムとカーボンブラックとを少なくとも含むゴム材料を評価する技術について記載されている。

しかし、特許文献１に記載された方法では、固体表面の境界に対して、ノンスリップ条件等の適切な境界条件を満足させるような計算が必要となり、ゴム・フィラー複合体のように境界が多い材料においては、計算コストが膨大になる問題がある。また、ＦＥＭ等においては、解析対象を小さな要素に分割するためのメッシュを作成するが、このメッシュをフィラーの形状に合わせて設定する必要があり、フィラーの位置が時々刻々変化する場合には、そのメッシュの再作成の計算コストが膨大となることも問題である。

また、特許文献２に記載された方法では、ＦＥＭにおいて必要となる境界条件を考慮する必要はないが、マクロな物性を評価するためには、膨大な計算コストがかかり、ゴム・フィラー複合体のような大きな系のマクロな物性を評価することは現実的には不可能である。

一方、非特許文献１には、液体の溶媒中にコロイド粒子が分散したコロイド分散系において、コロイド粒子を剛体ではなく粘性の高い液体ドメインと見なして、コロイド分散系の動的挙動を計算する、いわゆる流体粒子ダイナミクス（Ｆｌｕｉｄ Ｐａｒｔｉｃｌｅ Ｄｙｎａｍｉｃｓ，ＦＰＤ）法が提案されている。このＦＰＤ法においては、コロイド粒子の高粘性ドメインと周囲の粘性液体の物性とが滑らかに接続されるため、境界条件をあらわに解くことなく、流体効果を考慮したコロイド分散系の動的挙動を高効率にシミュレーションすることが可能である。

さらに、特許文献３には、非特許文献１に記載された方法を発展させて、荷電粒子が液体の溶媒中に分散した荷電コロイド分散系の動的挙動をシミュレーションする方法について記載されており、電気泳動ディスプレイにおける電極間の帯電顔料粒子のコロイド分散液において生じるコロイド電気泳動現象のシミュレーションが行われている。

このように、ＦＰＤ法により、コロイド分散系の動的挙動を高効率にシミュレーションできていることから、上記ＦＰＤ法をゴム・フィラー複合体の動的挙動のシミュレーションに適用して、高性能なゴム・フィラー複合体の設計に役立たせることが期待される。

特開２００５−４９３３３号公報 特開２００６−６４６５８号公報 特開２００５−３０８７８２号公報

Ｈ．Ｔａｎａｋａ ａｎｄ Ｔ．Ａｒａｋｉ，Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．８５（６），１３３８（２０００）

しかしながら、非特許文献１に記載されたＦＰＤ法においては、溶媒として粘性の液体を用いることが想定されているため、ＦＰＤ法の、ゴム・フィラー複合体のような、溶媒が粘弾性体の場合に対する適用は一切検討されて来なかった。そのため、どのような条件の下で粘弾性体・フィラー複合体の動的挙動のシミュレーションを行うことができるのか不明となっていた。
そこで、本発明の目的は、ゴムをはじめとする粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動を簡便にシミュレーションできるシミュレーション方法、シミュレーション装置およびシミュレーションプログラムを提供することにある。

発明者は、上記課題を解決する方途について鋭意検討した。上記ＦＰＤ法において、溶媒が従う流体の運動方程式であるナビエ・ストークス方程式は、液体の溶媒および固体であるコロイド粒子の粘度をパラメータとして含み、コロイド粒子の粘度を液体である溶媒の粘度よりも高く設定することにより、粒子内に液体の流れが入り込まず、コロイド粒子の境界条件が加味される。しかしながら、上記ＦＰＤ法においては、液体の溶媒中に剛体粒子であるフィラーを充填させても、フィラーによるゴム弾性の向上効果を評価できないことが判明した。そこで、この問題を解決するためには、溶媒である粘弾性体の弾性的な特性を加味した式を用いてシミュレーションを行う必要があるとの認識し、本発明を完成させるに至った。

すなわち、本発明の要旨構成は以下の通りである。
（１）粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動のシミュレーション方法であって、ａ）前記粘弾性体、前記フィラー、および前記粘弾性体と前記フィラーとを接続する拘束層からなる前記粘弾性体・フィラー複合体のモデルに基づいて前記シミュレーションを行うための、前記粘弾性体、前記フィラー、および前記拘束層の物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力ステップと、ｂ）入力された前記計算条件の下で前記粘弾性体および前記拘束層が従う流体の運動方程式を解く演算ステップと、
ｃ）前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出ステップと、ｄ）算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新ステップと、ｅ）前記ステップｂ）〜ｄ）の繰り返し処理により求めた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップとを含み、前記粘弾性体および前記拘束層が従う流体の運動方程式は、前記粘弾性体および拘束層の弾性的性質を加味したものであることを特徴とするシミュレーション方法。

（２）前記粘弾性体および前記拘束層が従う流体の運動方程式は以下の式（Ａ）で表される、前記（１）に記載のシミュレーション方法。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆ_ｐは外場、Πは連続体の応力テンソル、Π_ｊ ^（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Π_ｊ ^（ｆ）は連続体の応力テンソルのフィラー部分の成分、Π_ｊ ^（Ｉ）は連続体の応力テンソルの拘束層部分の成分、Ｎ_ｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_ｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_Ｉは前記拘束層における粘弾性Ｍａｘｗｅｌｌ要素の数、τ_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素の緩和時間、η_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φ_ｂは前記粘弾性体の位置にて１を返し前記粘弾性体以外の位置で０を返す関数、φ_ｆは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、φ_Ｉは前記拘束層の位置にて１を返し前記拘束層以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。

（３）前記式（Ａ）における粘度η_ｊ ^（ｂ）、η_ｊ ^（ｆ）およびη_ｊ ^（Ｉ）、緩和時間τ_ｊ ^（ｂ）、τ_ｊ ^（ｆ）およびτ_ｊ ^（Ｉ）は下記の式（Ｂ）のいずれかを満足する、前記（２）に記載のシミュレーション方法。

（４）前記データ出力ステップは、前記求められたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに含む、前記（１）〜（３）のいずれか一項に記載のシミュレーション方法。

（５）粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動をシミュレーションする装置であって、前記粘弾性体、前記フィラー、および前記粘弾性体と前記フィラーとを接続する拘束層からなる前記粘弾性体・フィラー複合体のモデルに基づいて前記シミュレーションを行うための、前記粘弾性体、前記フィラーおよび前記拘束層の物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力部と、入力された前記計算条件の下で前記粘弾性体および前記拘束層が従う流体の運動方程式を解く演算器と、前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出器と、算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新器と、前記演算器、前記移動量算出器および前記配置更新器により求められた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力器とを有する解析部とを備え、前記粘弾性体および前記拘束層が従う流体の運動方程式は、前記粘弾性体および拘束層の弾性的性質を加味したものであることを特徴とするシミュレーション装置。

（６）前記粘弾性体および前記拘束層が従う流体の運動方程式は以下の式（Ａ）で表される、前記（５）に記載のシミュレーション装置。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆ_ｐは外場、Πは連続体の応力テンソル、Π_ｊ ^（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Π_ｊ ^（ｆ）は連続体の応力テンソルのフィラー部分の成分、Π_ｊ ^（Ｉ）は連続体の応力テンソルの拘束層部分の成分、Ｎ_ｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_ｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_Ｉは前記拘束層における粘弾性Ｍａｘｗｅｌｌ要素の数、τ_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素の緩和時間、η_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の粘性要素の粘度、η_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φ_ｂは前記粘弾性体の位置にて１を返し前記粘弾性体以外の位置で０を返す関数、φ_ｆは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、φ_Ｉは前記拘束層の位置にて１を返し前記拘束層以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。

（７）前記式（Ａ）における粘度η_ｊ ^（ｂ）、η_ｊ ^（ｆ）およびη_ｊ ^（Ｉ）、緩和時間τ_ｊ ^（ｂ）、τ_ｊ ^（ｆ）およびτ_ｊ ^（Ｉ）は下記の式（Ｂ）のいずれかを満足する、前記（６）に記載のシミュレーション装置。

（８）前記解析部は、前記得られたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析器をさらに有する、前記（５）〜（７）のいずれか一項に記載のシミュレーション装置。

（９）粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動をシミュレーションする装置として構成するコンピュータに、ａ）前記粘弾性体、前記フィラー、および前記粘弾性体と前記フィラーとを接続する拘束層からなる前記粘弾性体・フィラー複合体のモデルに基づいて前記シミュレーションを行うための、前記粘弾性体、前記フィラー、および前記拘束層の物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力ステップと、ｂ）入力された前記計算条件の下で前記粘弾性体および前記拘束層が従う流体の運動方程式を解く演算ステップと、ｃ）前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出ステップと、ｄ）算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新ステップと、ｅ）前記ステップｂ）〜ｄ）の繰り返し処理により求めた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップとを実行させ、前記流体の運動方程式は、前記粘弾性体および前記拘束層の弾性的性質を加味したものであることを特徴とするシミュレーションプログラム。

（１０）前記粘弾性体および前記拘束層が従う流体の運動方程式は以下の式（Ａ）で表される、前記（９）に記載のシミュレーションプログラム。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆ_ｐは外場、Πは連続体の応力テンソル、Π_ｊ ^（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Π_ｊ ^（ｆ）は連続体の応力テンソルのフィラー部分の成分、Π_ｊ ^（Ｉ）は連続体の応力テンソルの拘束層部分の成分、Ｎ_ｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_ｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_Ｉは前記拘束層における粘弾性Ｍａｘｗｅｌｌ要素の数、τ_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素の緩和時間、η_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φ_ｂは前記粘弾性体の位置にて１を返し前記粘弾性体以外の位置で０を返す関数、φ_ｆは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、φ_Ｉは前記拘束層の位置にて１を返し前記拘束層以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。

（１１）前記式（Ａ）における粘度η_ｊ ^（ｂ）、η_ｊ ^（ｆ）およびη_ｊ ^（Ｉ）、緩和時間τ_ｊ ^（ｂ）、τ_ｊ ^（ｆ）およびτ_ｊ ^（Ｉ）は下記の式（Ｂ）のいずれかを満足する、前記（１０）に記載のシミュレーションプログラム。

（１２）前記データ出力ステップは、前記求められたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに含む、前記（９）〜（１１）のいずれか一項に記載のシミュレーションプログラム。

本発明によれば、粘弾性体・フィラー複合体の動的挙動を簡便にシミュレーションできるシミュレーション方法、シミュレーション装置およびシミュレーションプログラムを提供することができる。

非特許文献１に記載のＦＰＤ法により求めたゴム・フィラー複合体の物性を示す図である。 本発明のシミュレーション方法の一実施形態を示すフローチャートである。 （ａ）および（ｂ）は粘度および緩和時間とフィラーの流体効果との関係を示す図であり、（ｃ）は本発明の方法により得られた弾性率と理論から得られた弾性率との関係を示す図であり、（ｄ）はφ_ｆ（ｘ）および（η_Ｉ／η_ｆ）φ_Ｉ（ｘ）の一例を示す図である。 本発明のシミュレーション装置の一実施形態を示す図である。 フィラーの充填率とＧｕｔｈ−Ｇｏｌｄ因子との関係を示す図である。 （ａ）拘束層を含むフィラーの表面形状を示す図であり、（ｂ）は、拘束層の厚さと、本発明のシミュレーション方法により得られた複素粘弾性係数との関係を示す図である。

（シミュレーション方法）
以下、図面を参照して、本発明の実施形態について説明する。
本発明の粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動のシミュレーション方法は、ａ）予め作成しておいた、粘弾性体、フィラー、および粘弾性体とフィラーとを接続する拘束層からなる粘弾性体・フィラー複合体のモデルに基づいてシミュレーションを行うための、粘弾性体、フィラー、および拘束層の物性値、並びに粘弾性体中におけるフィラーの初期配置を含む、動的挙動をシミュレーションするための計算条件を入力する入力ステップと、ｂ）入力された計算条件の下で粘弾性体および拘束層が従う流体の運動方程式を解く演算ステップと、ｃ）流体の運動方程式の解に基づいてフィラーの移動量を算出する移動量算出ステップと、ｄ）算出されたフィラーの移動量を用いてフィラーの配置を更新する配置更新ステップと、ｅ）ステップｂ）〜ｄ）の繰り返し処理により求めた粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップとを含む。ここで、粘弾性体および拘束層が従う流体の運動方程式は、粘弾性体および拘束層の弾性的性質を加味したものであることが肝要である。本発明において、「フィラー」とは、カーボンブラックやシリカ等、粘弾性体の補強効果を図るために粘弾性体に配合される充填剤であり、粒子状のものを意味している。ここで、フィラーの形状は球形に限定されない。また、本発明において、「粘弾性体・フィラー複合体」とは、粘弾性体中にフィラーが分散された材料を意味している。以下の説明においては、粘弾性体がゴムの場合を例に説明するが、それに限定されない。

本発明により粘弾性体・フィラー複合体の動的挙動をシミュレーションするに当たって、粘弾性体とフィラーとを接続し、粘弾性体やフィラーとは異なる物性を有する拘束層を導入した粘弾性体・フィラー複合体モデルを採用する。すなわち、粘弾性体、例えばポリマーにフィラーを分散させたポリマー・フィラー複合体において、フィラーとポリマーとの界面近傍の領域に存在する分子は、界面との相互作用により、バルクとは異なった物性を有すると考えられている（例えば、特開２０１０−１８１１９４号公報、G. Tsagaropoulos, and A. Eisenberg, Macromolecules vol.28, p.6067(1995)参照）。そこで、上記界面近傍の領域を拘束層として粘弾性体・フィラー複合体モデルに導入することにより、より現実に近いシミュレーションを行うことができ、シミュレーション精度の向上が期待できる。

ここで、拘束層は、フィラーの周囲に設けられ、具体的には、フィラー周部に所定の厚みを有する領域として設ける。また、拘束層の物性は、パラメータηやτの値を変更することにより調整でき、例えば粘弾性体がゴムの場合、低温で測定されたゴムの物性を再現するように、パラメータηやτを設定する。

さて、上述のように、非特許文献１に記載されたＦＰＤ法は、特許文献１に記載されたＦＥＭのような境界条件をあらわに解く必要なしに、粘性液体溶媒の流体効果を考慮したコロイド分散系の動的挙動を高効率にシミュレーションすることを可能にする。このＦＰＤ法において、溶媒である粘性液体の運動は、下記のナビエ・ストークス方程式で表される。

ここで、Ｄ／Ｄｔはラグランジュ微分（＝∂／∂ｔ＋ｖ・∇）、Ｆ_ｐは力場、ρは質量密度、ｖは速度、ｐは圧力、ηは粘度、ζは熱遥動力である。

ＦＰＤ法を、粘弾性体・フィラー複合体の動的挙動のシミュレーションに適用する場合、上記式（１）において、粘度ηの値を、溶媒として粘性液体を用いる場合よりも大きな値に設定することにより、粘弾性体のフィラーによる補強性を近似的に取り込むことができると考えられる。そこで、式（１）を用いて、粘弾性体・フィラー複合体の動的挙動をシミュレーションしてその粘弾性特性を調べた。具体的には、粘弾性体の粘度に対して熱エネルギーは非常に小さいと考えられるためゼロとし、２次元のモデルに剛体フィラーの充填率（面積分率）を０％〜７％まで変化させ、周期せん断条件の下でシミュレーションを行った。無次元化は、密度ρ=１、粘度η＝０．５、フィラー半径Ｒ＝３．８となる条件で行った。得られた結果を図１に示す。ここで、図１（ａ）はフィラー充填率と貯蔵弾性率Ｇ’との関係を、図１（ｂ）はフィラー充填率と損失弾性率Ｇ’’との関係をそれぞれ示している。ここで、貯蔵弾性率Ｇ’および損失弾性率Ｇ’’は、後述するように、動的粘弾性を特徴付ける量であり、貯蔵弾性率Ｇ’および損失弾性率Ｇ’’は、粘弾性体の弾性的性質、粘性的性質をそれぞれ特徴付ける量である。これらの図から明らかなように、フィラー充填率の増加とともに損失弾性率Ｇ’’は増加する一方、貯蔵弾性率Ｇ’は、フィラー充填率にかかわらず０であり、ＦＰＤ法では粘弾性体の弾性的性質を取り込めていないことが分かる。すなわち、ＦＰＤ法によりシミュレーションを行うと、粘弾性体にフィラーを充填しても弾性率を求めることができず、適切な結果が得られない。

このように、発明者は、従来のＦＰＤ法では、粘弾性体の動的挙動のシミュレーションを正しく行えないことを受けて、溶媒としての粘弾性体の弾性的な性質を取り込むことができる流体の運動方程式を用いてシミュレーションを行う必要があるとの認識にたどり着いた。そこで、発明者は、試行錯誤の末、以下に示す、本発明の粘弾性体・フィラー複合体の動的挙動のシミュレーション方法に用いる流体の運動方程式を導出するに至ったのである。

すなわち、まず、非圧縮性粘性流体の運動量保存則は、以下の式（２）で与えられる。

ここで、Πは連続体の応力テンソルであり、

と表され、また、ｖは、連続の条件から以下の式（４）の関係を満たす。

本発明においては、粘弾性体の弾性的性質を加味すべく鋭意検討した結果、上述のモデルに従って、上記式（２）における応力テンソルΠの構成方程式を、粘弾性体部分、フィラー部分および拘束層部分に分けて考え、Ｍａｘｗｅｌｌモデルに基づいて下記の式（５）および（６）で書き換える。これにより、粘弾性体の弾性的性質を加味することができる。

ここで、Π_ｊ ^（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Π_ｊ ^（ｆ）は連続体の応力テンソルのフィラー部分の成分、Π_ｊ ^（Ｉ）は連続体の応力テンソルの拘束層部分の成分、Ｄは連続体の歪み速度テンソル、Ｎ_ｂは粘弾性体部分の粘弾性Ｍａｘｗｅｌｌ要素数、Ｎ_ｆはフィラー部分の粘弾性Ｍａｘｗｅｌｌ要素数、Ｎ_Ｉは拘束層部分の粘弾性Ｍａｘｗｅｌｌ要素数、τ_ｊ ^（ｂ）は粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（ｆ）はフィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（Ｉ）は拘束層におけるＭａｘｗｅｌｌ要素の緩和時間、η_ｊ ^（ｂ）は粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（ｆ）はフィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（Ｉ）は拘束層におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η／τはＭａｘｗｅｌｌ要素における弾性要素の弾性率を表す。また、∇は、上対流時間微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ Ｔｉｍｅ Ｄｅｒｉｖａｔｉｖｅ）であり、

である。よって、

である。また、Ｄは、

である。また、φ_ｂは、粘弾性体の位置で１を返し粘弾性体以外の位置で０を返す連続的な関数、φ_ｆは、フィラーの位置で１を返しフィラー以外の位置で０を返す連続的な関数、およびφ_Ｉは、拘束層の位置で１を返し拘束層以外の位置で０を返す連続的な関数である。

上記式（５）および（６）は、Ｍａｘｗｅｌｌモデルに基づいたものであるが、原理的には、非線形性を含むＭａｘｗｅｌｌモデル以外のモデルに基づいた粘弾性構成則を導入することも可能である。

ここで、式（Ａ）におけるφ_ｂ、φ_ｆおよびφ_Ｉとしては、例えば、

を用いることができ。ここで、Ｒはフィラーの半径、ａはフィラーと拘束層および拘束層とバルクを接続する界面厚み、ｂは拘束層自身の厚みである。

特に、η_ｊ ^（Ｉ）はη_ｊ ^（ｆ）/η_ｊ ^（ｂ）が一定で、τ_ｊ ^（ｆ）=τ_ｊ ^（Ｉ）と近似できる特別の場合には、

として、式（５）および（６）の計算においてφ_ｆおよびφ_Iをまとめて表現し、計算式を簡略化でき、計算工数を低減して計算効率を向上させることができる。例えば、

とまとめて表現することも可能である。この場合には、フィラー半径Ｒに対して、界面厚みａで接続した、厚みｂの部分が拘束層の役割を果たす。

また、フィラーが球形で原点にある場合、φとしては、例えば、

を用いることができ、ここで、δｘは、空間離散化のメッシュサイズである。あるいは、

と表すこともできる。

また、フィラーが球形でない場合には、上記式（１４）〜（１７）を用いる代わりに、フィラーの形状を適切に表現した関数を用いるか、あるいはフィラーの位置においては１を返しフィラー以外の位置で０を返すように構成することもできる。このようにして、フィラーが球形でない場合にも対応することができる。

こうして、溶媒としての粘弾性体中でフィラーを粘弾性特性の大きく異なるドメインとして表現した流体の運動方程式（６）を導出することができた。この式（６）と、式（２）、（４）および（５）（特許請求の範囲における式（Ａ）に対応する）を連立させて解くことにより、粘弾性体・フィラーの動的挙動をシミュレーションすることができる。式（６）においては、前記応力Πの構成方程式の弾性成分に対する粘性成分の比が、フィラーが存在する位置でフィラーの形状に合わせて変化するため、フィラーと粘弾性体との境界条件をあらわに扱わずに計算コストを下げ、粘弾性体・フィラー複合体の動的挙動のシミュレーションが可能となるのである。

また、粘弾性体からの応力とフィラー間の相互作用力によって、フィラーの重心の位置が変化するため、フィラー同士の相対位置が変化した場合や、流れるような大きな変形を受けた場合についてシミュレーションできるようになる。以下、図２のフローチャートを参照して、本発明のシミュレーション方法の各ステップを説明する。ここでも、粘弾性体がゴムの場合を例に説明するが、これに限定されない。

まず、ステップＳ１において、ゴム、フィラー、およびゴムとフィラーとを接続する拘束層からなるゴム・フィラー複合体のモデルに基づいてシミュレーションを行うための、ゴム、フィラー、および拘束層の物性値、並びにゴム中におけるフィラーの初期配置を含む、動的挙動をシミュレーションするための計算条件を入力する。ここで、ゴム、フィラーおよび拘束層の物性値とは、具体的には、ゴムの粘度、弾性率、その周波数特性、拘束層の厚さ、粘度、弾性率、その周波数特性、およびフィラーのサイズ、個数、粘度、弾性率、その周波数特性等である。

また、ゴム・フィラー複合体の動的挙動をシミュレーションするための計算条件とは、具体的には、フィラーの配置の更新回数や時間および空間のメッシュ、システムに与える周期せん断ひずみの大きさおよび周波数等である。フィラーの配置は、ランダムに重心の座標を指定することにより配置する。

次いで、ステップＳ２において、ステップＳ１にて入力された計算条件の下で、溶媒であるゴムおよび拘束層が従う流体の運動方程式（６）を解く。具体的には、まず、モデルポテンシャルＵ（ｒ）を用いて、各フィラーに働く力Ｆ_ｉを計算する。このモデルポテンシャルＵ（ｒ）としては、例えば、下記の式（１７）に示すレナード・ジョーンズ（ＬＪ）型のポテンシャルを用いることができる。

ここで、εおよびαはフィラー間の相互作用の大きさを決定するパラメータである。

こうしたモデルポテンシャルＵ（ｒ）を用いることにより、各粒子に働く力Ｆ_ｉは、Ｕ（ｒ）の勾配を計算することにより（Ｕ（ｒ）をフィラーの位置ベクトルｒ_ｉで偏微分することにより）求めることができる。

ここで、Ｎはフィラーの数であり、「力Ｆ_ｉ」は、ベクトル量である。こうして求められた各フィラーに働く力Ｆ_ｉに、i番目のフィラーのφ_i（ｒ）を掛け合わせて全てのフィラーについて足し合わせることにより、流体の運動方程式（６）における力場Ｆ_pが得られる。

ここで、Ｖ_ｉは、要素ｉの体積である。この力場Ｆ（ｒ）を式（６）に代入して式（６）を解くことにより、各フィラーの速度ｖ_ｉが求められる。また、「速度ｖ_ｉ」は、ベクトル量である。

続いて、ステップＳ３において、上記流体の運動方程式の解に基づいてフィラーの移動量を算出する。具体的には、ステップＳ２において運動方程式（６）を解いて算出された各フィラーの速度ｖ_ｉを下記式（２０）に代入することにより、Δｔ後（例えば、０．０１秒後）の各フィラーの移動量Δｒ_ｉを算出する。ここで、「移動量Δｒ_ｉ」はベクトル量である。

その後、ステップＳ４において、得られたフィラーの移動量を用いてフィラーの配置を更新する。すなわち、ステップＳ３において得られた、Δｔ後の各フィラーの移動量Δｒ_ｉを、各フィラーの座標に足し合わせることにより、各フィラーの配置を更新する。上記ステップＳ２〜Ｓ４を繰り返し行うことにより、ゴム・フィラー複合体の動的挙動をシミュレーションできる。

粘弾性を評価する際には、システムに指定の周期せん断の境界条件を設定して上記連立方程式（２）、（４）〜（６）を解く。境界条件としては、実際の速度ｖからせん断速度を引いた量に、例えば、Ｌｅｅｓ−Ｅｄｗａｒｄｓ法による周期境界条件を課す。フィラーが境界を跨ぐ際にも、その位置、速度をせん断速度に合わせて周期境界条件を課して計算を行う。

続いて、ステップＳ５において、ステップＳ２〜Ｓ４の処理を所定の繰り返し回数だけ行ったかを判定し（例えば、Δｔ＝０．０１秒で１００秒のシミュレーションを行う場合には１００００回）、行っている場合には、ステップＳ６において、得られたゴム・フィラー複合体の動的挙動に関するデータを出力して終了し、行っていない場合にはステップＳ２に戻る。ここで、ゴム・フィラー複合体の動的挙動に関するデータは、各繰り返し回におけるフィラーの配置や速度等のデータを含んでいる。
こうして、ステップＳ２〜Ｓ４の処理を繰り返し行うことにより、ゴム・フィラー複合体の動的挙動をシミュレーションすることができる。

ここで、式（６）におけるゴムの粘度ηおよび緩和時間τとフィラーの流体効果との関係について調べる。図３は、ゴム・フィラー複合体の粘度および弾性率とフィラーの流体効果との関係を示している。図３（ａ）および（ｂ）ともに、拘束層のない場合について、Ｎ_ｆ＝Ｎ_ｂ＝１とした時のフィラーの流体効果を示しており、図３（ａ）は、τ_ｆ／τ_ｂ＝１の条件下で、Ｒ_１＝（η_ｆ／τ_ｆ）／（η_ｂ／τ_ｂ）とフィラーの流体効果との関係、図３（ｂ）は、η_ｆ／η_b＝１の条件下でＲ_２=τ_ｆ／τ_ｂとフィラーの流体効果との関係をそれぞれ示している。これらの図に示された結果は、ゴム・フィラー複合体に５％の歪み、周波数０．０１の振動が与えられた条件下で得られたものである。また、シミュレーションの数値は、密度ρ=１、粘度η_ｂ＝０．５、フィラー半径Ｒ＝３．８となる条件で無次元化して計算することにより得られたものである。

これらの図から、以下の式（２１）のいずれかを満たすことによって、フィラーによる剛体としての流体効果が適切に取り込まれていることが分かる。

ここで、「フィラーの流体効果」とは、貯蔵弾性率および損失弾性率の増分に対する、貯蔵弾性率および損失弾性率の実際の増分（シミュレーションにより得られた増分）を意味している（理論値＝１＋２．５ｇ、ｇは充填率）。従って、上記式（２１）のいずれかを満たす条件において、ゴム・フィラー複合体の粘弾性を適切に計算できていることが分かる。このように、粘度η_ｊ ^（ｂ）、η_ｊ ^（ｆ）、緩和時間τ_ｊ ^（ｂ）、τ_ｊ ^（ｆ）は、フィラー部の剛体球としての流体効果に影響を与える。一方、η_ｊ ^（I）τ_ｊ ^（I）の値はフィラー部の剛体球としての流体効果に影響を与えず、特に限定されない。

また、図３（ｃ）は、式（２１）を満たす条件下で、拘束層が適切に導入できているか検討した結果を示す図である。拘束層が適切に導入できているか否かは、以下のようにして評価した。すなわち、バルク部（粘弾性体）が弾性体の場合には、拘束層のような中間物性を有するフィラー分散系の弾性率を、例えば、G.Huber, and T. Vilgis, Macromolecules, vol. 35, p. 9204 (2002)に記載された理論式によって予測することができる。そこで、理論式から得られた弾性率と、本発明のシミュレーション方法から得られた弾性率とを比較した。図３（ｃ）から明らかなように、シミュレーション結果は理論値と良く一致していることが分かる。

同様に、バルク部（粘弾性体）が粘性体の場合にも同じ理論式で評価した結果、シミュレーション結果は理論値と良く一致しており、本発明のシミュレーション方法におけるモデル化によって、拘束層が適切に導入できていることが分かる。なお、これらの計算は、６４×６４の２次元系で次の条件の下で行った。すなわち、粘性１および粘性２は、密度ρ＝１、粘度η_ｂ＝０．５、η_ｆ＝５０、η_Ｉ＝５．０、緩和時間τ_ｂ＝τ_ｆ＝τ_Ｉ＝０．０の条件でフィラーを１つ置き、周波数０．００６で歪み１０％を与えて行った。また、粘性１は、ａ＝１．６、ｂ＝１．３３３とし、Ｒ＝４．０、４．６６６、５．３３３、６．０とした場合の結果である。また、粘性２は、ａ＝１．６、（Ｒ，ｂ）＝（４．０，１．３３３）、（４．６６６，２．６６６）、（５．３３３，４．０）、（６．０，５．３３３）とした場合の結果である。

一方、弾性１および弾性２は、密度ρ＝１、粘度η_ｂ／τ_ｂ＝０．５、η_ｆ／τ_ｆ＝５０、η_Ｉ／τ_Ｉ＝５．０、１／τ_ｂ＝１／τ_ｆ＝１／τ_Ｉ＝０．０となる条件でフィラーを１つ置き、周波数０．００６で歪み１０％を与えて行った。また、パラメータＲ，ａおよびｂは、弾性１が粘性１と同じにし、弾性２が粘性２と同じにした。なお、｜Ｇ^＊｜／|Ｇ^＊ _０|は、求められた各貯蔵弾性率をフィラーがない場合の貯蔵弾性率の理論値で割った値である。図３（ｄ）に、例として、（Ｒ，ａ，ｂ）＝（４．０，１．６，１．３３３）の場合について、φ_ｆ（ｘ）および（η_Ｉ／η_ｆ）φ_Ｉ（ｘ）を示す。ここで、横軸はフィラーの中心からの距離ｘを表す。

以上、図３から明らかなように、式（２１）のいずれかの式を満足させることにより、フィラーを剛体球として模擬し、本発明のシミュレーション方法におけるモデル化によって、拘束層が適切に導入できることが分かる。

以上の本発明の粘弾性体・フィラー複合体の動的挙動のシミュレーション方法において、データ出力ステップにて、求められたデータに基づいて粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに行うことにより、粘弾性体・フィラー複合体の粘度や弾性率等の物性値を求めることができる。粘弾性体としてゴムを例に説明すると、具体的には、ゴム・フィラー複合体の物性値は、動的粘弾性の貯蔵弾性率Ｇ’および損失弾性率Ｇ’’として求める。このＧ’およびＧ’’は以下のように求める。すなわち、大きさ（振幅）γ_０、周波数ωである周期せん断ひずみγ_０ｓｉｎ（ωt）を、システムのある方向（ｘ方向）にせん断流れを、その方向に直交する方向（ｙ方向）に勾配を付けて作用させて、その時の平均の応力σを次式（２２）から計算することにより求めることができる。

この応力が、

となった時の、Ａ／γ_０、Ｂ／γ_０がそれぞれ、Ｇ’、Ｇ’’である。これらＧ’およびＧ’’は、完全弾性体ではＧ’≠０、Ｇ’’＝０であり、完全粘性体ではＧ’＝０、Ｇ’’≠０である。また、ｔａｎδ＝Ｇ’’／Ｇ’が小さいほど、弾性体に近い粘弾性材料を意味する。
こうして、従来のＦＰＤ法では実質的に不可能だった、粘弾性体・フィラー複合体のマクロな物性を求めて評価することができ、ゴム等の粘弾性材料の設計に有効利用することができる。

（シミュレーション装置）
次に、本発明の粘弾性体・フィラー複合体の動的挙動のシミュレーション装置について説明する。図４は、本発明の粘弾性体・フィラー複合体の動的挙動のシミュレーション装置の一実施形態を示す図である。この装置１は、入力部１１と、解析部１２と、表示部１３とを備える。また、解析部１２は、演算器２１と、移動量算出器２２と、配置更新器２３と、データ出力器２４とを有する。以下、粘弾性体がゴムの場合を例に説明するが、これに限定されない。

入力部１１は、ゴム、フィラー、およびゴムとフィラーとを接続する拘束層からなるゴム・フィラー複合体のモデルに基づいてシミュレーションを行うための、ゴム、フィラーおよび拘束層の物性値、並びにゴム中におけるフィラーの初期配置を含む、動的挙動をシミュレーションするための計算条件を入力する。この入力部１１としては、例えばキーボードやマウス等を用いることができる。

解析部１２は、入力部１１により入力されたゴム、フィラーおよび拘束層の物性値、並びに動的挙動を求めるための計算条件下でゴムが従う流体の運動方程式である式（６）を解いて、フィラーの移動量を算出することにより、ゴム・フィラー複合体の動的挙動をシミュレーションする。以下、解析部１２の各構成について説明する。

演算器２１は、入力部１１により入力されたフィラーの配置の下で、ゴムおよび拘束層が従う流体の運動方程式を解く。ここで、ゴムおよび拘束層の弾性的性質を加味された流体の運動方程式である式（６）は、記憶装置（図示せず）に予め格納しておき、流体の運動方程式を解く際に記憶装置から読み出されるように構成されている。

移動量算出器２２は、演算器２１により得られた流体の運動方程式の解に基づいて、フィラーの移動量を算出する。

配置更新器２３は、移動量算出器２２により算出されたフィラーの移動量を用いて前記フィラーの配置を更新する。

データ出力器２４は、演算器２１、移動量算出器２２および配置更新器２３により求められたゴム・フィラー複合体の動的挙動に関するデータを出力する。

上記演算器２１、移動量算出器２２、配置更新器２３およびデータ出力器２４からなる解析部１２は、パーソナルコンピュータ中のＣＰＵやメモリ、記憶装置等により構成することができる。

表示部１３は、データ出力器２４により出力されたゴム・フィラー複合体に関するデータを表示する。この表示部１３としては、例えばディスプレイを用いることができる。

また、ゴムが従う流体の運動方程式（６）において、図３に示した、ゴムの粘度および緩和時間とフィラー効果との関係から、式（Ｂ）のいずれかの式を満足させることにより、フィラー部が剛体球としての適切な流体効果を与えることができるのは上述の通りである。
こうして、本発明のシミュレーション装置により、粘弾性体・フィラー複合体の動的挙動を簡便にシミュレーションすることができる。

以上の本発明のシミュレーション装置において、解析部１２は、上記演算器２１、移動量算出器２２、配置更新器２３により求められた粘弾性体・フィラー複合体の動的挙動に関するデータに基づいて、粘弾性体・フィラー複合体の粘弾性特性を解析する粘弾性解析器をさらに有することができる。こうして、従来のシミュレーション装置では実質的に不可能だった、粘弾性体・フィラー複合体のマクロな物性を求めて評価することができ、ゴム等の粘弾性材料の設計に有効利用することができる。

（シミュレーションプログラム）
上述した本発明のシミュレーション装置１として、コンピュータを好適に用いることができ、そのようなコンピュータは、図４に示した手段の各機能を実現する処理内容を記述したプログラムを、当該コンピュータの記憶部に格納しておき、当該コンピュータの中央演算処理装置（ＣＰＵ）によってこのプログラムを読み出して実行させることで実現することができる。具体的には、本発明に係るシミュレーションプログラムは、粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動をシミュレーションする装置として構成するコンピュータに、ａ）粘弾性体、フィラー、および粘弾性体とフィラーとを接続する拘束層からなる粘弾性体・フィラー複合体のモデルに基づいてシミュレーションを行うための、粘弾性体、フィラー、および拘束層の物性値、並びに粘弾性体中におけるフィラーの初期配置を含む、動的挙動をシミュレーションするための計算条件を入力する入力ステップと、ｂ）入力された計算条件の下で粘弾性体および拘束層が従う流体の運動方程式を解く演算ステップと、ｃ）流体の運動方程式の解に基づいてフィラーの移動量を算出する移動量算出ステップと、ｄ）算出されたフィラーの移動量を用いてフィラーの配置を更新する配置更新ステップと、ｅ）ステップｂ）〜ｄ）の繰り返し処理により求めた粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップとを実行させる。ここで、流体の運動方程式は、粘弾性体および拘束層の弾性的性質を加味したものである。こうして、本発明のシミュレーションプログラムをコンピュータに実行させることにより、粘弾性体・フィラー複合体の動的挙動を簡便にシミュレーションすることができる。

＜シミュレーション方法の妥当性の検証＞
まず、本発明のシミュレーション方法の妥当性を検証した。ここで、粘弾性体としてゴムを想定した。液体中や弾性体中に剛体フィラーを希薄に分散させた場合、その流体効果として、全体の粘度や弾性率がフィラーを入れた数に応じて変化することがＥｉｎｓｔｅｉｎ則やＧｕｔｈ−Ｇｏｌｄ則として知られている。理論的には小さいフィラーを分散させた場合には、通常レイノルズ数が十分小さく、慣性の影響を無視したＳｔｏｋｅｓ近似で挙動を解析することが可能であり、希薄な極限として、フィラーを１つ考慮した際の結果と一致する。
粘弾性媒質中のフィラーの流体効果についても同様のことが成立すると考えられる。そこで、ゴムにフィラーを配合した際に、弾性率が上昇する効果の１つであるＧｕｔｈ−Ｇｏｌｄ則が成立していることを検証した。すなわち、上記法則が成立している場合には、フィラーの面積分率ｇに対して、

となることが期待される。ここで、ωは周波数、Ｇ’およびＧ’’はゴム・フィラー複合体の粘弾性体の粘弾性係数であり、Ｇ’は貯蔵弾性率、Ｇ’’は損失弾性率である。また、Ｇ’_ｂｕｌｋおよびＧ’’_ｂｕｌｋはフィラーがない場合のゴムの粘弾性係数である。

図５は、本発明のシミュレーション方法による、フィラーの充填率とＧｕｔｈ−Ｇｏｌｄ因子Ａとの関係を示す図である。ここで、この関係は、ゴム・フィラー複合体に対して１０％のひずみ、周波数０．０１の振動を与えることにより得られたものである。数値は、密度ρ＝１、粘度η＝０．５、フィラー半径Ｒ＝３．８となる条件で無次元化した。この図から明らかなように、フィラーの充填率の増大とともに、Ｇｕｔｈ−Ｇｏｌｄ因子が増大しており、ゴム・フィラー複合体の応力が増大していることが分かる。

図６（ａ）は、拘束層を含むフィラーの表面形状を示す図である。この図は、式（１４）を用いて、Ｒ＝３、ａ＝１．０、ｂ＝３．０とした場合のφ_Iおよびφ_ｆを示しており、横軸はフィラー中心からの距離、縦軸はφ_＋を表している。この図において、ｘ＝０〜３までがフィラーの領域であり、ｘ＝３〜６まで拘束層の領域あることを示している。

また、図６（ｂ）は、拘束層の厚さと本発明のシミュレーション方法により得られた複素粘弾性係数Ｇ^＊＝Ｇ’＋ｉＧ’’（ｉは虚数）との関係を示す図である。この関係は、ゴム・フィラー複合体に対して１０％の歪み、周波数０．００６の振動を与えることにより得られたものである。また、フィラー部および拘束層は、式（１４）を用いて規定した。図６（ｂ）における数値は、密度ρ＝１．０、粘度η_ｂ＝０．５、η_ｆ＝５０、フィラー半径Ｒ＝３となる条件で無次元化した。緩和時間としてはτ_ｂ＝τ_ｆ＝０，１０および３００の３つの場合について計算した。この図から明らかなように、フィラーがない時の複素粘弾性係数Ｇ^＊ _０に対して、拘束層の厚さの増加とともに｜Ｇ^＊｜／|Ｇ^＊ _０|が増加し、ゴム・フィラー複合体の応力が増加していることが分かる。これらの結果から、τ≠０の粘弾性条件においても、拘束層によって発生する応力が増加し、拘束層を導入した効果が得られたことが分かる。

１ シミュレーション装置
１１ 入力部
１２ 解析部
１３ 表示部
２１ 演算器
２２ 移動量算出器
２３ 配置更新器
２４ データ出力器
Ｓ１ 入力ステップ
Ｓ２ 演算ステップ
Ｓ３ 移動量算出ステップ
Ｓ４ 配置更新ステップ
Ｓ５ 繰り返し回数判定ステップ
Ｓ６ データ出力ステップ

Claims (12)

1. 粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動のシミュレーション方法であって、
   ａ）前記粘弾性体、前記フィラー、および前記粘弾性体と前記フィラーとを接続する拘束層からなる前記粘弾性体・フィラー複合体のモデルに基づいて前記シミュレーションを行うための、前記粘弾性体、前記フィラー、および前記拘束層の物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力ステップと、
   ｂ）入力された前記計算条件の下で前記粘弾性体および前記拘束層が従う流体の運動方程式を解く演算ステップと、
   ｃ）前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出ステップと、
   ｄ）算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新ステップと、
   ｅ）前記ステップｂ）〜ｄ）の繰り返し処理により求めた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップと、
   を含み、前記粘弾性体および前記拘束層が従う流体の運動方程式は、前記粘弾性体および拘束層の弾性的性質を加味したものであることを特徴とするシミュレーション方法。
2. 前記粘弾性体および前記拘束層が従う流体の運動方程式は以下の式（Ａ）で表される、請求項１に記載のシミュレーション方法。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆ_ｐは外場、Πは連続体の応力テンソル、Π_ｊ ^（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Π_ｊ ^（ｆ）は連続体の応力テンソルのフィラー部分の成分、Π_ｊ ^（Ｉ）は連続体の応力テンソルの拘束層部分の成分、Ｎ_ｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_ｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_Ｉは前記拘束層における粘弾性Ｍａｘｗｅｌｌ要素の数、τ_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素の緩和時間、η_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φ_ｂは前記粘弾性体の位置にて１を返し前記粘弾性体以外の位置で０を返す関数、φ_ｆは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、φ_Ｉは前記拘束層の位置にて１を返し前記拘束層以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。
   

   
3. 前記式（Ａ）における粘度η_ｊ ^（ｂ）、η_ｊ ^（ｆ）、およびη_ｊ ^（Ｉ）、緩和時間τ_ｊ ^（ｂ）、τ_ｊ ^（ｆ）、およびτ_ｊ ^（Ｉ）は下記の式（Ｂ）のいずれかを満足する、請求項２に記載のシミュレーション方法。
   

   
4. 前記データ出力ステップは、前記求められたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに含む、請求項１〜３のいずれか一項に記載のシミュレーション方法。
5. 粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動をシミュレーションする装置であって、
   前記粘弾性体、前記フィラー、および前記粘弾性体と前記フィラーとを接続する拘束層からなる前記粘弾性体・フィラー複合体のモデルに基づいて前記シミュレーションを行うための、前記粘弾性体、前記フィラーおよび前記拘束層の物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力部と、
   入力された前記計算条件の下で前記粘弾性体および前記拘束層が従う流体の運動方程式を解く演算器と、前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出器と、算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新器と、前記演算器、前記移動量算出器および前記配置更新器により求められた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力器とを有する解析部と、
   を備え、前記粘弾性体および前記拘束層が従う流体の運動方程式は、前記粘弾性体および拘束層の弾性的性質を加味したものであることを特徴とするシミュレーション装置。
6. 前記粘弾性体および前記拘束層が従う流体の運動方程式は以下の式（Ａ）で表される、請求項５に記載のシミュレーション装置。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆ_ｐは外場、Πは連続体の応力テンソル、Π_ｊ ^（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Π_ｊ ^（ｆ）は連続体の応力テンソルのフィラー部分の成分、Π_ｊ ^（Ｉ）は連続体の応力テンソルの拘束層部分の成分、Ｎ_ｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_ｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_Ｉは前記拘束層における粘弾性Ｍａｘｗｅｌｌ要素の数、τ_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素の緩和時間、η_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φ_ｂは前記粘弾性体の位置にて１を返し前記粘弾性体以外の位置で０を返す関数、φ_ｆは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、φ_Ｉは前記拘束層の位置にて１を返し前記拘束層以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。
   

   
7. 前記式（Ａ）における粘度η_ｊ ^（ｂ）、η_ｊ ^（ｆ）およびη_ｊ ^（Ｉ）、緩和時間τ_ｊ ^（ｂ）、τ_ｊ ^（ｆ）およびτ_ｊ ^（Ｉ）は下記の式（Ｂ）のいずれかを満足する、請求項６に記載のシミュレーション装置。
   

   
8. 前記解析部は、前記得られたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析器をさらに有する、請求項５〜７のいずれか一項に記載のシミュレーション装置。
9. 粘弾性体中にフィラーが分散した粘弾性体・フィラー複合体の動的挙動をシミュレーションする装置として構成するコンピュータに、
   ａ）前記粘弾性体、前記フィラー、および前記粘弾性体と前記フィラーとを接続する拘束層からなる前記粘弾性体・フィラー複合体のモデルに基づいて前記シミュレーションを行うための、前記粘弾性体、前記フィラー、および前記拘束層の物性値、並びに前記粘弾性体中における前記フィラーの初期配置を含む、前記動的挙動をシミュレーションするための計算条件を入力する入力ステップと、
   ｂ）入力された前記計算条件の下で前記粘弾性体および前記拘束層が従う流体の運動方程式を解く演算ステップと、
   ｃ）前記流体の運動方程式の解に基づいて前記フィラーの移動量を算出する移動量算出ステップと、
   ｄ）算出された前記フィラーの移動量を用いて前記フィラーの配置を更新する配置更新ステップと、
   ｅ）前記ステップｂ）〜ｄ）の繰り返し処理により求めた前記粘弾性体・フィラー複合体の動的挙動に関するデータを出力するデータ出力ステップと、
   を実行させ、前記流体の運動方程式は、前記粘弾性体および前記拘束層の弾性的性質を加味したものであることを特徴とするシミュレーションプログラム。
10. 前記粘弾性体および前記拘束層が従う流体の運動方程式は以下の式（Ａ）で表される、請求項９に記載のシミュレーションプログラム。ただし、ρは質量密度、Ｄは連続体の歪み速度テンソル、ｖは速度、ｐは圧力、Ｆ_ｐは外場、Πは連続体の応力テンソル、Π_ｊ ^（ｂ）は連続体の応力テンソルの粘弾性体部分の成分、Π_ｊ ^（ｆ）は連続体の応力テンソルのフィラー部分の成分、Π_ｊ ^（Ｉ）は連続体の応力テンソルの拘束層部分の成分、Ｎ_ｂは前記粘弾性体における粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_ｆは前記フィラーにおける粘弾性Ｍａｘｗｅｌｌ要素の数、Ｎ_Ｉは前記拘束層における粘弾性Ｍａｘｗｅｌｌ要素の数、τ_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素の緩和時間、τ_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素の緩和時間、η_ｊ ^（ｂ）は前記粘弾性体におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（ｆ）は前記フィラーにおけるＭａｘｗｅｌｌ要素における粘性要素の粘度、η_ｊ ^（Ｉ）は前記拘束層におけるＭａｘｗｅｌｌ要素における粘性要素の粘度、φ_ｂは前記粘弾性体の位置にて１を返し前記粘弾性体以外の位置で０を返す関数、φ_ｆは前記フィラーの位置にて１を返し前記フィラー以外の位置にて０を返す関数、φ_Ｉは前記拘束層の位置にて１を返し前記拘束層以外の位置にて０を返す関数、∇は上対流微分（Ｕｐｐｅｒ Ｃｏｎｖｅｃｔｉｖｅ ｔｉｍｅ ｄｅｒｉｖａｔｉｖｅ）である。
    

    
11. 前記式（Ａ）における粘度η_ｊ ^（ｂ）、η_ｊ ^（ｆ）およびη_ｊ ^（Ｉ）、緩和時間τ_ｊ ^（ｂ）、τ_ｊ ^（ｆ）およびτ_ｊ ^（Ｉ）は下記の式（Ｂ）のいずれかを満足する、請求項１０に記載のシミュレーションプログラム。
    

    
12. 前記データ出力ステップは、前記求められたデータに基づいて前記粘弾性体・フィラー複合体の粘弾性を解析する粘弾性解析ステップをさらに含む、請求項９〜１１のいずれか一項に記載のシミュレーションプログラム。

Images