JP6030986B2 - ゴム材料の接触シミュレーション方法 - Google Patents

Description

本発明は、ゴム材料と、任意の材料の接触面との接触状態を解析しうるゴム材料の接触シミュレーション方法に関する。

従来、ゴム材料と、任意の材料の接触面との接触状態を解析するには、例えば、有限要素法に基づいて、ゴム材料及び接触面を有する接触部それぞれを有限個の要素でモデル化するシミュレーション方法が提案されている。このようなシミュレーション方法では、各要素に定義される材料特性等に基づいて、ゴム材料をモデル化したゴムモデルと、接触部をモデル化した接触部モデルとの接触状態を解析することができる。

しかしながら、各要素は、隣り合う要素で拘束されているため、ゴムモデルの引張変形が１００％以上となるような大きな変形計算が行われると、要素潰れ等による計算落ちが生じやすい。従って、有限要素法に基づくシミュレーション方法では、大変形を伴うゴム材料と接触面との接触状態を十分に解析することができないという問題があった。

また、下記非特許文献１では、分子モデル間に相互作用を定義して、分子間の接触状態を解析することができる分子動力学に基づくシミュレーション方法が提案されている。このようなシミュレーション方法では、上記のような計算落ちが発生することはない。

しかしながら、分子動力学に基づくシミュレーション方法では、分子間のミクロ・ナノレベルの接触状態のみを解析するものであるため、ゴム材料と接触面とのマクロレベルの接触状態を解析するには現実的に不可能であった。

さらに、下記非特許文献２及び非特許文献３では、解析対象物を、有限個の粒子でモデル化する粒子法に基づくシミュレーション方法が提案されている。各粒子には、上記のような分子間の相互作用ではなく、ニュートンの運動方程式がそれぞれ定義されている。

H. Morita, T. Ikehara, T. Nishi, M. Doi, Polymer J. 36, 265 (2004) Yoichi KAWASHIMA and Yuzuru SAKAI 著「Large Deformation Analysis of Hyperelastic Materials Using SPH Method」 e-Journal of Soft Materials Vol.3, pp.22-28 (2007) 酒井 譲、"ＳＰＨ粒子法の基礎と応用"、［online］、ＣＡＥ懇話会、［平成２４年１０月１日検索］、インターネット<URL:http://www.cae21.org/kaisekijuku2006/SPH.pdf>

しかしながら、各粒子の運動方程式には、ゴム材料に有する復元力が定義されていない。このため、ゴムモデルは、大きな変形が計算されると、変形したままとなり、元の形状に復元することができない。従って、このような粒子法に基づくシミュレーション方法では、大変形を伴うゴム材料と接触面との接触状態を解析できないという問題があった。

本発明は、以上のような実状に鑑み案出されたもので、ゴムモデルの運動方程式に、引張変形に対するゴム材料の復元力を定義することを基本として、大きな変形を伴うゴム材料と接触面との接触状態を精度良く解析しうるゴム材料の接触シミュレーション方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、コンピュータを用いて、ゴム材料と、該ゴム材料が接触する接触面との接触状態を、粒子法に基づいて解析するゴム材料の接触シミュレーション方法であって、前記コンピュータに、前記ゴム材料を、有限個の粒子でモデル化し、かつ、該粒子に運動方程式を定義したゴムモデルを入力するステップと、前記コンピュータに、前記接触面を含む接触部を、有限個の粒子でモデル化した接触部モデルを入力する接触部モデル入力ステップと、前記コンピュータが、前記ゴムモデルを前記接触部モデルに接触させて、該ゴムモデルの変形計算を行う変形計算ステップとを含み、前記ゴムモデルの前記運動方程式には、引張変形に対する前記ゴム材料の復元力が定義されることを特徴とする。

また、請求項２記載の発明は、前記変形計算ステップでは、前記ゴムモデルの前記粒子と前記接触部モデルの前記粒子との間に、前記ゴムモデルと前記接触部モデルとの接触面内において、前記ゴムモデルと前記接触部モデルとの接触を表現する接触バネを設定するステップを含む請求項１に記載のゴム材料の接触シミュレーション方法である。

また、請求項３記載の発明は、前記接触部モデルは、前記接触面に、凹凸を有する請求項１又は２に記載のゴム材料の接触シミュレーション方法である。

また、請求項４記載の発明は、前記変形計算ステップは、前記接触部モデル側への荷重を前記ゴムモデルに定義し、かつ該ゴムモデルを、該接触部モデルの前記接触面に沿ってスライドさせた状態を計算するステップを含む請求項１乃至３のいずれかに記載のゴム材料の接触シミュレーション方法である。

