JP4697870B2

JP4697870B2 - 粘弾性材料のシミュレーション方法

Info

Publication number: JP4697870B2
Application number: JP2005295316A
Authority: JP
Inventors: 正登内藤
Original assignee: Sumitomo Rubber Industries Ltd
Current assignee: Sumitomo Rubber Industries Ltd
Priority date: 2005-10-07
Filing date: 2005-10-07
Publication date: 2011-06-08
Anticipated expiration: 2025-10-07
Also published as: JP2007101499A

Description

本発明は、ゴム等の粘弾性材料の変形状態等を精度良く解析するのに役立つ粘弾性材料のシミュレーション方法に関する。

ゴム材料は、タイヤ、スポーツ用品、その他各種の工業製品に広く使用される。これらの試作の手間とコストとを減じるために、ゴム材料の変形状態などを予めコンピュータを用いて計算しかつ予測することが行われている。従来、このような方法として、下記の非特許文献１ないし２が知られている。

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までさらに引っ張ると、経路０、Ｒ３、Ｒ４、Ｒ５、０を通る。マリンス効果は、以上のように、負荷及び除荷という変形の１サイクルを繰り返すと、２サイクル目以降はゴムが軟化する現象である。

しかしながら、前記非特許文献１ないし２は、分子鎖網目理論（後で詳しく述べる。）等を用いることにより、ゴム材料を解析することが記載されているが、これらの文献に記載された方法では、上述のマリンス効果を再現できない。

特に、空気入りタイヤ等のゴム製品では、変形サイクルが周期的に繰り返される。従って、それらに用いる粘弾性材料のシミュレーションを精度良く行うためには、このようなマリンス効果を無視することは妥当ではない。

本発明は、以上のような問題点に鑑み案出なされたもので、粘弾性材料モデルに、ひずみを増加させる負荷変形及び該ひずみを減少させる除荷変形からなる変形１サイクルにおいて、前記ひずみ量に基づいてエネルギーロスを生成させる第１のエネルギーロス特性と、ひずみ速度に基づいてエネルギーロスを生成させる第２のエネルギーロス特性とを定義するとともに、２サイクルの変形計算を行う場合に、これらのエネルギーロスを制御することを基本として、いわゆるマリンス効果を再現しうる粘弾性材料のシミュレーション方法を提供することを目的としている。

本発明のうち請求項１記載の発明は、粘弾性材料の変形をシミュレートする方法であって、前記粘弾性材料を数値解析が可能な要素で分割した粘弾性材料モデルを設定するステップと、前記粘弾性材料モデルに条件を設定して変形計算を行うステップと、前記変形計算から必要な物理量を取得するステップとを含み、前記粘弾性材料モデルには、ひずみを増加させる負荷変形及び該ひずみを減少させる除荷変形からなる変形１サイクルにおいて、前記ひずみ量に基づいた第１のエネルギーロスを生成させる第１のエネルギーロス特性と、ひずみ速度に基づいた第２のエネルギーロスを生成させる第２のエネルギーロス特性とが定義されるとともに、少なくとも２サイクルの変形計算を行う場合、１サイクル目において、前記第１及び第２のエネルギーロスを生成させる処理と、２サイクル目において、ひずみが１サイクル目の最大ひずみを超えない場合、前記第２のエネルギーロスのみを生成させる処理とを含むことを特徴とする。

また請求項２記載の発明は、前記粘弾性材料モデルには、下記式（１）で示される応力とひずみとの関係が定義されるとともに、前記第１のエネルギーロスは、下記式（１）における１本の分子鎖当たりの平均セグメント数Ｎを、負荷変形時ではひずみ量に基づいて増加させ、かつ、除荷変形時では除荷を開始したときの値を維持させることにより生成されることを特徴とする請求項１記載の粘弾性材料のシミュレーション方法である。

