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

Abstract

【課題】 本発明は、ゴム材料の変形状態等を精度良く解析するのに役立つゴム材料のシミュレーション方法を提供する。
【解決手段】 ゴムとフィラーとを含むフィラー充填ゴムの変形をシミュレーションするゴム材料のシミュレーション方法であって、前記フィラー充填ゴムを数値解析が可能な要素で分割した材料モデルを設定するステップＳ１と、前記材料モデルに条件を設定して変形計算を行うステップＳ２、Ｓ３と、前記変形計算から必要な物理量を取得するステップＳ４とを含み、かつ前記材料モデルは、前記フィラーを有限個の要素で分割したフィラーモデルと、前記フィラーモデルを取り囲むとともに前記ゴムを有限個の要素で分割したゴムモデルとを含むとともに、前記ゴムモデルには、変形速度に応じて発生する応力が異なる変形速度依存性が定義されることを特徴とする。
【選択図】 図１

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)

非特許文献１は、分子鎖網目理論（後で詳しく述べる。）を用いることにより、ゴム材料を解析することが記載されている。しかしながら、この文献に記載された方法では、ゴム材料が一般に有する変形速度依存性を再現することができない。

変形速度依存性とは、ゴム材料が、その変形速度によって異なる粘弾性特性を示す性質である。簡単に述べると、同じゴム試験片に異なる周波数の振幅歪（例えば１０Hzと１０００Hz）を与えた場合、それぞれの周波数でのエネルギーロスは異なった値を示す。上記非特許文献１では、このようなゴム材料の変形速度依存性が考慮に入れられていないため、一つのモデルから、互いに異なる周波数での変形性能を正しく評価することができない。

ここで、ゴム材料が主体的に用いられている空気入りタイヤを例に挙げ、その燃費性能とグリップ性能（加減速時等でのタイヤと路面との密着性の指標）という２つの性能について考える。タイヤの燃費性能を向上させるためには、一般的な定常走行速度においてエネルギーロスができるだけ小さいゴムをトレッドに用いることが有効である。つまり、燃費性能の評価については、周波数で約１０〜１００Hz程度の低い変形速度におけるゴムの粘弾性（エネルギーロス）を調べる必要がある。

一方、グリップ性能の向上には、トレッドゴムと路面との接地瞬間において、該路面の微小な凹凸にトレッドゴムを柔軟に変形させて接触させる必要がある。このため、グリップ性能は、高速変形時において、エネルギーロスが大きいほど良い。つまり、グリップ性能の評価については、周波数で約１００００〜１００００００Hz程度の大きな変形速度でのゴムの粘弾性（エネルギーロス）を調べる必要がある。

図２１には、タイヤ用ゴム組成物について、周波数・温度換算則を用いて計算された周波数−エネルギーロスのグラフを実線にて示す。図から明らかなように、周波数が変化すると、エネルギーロスも大きく変わる。然るに、前記非特許文献１では、ゴムモデルの変形速度を異ならせても、二点鎖線で示されるような同じ値しか得られない。これでは、ゴム開発に有用な情報は得られない。なお変形速度は、変形周波数とひずみとの積とする。

また、前記非特許文献２は、フィラーが充填されたゴム材料からゴムモデルを設定し、そこに変形速度依存性を考慮することが記載されている。しかし、フィラーを配合したゴム材料の変形計算を行う場合、フィラーとマトリックスゴムとの界面をどのように取り扱うかは重要な問題である。種々の研究の結果、前記界面には、様々な現象が生じており、とりわけ界面におけるゴムとフィラーとの滑りないし摩擦現象は、比較的大きなエネルギーロスを生じさせることが昨今分かってきた。従って、精度の良いゴム材料のシミュレーションを行うためには、フィラーとゴムとを分けてモデル化することが重要である。

しかしながら、前記非特許文献２に記載されている技術は、フィラーとゴムとを区別することなく一つにまとめてモデル化している。このような方法では、ゴムのフィラーとの界面の様子や、変形時の応力の分布状態などの情報を知ることができない。

本発明は、以上のような問題点に鑑み案出なされたもので、材料モデルを、フィラーを有限個の要素で分割したフィラーモデルと、このフィラーモデルを囲むとともにゴムを有限個の要素で分割したゴムモデルとを含んで設定し、かつ、前記ゴムモデルには、変形速度に応じて発生する応力が異なる変形速度依存性を定義することを基本として、フィラーとゴムとの界面を含めて、精度良くゴム材料の変形状態等をシミュレーションしうるゴム材料のシミュレーション方法を提供することを目的としている。