また、請求項５記載の発明は、前記ゴムモデルは、前記ゴム材料の二次元の断面形状に形成され、前記ゴムモデルの前記接触面と反対側の面には、変形不能に定義された剛部材モデルが定義され、 前記変形計算ステップでは、前記剛部材モデルの上下方向の座標を一定に保持する請求項４に記載のゴム材料の接触シミュレーション方法である。

また、請求項６記載の発明は、前記ゴムモデルには、Ａｒｒｕｄａ−Ｂｏｙｃｅ型構成方程式が定義され、前記変形計算ステップは、前記運動方程式に基づいて、前記ゴムモデルの前記各粒子の加速度を計算する第１加速度計算ステップ、前記加速度で移動した前記ゴムモデルの歪を計算する歪計算ステップ、前記ゴムモデルの歪を、前記Ａｒｒｕｄａ−Ｂｏｙｃｅ型構成方程式に代入して、前記ゴムモデルが前記接触部モデルから受ける応力を計算する応力計算ステップ、前記応力を前記運動方程式に代入して、前記ゴムモデルの前記各粒子の加速度を計算する第２加速度計算ステップを含み、前記変形計算ステップが終了するまで、前記歪計算ステップ、前記応力計算ステップ、及び、前記第２加速度計算ステップが繰り返し実行される請求項１乃至５のいずれかに記載のゴム材料の接触シミュレーション方法である。

また、請求項７記載の発明は、前記応力計算ステップでは、前記ゴムモデルの前記応力から、前記ゴムモデルと前記接触部モデルとの摩擦力を求める請求項６に記載のゴム材料の接触シミュレーション方法である。

本発明のゴム材料の接触シミュレーション方法は、コンピュータに、ゴム材料を、有限個の粒子でモデル化し、かつ、該粒子に運動方程式を定義したゴムモデルを入力するステップと、コンピュータに、接触面を含む接触部を、有限個の粒子でモデル化した接触部モデルを入力する接触部モデル入力ステップと、コンピュータが、ゴムモデルを接触部モデルに接触させて、該ゴムモデルの変形計算を行う変形計算ステップとを含んでいる。

このようなシミュレーション方法では、ゴム材料が、有限個の粒子によってモデル化されるため、例えば、引張変形が１００％以上となるような大きな変形計算を柔軟に行うことができる。従って、本発明のゴム材料の接触シミュレーション方法では、従来のような、要素潰れによる計算落ちを防ぐことができる。

また、ゴムモデルの運動方程式には、引張変形に対するゴム材料の復元力が定義される。従って、ゴムモデルに大きな変形が与えられるシミュレーションであっても、ゴムモデルは、元の形状に復元することができる。従って、本発明のシミュレーション方法では、ゴム材料と接触面との接触状態を精度良く解析することができる。

本実施形態のゴム材料及び接触面の一例を示す断面図である。 コンピュータ装置の斜視図である。 本実施形態のシミュレーション方法を示すフローチャートである。 ゴムモデル及び接触部モデルを視覚化した側面図である。 連続体を離散化した粒子を視覚化して示す平面図である。 変形計算ステップを示すフローチャートである。 接触ステップを示すフローチャートである。 （ａ）は、接触面に接触したゴムモデルを視覚化して示す側面図、（ｂ）は、荷重によって変形したゴムモデルを視覚化して示す側面図である。 ゴムモデルの粒子と接触部モデルの粒子との間に定義される接触バネを視覚化して示す側面図である。 変形ステップでのゴムモデル及び接触部モデルを視覚化した側面図である。 スライドステップを示すフローチャートである。 応力と時間との関係を示すグラフである。 （ａ）は、実施例１のゴムモデルを視覚化して示す側面図、（ｂ）は、実施例３のゴムモデルを視覚化して示す側面図である。

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態のゴム材料１の接触シミュレーション方法（以下、単に「接触シミュレーション方法」ということがある）は、図１に示されるように、コンピュータを用いて、ゴム材料１と、該ゴム材料１が接触する接触面２Ｓとの接触状態を、粒子法に基づいて解析するための方法である。

図２に示されるように、コンピュータ３は、本体３ａ、キーボード３ｂ、マウス３ｃ、及び、ディスプレイ装置３ｄを含む。この本体３ａには、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリー、磁気ディスクなどの記憶装置、及び、ディスクドライブ装置３ａ１、３ａ２などが設けられる。また、記憶装置には、本実施形態の接触シミュレーション方法を実行するためのソフトウェアや、処理手順が予め記憶されている。

図３には、本実施形態の接触シミュレーション方法の具体的な処理手順が示される。
本実施形態では、先ず、図４に示されるように、コンピュータ３に、ゴム材料１（図１に示す）をモデル化したゴムモデル５を入力する（ステップＳ１）。ゴムモデル５は、粒子法に基づいて、例えば、ゴム材料１の二次元の断面形状（例えば、断面略矩形状）を、有限個の粒子６でモデル化して設定される。また、各粒子６は、コンピュータ３で取り扱うための数値データである。これらの数値データは、コンピュータ３に記憶される。なお、数値データとしては、例えば、質量、体積、直径、座標及び材料特性等が含まれる。

