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

Description

本発明は、ゴム材料の大変形を精度良く再現するのに役立つシミュレーション方法に関する。

本件出願人は、下記特許文献１において、ゴム材料の大変形を精度良く再現するのに役立つシミュレーション方法を提案した。このシミュレーション方法では、ゴム解析モデルの最大主歪が予め定められた閾値を超えた場合、ゴム解析モデルに定義されているポアソン比を小さくする（変化させる）ステップを含むことを特徴としている。このようなステップは、ゴム材料の大変形時に生じる見かけ上の体積の増加を再現し、ひいては精度の良い変形の解析を可能としている。

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

特許文献１では、最大主歪とポアソン比との関係を規定したテーブルが、予めコンピュータに記憶されており、ポアソン比は、このテーブルに基づいて変更される。

ところで、最大主歪とポアソン比との関係については、ゴム材料のポリマーの分子構造によって異なることが予測される。しかしながら、特許文献１には、これらの関係をどのように決定するのかについては、具体的に記載されていない。

本発明は、以上のような実情に鑑み案出なされたもので、破壊直前のような大変形時の挙動をも精度良く再現するのに役立つゴム材料のシミュレーション方法を提供することを主たる目的としている。

本発明は、ゴム材料の変形計算を微小な時間増分で行うシミュレーション方法であって、 数値解析が可能な複数個の要素でゴム材料をモデル化したゴム解析モデルを設定する処理と、前記ゴム解析モデルが実質的に非圧縮特性を示すように、前記ゴム解析モデルの体積変化を規定する第１の材料パラメータの値を設定する処理と、少なくとも前記第１の材料パラメータと予め定めた解析条件とを用いて、前記ゴム解析モデルの第１の変形計算を行う処理と、前記第１の変形計算後の各要素の歪又は応力を計算する処理と、前記歪又は応力が予め定めた値を超えているときに、当該要素に実質的な体積変化が生じるように前記第１の材料パラメータの値を変更するパラメータ変更処理とを含み、前記第１の材料パラメータがポアソン比であり、前記パラメータ変更処理は、前記ポアソン比と圧力との関係を、前記ゴム材料に基づいて予め行われた分子動力学シミュレーションの結果に基づいて計算するステップと、前記関係を用いて、前記ゴム解析モデルの歪又は応力の増加に伴って前記ポアソン比を減少させるステップとを含むことを特徴とする。

本発明のゴム材料のシミュレーション方法において、ゴム解析モデルには、当初、実質的に非圧縮性を示すように（即ち、実質的に密度変化が生じないように）第１の材料パラメータの値が設定される。一方、変形後のゴム解析モデルの要素の歪又は応力が予め定めた値を超えているときに、当該要素に実質的な体積変化が生じるように該体積変化を規定する第１の材料パラメータの値を変更する処理が行われる。このため、ゴム解析モデルは、その変形計算において、変形が小さいときには非圧縮性特性を示す一方、変形が大きいときには、体積変化が生じさせることができる。従って、本発明のシミュレーション方法によれば、ゴム解析モデルは、現実のゴム材料が示す大変形時の見かけ上の体積変化を伴った形状変化を再現できる。従って、コンピュータシミュレーションであっても、ゴム材料の実物の変形挙動を精度良く再現することが可能になる。

また、本発明では、ゴム材料に基づいて予め行われた分子動力学シミュレーションの結果に基づいて、前記第１の材料パラメータの値が変更される。従って、ゴム材料のポリマーの分子構造や配合等に関係なく、上述の大変形をより精度良く再現しうる。

コンピュータ装置の一例を示す斜視図である。 本実施形態の処理手順の一例を示すフローチャートである。 （ａ）は解析対象のゴム材料の斜視図、（ｂ）はそれを対象としたゴム解析モデルの側面図である。 本実施形態の分子動力学シミュレーションの処理手順の一例を示すフローチャートである。 分子鎖モデルが視覚化された線図である。 ゴム材料モデルが視覚化された線図である。 ゴム材料モデルの時刻と密度との関係を示すグラフである。 ゴム材料モデルの圧力と体積との関係を示すグラフである。 ゴム材料モデルの圧力とポアソン比との関係を示すグラフである。 ゴム解析モデルのシミュレーション結果を示す側面図である。 実物のゴム材料の引張試験結果を示す側面図である。 比較例のゴム解析モデルのシミュレーション結果を示す側面図である。