本発明のうち請求項１記載の発明は、ゴムとフィラーとを含むフィラー充填ゴムの変形をシミュレーションするゴム材料のシミュレーション方法であって、前記フィラー充填ゴムを数値解析が可能な要素で分割した材料モデルを設定するステップと、前記材料モデルに条件を設定して変形計算を行うステップと、前記変形計算から必要な物理量を取得するステップとを含み、かつ前記材料モデルは、前記フィラーを有限個の要素で分割したフィラーモデルと、前記フィラーモデルを囲むとともに前記ゴムを有限個の要素で分割したゴムモデルとを含み、前記ゴムモデルには、変形速度に応じて発生する応力が異なる変形速度依存性が定義されたことを特徴とする。

また請求項２記載の発明は、前記物理量は、前記ゴムモデルの各要素に生じた応力又はエネルギーロスを含み、かつ、前記物理量を視覚化しその分布図を得るステップをさらに含むことを特徴とする請求項１又は２記載のゴム材料のシミュレーション方法である。

また請求項３記載の発明は、前記ゴムモデルは、前記フィラーモデルの周囲を囲む厚さが小さい界面モデルと、この界面モデルを取り囲むマトリックスモデルとを含むとともに、その少なくとも一方に、前記変形速度依存性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法である。

また請求項４記載の発明は、前記界面モデルは、任意の変形速度において、前記マトリックスモデルよりも軟い粘弾性特性が定義されることを特徴とする請求項３記載のゴム材料のシミュレーション方法である。

本発明のシミュレーション方法では、変形計算に供される材料モデルが、フィラーモデル及びそれを囲むゴムモデルを含む。このため、ゴムモデルについては、フィラー周辺での物理量のみならず、フィラーから離れた位置での物理量を得ることができる。言い換えると、フィラーの周りの物理量（例えば応力やエネルギーロス等）の分布を知り得る。従って、より詳細かつ具体的な解析が可能となる。これらの物理量は、請求項２に記載の発明のように、分布図として視覚化されることによって、より一層、ゴム及びフィラー材料の開発において有用な情報となる。

また本発明のシミュレーション方法では、ゴムモデルに変形速度依存性が定義される。ゴムに周期的な歪を与えて粘弾性を測定する物理的な試験器は存在する。しかし、現在の測定技術では歪の周波数はせいぜい１０００Hz程度である。これは、例えば、タイヤトレッドゴムのグリップ性能を評価するためには十分ではない。本発明のようにゴムモデルに変形速度依存性を与えることによって、試験器では得ることができない次元での変形速度シミュレーションができる。

また請求項３に記載された発明のように、ゴムモデルは、界面モデルとマトリックスモデルとを含むことが望ましい。フィラー充填ゴムは、フィラーとゴムとの間に物理的な界面がある。フィラーがカーボンブラックである場合、ゴムとの界面は物理的な結合である。またフィラーがシリカである場合、ゴムとの界面はカップリング剤による化学的結合である。これらの結合は、界面モデルのゴムの粘弾性特性の調節により、擬似的かつ比較的正確に再現できる。したがって、カップリング剤の検討及び評価を行うのに役立つ。

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

図２には、本実施形態のシミュレーション方法の処理手順の一例が示される。本実施形態では、先ずゴムとフィラーとを含むフィラー充填ゴムの材料モデルが設定される（ステップＳ１）。図３には、微視構造（「ユニットセル」とも呼ばれる）Ｕとしての材料モデル２の一例が視覚化して示されている。

該材料モデル２は、解析しようとするフィラー充填ゴムの微小領域を、有限個の小さな要素２ａ、２ｂ、２ｃ…に分割（置換）することにより得られる。解析対象のフィラー充填ゴムは、実在するか否かは問われない。各要素２ａ、２ｂ、２ｃ…は、数値解析が可能に定義される。数値解析が可能とは、例えば有限要素法、有限体積法、差分法又は境界要素法といった数値解析法により各要素ないし系全体についての変形計算が可能なことを意味する。

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

この実施形態の材料モデル２は、後述する変形シミュレーションにおいて平面ひずみ状態で解析が行われる。したがってＺ方向（紙面と直角方向）にはひずみを持たない。また、この例の材料モデル２は、例えば縦横それぞれ３００×３００ナノメータで設定されてる。

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