また請求項３記載の発明は、前記粘弾性材料モデルは、ゴムをモデル化したゴムモデルと、前記ゴムに囲まれるフィラーを前記ゴムモデルで囲まれるようにモデル化したフィラーモデルとを含み、かつ前記ゴムモデルは、前記フィラーモデルの周囲を囲む厚さが小さい界面モデルと、この界面モデルを取り囲むマトリックスモデルとからなり、しかもその少なくとも一方に、前記第１及び第２のエネルギーロス特性が定義されることを特徴とする請求項１又は２記載の粘弾性材料のシミュレーション方法である。

本発明の粘弾性材料のシミュレーション方法によれば、２サイクル目において、マリンス効果が再現される。従って、現実の加硫ゴムの挙動に近い精度の良いシミュレーション結果を得ることができる。

以下、本発明の実施の一形態を図面に基づき説明する。
図１には、本発明のシミュレーション方法を実施するためのコンピュータ装置１が示される。該コンピュータ装置１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含む。本体１ａは、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリー及び磁気ディスク（いずれも図示せず）が設けられる他、ＣＤ−ＲＯＭやフレキシブルディスクなどのドライブ装置１ａ１、１ａ２が設けられる。そして、前記磁気ディスクには、本発明のシミュレーション方法を実行するために必要なプログラムが記憶される。

図２には、本実施形態のシミュレーション方法の処理手順の一例が示される。本実施形態では、先ず粘弾性材料モデル２が設定される（ステップＳ１）。

前記粘弾性材料としては特に限定されないが、本実施形態では、ゴムと、その中に配合されたフィラーとを含むフィラーが充填された加硫ゴムの場合が示される。ただし、解析対象のゴムないし他の粘弾性材料は、実在するか否かは問わない。粘弾性材料モデルは、このゴムに基づいて設定される。図３には、粘弾性材料モデル２の微視構造であるユニットセルＵの一例が視覚化して示される。

前記粘弾性材料モデル２のユニットセルＵは、解析しようとするゴムの微小領域を、有限個の小さな要素２ａ、２ｂ、２ｃ…に分割（置換）することにより設定される。微小領域としては、解析を行いうるサイズであれば特に限定されないが、この例では、縦横それぞれ３００×３００ナノメータの矩形領域として設定される。

各要素２ａ、２ｂ、２ｃ…は、数値解析が可能に定義される。ここで、数値解析が可能とは、例えば有限要素法、有限体積法、差分法又は境界要素法といった数値解析法により各要素ないし系全体についての変形計算が可能なことを意味する。

具体的には、各要素２ａ、２ｂ、２ｃ…について、座標系における節点の座標値、要素の形状、材料特性などが数値データにて定義される。各要素２ａ、２ｂ、２ｃ…には、例えば２次元平面要素として三ないし四辺形の要素、３次元要素としては、例えば４ないし６面体の要素が用いられる。本実施形態では、前者が用いられる。これらにより、粘弾性材料モデル２は、前記コンピュータ装置１にて変形計算が可能なデータをなす。

この実施形態の粘弾性材料モデル２は、後述する変形シミュレーションにおいて平面ひずみ状態で解析が行われる。したがって、Ｚ軸方向（図３において紙面と直角方向）にはひずみを持たないものとする。

また、本実施形態において、粘弾性材料モデル２は、フィラーがモデル化されたフィラーモデル３と、ゴムがモデル化されたゴムモデル４とを含む。さらに、本実施形態のゴムモデル４は、前記フィラーモデル３の周囲に配されかつ厚さが小さい界面モデル５と、この界面モデル５を取り囲むマトリックスモデル６とを含んでいる。図３では、理解しやすいように、フィラーモデル３は白色、界面モデル５はやや薄い灰色及びマトリックスモデル６は最も濃い灰色でそれぞれ視覚化されている。

前記フィラーモデル３は、フィラーとして、例えばカーボンブラックに基づいて設定される。但しフィラーは、上記の具体例に限定されるものではなく、例えばシリカ等の他のフィラーであっても良いのは言うまでもない。また、フィラーモデル３の形状は、例えば、電子顕微鏡にて撮像した実際のゴム中のカーボンブラックの形状に基づいて定められるのが望ましい。