以下、本発明の実施の一形態が図面に基づき説明される。
本発明では、ゴム材料の変形挙動を調べるためのゴム材料のシミュレーション方法が提供される。このようなシミュレーション方法は、図１に示されるようなコンピュータ装置１を用いて行われる。コンピュータ装置１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含んで構成される。本体１ａには、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリー、磁気ディスクなどの記憶装置及びディスクドライブ装置１ａ１、１ａ２などが適宜設けられる。記憶装置には、後述する方法を実行するための処理手順（プログラム）が予め記憶される。

図２には、本実施形態のシミュレーション方法のフローチャートを示す。
本実施形態では、先ずゴム解析モデルが設定される（ステップＳ１）。

ゴム解析モデルは、解析対象となるゴム材料（実在するか否かは問わない。）を、数値解析が可能な要素でモデル化することにより設定される。ゴム材料は、ゴムを主要部に含むものであればいかなる形状のものでも良い。例えば、ゴム材料は、全てがゴムからなる塊状のブロックやシートでも良いし、ゴムとコード材料とを複合させた複合体であっても良い。さらに、ゴム材料は、例えばゴルフボールや空気入りタイヤなどのより具体的なゴム製品（ともに図示せず）の形状を有するものでも良い。

図３（ａ）には、本実施形態で解析対象となるゴム材料２の斜視図を示す。ゴム材料２は、全体が薄肉六面体をなし、その中央に幅が小さい切り込みからなるスリット２ａが形成されている。本実施形態のシミュレーション方法では、このようなゴム材料２にそのスリット２ａを開くような引張荷重Ｐを与えたときの変形挙動が調べられる。ただし、ゴム材料２の形状等は、例示の態様に限定されるものではなく自由に定めることができるのは言うまでもない。

図３（ｂ）には、図３（ａ）のゴム材料２をモデル化したゴム解析モデル３の側面図を視覚化して示す。ゴム解析モデル３は、ゴム材料２の形状を表すように、有限個の小さい要素ｅを用いて分割（メッシュ分割）することにより設定される。従って、ゴム解析モデル３もスリット３ａを有している。メッシュ分割処理は、コンピュータ装置１上でメッシュ作成ソフトウエア等を使用して任意に行うことができる。

ゴム解析モデル３は、コンピュータ装置１で取り扱いが可能な数値データである。数値データとして、ゴム解析モデル３は、例えば、各要素ｅの節点座標値、節点番号及び要素番号などの情報を含んでいる。これらの数値データは、コンピュータ装置１での計算に用いられるとともに、ソフトウエアを用いて図３（ｂ）のように容易に視覚化され得る。

要素ｅには、例えば、四つの節点を有する四面体要素又は六つの節点を有する六面体要素といった三次元のソリッド要素が好適に用いられる。本実施形態のゴム解析モデル３は、全てが正六面体のソリッド要素でモデル化されている。切り欠き効果が現れるスリット３ａ周辺の変形をより詳細に調べるために、スリット３ａの周辺部は、他の部分に比して相対的に小さい要素で構成される（節点を密とする）ことが望ましい。また、ゴム材料２がコード層等の複合材を含む場合、要素ｅには、例えば直交異方性が定義されたシェル要素などの二次元平面要素が用いられても良い。

次に、図２に示したように、本実施形態では、ゴム解析モデル３に、材料特性と全体剛性マトリックスとが設定される（ステップＳ２）。これらは、コンピュータ装置１を使用してユーザーにより入力、設定される。

材料特性として、例えばゴム解析モデル３の各要素ｅの弾性率Ｅ及び密度ρが定義される。材料特性の主な値は、解析対象となるゴム材料２に即して定められ、その値が各要素ｅに割り当てられる。本実施形態のゴム解析モデル３は、均質な材料からなるため、初期値として、全ての要素ｅに同一の材料特性が定義される。