本実施形態の粒子法としては、例えば、本実施形態では上記非特許文献２及び非特許文献３に示したＳＰＨ粒子法が採用される。ここで、ＳＰＨ粒子法は、図５に示されるように、任意の連続体１１を、複数の粒子１２に離散化し、これらの各粒子１２の集合を連続体１１として取り扱う手法である。各粒子１２の物理量ｆ（ｘ）は、非特許文献２及び非特許文献３の記載に従って、下記式（１）で定義される。

ここで、符号は次のとおりである。
Ｗ（ｘ−ｘ’，ｈ）：重み関数
ｈ：影響半径
ｘ：評価点の中心座標値
ｘ’：任意粒子の座標値

また、重み関数Ｗ（ｘ−ｘ’，ｈ）は、影響半径ｈの内部に分布する粒子１２の重み付けを行う関数である。このため、上記式（１）で示した物理量ｆ（ｘ）を求めるには、下記式（２）乃至下記式（４）で示される条件を満たす必要がある。

上記式（２）では、重み関数（ｘ−ｘ’，ｈ）の極限（ｈ→０）がデルタ関数であることを示している。また、上記式（３）では、重み関数Ｗ（ｘ−ｘ’，ｈ）が規格化されることを示している。さらに、上記式（４）では、影響半径ｈの範囲以外で、重み関数Ｗ（ｘ−ｘ’，ｈ）がゼロになるコンパクト化が可能であることを示している。なお、上記式（１）は、上記非特許文献２及び非特許文献３に基づいて、適宜変形して利用することができる。

また、ＳＰＨ粒子法では、各粒子１２に、下記式（５）に示される運動方程式が定義される。この下記式（５）では、各粒子１２の加速度ａを定義している。この加速度ａに基づいて、各粒子１２の運動が定義される。なお、下記式（５）は、上記非特許文献２及び非特許文献３に基づいて、適宜変形して利用することができる。

ここで、
ａ：加速度
ρ：粒子の密度
σ：応力テンソル

さらに、ＳＰＨ法では、下記式（６）に示される歪速度、及び、粒子速度の関係式を用いて、二次元弾性場、即ち、連続体１１の歪ε_ijを求めることができる。

ここで：
Ｖ_i：評価粒子の速度
Ｖ_j：影響半径内にある他の粒子の速度
ｘ_i：評価粒子の座標
ｘ_j：影響半径内にある他の粒子の座標

図４に示されるように、本実施形態のゴムモデル５には、上記式（１）に示した物理量ｆ（ｘ）が設定される。なお、本実施形態の重み付け関数Ｗ（ｘ−ｘ’，ｈ）は、３次のスプライン関数が使用されている。また、本実施形態の影響半径ｈには、変形計算が行わる前の初期状態において、隣接する粒子６、６間の距離Ｄ１（図９に示す）の２倍に設定されている。

また、本実施形態のゴムモデル５の粒子６には、上記式（５）の運動方程式に代えて、下記式（７）、（８）で示される運動方程式が設定される。

ここで、
ｆ：粒子間の粘弾性力
ｋ：バネ係数
η：せん断バネ係数
ΔＬ_i：バネの長さ
ΔC_i：伸縮量
Δｔ：時間増分
Ｗ_i：重み係数

上記式（７）は、上記式（５）に、粒子６、６間の粘弾性力ｆが加えられている。この粘弾性力ｆは、加速度ａによる粒子６、６間の広がり（例えば、引張変形）を抑える働きを示す関数である。これにより、ゴムモデル５の運動方程式には、引張変形に対するゴム材料１（図１に示す）の復元力が定義される。従って、ゴムモデル５は、大きな変形が与えられるようなシミュレーションであっても、元の形状に復元することができる。

また、上記式（８）のバネ係数ｋは、実際のゴム材料１（図１に示す）で測定されたヤング率が設定される。また、せん断バネ係数ηは、ゴム材料１で測定された粘度（ポアズ）が設定される。このように、各係数ｋ、ηには、ゴム材料１の実際の測定値が設定されるため、シミュレーション精度を向上しうる。なお、せん断バネ係数ηは、図１に示したゴム材料１と接触面２Ｓとの摩擦を示すパラメータである。このため、本実施形態のように、ゴム材料１と接触面２Ｓとの接触状態を解析する場合には、上記式（８）において、せん断バネ係数ηを省略することができる。