図４には、電子顕微鏡にて撮像されたカーボンブラックの斜視図が示される。図のように、カーボンブラックＣＢは、炭素原子からなる直径数１０ナノメータ程度の球状の一次粒子７が不規則に３次元的に結合した葡萄の房のような構造を持つ。本実施形態のフィラーモデル３は、この構造の平面視を基準として設定されている。

またカーボンブラックは、ゴム等の粘弾性材料に比べると数百倍の硬さを持つ。このため、本実施形態において、フィラーモデル３は弾性体として、エネルギーロスは生じないものとして取り扱われる。また、セルユニットＵの中で占めるフィラーモデル３の面積ないし体積等は、解析対象のゴムのフィラー充填量などに基づいて定めるのが望ましい。

本実施形態において、前記界面モデル５は、フィラーモデル３を連続して取り囲み、かつ、小さい厚さで設定される。ゴムとフィラーとは界面で物理的ないし化学的に結合しており、他の部分では見られない様々な現象、例えば大きな滑りやこれに伴う摩擦が生じていると推察されている。従って、このようなモデル５を設定し、そこに適切なパラメータを与えることは、正確なシミュレーションを行う上で役に立つ。また、界面モデル５は必ずしも連続してフィラーモデル３を連続して取り囲む必要はない。さらに、界面モデル５の厚さｔは特に限定はされないが、種々の実験の結果、好ましくは１〜２０ナノメートル、より好ましくは５〜１０ナノメートルが望ましい。

前記マトリックスモデル６は、界面モデル５の外側で該界面モデル５を取り囲むように設定される。ユニットセルＵにおいて、マトリックスモデル６は、フィラーモデル３及び界面モデル５を除いた部分であり、ゴムモデル４の主要部を構成している。

本実施形態では、ゴムモデル４に、第１のエネルギーロス特性及び第２のエネルギーロス特性が定義される。

図５には、このようなエネルギーロス特性を有するゴムモデル４の変形原理を説明する概念図が示される。ゴムモデル４の個々の要素は、第１の変形ネットワークＡと、第２の変形ネットワークＢとが並列に結合されたものと等価なもとして考えられる。系全体に発生する応力Ｓは、第１の変形ネットワークＡに生じた応力と、第２の変形ネットワークＢに生じた応力との和である。従って、個々のネットワークＡ、Ｂの応力を求めることにより、系全体の応力が計算できる。なお、両ネットワークＡ、Ｂは並列に結合されているため、各々のネットワークＡ及びＢに生じるひずみ（伸び）はともに同一になる。以下、各変形ネットワークＡ、Ｂに基づいたゴムモデル４の応力計算手順について述べる。

＜第１の変形ネットワークＡ＞
前記第１の変形ネットワークＡは、バネと等価なものとして考えられ、その応力は、バネの伸び（ひずみ）に応じた値になる。このバネの伸びは、ひずみ速度の影響を受けない。従って、第１の変形ネットワークＡは、ひずみ速度に依存しない応力を示す。

第１の変形ネットワークＡの応力の計算には、分子鎖網目理論に基づく８鎖モデルが用いられる。その応力は、速度形式において例えば下記式（１）のものが適用できる。

図６（Ａ）、（Ｂ）に示されるように、分子鎖網目理論は、連続体としてのゴム材料ａは微視構造ｄａとして、無秩序に配向された分子鎖ｃが接合点ｂで連結された網目構造を持つ、と考える理論である。前記接合点ｂは、例えば分子間の化学的結合等であってそれには架橋点等が含まれる。

図６（Ｃ）に示されるように、１本の分子鎖ｃは、複数のセグメントｄから構成されるものとする。一つのセグメントｄは、分子鎖網目理論においては繰り返しの最小構成単位である。また一つのセグメントｄは、化学的には同図（Ｄ）に示されるように、炭素原子が共有結合によって連結した複数個のモノマーｆが連結したものと等価と考える。個々の炭素原子は、原子同士の結合軸の周りで互いに自由に回転しうるため、セグメントｄは全体として曲がりくねるなど様々な形態をとることができる。