また、ゴム解析モデル３の材料特性として、体積変化に関する第１の材料パラメータが定義される。第１の材料パラメータとしては、例えば、ポアソン比が定義される。ポアソン比は、変形時のゴム解析モデル３の体積変化を規定することができる。ポアソン比の初期値として、ゴム解析モデル３が実質的に非圧縮特性を示す値に設定される。体積弾性率ｋ、弾性率Ｅ及びポアソン比νは、下式（１）の関係を満たす。
ｋ＝Ｅ／｛３（１−２ν）｝ …（１）

上記式（１）において、ポアソン比νを０．５とおくと、体積弾性率ｋが無限大になり、コンピュータ計算時にエラーが生じる。このような不具合を防止するために、ポアソン比の初期値として、０．５にきわめて近似した値、例えば０．４９以上かつ０．５０未満の範囲で設定されるのが望ましい。従って、「ゴム解析モデル３が実質的に非圧縮特性を示す」とは、ポアソン比が０．５である完全な非圧縮特性を意味するのではなく、上述のように０．５ではないがそれに近似したポアソン比を具える物体が示す特性で足りるものとする。

ゴム解析モデル３の全体剛性マトリックス「Ｋ］は、各要素ｅに設定される要素剛性マトリックス［Ｄ］をモデル全体について足し合わせたもので、要素剛性マトリックス［Ｋ］e を用いて下記式（２）により得ることができる。
［Ｋ］＝Σ［Ｋ］e …（２）

ゴム解析モデル３の要素剛性マトリックス［Ｋ］eは、下記式（３）のように、剛性マトリックス［Ｄ］とＢマトリックスとを用いて計算される。剛性マトリックス［Ｄ］は、要素の荷重に対する変形し難さ、即ち下記式（４）のように、応力σと歪εとの関係を示し前記材料特性に応じて決めることができる。さらに、Ｂマトリックスは下記式（５）のように、変位δと歪εとの関係を示すものである。
［Ｋ］e ＝∫［Ｂ］^-1［Ｄ］［Ｂ］ｄｖ …（３）
［σ］＝［Ｄ］｛ε｝ …（４）
｛ε｝＝［Ｂ］｛δ｝ …（５）

本実施形態では、剛性マトリックス［Ｄ］は、次式（６）のように第１の材料パラメータとしてのポアソン比を含むものとなる。

次に、図２に示したように、解析終了か否かがコンピュータ装置１によって判断される（ステップＳ３）。解析を終了させる条件は、必要に応じて種々定めることができ、例えば、荷重負荷時からの経過時間、ゴム解析モデル３の破断の有無等を挙げることができる。本実施形態では後者が採用される。そして、解析終了条件を満たす場合、処理を終える（ステップＳ３でＹ）。

一方、解析終了ではないと判断された場合（ステップＳ３でＮ）、本実施形態のシミュレーション方法では、境界条件が更新される（ステップＳ４）。更新される境界条件としては、引張荷重Ｐの大きさ、向き、荷重速度及び／又は変形計算のための時間増分などが挙げられる。

次に、ゴム解析モデル３の第１の変形計算が、例えば、有限要素法により行われる（ステップＳ５）。第１の変形計算では、ゴム解析モデル３の各要素ｅの節点変位、即ち、変形後のゴム解析モデル３の変形形状が計算される。全体剛性マトリックス［Ｋ］、各節点の変位を示す節点変位ベクトル｛δ｝及び外力ベクトル｛Ｐ｝は、下記式（７）の関係があるので、全体剛性マトリックスの逆マトリックス［Ｋ］^-1を求めることにより、節点変位ベクトル｛δ｝を得ることができる。
｛Ｐ｝＝［Ｋ］｛δ｝ …（７）

次に、ゴム解析モデル３の各要素ｅについて最大主歪が計算される（ステップＳ６）。最大主歪は、各要素ｅの様々な方向に生じる歪の大きさを代表的に示す。従って、最大主歪の値を調べることにより、各要素ｅの変形の大きさを定量的に把握することができる。しかし、このようパラメータに代えて、各要素ｅの最大応力や平均応力などの物理量が用いられても良い。