さらに、本実施形態のゴムモデル５には、分子鎖網目理論に基づくＡｒｒｕｄａ−Ｂｏｙｃｅ型構成方程式が定義される。なお、Ａｒｒｕｄａ−Ｂｏｙｃｅ型構成方程式は、「Ellen M. Arruda and Marry C. Boyce著「A THREE-DIMENSIONAL CONSTITUTIVE MODEL FOR THE LARGE STRECH BEHAVIOR OF RUBBER ELASTIC MATERIALSS」 Journal of the Mechanics and Physics of Solids Volume 41, Issue 2, Pages 389-412 (February 1993)」に基づいて、例えば、特開２０１１−２４２３３６号公報に記載されるように、下記式（９）に変形して使用することができる。

ここで、
σ：応力
Ｋ_B：Ｂｏｌｔｚｍａｎｎ定数
Ｔ：絶対温度
λ：ストレッチ
Ｎ：分子鎖一本あたりの平均セグメント数
Ｌ（ｘ）：下記式（１０）で定義されるＬａｎｇｅｖｉｎ関数である。

上記式（９）では、上記式（６）で求められるゴムモデル５の歪ε_ijが、ストレッチλに代入されると、ゴムモデル５が接触面９Ｓから受ける応力σを計算することができる。このような応力σは、ゴムモデル５と、接触部２（図１に示す）をモデル化した接触部モデル９との摩擦力とみなすことができる。

本実施形態のゴムモデル５には、接触部モデル９の接触面９Ｓと反対側の面、即ち、本実施形態では上面５ｕに、変形不能に定義された剛部材モデル７が定義される。この剛部材モデル７は、粒子法に基づいて、有限個の粒子８でモデル化される。これらの粒子８も、コンピュータ３で取り扱うための数値データであり、コンピュータ３にそれぞれ記憶される。

また、剛部材モデル７の粒子８には、上記式（１）で示した物理量ｆ（ｘ）や、上記式（５）で示される運動方程式が定義されている。さらに、剛部材モデル７は、例えば、下記式（１１）に示される金属の弾性構成方程式が定義される。

ここで、
{σ}：応力マトリクス
{ε}：ひずみマトリクス
[D]：応力・ひずみマトリクス

次に、コンピュータ３に、接触面２Ｓを含む接触部２（図１に示す）をモデル化した接触部モデル９が入力される（ステップＳ２）。この接触部モデル９も、粒子法に基づいて、例えば、接触部２を、有限個の粒子１０でモデル化される。これらの粒子１０も、コンピュータ３で取り扱うための数値データであり、コンピュータ３にそれぞれ記憶される。

また、本実施形態の接触部モデル９の接触面９Ｓには、凹凸が形成される。この凹凸は、上方に突出する凸部１５と、下方に凹む凹部１６とが、水平方向で交互に形成されている。また、凸部１５の上端１５ｕと、凹部１６の下端１６ｄとの間の凹凸高さＨ１は、任意に設定することができる。

なお、本実施形態の接触部モデル９の粒子１０には、上記式（１）で示した物理量ｆ（ｘ）や、上記式（５）で示される運動方程式が定義されている。さらに、本実施形態の接触部モデル９は、例えば、鋼鉄（ヤング率２０８ＧＰａ、ポアソン比０．３３））として定義される。このような接触部モデル９は、ゴムモデル５（例えば、ヤング率１０ＭＰａ程度）に比べると、完全剛体とみなすことができる。なお、接触部モデル９を鋼鉄（完全剛体）として定義するには、上記式（１１）を用いることができる。なお、接触部モデル９は、軟体として定義されてもよい。

次に、コンピュータ３が、ゴムモデル５を接触部モデル９に接触させて、該ゴムモデル５の変形計算を行う（変形計算ステップＳ３）。本実施形態では、図１に示されるように、ゴムモデル５を接触部モデル９の接触面９Ｓに接触させて、接触面９Ｓに沿ってスライドさせた状態を計算する。図６には、本実施形態の変形計算ステップＳ３の具体的な処理手順が示される。

変形計算ステップＳ３では、先ず、ゴムモデル５を、接触部モデル９に接触させる（接触ステップＳ４）。図７には、本実施形態の接触ステップＳ４の具体的な処理手順が示される。

図８（ａ）に示されるように、接触ステップＳ４では、ゴムモデル５を接触面９Ｓに接触させて、ゴムモデル５の剛部材モデル７に、接触部モデル９側への荷重Ｇを定義する（ステップＳ４１）。本実施形態のステップＳ４１では、有限要素法に基づいて、剛部材モデル７の粒子８の運動方程式に、荷重Ｇが定義される。これにより、剛部材モデル７の移動が計算される。また、剛部材モデル７は、実質的に変形不能に設定されるため、ゴムモデル５の上面５ｕの全体に、荷重Ｇが連続的に設定される。