また分子鎖網目理論では、接合点ｂが原子の揺らぎ周期に対して長時間的には平均位置が変化しないものと仮定し、接合点ｂの回りの摂動は無視される。さらに二つの接合点ｂ、ｂを両端に持つ分子鎖ｃの端−端ベクトル（end-to-end vector）は、それが埋め込まれているゴム材料の連続体と共に変形するものとする。

また図７に示されるように、分子鎖網目理論は、粘弾性材料を巨視的には微小な８鎖モデルｇが集合した立方体状の網目構造体ｈとして考える。一つの８鎖モデルｇは、右側に拡大して示されるように、立方体の中心に定められた一つの接合点ｂ１から、各頂点に設けられた８つの接合点ｂ２にそれぞれのびる分子鎖ｃを有するものとし、これらをベースとして応力等の計算が行われる。この８鎖モデルについては、前記非特許文献１において述べられている。

しかし、本実施形態では、第１の変形ネットワークＡにおいて、非特許文献１などが提案する一般的な８鎖モデルを、ひずみ量に基づいた第１のエネルギーロスが生成されるように、以下のような修正を加えて用いられる。

ゴム材料は、荷重負荷における変形過程において、数百％にも達し得る大きな歪が許容される。このメカニズムは、複雑に絡み合った前記分子鎖ｃの絡み部分（接合点ｂ）がほどける（消滅する）ことにより得られると考えられる。具体的に述べると、図８（Ａ）に示されるように、一つの接合点ｂで絡み合っている分子鎖ｃ１ないしｃ４に矢印方向の引張応力が作用すると、各分子鎖ｃ１ないしｃ４は伸び、接合点ｂは大きな歪を受けて破損（消滅）すると考えられる。そして、同図（Ｂ）に示されるように、これまで２本であった分子鎖ｃ１及びｃ２は、あたかも１本の長い分子鎖ｃ５のように振る舞う。分子鎖ｃ３及びｃ４についても同様である。そして、この接合点ｂの破損によってエネルギーロスが生じると考えられる。本実施形態のシミュレーションでは、このような分子鎖ｃの接合点ｂの破損（消滅）による影響が前記式（１）に取り入れられる。

ここで、８鎖モデルにおける分子鎖ｃのからみの減少と前記セグメントｄとの関係について考える。図７に示した網目構造体ｈは、幅方向、高さ方向及び奥行き方向にそれぞれ８鎖モデルｇがｋ個結合したものとする（ただし、ｋは、十分に大きい数とする。）。いま、網目構造体ｈに含まれる接合点ｂの総数を「からみ数」として符号ｍで表すと、からみ数ｍは、式（２）で表される。同様に、網目構造体ｈに含まれる分子鎖ｃの数、即ち、ゴムモデル４の単位体積中に含まれる分子鎖の数ｎは、式（３）で表し得る。
ｍ＝（ｋ＋１）³ ＋ｋ³…式（２）
ｎ＝８ｋ³ …式（３）

ここで、ｋは十分に大きい数のため、その３次項以外を省略すると、上式（２）は次式（４）で表し得る。
ｍ＝２ｋ³ …式（４）
さらに上記式（３）及び（４）の関係から、からみ数ｍは、ｎを用いて、式（５）で表し得る。
ｎ＝４ｍ …式（５）

ここで、ゴムモデル４は、変形してもその単位体積中に含まれる分子鎖のセグメントの総数Ｎ_A には変化がないから、１本の分子鎖に含まれる平均セグメント数をＮとすると、式（６）が成り立つ。
Ｎ_A ＝ｎ・Ｎ …式（６）

式（６）をＮについて解くと、式（７）が得られる。
Ｎ＝Ｎ_A ／ｎ …式（７）