本実施形態のシミュレーション方法は、変形後のゴム解析モデル３において、ある要素ｅの最大主歪が予め定めた値を超えている場合、当該要素ｅに実質的な体積変化が生じるように第１の材料パラメータの値を変更する処理を行う。このために、本実施形態では、各要素ｅの最大主歪の大きさを予め定められた閾値と比較し、第１の材料パラメータ（ポアソン比）の値を変更する必要がある要素の有無が調べられる（ステップＳ７）。最大主歪の閾値については、任意に定めうるが、これまでの実験結果等に鑑み、例えば、３．０程度が設定されるのが望ましい

次に、ステップＳ７において、第１の材料パラメータを変更しなければならない要素（以下、このような要素を「対象要素」と呼ぶ。）が存在しないと判断された場合（ステップＳ７でＮ）、解析終了条件を満たすまで、ステップＳ４〜６が繰り返し行われる。これにより、引張荷重が作用した瞬間から、時間増分きざみで刻々と変化するゴム解析モデル３の形状が実質的な体積変化なしで計算されかつ出力される。

一方、ステップＳ７で対象要素が存在すると判断された場合（ステップＳ７でＹ）、該対象要素の第１の材料パラメータ（ポアソン比）が小さい値変更され（パラメータ変更処理）、他の境界条件を用いて全体剛性マトリックスが再度設定される（ステップＳ８）。そして、この新たに設定された全体剛性マトリックスを用いて変形形状の再計算、即ち第２の変形計算が行われる（ステップＳ９）。

式（６）で示したように、剛性マトリックス［Ｄ］には、ポアソン比が含まれている。従って、変形度合いが大きい要素については、ポアソン比を小さく変更することで、第２の変形計算時、第１の変形計算に比べて体積変化を大きく発現させることができる。即ち、現実のゴム材料が大変形するときに示すような見かけ上の体積変化を伴った形状変化を、本来変化しない第１の材料パラメータを可変とすることにより再現できる。従って、実物の変形挙動を精度良く再現することが可能になる。

一方、ゴム材料は、本質的に非圧縮性であり、体積弾性率も高く、実験的に変形と体積変化との関係を調べるのは容易ではない。そこで、本発明で行われる上記パラメータ変更処理（ステップＳ８）は、ゴム材料に基づいて予め行われた分子動力学シミュレーションの結果に基づいて、第１の材料パラメータの値を変更することとした。

このような分子動力学シミュレーションの一例については、図４に示されている。分子シミュレーションでは、先ず、コンピュータ装置１に、該コンピュータ装置１で取り扱い可能な分子鎖モデルが設定（入力）される（ステップＳ１００）。

図５には、分子鎖モデル４の一実施形態が視覚化されている。本実施形態の分子鎖モデル４は、解析対象のコム材料（例えば、cis−BR）の高分子をモデル化したものであり、複数の粒子モデル５と、これらの粒子モデル５、５間を接続する結合鎖モデル６とから構成された直鎖状の三次元構造を有する。

本実施形態の分子鎖モデル４は、一つの粒子モデル５が一つの原子に対応したフルアトムモデルとして、cis−BRの重合度３の１分子が示されている。ただし、分子鎖モデル４の重合度は任意に定めうる。また、分子鎖モデル４は、例えば、計算時間を短縮するために、数個の原子からなる集団が一つの粒子モデル５として粗視化された粗視化モデルで設定されても良い。また、結合鎖モデル６は、各粒子モデル５、５間を拘束するもので、例えば、平衡長とバネ定数とがコンピュータ装置１に入力される。

分子鎖モデル４の粒子モデル５は、それぞれ分子動力学法によるシミュレーションにおいて、運動方程式の質点として取り扱われる。従って、コンピュータ装置１には、各粒子モデル５の質量、体積、直径及び初期座標などが少なくとも入力される。