次に、上記式（７）に示した運動方程式に基づいて、荷重Ｇによるゴムモデル５の各粒子６（図４に示す）の加速度ａを計算する（第１加速度計算ステップＳ４２）。本実施形態の第１加速度計算ステップＳ４２では、剛部材モデル７の移動に基づいて、ゴムモデル５の各粒子６の加速度ａが計算される。各粒子６は、加速度ａに基づいてそれぞれ移動し、図８（ｂ）に示されるように、荷重Ｇによるゴムモデル５の変形が計算される。

次に、図９に示されるように、ゴムモデル５の各粒子６と接触部モデル９の各粒子１０との間に、接触バネ２０を定義する（接触バネ設定ステップＳ４３）。接触バネ２０は、ゴムモデル５の各粒子６と、接触部モデル９の粒子１０との接触を妨げる方向に復元力を作用させるものである。本実施形態の接触バネ２０は、上記式（８）に示される粘弾性力ｆによって定義される。なお、接触バネ２０が定義された粒子６の運動は、粘弾性力ｆを考慮した上記式（７）の運動方程式によって定義することができる。

また、接触バネ２０は、接触部モデル９の各粒子１０との距離Ｌ３が、影響半径ｈ以下となるゴムモデル５の各粒子６のみに定義される。これにより、接触バネ２０は、ゴムモデル５と接触部モデル９との接触面９Ｓ内において、ゴムモデル５の粒子６が、接触部モデル９の粒子１０への浸透を抑制することができ、ゴムモデル５と接触部モデル９との接触状態を維持することができる。

次に、ゴムモデル５の歪ε_ijを計算する（歪計算ステップＳ４４）。この歪計算ステップＳ４４では、上記式（６）に基づいて、加速度ａで移動した後に生じるゴムモデル５の歪ε_ijが計算される。本実施形態では、接触面９Ｓ内において、接触バネ２０が定義されているため、ゴムモデル５の歪ε_ijを精度良く計算することができる。

次に、ゴムモデル５の歪ε_ijに基づいて、ゴムモデル５が接触部モデル９から受ける応力σを計算する（応力計算ステップＳ４５）。この応力計算ステップＳ４５では、上記式（９）に示したＡｒｒｕｄａ−Ｂｏｙｃｅ型構成方程式に、歪ε_ijを代入して、ゴムモデル５が接触部モデル９から受ける応力σを計算する。これにより、ゴムモデル５が接触部モデル９から受ける応力σ、即ち、ゴムモデル５と接触部モデル９との摩擦力を求めることができる。

次に、応力σに基づいて、ゴムモデル５の各粒子６（図４に示す）の新たな加速度ａを計算する（第２加速度計算ステップＳ４６）。第２加速度計算ステップＳ４６では、ゴムモデル５の各粒子６に設定されている上記式（７）の運動方程式に、応力σを代入して、各粒子６の加速度ａを求めている。これにより、各粒子６には、ゴムモデル５が接触部モデル９から受ける応力σに基づく新たな加速度ａが設定される。各粒子６は、この新たな加速度ａに基づいてそれぞれ移動し、ゴムモデル５の変形が計算される。

次に、コンピュータ３が、ゴムモデル５の変形計算が終了したか否かを判断する（終了判断ステップＳ４７）。この終了判断ステップＳ４７では、ゴムモデル５の各粒子６（図４に示す）の振動が収束した場合に、ゴムモデル５の変形計算が終了したと判断している。

本実施形態のステップＳ４７では、コンピュータ３がゴムモデル５の変形計算が終了したと判断した場合、接触ステップＳ４を終了し、次のスライドステップＳ５を実施する。一方、コンピュータ３がゴムモデル５の変形が終了していないと判断した場合には、単位時間Ｔｘを一つ進めて（ステップＳ４８）、接触バネ設定ステップＳ４３、歪計算ステップＳ４４、応力計算ステップＳ４５、第２加速度計算ステップＳ４６、及び、終了判断ステップＳ４７の各ステップが順次行われる。これにより、コンピュータ３は、荷重Ｇによるゴムモデル５の変形計算を確実に終了させることができ、次のスライドステップＳ５でのシミュレーション精度を向上させることができる。

次に、図１０に示されるように、荷重Ｇをゴムモデル５に定義したまま、ゴムモデル５を、接触部モデル９の接触面９Ｓに沿ってスライドさせた状態を計算する（スライドステップＳ５）。図１１には、本実施形態のスライドステップＳ５の具体的な処理手順が示される。

図１０に示されるように、本実施形態のスライドステップＳ５は、先ず、ゴムモデル５の各粒子６の加速度ａを計算する（第１加速度計算ステップＳ５１）。このステップＳ５１では、上記式（７）で示した運動方程式に基づいて、ゴムモデル５のスライド速度Ｖ１による各粒子６の加速度ａが計算される。これにより、各粒子６は、加速度ａに基づいてそれぞれ移動し、ゴムモデル５のスライドが計算される。