また式（７）及び式（５）から式（８）が得られる。
Ｎ＝Ｎ_A ／４ｍ …式（８）

式（８）から明らかなように、負荷変形によって、ゴムの分子鎖ｃのからみ数ｍが減少していくと、１本の分子鎖ｃに含まれる平均セグメント数Ｎは増加することになる。従って、負荷、言い換えるとひずみ量に応じて、式（１）の平均セグメント数Ｎを増加させることにより、分子鎖ｃのからみの減少を式（１）の中に取り込むことができる。これは、ゴムの変形メカニズムをより精度良く擬似化でき、計算精度を向上させる。ここで、負荷変形とは、微小時間の間で粘弾性材料モデル２のひずみが増大する変形であり、逆にひずみが減少する変形を除荷変形とする。

前記平均セグメント数Ｎは、ひずみ量に基づいて様々な方法で増加させることができる。例えば、平均セグメント数Ｎは、負荷変形時におけるひずみ量自体又はそれに関連したパラメータに基づいて増大するものが望ましい。具体的なパラメータとしては、ひずみの１次の不変量Ｉ1 などを挙げることができる。本実施形態では、前記平均セグメント数Ｎを、下記式（９）で定める。式（９）は、平均セグメント数Ｎが、各要素個々において、それぞれ歪の１次の不変量Ｉ1 （より詳しくはその平方根であるパラメータλc ）の関数であることを示す。

式（９）の上記ＡないしＥは、いずれも定数である。これは、ゴム試験片の単純な１軸引張試験などの実測結果から容易に定めることができる。例えば、先ず解析対象となるゴム材料の応力−歪曲線を得る。そして、その荷重除荷時の曲線に沿うように前記ｎ、Ｎを定める。これにより、分子鎖のセグメントの総数Ｎ_A （＝ｎ・Ｎ）が決まる。次に、荷重負荷時の曲線に整合するよう、各ひずみにおける平均セグメント数Ｎを求める。そして、決定された負荷時の平均セグメントＮに一致するよう、式（９）のパラメータＡないしＥが決定される。

図９には、前記平均セグメント数Ｎとパラメータλc との関係が示される。この例では、上記定数が以下のように設定されている。
Ａ＝＋２．９４９３
Ｂ＝−５．８０２９
Ｃ＝＋５．５２２０
Ｄ＝−１．３５８２
Ｅ＝＋０．１３２５
ひずみ量に関するパラメータであるλ_c が大きくなると、平均セグメント数Ｎは滑らかに増大する。この例では、λ_c の上限は約２．５である。後述する変形シミュレーションにおいては、ゴムモデル４の各要素について、負荷変形時ではλ_cが常に計算される。計算されたλ_c は、式（９）に代入され、当該要素の当該歪状態における平均セグメント数Ｎが計算される。平均セグメント数Ｎは、所定のタイミングで計算され、式（１）に取り込まれる。

また、ゴムモデル４の除荷変形時では、１本の分子鎖ｃに含まれる平均セグメント数Ｎは、その除荷を開始した時の平均セグメント数Ｎが、ひずみが零になるまで維持される。従って、除荷時には、分子鎖のからみが解ける（又は消滅する）ことがないため、ゴムモデル４は、負荷時に比べて低い応力で変形できる。これにより、除荷時の軟化現象がシミュレーションに再現される。また、このような第１の変形ネットワークＡでは、負荷時と除荷時とでは、同一のひずみ状態において応力が異なる。これにより、応力−ひずみ曲線において、前記第１のエネルギーロスが生成され、これがヒステリシスループを描く。

このように、前記式（１）における１本の分子鎖当たりの平均セグメント数Ｎを、負荷変形時ではひずみ量に基づいて増加させ、かつ、除荷変形時では除荷を開始したときの値を維持させることにより、第１の変形ネットワークＡに、ひずみ量に基づいた第１のエネルギーロスが生成される第１のエネルギーロス特性が定義される。