また、分子鎖モデル４には、相互作用（相互作用ポテンシャルエネルギーで、以下、単に「ポテンシャル」ということがある。）が定義される。ポテンシャルは、粒子モデル５，５間の距離の関数であって、２つの粒子モデル５、５の間に作用する力を計算する際に用いられる。例えば、粒子モデル５，５間の距離が予め定められたカットオフ距離よりも大きい場合、当該２つの粒子モデル５，５間には相互作用が働かないように定義される。一方、粒子モデル５，５間の距離が予め定められたカットオフ距離よりも小さい場合、当該２つの粒子モデル間５，５にはポテンシャルに基づいた力（本実施形態では斥力）が働くように定義されている。

次に、図４に示したように、分子鎖モデル４をセル内に配置したゴム材料モデル７がコンピュータ装置１に設定（入力）される（ステップＳ１０１）。

図６には、ゴム材料モデル７の一実施形態が視覚化されている。セルＳは、解析対象のゴム材料の微小構造部分に相当するもので、本実施形態では、微小な立方体として定義される。一つのセルＳの中には、上記分子鎖モデル４が複数本、ランダムに初期配置される。また、セルＳには、周期境界条件が適用されている。周期境界条件は、一つのセルＳの周囲に、それと同じ分子鎖モデル４の配置及び動きを持つセルＳが仮想的に連続するように配置されているものとして取り扱われる。これにより、計算上、連続したゴム材料の空間が提供され得る。

次に、図４に示したように、ゴム材料モデル７についての構造緩和計算が行われる（ステップＳ１０２）。構造緩和計算では、所定の時間、配置した全ての分子鎖モデル４が古典力学に従うものとしてニュートンの運動方程式が適用される。そして、各時刻における全ての分子鎖モデル４の動きが追跡される。これにより、ステップＳ１０１で決められた分子鎖モデル４の人為的な初期配置の影響を減らすことができるように、ゴム材料モデル７の構造が緩和される。構造緩和計算では、系内の圧力及び温度は一定で行われる。

ところで、ゴム材料モデル７の密度は、セルＳの容積と、その中に配置される分子鎖モデル４の質量によって決まる。従って、セルＳ内に配置されるべき分子鎖モデル４の数は、解析対象となるゴム材料の密度に基づいて決定される。しかし、ゴム材料モデル７の設定当初から、本来のゴム材料の密度となるように分子鎖モデル４が配置されていると、例えば、分子鎖モデル４同士の重なり等によって、計算途中で異常な力が計算され、ひいては、計算落ちが生じるおそれがある。このような不具合を無くすために、好ましい態様では、ゴム材料モデル７の初期の設定時には、ゴム材料モデル７の密度がゴム材料の密度よりも小さくなるように、分子鎖モデル４が配置されることが望ましい。そして、その後、セルＳの容積を少しずつ縮めながらゴム材料モデル７の密度を、解析対象のゴム材料の密度に調整することもできる。

図７には、密度が調整されたゴム材料モデル７の密度と時刻との関係を示すグラフが示されている。縦軸は、ゴム材料モデル７の密度、横軸は、時刻である。ゴム材料モデル７（系）の密度は、時刻とともに単調に増加した後、目的とする値に安定する（単調減少及び単調増加がなくなる）。ステップＳ１０２の構造緩和計算は、このようなゴム材料モデル７の密度の安定を確認した後、終了されるのが望ましい。他の実施形態として、構造緩和計算は、このような安定的な密度が得られるであろう一定の時間（例えば１０ps以上）を終えた時点で計算終了とされても良い。

次に、図４に示したように、ゴム材料モデル７について、圧力と体積との関係を取得するための計算が行われる（ステップＳ１０３）。上記ステップにおいて、設定されたゴム材料モデル７に対してある圧力が設定される。この圧力は、セルＳに作用する圧力であり、例えば、セルＳにｘ、ｙ及びｚ軸の各面に等しい力を加えることにより、調節され得る。そして、目標とする圧力の作用下でのゴム材料モデル７の体積が計算される。このような処理を、圧力を様々に変化させて行うことにより、図８に示されるような、ゴム材料モデル７の圧力と体積との関係を取得することができる。本実施形態では、セルＳを膨張させるように、セルＳに、標準大気圧よりも低い負圧を作用させている。