次に、接触部モデル９の各粒子１０の加速度ａを計算する（ステップＳ５２）。このステップＳ５２では、上記式（５）で示した運動方程式に基づいて、スライド速度Ｖ１に対応する並進速度Ｖ２による各粒子１０の加速度ａが計算される。これにより、各粒子１０は、加速度ａに基づいてそれぞれ移動し、接触部モデル９の移動が計算される。従って、ゴムモデル５が、接触部モデル９の接触面９Ｓに沿ってスライドした状態を計算することができる。

次に、図９に示されるように、ゴムモデル５各粒子６と、接触部モデル９の各粒子１０との間に、接触バネ２０を定義する（接触バネ設定ステップＳ５３）。この接触バネ設定ステップＳ５３では、接触ステップＳ４の接触バネ設定ステップＳ４３と同様に、接触部モデル９の各粒子１０との距離Ｌ３が影響半径ｈ以下となるゴムモデル５の各粒子６にのみ、接触バネ２０を設定している。これにより、ゴムモデル５は、接触部モデル９の接触面９Ｓに沿ってスライドした状態において、接触部モデル９との接触状態を維持することができる。

次に、図１０に示されるように、ゴムモデル５の歪ε_ijを計算する（歪計算ステップＳ５４）。この歪計算ステップＳ５４では、上記式（６）に基づいて、加速度ａで移動した後に生じるゴムモデル５の歪ε_ijが計算される。本実施形態のスライドステップＳ５では、接触ステップＳ４と同様に、歪計算ステップＳ５４に先立ち、接触面９Ｓ内において、接触バネ２０が定義されるため、ゴムモデル５の歪ε_ijを精度良く計算することができる。

次に、ゴムモデル５の歪ε_ijに基づいて、ゴムモデル５が接触部モデル９から受ける応力σを計算する（応力計算ステップＳ５５）。この応力計算ステップＳ４５では、上記式（９）に示すＡｒｒｕｄａ−Ｂｏｙｃｅ型構成方程式に、歪ε_ijを代入して、ゴムモデル５が接触部モデル９から受ける応力σを計算する。このような応力σを計算することにより、ゴムモデル５が接触部モデル９の接触面９Ｓをスライドする際に生じる摩擦力を求めることができる。

次に、応力σに基づいて、ゴムモデル５の各粒子６の新たな加速度ａを計算する（第２加速度計算ステップＳ５６）。この第２加速度計算ステップＳ５６では、ゴムモデル５の各粒子６に設定されている上記式（７）の運動方程式に、応力σを代入して、各粒子６の加速度ａを求めている。これにより、各粒子６には、ゴムモデル５が接触部モデル９から受ける応力σに基づく新たな加速度ａが設定される。各粒子６は、この新たな加速度ａに基づいてそれぞれ移動し、ゴムモデル５のスライドが計算される。

次に、コンピュータ３が、予め定められた終了時間が経過したか否かを判断する（終了判断ステップＳ５７）。この終了判断ステップＳ５７では、コンピュータ３が、終了時間が経過したと判断した場合、スライドステップＳ５を終了させる。一方、コンピュータ３が、終了時間が経過していないと判断した場合、単位時間Ｔｘを一つ進めて（ステップＳ５８）、接触バネ設定ステップＳ５３、歪計算ステップＳ５４、応力計算ステップＳ５５、第２加速度計算ステップＳ５６、及び、終了判断ステップＳ５７の各ステップが順次行われる。

これにより、コンピュータ３は、スライド開始から終了までのゴムモデル５の歪ε_ijや応力σを、単位時間Ｔｘごとの時系列データとして求めることができる。これらの時系列データは、コンピュータ３に記憶される。なお、終了時間は、実行するシミュレーションに応じて適宜定められる。

このように、本実施形態の接触シミュレーション方法では、ゴム材料１が、有限個の粒子６によってモデル化されるため、変形計算ステップＳ３において、例えば、引張変形が１００％以上となるような大きな変形計算を柔軟に行うことができる。

さらに、ゴムモデル５の粒子６は、隣り合う粒子６に拘束されないため、本実施形態のように接触面９Ｓに凹凸が形成されても、ゴムモデル５を接触面９Ｓに追従させて変形計算することができる。しかも、上記式（７）に示される運動方程式には、引張変形に対するゴム材料１の復元力が定義されるため、ゴムモデル５に大きな変形が与えられても、元の形状に復元することができる。従って、本実施形態の接触シミュレーション方法では、従来のような、要素潰れによる計算落ちを防ぎつつ、図１に示したゴム材料１と接触面２Ｓとの接触状態を精度良く解析することができる。