＜第２の変形ネットワークＢ＞
前記第２の変形ネットワークＢは、バネｅとダッシュポットｐとの直列接続体と等価なものとして考えられる。ダッシュポットｐは、ニュートンの粘性法則に従う流体（油等）が封止されたシリンダ筒の中を、オリフィス等が設けられたピストンが移動する減衰装置であり、ピストンの移動速度に比例した抵抗が生じる。また、その抵抗はエネルギーロスになる。

第２の変形ネットワークＢの全体の伸び（ひずみ）は、バネｅの伸びとダッシュポットの伸び（ひずみ）との和になる。また、第２の変形ネットワークＢに生じる応力も、バネｅの伸び（ひずみ）に応じた値になる。しかし、バネｅの伸びは、荷重によって単純に決定されず、ひずみ速度に依存したダッシュポットｐの抵抗の影響を受ける。つまり、第２の変形ネットワークＢは、ひずみ速度に依存した応力を示す。

ここで、分子鎖網目理論に基づく８鎖モデルは、第２の変形ネットワークＢのバネｅにも適用される。従って、前記バネｅに生じる応力も、前記式（１）を用いて計算される。この計算において、前記１本の分子鎖ｃに含まれる平均セグメント数Ｎは、第１の変形ネットワークＡの場合のように変化するものでも良く、また固定値であっても良い。

また、第２の変形ネットワークＢには、第２のエネルギーロス特性が定義される。該第２のエネルギーロス特性は、ひずみ速度に基づいて第２のエネルギーロスを生成させる。即ち、第２の変形ネットワークＢでは、ひずみ速度に依存した抵抗がダッシュポットｐに生じ、これがバネｅの伸長比λ_c に影響する。伸長比λ_c は、式（１）のただし書きのように、分子鎖ｃの各方向の主ストレッチから計算できる。

次に、第２の変形ネットワークＢにおいて、各主ストレッチλ_i ^(Be)は、次の手順により計算される。なおλの添字ｉは各方向を表し（ｉ＝１、２…）、二次元であればｘ，ｙを意味する。また右肩の添字（Be）は第１の変形ネットワークＡの主ストレッチと区別するための符号である。

先ず第２の変形ネットワークＢにおける相当せん断応力τ^(B) が、式（１０）で計算される。式（１０）のσ^(B)' は、前ステップの計算で得られた第２の変形ネットワークＢの相当応力である。また、ダッシュポットｐのひずみ速度が、相当せん断応力τ^(B)及びλ_c ^(Bp)の値を用いて式（１１）から計算される。なお式（１１）のλ_c ^(Bp)は、ダッシュポットの伸長比であり、この値は、前回の計算値λ_i ^(Bp)（後述の式（１４）で計算される。）から式（１）のただし書きの方法に準じて計算される。

また前記ひずみ速度から、ダッシュポットｐの変形速度Ｄ^(Bp)が、式（１２）より計算され、該変形速度Ｄ^(Bp)からダッシュポットｐの伸びが式（１３）及び（１４）で計算される。ここで、添字ｔは計算の時間ステップを示し、Δｔはその増分である。

そして、第２の変形ネットワークＢのバネｅの各ストレッチλ_i ^(Be)は、下記式（１５）から得られる。なお式（１５）中、λ_i は、第２の変形ネットワークＢの全体の伸長比である。

前記λ_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 ＝ｎ・Ｎ＝４．３２８×１０²⁶
フィラーモデルの体積含有率３０％
フィラーモデルの縦弾性係数Ｅ：１００ＭＰａ
フィラーモデルのポアソン比ν：０．３

さらに変形シミュレーションは、材料全体Ｍの解析領域に一軸引張変形を発生させた。引張条件は、図１２のｘ₂方向にひずみ速度＝０．０１（／ｓ）とし、４サイクルの変形が順番に行われた。各サイクルの引張り量は次の通りとした。
１サイクル目：１００％（伸長比２．０）
２サイクル目：１００％（伸長比２．０）
３サイクル目：２００％（伸長比３．０）
４サイクル目：２００％（伸長比３．０）