図８の横軸の圧力は、右側に向かうほど圧力が低くなる（真空度が良くなる）ものとして表示されている。従って、圧力が右側に向かうほど、ゴム材料モデル７は、より大きく３軸引張変形（膨張変形）させられる。図８から明らかなように、圧力の低下に従って、ゴム材料モデル７の体積が増加しており、ある圧力を境に、ゴム材料モデル７の体積が顕著に増加するという関係が得られた。このように、ゴム材料モデル７の体積の変化が急激に大きくなるのは、引張変形の過程において、ゴム材料モデル７にいわゆる空隙（ボイド）が生成され見かけの体積が膨張したためと考えられる。つまり、分子動力学シミュレーションの中でも、このようなゴム材料の特有の現象が正しく再現されている。このようなゴム材料モデル７の圧力と体積との関係は、ゴム材料の分子構造等によって異なる。

次に、図４に示したように、ゴム材料モデル７の圧力と体積弾性率との関係が計算される（ステップＳ１０４）。ゴム材料モデル７の圧力、体積弾性率及び体積歪との関係は、下記式（８）で表される（「体積歪」は、ゴム材料モデル７の初期体積に対する体積変化量の比である。）。
圧力 ＝ 体積弾性率 × 体積歪 … （８）

従って、ゴム材料モデル７の体積弾性率は、式（８）から、圧力をそのときの体積歪で除すことにより計算される。さらに、図８に示したように、予め取得されたゴム材料モデル７の圧力と体積との関係に基づいて、コンピュータ装置１は、ゴム材料モデル７の各圧力での体積弾性率を計算し、これを記憶する。

次に、図２のステップＳ７で「対象要素」と判定された要素については、上記の手順で予め得られているゴム材料モデル７の各圧力（３軸応力）と体積弾性率との関係に基づいて、その第１のパラメータであるポアソン比が計算される。ポアソン比は、上記式（１）の変形式から直ちに得られる。図９には、このようにして得られたゴム材料モデル７の圧力とポアソン比との関係の一例が実線で示されている。一方、仮想線は、分子動力学シミュレーションではなく、直感的に、ポアソン比が線形に小さくなるように設定された比較例が示されている。

図２に示したように、本実施形態のシミュレーション方法では、第２の変形計算の後、第２の変形計算結果と前記第１の変形計算の結果との差が計算され、その差が予め定めた許容範囲内か否かが判断される（ステップＳ１０）。この判断基準は、シミュレーションの内容に応じて適宜設定できる。一例として、ゴム解析モデル３の全節点数をＮ、第ｉ節点（ｉは１〜Ｎまでの整数）の第１の変形計算での変位をｄ１ｉ、第２の変形計算での変位をｄ２ｉとしたとき、変位の差｜ｄ１ｉ−ｄ２ｉ｜の最大値が、第１の変形計算での全要素の平均変位（Σｄ１ｉ／Ｎ）の所定パーセント以下の場合に許容範囲内とすることが望ましい。なお、数値範囲などは、材料に応じて変えることができるが、一例として５％程度が好適である。

そして、第２の変形計算結果と第１の変形計算の結果との差が許容範囲内の場合（ステップＳ１０でＹ）、十分な解析終了条件を満たすまでステップＳ４〜６が繰り返され、ゴム解析モデル３の変形形状が再度計算される。他方、計算結果の差が許容範囲外の場合（ステップＳ１０でＮ）、時間増分を小さくし（ステップＳ１１）、該対象要素の第１の材料パラメータ（ポアソン比）を変更前の値に戻し、他の境界条件を用いて全体剛性マトリックスが再度設定される（ステップＳ１２）。そして、ステップＳ４以降が実行される。