本実施形態のフィラーモデル３は、カーボンブラックをモデル化することにより設定されている。但しフィラーは、上記のものに限定されるものではなく、例えばシリカ等の他のフィラーであっても良い。フィラーモデル３の形状は、例えば、電子顕微鏡にて撮像した実際のゴム中のフィラー（カーボンブラック）の形状に基づいて定め得る。図４に示されるように、カーボンブラックＣＢは、具体的には、炭素原子からなる直径数１０ナノメータ程度の球状の一次粒子７が不規則に３次元的に結合した葡萄の房状構造を持つ。本実施形態のフィラーモデル３は、この構造の平面視を基準に設定されている。

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

次に、前記界面モデル５について述べる。ゴムとフィラーとは界面で物理的ないし化学的に結合しており、他の部分では見られない様々な現象、例えばゴムとフィラーとの間の大きな滑りやこれに伴う摩擦が生じていると推察されている。

このため、本実施形態では、フィラーモデル３を連続して取り囲み、かつ、薄い厚さで界面モデル５を設定する。界面モデル５は必ずしも連続してフィラーモデル３を取り囲む必要はない。また界面モデル５の厚さｔも特に限定はされないが、種々の実験の結果、前記厚さｔは、例えば１〜２０ナノメートル、より好ましくは５〜１０ナノメートルが望ましい。

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

本発明では、ゴムモデル４に、変形速度に応じて発生する応力が異なる変形速度依存性が定義される。図５には、変形速度依存性を有するゴムモデル４の変形の原理を説明する概念図が示される。ゴムモデル４の個々の要素は、第１の変形ネットワークＡと、第２の変形ネットワークＢとが並列に結合されたものと考えることができる。系全体に発生する応力Ｓは、第１の変形ネットワークＡに生じた応力と、第２の変形ネットワークＢに生じた応力との和である。また、両ネットワークＡ、Ｂは並列結合であるため、各々のネットワークＡないしＢに生じる歪（伸び）はともに同一である。

第１の変形ネットワークＡは、バネと等価なものとして考える。第１の変形ネットワークＡに生じる応力は、バネの伸び（歪）に応じた値になる。またバネの伸びは、速度に関係なく単に荷重に基づいて生じる。従って、第１の変形ネットワークＡは、変形速度に依存しない応力を表す。

第２の変形ネットワークＢは、バネｅとダッシュポットｐとの直列接続体と等価なものとして考え得る。ダッシュポットｐは、ニュートンの粘性法則に従う流体（油等）が封止されたシリンダ筒の中を、オリフィス等が設けられたピストンが移動する減衰装置であって、ピストンの移動速度に比例した抵抗が生じるものとする。第２の変形ネットワークＢの全体の伸び（歪）は、バネｅの伸びとダッシュポットの伸びとの和である。さらに第２の変形ネットワークＢに生じる応力も、バネｅの伸び（歪）に応じた値になる。しかし、バネｅの伸びは、荷重によって単純に決定されるものではなく、ゴムモデルの変形速度に依存したダッシュポットｐの抵抗の影響を受ける。つまり、第２の変形ネットワークＢは、変形速度に依存した応力を表す。

次に、各変形ネットワークＡ、Ｂに基づいたゴムモデル４の応力計算手順について述べる。

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

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

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

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

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

本実施形態では、第１の変形ネットワークＡにおいて、非特許文献１などが提案する一般的な８鎖モデルを以下のように修正して用いる。

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

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

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

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

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

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

式（８）から明らかなように、負荷変形によって、ゴムの分子鎖ｃのからみ数ｍが減少していくと、１本の分子鎖ｃに含まれる平均セグメント数Ｎは増加する。従って、負荷変形に応じて、式（１）の平均セグメント数Ｎを増加させることにより、分子鎖ｃのからみの減少を式（１）の中に取り込み得る。これは、ゴムの変形メカニズムをより精度良く擬似化でき、計算精度の向上に役立つ。ここで、負荷変形とは、微小時間の間でマトリックスモデル４の歪が増大する変形である。逆に、歪が減少する変形を除荷変形とする。

前記平均セグメント数Ｎは、様々な方法で増加させ得る。例えば、平均セグメント数Ｎは、負荷変形時における歪に関連したパラメータに基づいて増大するものが望ましい。歪に関連したパラメータとしては、特に限定されるものではないが、例えば、歪、歪速度又は歪の１次の不変量Ｉ1 などがある。本実施形態では、前記平均セグメント数Ｎを、下記式（９）で定める。式（９）は、平均セグメント数Ｎが、当該マトリックスモデル４の各要素個々において、それぞれ歪の１次の不変量Ｉ₁ （より詳しくはその平方根であるパラメータλc ）の関数であることを示す。

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