また変形シミュレーションでは、各ひずみ状態において、ゴムモデル４の平均セグメント数Ｎが式（９）により計算され、その値を式（１）へ逐次代入した。また材料モデルは、厚さ方向（図３のＺ軸方向）に変化しないように拘束条件が与えられた。さらに、界面モデル５及びマトリックスモデル６の平均セグメント数Ｎは、次のように設定した。

＜マトリックスモデル＞
・負荷変形時の平均セグメント数
Ｎ＝-3.2368+20.6175λ_c-21.8168λ_c ²+10.8227λ_c ³ -1.9003λ_c ⁴
・除荷変形時の平均セグメント数
Ｎ＝６．６
・分子鎖のセグメントの総数（一定）
Ｎ_A ＝４．３２８１×１０²⁶

＜界面モデル＞
・負荷変形時の平均セグメント数
Ｎ＝-5.9286+20.6175λ_c-21.8168λ_c ²+10.8227λ_c ³ -1.9003λ_c ⁴
・除荷変形時の平均セグメント数
Ｎ＝３．９１
・分子鎖のセグメントの総数（一定）
Ｎ_A ＝３．２０３×１０²⁵

前記変形計算が行われると、その結果から、各計算ステップのタイミング毎に必要な物理量を得る（ステップＳ４）。図１３（Ａ）、（Ｂ）には、粘弾性材料モデルの１ないし２サイクル目の応力−伸長比の関係を示すグラフが、また、図１４（Ａ）、（Ｂ）には、粘弾性材料モデルの３ないし４サイクル目の応力−伸長比の関係を示すグラフが示される。

図から明らかなように、粘弾性材料モデル２は、２サイクル目以降、マリンス効果が再現されていることが確認できる。

即ち、１サイクル目の応力−伸長比のグラフが描くループは、第１のエネルギーロスＥ１及び第２のエネルギーロスＥ２をともに含んでいる。即ち、負荷変形では、概念的には、ひずみ量に基づいて分子鎖ｃの結合点ｂの破損が生じるので、１本の分子鎖当たりの平均セグメント数が増加し、第１のエネルギーロスが生成される。また、変形ネットワークＢでは、ひずみ速度に基づいたダッシュポットの抵抗（エネルギーロス）が生じるので、これに基づいた第２のエネルギーロスが生成される。

また、図１３において、１サイクル目と同じひずみ量の２サイクル目では、第１のエネルギーロスが生成されず、第２のエネルギーロスだけが生成されていることがわかる。これは、２サイクル目のひずみ量が１サイクル目の最大ひずみを超えていないので、分子鎖ｃの結合点ｂが新たに破損しないことに基づいている、従って、この実施形態では、２サイクル目の負荷変形は、１サイクル目の除荷変形と同じ経路を通ることになる。他方、２サイクル目においても、そのひずみ速度に基づいた第２のエネルギーロスは生じるので、図１３（Ｂ）のように、応力−伸長比曲線は、１サイクル目の除荷変形と同じ経路を通ることになる。

また、３サイクル目の曲線は、負荷変形時、伸長比が１ないし２サイクル目と等しい２．０までの間は、２サイクル目と同じ経路を通る、しかし、伸長比が２．０を超えると、新たな分子鎖ｃの結合点ｂが破損するため、応力が上昇し、ひずみ量に基づいた第１のエネルギーロスが再び生成されることになる。また、１ないし２サイクル目と同様、ひずみ速度に応じて第２のエネルギーロスも生じ得る。

また、４サイクル目では、２サイクル目と同様、ひずみが３サイクル目の最大ひずみを超えないので、分子鎖ｃの結合点ｂが新たに破損して消滅することはない。このため、４サイクル目の負荷変形は、３サイクル目の除荷変形と同じ経路を通る。このため、ひずみ量に基づいた第１のエネルギーロスは生成されず、ひずみ速度に応じた前記第２のエネルギーロスのみが生じることになる。