第１の材料パラメータを変更したことにより計算結果が急激に変化する場合、この変化の過程を細かく追跡しなければ正確な計算結果（正解）には到達できないおそれがある。従って、本実施形態のように、ステップＳ１０において、第１の変形計算の結果と第２の変形計算結果との差が許容範囲ではないと判断された場合、境界条件の更新増分量を小さくして剛性マトリックス等を再設定することにより、計算精度を高めることができる。なお、境界条件更新増分量を小さくする態様としては、例えば計算の時間増分を小さくすることが挙げられる。また、静的な釣合計算の場合には、例えば本来１mm変形させるところを０．５mmとし、荷重の負荷量を減らすことなどを挙げることができる。なお、本実施形態によれば、ステップＳ１０でＮの場合、先の境界条件に基づいて計算された第１の変形計算結果及び第２の変形計算結果はいずれも破棄されることになる。

以上説明したシミュレーション方法の具体例を示す。図１０には、図２の手順に従って行ったゴム解析モデルの引張試験のシミュレーション結果として、ゴム解析モデルの変形形状を示す側面図である。また、図１１には、実物のゴム材料にシミュレーションと同様の引張試験を行ったときの側面図を示す。いずれも、引張荷重は、できるだけゆっくりかつ徐々に増大するように作用させ、モデルが破断することにより、シミュレーションを停止させた、他の条件は次の通りである。
ゴム解析モデルの形状：高さ３０mm×水平長さ３０mm×厚さ１mmの薄肉六面体
スリット：切り込み深さ１０mm
要素の弾性率：５ＭＰａ
密度：１．１ｇ／cm³

図１０及び図１１から明らかなように、本実施形態のシミュレーションでは、実物のゴム材料と非常に近似した変形形状が得られていることが確認できる。他方、図１２には、ポアソン比を図９の比較例の関係で変化させた同一時刻の比較例のシミュレーション結果を示す。このものでは、スリットの変形形状が、実物のそれと大きく異なっていることが分かる。

以上本発明のシミュレーション方法について説明したが、本発明はその要旨を逸脱しない範囲において種々の変更が可能である。例えば、上記実施形態では、第１の材料パラメータとして、ポアソン比を用いたが、その逆数であるポアソン数を採用しても良いのは言うまでもない。この場合、ポアソン数は、要素の変形の増大に伴って増加するように変更される。また、第１の材料パラメータとして、体積弾性率を採用することもできる。この場合、体積弾性率の初期の値は、ポアソン比０．４９に相当する程度に定められ、好ましくはゴム解析モデル３の要素ｅに定義される弾性率（縦弾性率）の２０倍以上かつ１０００倍以下に設定されるのが望ましい。そして、パラメータ変更処理では、ポアソン比の場合と同様に、ゴム解析モデル３の歪又は応力の増大に伴って、前記体積弾性率を減少させることができる。

１ コンピュータ装置
２ ゴム材料
３ ゴム解析モデル
４ 分子鎖モデル
５ 粒子モデル
６ 結合鎖モデル
７ ゴム材料モデル
Ｓ セル

Claims (1)

1. ゴム材料の変形計算を微小な時間増分で行うシミュレーション方法であって、
   数値解析が可能な複数個の要素でゴム材料をモデル化したゴム解析モデルを設定する処理と、
   前記ゴム解析モデルが実質的に非圧縮特性を示すように、前記ゴム解析モデルの体積変化を規定する第１の材料パラメータの値を設定する処理と、
   少なくとも前記第１の材料パラメータと予め定めた解析条件とを用いて、前記ゴム解析モデルの第１の変形計算を行う処理と、
   前記第１の変形計算後の各要素の歪又は応力を計算する処理と、
   前記歪又は応力が予め定めた値を超えているときに、当該要素に実質的な体積変化が生じるように前記第１の材料パラメータの値を変更するパラメータ変更処理とを含み、
   前記第１の材料パラメータがポアソン比であり、
   前記パラメータ変更処理は、前記ポアソン比と圧力との関係を、前記ゴム材料に基づいて予め行われた分子動力学シミュレーションの結果に基づいて計算するステップと、
   前記関係を用いて、前記ゴム解析モデルの歪又は応力の増加に伴って前記ポアソン比を減少させるステップとを含むことを特徴とするゴム材料のシミュレーション方法。