図９には、前記平均セグメント数Ｎとλc との関係が示される。歪に関するパラメータであるλc が大きくなると、平均セグメント数Ｎは滑らかに増大する。この例では、λc の上限は約２．５である。後述する変形シミュレーションにおいては、ゴムモデル４の各要素について、負荷変形時ではλc が常に計算される。計算されたλc は、式（９）に代入され、当該要素の当該歪状態における平均セグメント数Ｎが計算される。平均セグメント数Ｎは、適宜のタイミングで計算されて、式（１）に取り込まれる。なお本実施形態では、ゴムモデル４の除荷変形においては、その除荷開始時の平均セグメント数Ｎを、ひずみが零になるまで保持とする。

＜第２の変形ネットワークＢ＞
分子鎖網目理論に基づく８鎖モデルは、第２の変形ネットワークＢのバネｅにも適用される。従って、第２の変形ネットワークＢに生じる応力も、前記式（１）を用いて計算しうる。ただし、第２の変形ネットワークＢでは、変形速度に依存した抵抗が生じ、これは伸長比λc に影響する。

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

先ず第２の変形ネットワークＢにおける相当せん断応力τ^(B) を式（１０）で計算する。式（１０）のσ^{(B) '} は、前ステップの計算で得られた第２の変形ネットワークＢの相当応力である。また、相当せん断応力τ^(B) 及びλ_c ^(Bp)の値を用いて式（１１）からダッシュポットｐのひずみ速度を計算する。なお式（１１）のλ_c ^(Bp)は、ダッシュポットの伸長比であり、この値は、前回の計算値λ_i ^(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 ＝ｎ・Ｎ＝４．３２８×１０²⁶
フィラーモデルの体積含有率 ３０％
フィラーモデルの縦弾性係数Ｅ：１００ＭＰａ
フィラーモデルのポアソン比ν：０．３

さらに変形シミュレーションは、材料全体Ｍの解析領域に一様な一軸引張変形を発生させた。引張条件は、図１２のｘ２方向にひずみ速度＝１（／ｓ）の高速変形と、ひずみ速度＝０．０１（／ｓ）の低速変形との２種の条件を用い、いずれもひずみを０．１５とした後、それぞれの速度で除荷変形させた。

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

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

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

前記変形計算が行われると、その結果から、各計算ステップのタイミング毎に必要な物理量を得ることができる（ステップＳ４）。図１３は、その計算結果を利用して作られた材料モデルの真応力とひずみとの関係を示すグラフである。図から明らかなように、変形速度によって、材料モデル２は異なる粘弾性特性を示す。本シミュレーションの結果、実際のゴム材料と同様、高速変形では、低速変形（一点鎖線）に比して材料の硬化現象が再現されている。また高速変形の方が、低速変形よりもエネルギーロスが大きくなっていることが確認された。

図１４には、本発明の他の実施形態を示す。この実施形態では、変形速度依存性が、界面モデルのみに定義され、マトリックスモデルには定義されていない。フィラーの代表例として、カーボンブラックとシリカとが挙げられる。両者の最も大きな違いは、フィラー自体よりも、むしろフィラー表面とゴムとの境界層、すなわち界面にあると考えられており、この界面の付近のゴムの方が、マトリックスゴムよりも、フィラー充填ゴムの特性に対して敏感に影響を及ぼし、変形速度に対する依存性も大きいと考えられる。従って、界面モデルとマトリックスモデルとの変形速度依存性を変えて計算した結果は、実際にフィラーとポリマーとをどのように結合させれば良いかという重要な知見を得るのに役立つ。

各界面モデルには、図１４に示されるような変形速度依存性が定義される。この例では、カーボンブラックを例にした界面モデル１は変形速度依存性が大きくかつエネルギーロスも大きく設定されている。シリカを例にした界面モデル２は、変形速度依存性が小さくかつエネルギーロスも小さく設定されている。なお図中、高速変形は、変形速度１０００（／ｓ）であり、低速変形は０．１（／ｓ）であることを示す。

図１５には、前記界面モデル１ないし２を用いた２種の材料モデル１ないし２を使用し、変形速度１０００（／ｓ）で変形シミュレーションを行った結果を示す。縦軸はひずみ、横軸は真応力である。この結果より、厚さが小さい界面モデルの速度依存性を変えただけであるにも拘わらず、材料全体の粘弾性特性に大きな差異が生じることが分かる。図１６には、横軸に弾性率（図１５のグラフの横軸に対する傾き）、縦軸にエネルギーロス（ヒステリシスループ）をとったグラフである。前述のタイヤのグリップ性能を例に挙げると、グラフの左斜め上方に向かう粘弾性特性が好ましいものとなる。

また、本実施形態のシミュレーションでは、得られた物理量から分布図を作るステップを含む（ステップＳ５）。図１７には、高速変形時かつ最大歪状態のユニットセルＵ（なおユニットセルＵは、図１２に示すもので、図３のものとは異なる。）の各要素に生じた応力の分布図の一例を示す。併せて、ユニットセルＵの輪郭形状の変化も視覚化されている。この分布図を作るステップは、各要素それぞれについて得られた応力値を色情報（彩度、色相又は明度等）等、視覚にて識別しうる情報に置き換えることで行われる。図１７では、色彩が濃いほど大きなひずみであることを示している。この結果より、フィラーモデル３、３の間の距離が小さい部分に、より大きな応力が集中していることが分かる。また、応力が高い部分は、フィラーモデル３を介して引張方向に沿って線状につながっていることも分かる。

また図１８には、高速変形時のユニットセルＵの各要素に生じたエネルギーロスの分布図を示す。エネルギーロスを得るために、各変形速度条件と等しい速度で除荷を行い１サイクルの応力ひずみ曲線を得た。またユニットセルＵの輪郭形状は、図１７のものを採用した。この分布図においても、各要素それぞれについて得られたエネルギーロスを色情報（彩度、色相又は明度等）に置き換えて視覚化したものである。色彩が濃いほど大きなエネルギーロスが生じていることを示す。この結果からも、フィラーモデル３、３の間の距離が小さい部分に、より大きなエネルギーロスが生じていることが分かる。なおフィラーモデル３は着色しているがロスはない。また分布図に用いられる物理量は、応力値やエネルギーロスに限られるものではなく、解析の目的応じて種々のものが採用される。

また図１９（Ａ）には、高速変形時かつ最大歪状態のユニットセルＵの応力分布の部分拡大図、同図（Ｂ）には低速変形時のそれが示される。歪が同一であるにも拘わらず、変形周波数の違いによって、応力の分布は異なった様相を示す。即ち、高速変形時の場合、高い応力が比較的広い範囲で生じており、特にフィラーモデル３、３間ではそれが顕著に表れている。

前記非特許文献２のように、フィラーとゴムとをまとめて一つにモデル化した場合、ユニットセルの応力分布は均一であり、かつ物理的な変形も単に長方形になるに過ぎない。従って、これらの分布図は、フィラー充填ゴムのフィラーとゴムとをそれぞれ分けてモデル化することによって初めて得ることができる。

前記分布図は、実際のゴム材料の開発に役立つ非常に有用な情報を提供する。例えば、フィラー充填ゴムにおいて、変形時、より均一な応力分布が得られるものが好ましい。このようなフィラー充填ゴム材料を開発するためには、先ず現状の応力分布状態を知ることが重要である。そして、その結果から、フィラー間の距離のバラツキが極力なくなるように、フィラーの分散性やその物理的形状についてさらなる改善を試みることができる。そして、改良されたたフィラーモデル等を用いて、再度ユニットセルを作成し、より好ましい均一な応力分布が得られるまでコンピュータ上で繰り返し変形シミュレーションを行うことができる。

そして、計算によって、より均一な応力分布が得られた場合、そのユニットセルの物理的構成、即ちフィラーの粒子形状や分散度等を実際の製品に反映させれば良い。特にフィラーの分散状態は、ゴムへの配合量、混練時間及び混練温度の３つのパラメータによってコントロールしうるため、理想的なフィラーモデル３の分散性が得られれば、グリップ力に優れたタイヤ用のゴム組成物或いは燃費性能に優れたゴム材料を比較的短期に開発することができる。

また、本発明のシミュレーション方法は、空気入りタイヤのトレッドゴムの開発に有効に用いることができる。先に述べたように、トレッドゴムの燃費性能は、ゴムの低速変形時の複素弾性率が大きくかつエネルギーロスが小さいほど良い。ここで、トレッドゴムの燃費性能の指標Ｘ１は、下記式（１６）で表すことができ、この指標Ｘ１を大きくするゴムが好適である。
Ｘ１＝Ｋ１・ｔａｎδ_L ＋Ｋ２・Ｅ_L ^{* -1} …式（１６）
ただし、Ｋ１及びＫ２は定数であり、ｔａｎδ_L は、周波数１０〜１００Hzの範囲内の特定周波数でのフィラー充填ゴムの損失正接、Ｅ_L ^* は、変形周波数１０〜１００Hzの範囲内の特定周波数でのフィラー充填ゴムの複素弾性率とする。

他方、トレッドゴムのグリップ性能は、ゴムの高速変形時のエネルギーロスが大きいほど良い。ここで、トレッドゴムのグリップ性能の指標Ｘ２は、下記式（１７）で表すことができ、この指標Ｘ２を大きくするゴムが好適である。
Ｘ２＝Ｋ３・ｔａｎδ_H …式（１７）
ただし、Ｋ３は定数であり、ｔａｎδ_H は変形周波数１０⁴ 〜１０⁶ Hzの範囲内の特定周波数でのフィラー充填ゴムの損失正接とする。

図２０には、本発明のシミュレーション方法を利用したトレッドゴムの開発方法の一例を示すフローチャートが示される。本実施形態では、先ず、材料モデル２を設定し（ステップＳａ１）、それに基づいて低速変形条件で変形シミュレーションを行い前記燃費性能に関する指標Ｘ１が計算される（ステップＳａ２）。また同じ材料モデル２に基づいて高速変形条件で変形シミュレーションを行い前記グリップ性能に関する指標Ｘ２が計算される（ステップＳａ３）。なおステップＳａ２及びＳａ３は、順序を逆としても良い。

燃費性能とグリップ性能とを両立させるためには、前記各指標Ｘ１及びＸ２の和が大きいほど良いため、本実施形態では、指標Ｘ１及びＸ２の和が最大になったか否かを判定する（ステップＳａ４）。そして、ステップＳａ４で真の場合、設定された材料モデル２（ユニットセルＵ）に基づいて実際のフィラー、ゴム配合等を試作、設計等すれば良い。また、ステップＳａ４で偽の場合、ユニットセルＵのパラメータ等を変更しながら材料モデル２を再度設定し（ステップＳａ１）、前記指標の和（Ｘ１＋Ｘ２）が最大となるまで変形シミュレーションを繰り返せば良い。なお変更するパラメータとしては、例えば応力−歪線図の傾き（これは、式（１）のＣ^R やＮのパラメータを変更することでなしうる。）、ヒステリシスループの大きさ（これは例えば図９のλ_c 、Ｎの傾きを変えることによってなしうる。）、前記変形速度依存性（これは、式１１の材料定数を変えることによってなしうる。）等が挙げられる。このように、本発明のシミュレーション方法は、具体的な各種ゴム製品の開発に大いに役立つ。

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

符号の説明

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

Claims (4)

1. ゴムとフィラーとを含むフィラー充填ゴムの変形をシミュレーションするゴム材料のシミュレーション方法であって、
   前記フィラー充填ゴムを数値解析が可能な要素で分割した材料モデルを設定するステップと、
   前記材料モデルに条件を設定して変形計算を行うステップと、
   前記変形計算から必要な物理量を取得するステップとを含み、
   かつ前記材料モデルは、前記フィラーを有限個の要素で分割したフィラーモデルと、
   前記フィラーモデルを囲むとともに前記ゴムを有限個の要素で分割したゴムモデルとを含み、
   前記ゴムモデルには、変形速度に応じて発生する応力が異なる変形速度依存性が定義されたことを特徴とするゴム材料のシミュレーション方法。
2. 前記物理量は、前記ゴムモデルの各要素に生じた応力又はエネルギーロスを含み、かつ、前記物理量を視覚化しその分布図を得るステップをさらに含むことを特徴とする請求項１記載のゴム材料のシミュレーション方法。
3. 前記ゴムモデルは、前記フィラーモデルの周囲を囲む厚さが小さい界面モデルと、この界面モデルを取り囲むマトリックスモデルとを含むとともに、その少なくとも一方に、前記変形速度依存性が定義されることを特徴とする請求項１又は２記載のゴム材料のシミュレーション方法。
4. 前記界面モデルは、任意の変形速度において、前記マトリックスモデルよりも軟い粘弾性特性が定義されることを特徴とする請求項３記載のゴム材料のシミュレーション方法。