以上のように、本実施形態のシミュレーション方法によれば、粘弾性材料、とりわけ加硫ゴムに見られるマリンス効果をシミュレーション上に正確に再現することができる。従って、繰り返し荷重を受けて周期的に変形するような空気入りタイヤのゴム材料の性能評価をより精度良く行うことができる。

本実施形態で用いたコンピュータ装置の一例を示す斜視図である。本実施形態の処理手順を示すフローチャートである。粘弾性材料モデルの微視構造の一実施形態を示す線図である。カーボンブラックの形状を示す線図である。ゴムモデルの変形原理を説明する概念図である。（Ａ）は粘弾性材料、（Ｂ）はその分子鎖１構造を説明する線図、（Ｃ）は１本の分子鎖の拡大図、（Ｄ）はセグメントの拡大図である。粘弾性材料の網目構造体及びそれを構成する８鎖モデルの一例を示す斜視図である。（Ａ）、（Ｂ）は分子鎖の接合点の破断を説明する線図である。パラメータλc と分子鎖１本当たりの平均セグメント数Ｎとの関係を示すグラフである。任意の変形速度における界面モデル及びマトリックスモデルの応力−ひずみ曲線である。変形シミュレーションの手順を示すフローチャートである。均質化法を説明する全体構造とそのユニットセルとの関係を示す。（Ａ）は変形シミュレーションの１サイクル目、（Ｂ）はその２サイクル目を示し、それぞれ応力−伸長比の関係を示すグラフである。（Ａ）は変形シミュレーションの１サイクル目、（Ｂ）はその２サイクル目を示し、それぞれ応力−伸長比の関係を示すグラフである。ゴム材料のマリンス効果を説明する応力ひずみ線図である。

符号の説明

１コンピュータ装置
２材料モデル
３フィラーモデル
４ゴムモデル
５界面モデル
６マトリックスモデル

Claims

粘弾性材料の変形をシミュレートする方法であって、
前記粘弾性材料を数値解析が可能な要素で分割した粘弾性材料モデルを設定するステップと、
前記粘弾性材料モデルに条件を設定して変形計算を行うステップと、
前記変形計算から必要な物理量を取得するステップとを含み、
前記粘弾性材料モデルには、ひずみを増加させる負荷変形及び該ひずみを減少させる除荷変形からなる変形１サイクルにおいて、前記ひずみ量に基づいた第１のエネルギーロスを生成させる第１のエネルギーロス特性と、ひずみ速度に基づいた第２のエネルギーロスを生成させる第２のエネルギーロス特性とが定義されるとともに、
少なくとも２サイクルの変形計算を行う場合、１サイクル目において、前記第１及び第２のエネルギーロスを生成させる処理と、
２サイクル目において、ひずみが１サイクル目の最大ひずみを超えない場合、前記第２のエネルギーロスのみを生成させる処理とを含むことを特徴とする粘弾性材料のシミュレーション方法。
前記粘弾性材料モデルには、下記式（１）で示される応力とひずみとの関係が定義されるとともに、
前記第１のエネルギーロスは、下記式（１）における１本の分子鎖当たりの平均セグメント数Ｎを、負荷変形時ではひずみ量に基づいて増加させ、かつ、除荷変形時では除荷を開始したときの値を維持させることにより生成されることを特徴とする請求項１記載の粘弾性材料のシミュレーション方法。
前記粘弾性材料モデルは、ゴムをモデル化したゴムモデルと、前記ゴムに囲まれるフィラーを前記ゴムモデルで囲まれるようにモデル化したフィラーモデルとを含み、
かつ前記ゴムモデルは、前記フィラーモデルの周囲を囲む厚さが小さい界面モデルと、この界面モデルを取り囲むマトリックスモデルとからなり、しかもその少なくとも一方に、前記第１及び第２のエネルギーロス特性が定義されることを特徴とする請求項１又は２記載の粘弾性材料のシミュレーション方法。