さらに、本実施形態では、図９に示されるように、接触バネ２０が定義されることにより、ゴムモデル５の粒子６が、接触部モデル９の粒子１０に浸透するのを防ぐことができる。従って、ゴムモデル５の各粒子６は、接触面９Ｓ内において、接触部モデル９の各粒子１０との接触及び離間をスムーズに行うことができ、ゴム材料１と接触面２Ｓとの接触状態を精度良く解析することができる。

図９に示されるように、ゴムモデル５の各粒子６、及び、接触部モデル９の各粒子１０の直径Ｄ２は、適宜設定することができるが、ＳＰＨ粒子法に従って、互いに同一にするのが望ましい。なお、解析精度を高めるには、各粒子６、１０の直径Ｄ２を小さくほどよい。また、上記のように、凹凸高さＨ１は適宜設定することができるが、ゴムモデル５との接触挙動を精度よく表現するために、凹部１６の下端１６ｄとの間に、５層以上の粒子１０が配置されるのが望ましい。

図１０に示されるように、本実施形態のスライドステップＳ５は、剛部材モデル７の上下方向（ｙ軸方向）の下面７ｄの座標Ｃを、一定に保持して行われてもよい。これにより、ゴムモデル５の上面５ｕは、凸部１５を乗り越える際に生じる上下動が、剛部材モデル７を介して拘束されるため、ゴムモデル５の変形量を大きくすることができる。なお、剛部材モデル７の座標Ｃは、シミュレーションの態様によって適宜変更することができる。

また、剛部材モデル７は、接触部モデル９と同様に、完全剛体（鋼鉄）として定義されてもよい。これにより、剛部材モデル７は、ゴムモデル５の上下動をより効果的に拘束することができるため、シミュレーション精度を向上することができる。

次に、ゴムモデル５と、接触部モデル９の接触面９Ｓとの接触状態が、良好か否かが判断される（ステップＳ６）。このステップＳ６では、スライドステップＳ５で単位時間Ｔｘ毎に記憶されたゴムモデル５の歪ε_ij、及び、接触部モデル９から受ける応力σ（摩擦力）に基づいて、ゴムモデル５の接触状態が評価される。

ゴムモデル５の接触状態が良好と判断された場合は、ゴムモデル５に設定された材料特性等に基づいて、ゴム材料１が製造される（ステップＳ７）。一方、ゴムモデル５の接触状態が良好ではないと判断された場合は、ゴムモデル５の材料特性等を変更して（ステップＳ８）、再度シミュレーションが行われる（Ｓ３〜Ｓ６）。このように、本実施形態では、ゴムモデル５の接触状態が良好と判断されるまでゴムモデル５が変更されるため、接触状態が良好なゴム材料１を設計することができる。

以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。

[実施例]
図３に示した処理手順に従って、図４に示したゴムモデル及び接触部モデルがコンピュータに設定された。そして、図１０に示されるように、ゴムモデルを接触部モデルに接触させ、かつ、スライド速度Ｖ１をそれぞれ異ならせて変形計算が実施された（実施例１〜実施例３）。実施例１〜３のゴムモデルには、下記に示すスライド速度Ｖ１が夫々設定された。

また、比較のために、図１で示したゴム材料及び接触部を、有限要素法に基づいて、有限個の要素でモデル化したゴムモデル及び接触部モデルが、コンピュータに設定された。そして、ゴムモデルを接触部モデルに接触させた変形計算が実施された（比較例）。なお、各パラメータは次の通りである。
実施例１のスライド速度Ｖ１：３m／秒
実施例２のスライド速度Ｖ１：１０m／秒
実施例３のスライド速度Ｖ１：３０m／秒
比較例のスライド速度Ｖ１：３m／秒
荷重Ｇ：０．１４ＭＰａ
隣接する粒子間の距離Ｄ１：１mm
粒子の密度ρ：０．９５ｇ／cm³
凹凸高さＨ１：５mm
影響半径：距離Ｄ１×２
バネ係数ｋ：１００Ｎ／mm
せん断バネ係数η：０ポアズ
バネの長さΔＬ_i：粒子１個分
伸縮量ΔC_i：解析の過程で計算される
時間増分Δｔ：１×１０^−４秒
重み係数Ｗ_i：非特許文献３に基づく３次Ｓｐｌｉｎｅ関数
Ｂｏｌｔｚｍａｎｎ定数Ｋ_B：１．３８０６６×１０^−２９
絶対温度Ｔ：２９６Ｋ
分子鎖一本あたりの平均セグメント数Ｎ：８個
応力・ひずみマトリクス[D]：材料定数から決定されるマトリクスである（例えば、非特許文献３参照）
なお、評価粒子の速度Ｖ_i、影響半径内にある他の粒子の速度Ｖ_j、評価粒子の座標ｘ_i、影響半径内にある他の粒子の座標ｘ_j、応力マトリクス{σ}、及びひずみマトリクス{ε}は、解析の過程で求められる値である。

