(57)【特許請求の範囲】
【請求項１】
粘弾性材料の変形をシミュレートする方法であって、
前記粘弾性材料を数値解析が可能な要素で分割した粘弾性材料モデルを設定するステップと、
前記粘弾性材料モデルに条件を設定して変形計算を行うステップと、
前記変形計算から必要な物理量を取得するステップとを含み、
前記粘弾性材料モデルには、ひずみを増加させる負荷変形及び該ひずみを減少させる除荷変形からなる変形１サイクルにおいて、前記ひずみ量に基づいた第１のエネルギーロスを生成させる第１のエネルギーロス特性と、ひずみ速度に基づいた第２のエネルギーロスを生成させる第２のエネルギーロス特性とが定義されるとともに、
少なくとも２サイクルの変形計算を行う場合、１サイクル目において、前記第１及び第２のエネルギーロスを生成させる処理と、
２サイクル目において、ひずみが１サイクル目の最大ひずみを超えない場合、前記第２のエネルギーロスのみを生成させる処理とを含むことを特徴とする粘弾性材料のシミュレーション方法。
【請求項２】
前記粘弾性材料モデルには、下記式（１）で示される応力とひずみとの関係が定義されるとともに、
前記第１のエネルギーロスは、下記式（１）における１本の分子鎖当たりの平均セグメント数Ｎを、負荷変形時ではひずみ量に基づいて増加させ、かつ、除荷変形時では除荷を開始したときの値を維持させることにより生成されることを特徴とする請求項１記載の粘弾性材料のシミュレーション方法。
【数１】
0004697870.tif 000008
【請求項３】
前記粘弾性材料モデルは、ゴムをモデル化したゴムモデルと、前記ゴムに囲まれるフィラーを前記ゴムモデルで囲まれるようにモデル化したフィラーモデルとを含み、
かつ前記ゴムモデルは、前記フィラーモデルの周囲を囲む厚さが小さい界面モデルと、この界面モデルを取り囲むマトリックスモデルとからなり、しかもその少なくとも一方に、前記第１及び第２のエネルギーロス特性が定義されることを特徴とする請求項１又は２記載の粘弾性材料のシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ゴム等の粘弾性材料の変形状態等を精度良く解析するのに役立つ粘弾性材料のシミュレーション方法に関する。
【背景技術】
【０００２】
ゴム材料は、タイヤ、スポーツ用品、その他各種の工業製品に広く使用される。これらの試作の手間とコストとを減じるために、ゴム材料の変形状態などを予めコンピュータを用いて計算しかつ予測することが行われている。従来、このような方法として、下記の非特許文献１ないし２が知られている。
【０００３】
【非特許文献１】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)
【０００４】
【非特許文献２】J. S. BERGSTROM and M. C. BOYCE 著「 CONSTITUTIVE MODELING OF THE LARGE STRAIN TIME-DEPENDENT BEHAVIOR OF ELASTOMERS」 Journal of the Mechanics and Physics of Solids Volume 46, No.5 PP. 931-954 (1998)
【０００５】
ところで、ゴム材料にはマリンス（Mullins）効果があることが知られている。ここで、マリンス効果について簡単に述べる。図１５に示されるように、先ずゴムをひずみεmまで引っ張った後、除荷する（予備伸張）。このとき、応力−ひずみ曲線は、経路０、Ｒ１、Ｒ２、０を通ってエネルギーロスが生じる。次に、この予備伸張されたゴムを再度引っ張ると、応カーひずみ曲線は、経路０、Ｒ３、Ｒ２、０を通り、予備伸張時のそれとは一致しない低い応力を示す。一方、そのゴムをひずみεmを超えてεnまでさらに引っ張ると、経路０、Ｒ３、Ｒ４、Ｒ５、０を通る。マリンス効果は、以上のように、負荷及び除荷という変形の１サイクルを繰り返すと、２サイクル目以降はゴムが軟化する現象である。
【０００６】
しかしながら、前記非特許文献１ないし２は、分子鎖網目理論（後で詳しく述べる。）等を用いることにより、ゴム材料を解析することが記載されているが、これらの文献に記載された方法では、上述のマリンス効果を再現できない。
【０００７】
特に、空気入りタイヤ等のゴム製品では、変形サイクルが周期的に繰り返される。従って、それらに用いる粘弾性材料のシミュレーションを精度良く行うためには、このようなマリンス効果を無視することは妥当ではない。
【発明の開示】
【発明が解決しようとする課題】
【０００８】
本発明は、以上のような問題点に鑑み案出なされたもので、粘弾性材料モデルに、ひずみを増加させる負荷変形及び該ひずみを減少させる除荷変形からなる変形１サイクルにおいて、前記ひずみ量に基づいてエネルギーロスを生成させる第１のエネルギーロス特性と、ひずみ速度に基づいてエネルギーロスを生成させる第２のエネルギーロス特性とを定義するとともに、２サイクルの変形計算を行う場合に、これらのエネルギーロスを制御することを基本として、いわゆるマリンス効果を再現しうる粘弾性材料のシミュレーション方法を提供することを目的としている。
【課題を解決するための手段】
【０００９】
本発明のうち請求項１記載の発明は、粘弾性材料の変形をシミュレートする方法であって、前記粘弾性材料を数値解析が可能な要素で分割した粘弾性材料モデルを設定するステップと、前記粘弾性材料モデルに条件を設定して変形計算を行うステップと、前記変形計算から必要な物理量を取得するステップとを含み、前記粘弾性材料モデルには、ひずみを増加させる負荷変形及び該ひずみを減少させる除荷変形からなる変形１サイクルにおいて、前記ひずみ量に基づいた第１のエネルギーロスを生成させる第１のエネルギーロス特性と、ひずみ速度に基づいた第２のエネルギーロスを生成させる第２のエネルギーロス特性とが定義されるとともに、少なくとも２サイクルの変形計算を行う場合、１サイクル目において、前記第１及び第２のエネルギーロスを生成させる処理と、２サイクル目において、ひずみが１サイクル目の最大ひずみを超えない場合、前記第２のエネルギーロスのみを生成させる処理とを含むことを特徴とする。
【００１０】
また請求項２記載の発明は、前記粘弾性材料モデルには、下記式（１）で示される応力とひずみとの関係が定義されるとともに、前記第１のエネルギーロスは、下記式（１）における１本の分子鎖当たりの平均セグメント数Ｎを、負荷変形時ではひずみ量に基づいて増加させ、かつ、除荷変形時では除荷を開始したときの値を維持させることにより生成されることを特徴とする請求項１記載の粘弾性材料のシミュレーション方法である。
【数１】
0004697870.tif 000002
【００１１】
また請求項３記載の発明は、前記粘弾性材料モデルは、ゴムをモデル化したゴムモデルと、前記ゴムに囲まれるフィラーを前記ゴムモデルで囲まれるようにモデル化したフィラーモデルとを含み、かつ前記ゴムモデルは、前記フィラーモデルの周囲を囲む厚さが小さい界面モデルと、この界面モデルを取り囲むマトリックスモデルとからなり、しかもその少なくとも一方に、前記第１及び第２のエネルギーロス特性が定義されることを特徴とする請求項１又は２記載の粘弾性材料のシミュレーション方法である。
【発明の効果】
【００１２】
本発明の粘弾性材料のシミュレーション方法によれば、２サイクル目において、マリンス効果が再現される。従って、現実の加硫ゴムの挙動に近い精度の良いシミュレーション結果を得ることができる。
【発明を実施するための最良の形態】
【００１３】
以下、本発明の実施の一形態を図面に基づき説明する。
図１には、本発明のシミュレーション方法を実施するためのコンピュータ装置１が示される。該コンピュータ装置１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含む。本体１ａは、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリー及び磁気ディスク（いずれも図示せず）が設けられる他、ＣＤ−ＲＯＭやフレキシブルディスクなどのドライブ装置１ａ１、１ａ２が設けられる。そして、前記磁気ディスクには、本発明のシミュレーション方法を実行するために必要なプログラムが記憶される。
【００１４】
図２には、本実施形態のシミュレーション方法の処理手順の一例が示される。本実施形態では、先ず粘弾性材料モデル２が設定される（ステップＳ１）。
【００１５】
前記粘弾性材料としては特に限定されないが、本実施形態では、ゴムと、その中に配合されたフィラーとを含むフィラーが充填された加硫ゴムの場合が示される。ただし、解析対象のゴムないし他の粘弾性材料は、実在するか否かは問わない。粘弾性材料モデルは、このゴムに基づいて設定される。図３には、粘弾性材料モデル２の微視構造であるユニットセルＵの一例が視覚化して示される。
【００１６】
前記粘弾性材料モデル２のユニットセルＵは、解析しようとするゴムの微小領域を、有限個の小さな要素２ａ、２ｂ、２ｃ…に分割（置換）することにより設定される。微小領域としては、解析を行いうるサイズであれば特に限定されないが、この例では、縦横それぞれ３００×３００ナノメータの矩形領域として設定される。
【００１７】
各要素２ａ、２ｂ、２ｃ…は、数値解析が可能に定義される。ここで、数値解析が可能とは、例えば有限要素法、有限体積法、差分法又は境界要素法といった数値解析法により各要素ないし系全体についての変形計算が可能なことを意味する。
【００１８】
具体的には、各要素２ａ、２ｂ、２ｃ…について、座標系における節点の座標値、要素の形状、材料特性などが数値データにて定義される。各要素２ａ、２ｂ、２ｃ…には、例えば２次元平面要素として三ないし四辺形の要素、３次元要素としては、例えば４ないし６面体の要素が用いられる。本実施形態では、前者が用いられる。これらにより、粘弾性材料モデル２は、前記コンピュータ装置１にて変形計算が可能なデータをなす。
【００１９】
この実施形態の粘弾性材料モデル２は、後述する変形シミュレーションにおいて平面ひずみ状態で解析が行われる。したがって、Ｚ軸方向（図３において紙面と直角方向）にはひずみを持たないものとする。
【００２０】
また、本実施形態において、粘弾性材料モデル２は、フィラーがモデル化されたフィラーモデル３と、ゴムがモデル化されたゴムモデル４とを含む。さらに、本実施形態のゴムモデル４は、前記フィラーモデル３の周囲に配されかつ厚さが小さい界面モデル５と、この界面モデル５を取り囲むマトリックスモデル６とを含んでいる。図３では、理解しやすいように、フィラーモデル３は白色、界面モデル５はやや薄い灰色及びマトリックスモデル６は最も濃い灰色でそれぞれ視覚化されている。
【００２１】
前記フィラーモデル３は、フィラーとして、例えばカーボンブラックに基づいて設定される。但しフィラーは、上記の具体例に限定されるものではなく、例えばシリカ等の他のフィラーであっても良いのは言うまでもない。また、フィラーモデル３の形状は、例えば、電子顕微鏡にて撮像した実際のゴム中のカーボンブラックの形状に基づいて定められるのが望ましい。
【００２２】
図４には、電子顕微鏡にて撮像されたカーボンブラックの斜視図が示される。図のように、カーボンブラックＣＢは、炭素原子からなる直径数１０ナノメータ程度の球状の一次粒子７が不規則に３次元的に結合した葡萄の房のような構造を持つ。本実施形態のフィラーモデル３は、この構造の平面視を基準として設定されている。
【００２３】
またカーボンブラックは、ゴム等の粘弾性材料に比べると数百倍の硬さを持つ。このため、本実施形態において、フィラーモデル３は弾性体として、エネルギーロスは生じないものとして取り扱われる。また、セルユニットＵの中で占めるフィラーモデル３の面積ないし体積等は、解析対象のゴムのフィラー充填量などに基づいて定めるのが望ましい。
【００２４】
本実施形態において、前記界面モデル５は、フィラーモデル３を連続して取り囲み、かつ、小さい厚さで設定される。ゴムとフィラーとは界面で物理的ないし化学的に結合しており、他の部分では見られない様々な現象、例えば大きな滑りやこれに伴う摩擦が生じていると推察されている。従って、このようなモデル５を設定し、そこに適切なパラメータを与えることは、正確なシミュレーションを行う上で役に立つ。また、界面モデル５は必ずしも連続してフィラーモデル３を連続して取り囲む必要はない。さらに、界面モデル５の厚さｔは特に限定はされないが、種々の実験の結果、好ましくは１〜２０ナノメートル、より好ましくは５〜１０ナノメートルが望ましい。
【００２５】
前記マトリックスモデル６は、界面モデル５の外側で該界面モデル５を取り囲むように設定される。ユニットセルＵにおいて、マトリックスモデル６は、フィラーモデル３及び界面モデル５を除いた部分であり、ゴムモデル４の主要部を構成している。
【００２６】
本実施形態では、ゴムモデル４に、第１のエネルギーロス特性及び第２のエネルギーロス特性が定義される。
【００２７】
図５には、このようなエネルギーロス特性を有するゴムモデル４の変形原理を説明する概念図が示される。ゴムモデル４の個々の要素は、第１の変形ネットワークＡと、第２の変形ネットワークＢとが並列に結合されたものと等価なもとして考えられる。系全体に発生する応力Ｓは、第１の変形ネットワークＡに生じた応力と、第２の変形ネットワークＢに生じた応力との和である。従って、個々のネットワークＡ、Ｂの応力を求めることにより、系全体の応力が計算できる。なお、両ネットワークＡ、Ｂは並列に結合されているため、各々のネットワークＡ及びＢに生じるひずみ（伸び）はともに同一になる。以下、各変形ネットワークＡ、Ｂに基づいたゴムモデル４の応力計算手順について述べる。
【００２８】
＜第１の変形ネットワークＡ＞
前記第１の変形ネットワークＡは、バネと等価なものとして考えられ、その応力は、バネの伸び（ひずみ）に応じた値になる。このバネの伸びは、ひずみ速度の影響を受けない。従って、第１の変形ネットワークＡは、ひずみ速度に依存しない応力を示す。
【００２９】
第１の変形ネットワークＡの応力の計算には、分子鎖網目理論に基づく８鎖モデルが用いられる。その応力は、速度形式において例えば下記式（１）のものが適用できる。
【００３０】
【数１】
0004697870.tif 000003
【００３１】
図６（Ａ）、（Ｂ）に示されるように、分子鎖網目理論は、連続体としてのゴム材料ａは微視構造ｄａとして、無秩序に配向された分子鎖ｃが接合点ｂで連結された網目構造を持つ、と考える理論である。前記接合点ｂは、例えば分子間の化学的結合等であってそれには架橋点等が含まれる。
【００３２】
図６（Ｃ）に示されるように、１本の分子鎖ｃは、複数のセグメントｄから構成されるものとする。一つのセグメントｄは、分子鎖網目理論においては繰り返しの最小構成単位である。また一つのセグメントｄは、化学的には同図（Ｄ）に示されるように、炭素原子が共有結合によって連結した複数個のモノマーｆが連結したものと等価と考える。個々の炭素原子は、原子同士の結合軸の周りで互いに自由に回転しうるため、セグメントｄは全体として曲がりくねるなど様々な形態をとることができる。
【００３３】
また分子鎖網目理論では、接合点ｂが原子の揺らぎ周期に対して長時間的には平均位置が変化しないものと仮定し、接合点ｂの回りの摂動は無視される。さらに二つの接合点ｂ、ｂを両端に持つ分子鎖ｃの端−端ベクトル（end-to-end vector）は、それが埋め込まれているゴム材料の連続体と共に変形するものとする。
【００３４】
また図７に示されるように、分子鎖網目理論は、粘弾性材料を巨視的には微小な８鎖モデルｇが集合した立方体状の網目構造体ｈとして考える。一つの８鎖モデルｇは、右側に拡大して示されるように、立方体の中心に定められた一つの接合点ｂ１から、各頂点に設けられた８つの接合点ｂ２にそれぞれのびる分子鎖ｃを有するものとし、これらをベースとして応力等の計算が行われる。この８鎖モデルについては、前記非特許文献１において述べられている。
【００３５】
しかし、本実施形態では、第１の変形ネットワークＡにおいて、非特許文献１などが提案する一般的な８鎖モデルを、ひずみ量に基づいた第１のエネルギーロスが生成されるように、以下のような修正を加えて用いられる。
【００３６】
ゴム材料は、荷重負荷における変形過程において、数百％にも達し得る大きな歪が許容される。このメカニズムは、複雑に絡み合った前記分子鎖ｃの絡み部分（接合点ｂ）がほどける（消滅する）ことにより得られると考えられる。具体的に述べると、図８（Ａ）に示されるように、一つの接合点ｂで絡み合っている分子鎖ｃ１ないしｃ４に矢印方向の引張応力が作用すると、各分子鎖ｃ１ないしｃ４は伸び、接合点ｂは大きな歪を受けて破損（消滅）すると考えられる。そして、同図（Ｂ）に示されるように、これまで２本であった分子鎖ｃ１及びｃ２は、あたかも１本の長い分子鎖ｃ５のように振る舞う。分子鎖ｃ３及びｃ４についても同様である。そして、この接合点ｂの破損によってエネルギーロスが生じると考えられる。本実施形態のシミュレーションでは、このような分子鎖ｃの接合点ｂの破損（消滅）による影響が前記式（１）に取り入れられる。
【００３７】
ここで、８鎖モデルにおける分子鎖ｃのからみの減少と前記セグメントｄとの関係について考える。図７に示した網目構造体ｈは、幅方向、高さ方向及び奥行き方向にそれぞれ８鎖モデルｇがｋ個結合したものとする（ただし、ｋは、十分に大きい数とする。）。いま、網目構造体ｈに含まれる接合点ｂの総数を「からみ数」として符号ｍで表すと、からみ数ｍは、式（２）で表される。同様に、網目構造体ｈに含まれる分子鎖ｃの数、即ち、ゴムモデル４の単位体積中に含まれる分子鎖の数ｎは、式（３）で表し得る。
ｍ＝（ｋ＋１）3 ＋ｋ3…式（２）
ｎ＝８ｋ3 …式（３）
【００３８】
ここで、ｋは十分に大きい数のため、その３次項以外を省略すると、上式（２）は次式（４）で表し得る。
ｍ＝２ｋ3 …式（４）
さらに上記式（３）及び（４）の関係から、からみ数ｍは、ｎを用いて、式（５）で表し得る。
ｎ＝４ｍ …式（５）
【００３９】
ここで、ゴムモデル４は、変形してもその単位体積中に含まれる分子鎖のセグメントの総数ＮA には変化がないから、１本の分子鎖に含まれる平均セグメント数をＮとすると、式（６）が成り立つ。
ＮA ＝ｎ・Ｎ …式（６）
【００４０】
式（６）をＮについて解くと、式（７）が得られる。
Ｎ＝ＮA ／ｎ …式（７）
【００４１】
また式（７）及び式（５）から式（８）が得られる。
Ｎ＝ＮA ／４ｍ …式（８）
【００４２】
式（８）から明らかなように、負荷変形によって、ゴムの分子鎖ｃのからみ数ｍが減少していくと、１本の分子鎖ｃに含まれる平均セグメント数Ｎは増加することになる。従って、負荷、言い換えるとひずみ量に応じて、式（１）の平均セグメント数Ｎを増加させることにより、分子鎖ｃのからみの減少を式（１）の中に取り込むことができる。これは、ゴムの変形メカニズムをより精度良く擬似化でき、計算精度を向上させる。ここで、負荷変形とは、微小時間の間で粘弾性材料モデル２のひずみが増大する変形であり、逆にひずみが減少する変形を除荷変形とする。
【００４３】
前記平均セグメント数Ｎは、ひずみ量に基づいて様々な方法で増加させることができる。例えば、平均セグメント数Ｎは、負荷変形時におけるひずみ量自体又はそれに関連したパラメータに基づいて増大するものが望ましい。具体的なパラメータとしては、ひずみの１次の不変量Ｉ1 などを挙げることができる。本実施形態では、前記平均セグメント数Ｎを、下記式（９）で定める。式（９）は、平均セグメント数Ｎが、各要素個々において、それぞれ歪の１次の不変量Ｉ1 （より詳しくはその平方根であるパラメータλc ）の関数であることを示す。
【００４４】
【数２】
0004697870.tif 000004
【００４５】
式（９）の上記ＡないしＥは、いずれも定数である。これは、ゴム試験片の単純な１軸引張試験などの実測結果から容易に定めることができる。例えば、先ず解析対象となるゴム材料の応力−歪曲線を得る。そして、その荷重除荷時の曲線に沿うように前記ｎ、Ｎを定める。これにより、分子鎖のセグメントの総数ＮA （＝ｎ・Ｎ）が決まる。次に、荷重負荷時の曲線に整合するよう、各ひずみにおける平均セグメント数Ｎを求める。そして、決定された負荷時の平均セグメントＮに一致するよう、式（９）のパラメータＡないしＥが決定される。
【００４６】
図９には、前記平均セグメント数Ｎとパラメータλc との関係が示される。この例では、上記定数が以下のように設定されている。
Ａ＝＋２．９４９３
Ｂ＝−５．８０２９
Ｃ＝＋５．５２２０
Ｄ＝−１．３５８２
Ｅ＝＋０．１３２５
ひずみ量に関するパラメータであるλc が大きくなると、平均セグメント数Ｎは滑らかに増大する。この例では、λc の上限は約２．５である。後述する変形シミュレーションにおいては、ゴムモデル４の各要素について、負荷変形時ではλcが常に計算される。計算されたλc は、式（９）に代入され、当該要素の当該歪状態における平均セグメント数Ｎが計算される。平均セグメント数Ｎは、所定のタイミングで計算され、式（１）に取り込まれる。
【００４７】
また、ゴムモデル４の除荷変形時では、１本の分子鎖ｃに含まれる平均セグメント数Ｎは、その除荷を開始した時の平均セグメント数Ｎが、ひずみが零になるまで維持される。従って、除荷時には、分子鎖のからみが解ける（又は消滅する）ことがないため、ゴムモデル４は、負荷時に比べて低い応力で変形できる。これにより、除荷時の軟化現象がシミュレーションに再現される。また、このような第１の変形ネットワークＡでは、負荷時と除荷時とでは、同一のひずみ状態において応力が異なる。これにより、応力−ひずみ曲線において、前記第１のエネルギーロスが生成され、これがヒステリシスループを描く。
【００４８】
このように、前記式（１）における１本の分子鎖当たりの平均セグメント数Ｎを、負荷変形時ではひずみ量に基づいて増加させ、かつ、除荷変形時では除荷を開始したときの値を維持させることにより、第１の変形ネットワークＡに、ひずみ量に基づいた第１のエネルギーロスが生成される第１のエネルギーロス特性が定義される。
【００４９】
＜第２の変形ネットワークＢ＞
前記第２の変形ネットワークＢは、バネｅとダッシュポットｐとの直列接続体と等価なものとして考えられる。ダッシュポットｐは、ニュートンの粘性法則に従う流体（油等）が封止されたシリンダ筒の中を、オリフィス等が設けられたピストンが移動する減衰装置であり、ピストンの移動速度に比例した抵抗が生じる。また、その抵抗はエネルギーロスになる。
【００５０】
第２の変形ネットワークＢの全体の伸び（ひずみ）は、バネｅの伸びとダッシュポットの伸び（ひずみ）との和になる。また、第２の変形ネットワークＢに生じる応力も、バネｅの伸び（ひずみ）に応じた値になる。しかし、バネｅの伸びは、荷重によって単純に決定されず、ひずみ速度に依存したダッシュポットｐの抵抗の影響を受ける。つまり、第２の変形ネットワークＢは、ひずみ速度に依存した応力を示す。
【００５１】
ここで、分子鎖網目理論に基づく８鎖モデルは、第２の変形ネットワークＢのバネｅにも適用される。従って、前記バネｅに生じる応力も、前記式（１）を用いて計算される。この計算において、前記１本の分子鎖ｃに含まれる平均セグメント数Ｎは、第１の変形ネットワークＡの場合のように変化するものでも良く、また固定値であっても良い。
【００５２】
また、第２の変形ネットワークＢには、第２のエネルギーロス特性が定義される。該第２のエネルギーロス特性は、ひずみ速度に基づいて第２のエネルギーロスを生成させる。即ち、第２の変形ネットワークＢでは、ひずみ速度に依存した抵抗がダッシュポットｐに生じ、これがバネｅの伸長比λc に影響する。伸長比λc は、式（１）のただし書きのように、分子鎖ｃの各方向の主ストレッチから計算できる。
【００５３】
次に、第２の変形ネットワークＢにおいて、各主ストレッチλi (Be)は、次の手順により計算される。なおλの添字ｉは各方向を表し（ｉ＝１、２…）、二次元であればｘ，ｙを意味する。また右肩の添字（Be）は第１の変形ネットワークＡの主ストレッチと区別するための符号である。
【００５４】
先ず第２の変形ネットワークＢにおける相当せん断応力τ(B) が、式（１０）で計算される。式（１０）のσ(B)' は、前ステップの計算で得られた第２の変形ネットワークＢの相当応力である。また、ダッシュポットｐのひずみ速度が、相当せん断応力τ(B)及びλc (Bp)の値を用いて式（１１）から計算される。なお式（１１）のλc (Bp)は、ダッシュポットの伸長比であり、この値は、前回の計算値λi(Bp)（後述の式（１４）で計算される。）から式（１）のただし書きの方法に準じて計算される。
【００５５】
【数３】
0004697870.tif 000005
【００５６】
また前記ひずみ速度から、ダッシュポットｐの変形速度Ｄ(Bp)が、式（１２）より計算され、該変形速度Ｄ(Bp)からダッシュポットｐの伸びが式（１３）及び（１４）で計算される。ここで、添字ｔは計算の時間ステップを示し、Δｔはその増分である。
【００５７】
【数４】
0004697870.tif 000006
【００５８】
そして、第２の変形ネットワークＢのバネｅの各ストレッチλi (Be)は、下記式（１５）から得られる。なお式（１５）中、λi は、第２の変形ネットワークＢの全体の伸長比である。
【００５９】
【数５】
0004697870.tif 000007
【００６０】
前記λc (Be)は、式（１）に代入され、これにより、第２の変形ネットワークＢの応力が計算される。この応力は、ダッシュポットの作用と等価なメカニズムにより、ひずみ速度に依存したものになる。また、この変形ネットワークＢの応力−ひずみ線図では、ひずみ速度に依存したエネルギーロスを持つように、ヒステリシスループを描く。ひずみ速度の依存性については、前記式（１１）の各材料定数を適宜決定することにより調節できる。
【００６１】
そして、第１の変形ネットワークＡの応力と第２の変形ネットワークＢの応力との和によってゴムモデル４の全体の応力が計算される。
【００６２】
このような第１及び第２のエネルギーロス特性は、ゴムモデル４のうち、界面モデル５又はマトリックスモデル６の少なくとも一方に与えられれば良いが、好ましくは両方に定義されることが望ましい。
【００６３】
また界面モデル５と前記マトリックスモデル６とは、異なる粘弾性特性が定義されるのが望ましい。ここで言う粘弾性特性とは、図１０に示されるように、任意の変形速度における応力−ひずみ曲線で表される特性を含む。例えば図１０のように、界面モデル５をマトリックスモデル６よりも軟らかく設定することが望ましい。これにより、界面付近でのひずみを相対的に増加させ、界面でのゴムとフィラーとの間の大きなすべりを疑似化するのに役立つ。またヒステリシスループの面積が異なるようにも構成できる。この態様では、同一応力の場合、界面モデル５のひずみは、マトリックスモデル６よりも大きくなるように粘弾性特性が定められている。ただし、このような態様に限定されるものではない。
【００６４】
次に本実施形態のシミュレーション方法では、材料モデル２を用いて変形シミュレーションが行われる（ステップＳ３）。変形シミュレーションの具体的な処理手順は、図１１に示される。変形シミュレーションでは、先ず各種データがコンピュータ装置１に入力される（ステップＳ３１）。入力されるデータには、各要素に定義された節点の位置や材料特性といった情報が含まれる。
【００６５】
コンピュータ装置１では、入力されたデータに基づいて各要素の剛性マトリックスが作成され（ステップＳ３２）、しかる後、全体構造の剛性マトリックスが組み立てられる（ステップＳ３３）。全体構造の剛性マトリックスには、既知節点の変位、節点力が導入され（ステップＳ３４）、剛性方程式の解析が行われる。そして、未知節点変位が決定され（ステップＳ３５）、前記式などを用いて各要素の応力、主応力及び／又はひずみといった物理量が計算されかつ出力される（ステップＳ３６ないし３７）。この物理量は、適宜コンピュータ装置１に記憶される。ステップＳ３８では、計算を終了させるか否かの判定がなされ、否定的である場合には、ステップＳ３２以降が繰り返される。
【００６６】
なお、シミュレーションは、例えば有限要素法を用いたエンジニアリング系の解析アプリケーションソフトウエア（例えば米国リバモア・ソフトウェア・テクノロジー社で開発・改良されたＬＳ−ＤＹＮＡ等）を用いて行うことができる。
【００６７】
また本シミュレーションは、均質化法（漸近展開均質化法）に基づいて行われる態様を例示する。均質化法は、解析対象領域が任意の微視構造の繰り返しによって構成され、その繰り返し度合いが非常に密なために直接有限要素法で領域を離散化出来ない場合、解析対象を均質な等価モデルで代用して全体を解析し、その解析結果を任意の点での微視構造に戻すことによって微視構造自身の変形が近似的に求められる。
【００６８】
具体例で述べると、均質化法では、図１２に示されるように、材料モデル２のユニットセルＵを周期的に持っている材料全体Ｍを表現するｘI と、前記微視構造を表現するｙI との独立した２変数が用いられる。微視的スケールと巨視的スケールという異なる尺度の場におけるそれぞれ独立した変数を漸近展開することにより、ユニットセルＵのモデル構造を反映させた材料全体Ｍの平均的な力学応答を求め得る。漸近展開均質化法については、例えば次の文献に詳細に述べられている。
【００６９】
Higa,Y.and Tomita,Y,,Computational Prediction of Mechanical Properties of Nickel-based superalloy with gamma Prime Phase Precipitates,Proceedings of ICM8(Victoria,B.C.,Canada),Advance Materials and Modeling of Mechanical Behavior,(Edited by Ellyin,F,and Proven,J.W.),III(1999),1061-1066,Fleming Printing Ltd..
比嘉吉一，冨田佳宏，粒子強化型複合材の均質化法による変形挙動のモデル化とシミュレーション，日本機械学会論文集，A66(2000),1441-1446.
【００７０】
本実施形態では、式（１）及び（１０）の定数等を次のように設定して変形シミュレーションが行われた。
ＣR ＝０．２６８
Ｎ＝６．６
Ｔ＝２９６
ｋB ＝１．３８０６６×１０-29
ＮA ＝ｎ・Ｎ＝４．３２８×１０26
フィラーモデルの体積含有率 ３０％
フィラーモデルの縦弾性係数Ｅ：１００ＭＰａ
フィラーモデルのポアソン比ν：０．３
【００７１】
さらに変形シミュレーションは、材料全体Ｍの解析領域に一軸引張変形を発生させた。引張条件は、図１２のｘ2方向にひずみ速度＝０．０１（／ｓ）とし、４サイクルの変形が順番に行われた。各サイクルの引張り量は次の通りとした。
１サイクル目：１００％（伸長比２．０）
２サイクル目：１００％（伸長比２．０）
３サイクル目：２００％（伸長比３．０）
４サイクル目：２００％（伸長比３．０）
【００７２】
また変形シミュレーションでは、各ひずみ状態において、ゴムモデル４の平均セグメント数Ｎが式（９）により計算され、その値を式（１）へ逐次代入した。また材料モデルは、厚さ方向（図３のＺ軸方向）に変化しないように拘束条件が与えられた。さらに、界面モデル５及びマトリックスモデル６の平均セグメント数Ｎは、次のように設定した。
【００７３】
＜マトリックスモデル＞
・負荷変形時の平均セグメント数
Ｎ＝-3.2368+20.6175λc-21.8168λc2+10.8227λc3 -1.9003λc4
・除荷変形時の平均セグメント数
Ｎ＝６．６
・分子鎖のセグメントの総数（一定）
ＮA ＝４．３２８１×１０26
【００７４】
＜界面モデル＞
・負荷変形時の平均セグメント数
Ｎ＝-5.9286+20.6175λc-21.8168λc2+10.8227λc3 -1.9003λc4
・除荷変形時の平均セグメント数
Ｎ＝３．９１
・分子鎖のセグメントの総数（一定）
ＮA ＝３．２０３×１０25
【００７５】
前記変形計算が行われると、その結果から、各計算ステップのタイミング毎に必要な物理量を得る（ステップＳ４）。図１３（Ａ）、（Ｂ）には、粘弾性材料モデルの１ないし２サイクル目の応力−伸長比の関係を示すグラフが、また、図１４（Ａ）、（Ｂ）には、粘弾性材料モデルの３ないし４サイクル目の応力−伸長比の関係を示すグラフが示される。
【００７６】
図から明らかなように、粘弾性材料モデル２は、２サイクル目以降、マリンス効果が再現されていることが確認できる。
【００７７】
即ち、１サイクル目の応力−伸長比のグラフが描くループは、第１のエネルギーロスＥ１及び第２のエネルギーロスＥ２をともに含んでいる。即ち、負荷変形では、概念的には、ひずみ量に基づいて分子鎖ｃの結合点ｂの破損が生じるので、１本の分子鎖当たりの平均セグメント数が増加し、第１のエネルギーロスが生成される。また、変形ネットワークＢでは、ひずみ速度に基づいたダッシュポットの抵抗（エネルギーロス）が生じるので、これに基づいた第２のエネルギーロスが生成される。
【００７８】
また、図１３において、１サイクル目と同じひずみ量の２サイクル目では、第１のエネルギーロスが生成されず、第２のエネルギーロスだけが生成されていることがわかる。これは、２サイクル目のひずみ量が１サイクル目の最大ひずみを超えていないので、分子鎖ｃの結合点ｂが新たに破損しないことに基づいている、従って、この実施形態では、２サイクル目の負荷変形は、１サイクル目の除荷変形と同じ経路を通ることになる。他方、２サイクル目においても、そのひずみ速度に基づいた第２のエネルギーロスは生じるので、図１３（Ｂ）のように、応力−伸長比曲線は、１サイクル目の除荷変形と同じ経路を通ることになる。
【００７９】
また、３サイクル目の曲線は、負荷変形時、伸長比が１ないし２サイクル目と等しい２．０までの間は、２サイクル目と同じ経路を通る、しかし、伸長比が２．０を超えると、新たな分子鎖ｃの結合点ｂが破損するため、応力が上昇し、ひずみ量に基づいた第１のエネルギーロスが再び生成されることになる。また、１ないし２サイクル目と同様、ひずみ速度に応じて第２のエネルギーロスも生じ得る。
【００８０】
また、４サイクル目では、２サイクル目と同様、ひずみが３サイクル目の最大ひずみを超えないので、分子鎖ｃの結合点ｂが新たに破損して消滅することはない。このため、４サイクル目の負荷変形は、３サイクル目の除荷変形と同じ経路を通る。このため、ひずみ量に基づいた第１のエネルギーロスは生成されず、ひずみ速度に応じた前記第２のエネルギーロスのみが生じることになる。
【００８１】
以上のように、本実施形態のシミュレーション方法によれば、粘弾性材料、とりわけ加硫ゴムに見られるマリンス効果をシミュレーション上に正確に再現することができる。従って、繰り返し荷重を受けて周期的に変形するような空気入りタイヤのゴム材料の性能評価をより精度良く行うことができる。
【図面の簡単な説明】
【００８２】
【図１】本実施形態で用いたコンピュータ装置の一例を示す斜視図である。
【図２】本実施形態の処理手順を示すフローチャートである。
【図３】粘弾性材料モデルの微視構造の一実施形態を示す線図である。
【図４】カーボンブラックの形状を示す線図である。
【図５】ゴムモデルの変形原理を説明する概念図である。
【図６】（Ａ）は粘弾性材料、（Ｂ）はその分子鎖１構造を説明する線図、（Ｃ）は１本の分子鎖の拡大図、（Ｄ）はセグメントの拡大図である。
【図７】粘弾性材料の網目構造体及びそれを構成する８鎖モデルの一例を示す斜視図である。
【図８】（Ａ）、（Ｂ）は分子鎖の接合点の破断を説明する線図である。
【図９】パラメータλc と分子鎖１本当たりの平均セグメント数Ｎとの関係を示すグラフである。
【図１０】任意の変形速度における界面モデル及びマトリックスモデルの応力−ひずみ曲線である。
【図１１】変形シミュレーションの手順を示すフローチャートである。
【図１２】均質化法を説明する全体構造とそのユニットセルとの関係を示す。
【図１３】（Ａ）は変形シミュレーションの１サイクル目、（Ｂ）はその２サイクル目を示し、それぞれ応力−伸長比の関係を示すグラフである。
【図１４】（Ａ）は変形シミュレーションの１サイクル目、（Ｂ）はその２サイクル目を示し、それぞれ応力−伸長比の関係を示すグラフである。
【図１５】ゴム材料のマリンス効果を説明する応力ひずみ線図である。
【符号の説明】
【００８３】
１ コンピュータ装置
２ 材料モデル
３ フィラーモデル
４ ゴムモデル
５ 界面モデル
６ マトリックスモデル
【図１】
0004697870.tif 000009
【図２】
0004697870.tif 000010
【図３】
0004697870.tif 000011
【図４】
0004697870.tif 000012
【図５】
0004697870.tif 000013
【図６】
0004697870.tif 000014
【図７】
0004697870.tif 000015
【図８】
0004697870.tif 000016
【図９】
0004697870.tif 000017
【図１０】
0004697870.tif 000018
【図１１】
0004697870.tif 000019
【図１２】
0004697870.tif 000020
【図１３】
0004697870.tif 000021
【図１４】
0004697870.tif 000022
【図１５】
0004697870.tif 000023