テストの結果、実施例１〜３は、要素潰れ等による計算落ちが発生することなく、変形計算が終了した。一方、比較例は、要素潰れによる計算落ちが発生し、処理が中断された。

図１２には、応力（摩擦力）と時間との関係を示すグラフが示される。
このグラフに示されるように、スライド初期（〜０．００１秒）では、スライド速度Ｖ１が速いほど、応力σ（摩擦力）が大きくなることが確認できた。また、スライド中期（０．００１秒〜）においては、スライド速度Ｖ１が最も速い実施例３の応力σ（摩擦力）が大幅に小さくなることが確認できた。これらの結果より、ゴムモデルの応力σ（摩擦力）は、実際のゴム材料の摩擦力に近似することが確認できた。

図１３（ａ）には、実施例１のゴムモデルを視覚化した側面図が示される。また、図１３（ｂ）には、実施例３のゴムモデルを視覚化した側面図が示される。これらの図に示されるゴムモデルは、いずれもスライド初期のものである。これらの図から明らかなように、スライド速度Ｖ１が最も速い実施例３のゴムモデルは、スライド速度が最も遅い実施例１のゴムモデルに比べて反りが大きくなることが確認できた。これにより、実施例３のゴムモデルには、実施例１のゴムモデルよりも大きな摩擦力が生じていることが確認できた。従って、本発明の接触シミュレーション方法では、ゴム材料と接触面との接触状態を精度良く解析しうることが確認できた。

１ ゴム材料
２ 接触部
５ ゴムモデル
９ 接触部モデル

Claims (7)

1. コンピュータを用いて、ゴム材料と、該ゴム材料が接触する接触面との接触状態を、粒子法に基づいて解析するゴム材料の接触シミュレーション方法であって、
   前記コンピュータに、前記ゴム材料を、有限個の粒子でモデル化し、かつ、該粒子に運動方程式を定義したゴムモデルを入力するステップと、
   前記コンピュータに、前記接触面を含む接触部を、有限個の粒子でモデル化した接触部モデルを入力する接触部モデル入力ステップと、
   前記コンピュータが、前記ゴムモデルを前記接触部モデルに接触させて、該ゴムモデルの変形計算を行う変形計算ステップとを含み、
   前記ゴムモデルの前記運動方程式には、引張変形に対する前記ゴム材料の復元力が定義されることを特徴とするゴム材料の接触シミュレーション方法。
2. 前記変形計算ステップでは、前記ゴムモデルの前記粒子と前記接触部モデルの前記粒子との間に、前記ゴムモデルと前記接触部モデルとの接触面内において、前記ゴムモデルと前記接触部モデルとの接触を表現する接触バネを設定するステップを含む請求項１に記載のゴム材料の接触シミュレーション方法。
3. 前記接触部モデルは、前記接触面に、凹凸を有する請求項１又は２に記載のゴム材料の接触シミュレーション方法。
4. 前記変形計算ステップは、前記接触部モデル側への荷重を前記ゴムモデルに定義し、かつ該ゴムモデルを、該接触部モデルの前記接触面に沿ってスライドさせた状態を計算するステップを含む請求項１乃至３のいずれかに記載のゴム材料の接触シミュレーション方法。
5. 前記ゴムモデルは、前記ゴム材料の二次元の断面形状に形成され、
   前記ゴムモデルの前記接触面と反対側の面には、変形不能に定義された剛部材モデルが定義され、
   前記変形計算ステップでは、前記剛部材モデルの上下方向の座標を一定に保持する請求項４に記載のゴム材料の接触シミュレーション方法。
6. 前記ゴムモデルには、Ａｒｒｕｄａ−Ｂｏｙｃｅ型構成方程式が定義され、
   前記変形計算ステップは、前記運動方程式に基づいて、前記ゴムモデルの前記各粒子の加速度を計算する第１加速度計算ステップ、
   前記加速度で移動した前記ゴムモデルの歪を計算する歪計算ステップ、
   前記ゴムモデルの歪を、前記Ａｒｒｕｄａ−Ｂｏｙｃｅ型構成方程式に代入して、前記ゴムモデルが前記接触部モデルから受ける応力を計算する応力計算ステップ、
   前記応力を前記運動方程式に代入して、前記ゴムモデルの前記各粒子の加速度を計算する第２加速度計算ステップを含み、
   前記変形計算ステップが終了するまで、前記歪計算ステップ、前記応力計算ステップ、及び、前記第２加速度計算ステップが繰り返し実行される請求項１乃至５のいずれかに記載のゴム材料の接触シミュレーション方法。
7. 前記応力計算ステップでは、前記ゴムモデルの前記応力から、前記ゴムモデルと前記接触部モデルとの摩擦力を求める請求項６に記載のゴム材料の接触シミュレーション方法。