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

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)

特許文献１は、ゴム材料（粘弾性材料）がひずみ速度に応じて異なった縦弾性係数を示す点に着目している。具体的には、予め粘弾性材料の実使用状態を想定した測定条件で、当該粘弾性材料に生じるひずみ、ひずみ速度、応力などを測定し、縦弾性係数とひずみないしひずみ速度との対応関係を導出する工程が行われる。そして、解析対象である粘弾性材料の製品モデルに対して、所定のひずみ速度を与えるとともに、上記対応関係から縦弾性係数を適宜計算して変形シミュレーションを行っている。特許文献１の方法は、粘弾性材料の縦弾性係数が、ひずみ速度に応じて変化するという試みを含むため、ヒステリシスロスをシミュレーションの中に取り込み得る点で評価できる。

またArrudaらが提案した非特許文献１は、分子鎖網目理論を用いることにより、ゴム材料を高分子レベルにまでモデル化して計算することが記載されている。ここで、分子鎖網目理論について簡単に述べる。

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

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

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

Aruudaらは、分子鎖網目理論に基づいてさらに８鎖モデルを提案している。図４（Ａ）に示されるように、粘弾性材料は、巨視的には、微小な８鎖モデルｇが集合した立方体状の網目構造体ｈとして考えることができる。一つの８鎖モデルｇは、図４（Ｂ）に拡大して示すように、分子鎖ｃが立方体の中心に定められた一つの接合点ｂ１から、各頂点に設けられた８つの各接合点ｂ２にそれぞれのびているものと仮定される。

ゴム材料は、シミュレーションにおいては超弾性体（体積変化を殆ど生じず除荷後も元の形状に戻る材料）として取り扱われる。超弾性体は、下記式（１）で示されるように、Green ひずみの成分Ｅ_ijによって微分されることにより、共役なkirchhoff 応力Ｓ_ijを生じるようなひずみエネルギー関数Ｗが存在する物質として定義される。換言すれば、ひずみエネルギー関数は、ゴム材料が変形したときに蓄えられたポテンシャルエネルギの存在を仮定的に示す。従って、ひずみエネルギー関数Ｗの微分勾配から応力とひずみ関係を得ることができる。

Aruudaらは、非ガウス鎖理論により、ゴム材料の変形が大きくなるとエントロピー変化が急激に大きくなる（分子鎖が伸びきって配向する）と考え、式（２）のゴム弾性体のひずみエネルギ関数Ｗを示した。そして、このひずみエネルギー関数Ｗを上記式に代入することにより、ゴム材料の応力とひずみとの関係を取り出すことができる。

Aruudaらの応力とひずみとの関係を用いて、粘弾性材料の変形シミュレーションを行うことにより、例えば図１７に示すように、１軸引張変形シミュレーションにおいて、非線形な応力とひずみとの関係を得ることができる。この結果は、荷重の負荷変形時における実測値と良い相関を示す。

工業製品として使用されているゴム材料には、通常、カーボンブラック等のフィラー（充填剤）が配合されている。このようなフィラーを配合したゴム材料の変形シミュレーションを行う場合、フィラーとマトリックスゴムとの界面をどのように取り扱うかは重要な問題である。種々の研究の結果、前記界面には、様々な現象が生じており、とりわけ界面におけるゴムマトリックスとフィラーとの滑りないし摩擦現象は、比較的大きなエネルギーロスを生じさせることが昨今分かってきた。従って、精度の良いゴム材料のシミュレーションを行うためには、このような界面をモデルとして取り込むことが重要である。

しかるに、特許文献１では、ゴムマトリックスとフィラーとの界面については考慮されていない。非特許文献２についても同様である。

本発明は、以上のような問題点に鑑み案出なされたもので、ゴム材料モデルを、ゴムマトリックスをモデル化したマトリックスモデルと、フィラーをモデル化したフィラーモデルと、前記マトリックスモデルと前記フィラーモデルとの間の界面を形成する界面モデルとを含んで設定することを基本として、界面に独自の特性を与えて精度良くゴム材料を変形シミュレーションしうるゴム材料のシミュレーション方法を提供することを目的としている。

本発明のうち請求項１記載の発明は、フィラーが配合されたゴム材料の変形をシミュレーションするゴム材料のシミュレーション方法であって、前記ゴム材料を数値解析が可能な要素でモデル化したゴム材料モデルを設定するステップと、前記ゴム材料モデルに条件を設定して変形計算を行うステップと、前記変形計算から必要な物理量を取得するステップとを含むとともに、前記ゴム材料モデルは、ゴムマトリックスをモデル化したマトリックスモデルと、フィラーをモデル化したフィラーモデルと、前記マトリックスモデルと前記フィラーモデルとの間の界面を形成する界面モデルとを含み、前記界面モデルは、前記フィラーモデルを連続して取り囲みかつ厚さを有するとともに前記マトリックスモデルとは異なる粘弾性特性が定義されることを特徴とするゴム材料のシミュレーション方法である。

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

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

また請求項４記載の発明は、前記界面モデルは、前記マトリックスモデルよりもヒステリシスロスが大きい粘弾性特性が定義されることを特徴とする請求項１乃至３のいずれかに記載のゴム材料のシミュレーション方法である。

また請求項５記載の発明は、前記界面モデルの厚さは１〜２０nmである請求項１記載のゴム材料のシミュレーション方法であり、請求項６記載の発明は、前記フィラーモデルは、ゴムよりも硬い弾性率が定義される請求項１記載のゴム材料のシミュレーション方法である。

本発明のゴム材料のシミュレーション方法では、解析の対象となるゴム材料モデルが、ゴムマトリックスをモデル化したマトリックスモデルと、フィラーをモデル化したフィラーモデルと、前記マトリックスモデルと前記フィラーモデルとの間の界面を形成する界面モデルとを含む。従って、界面の影響を計算の結果に取り込むことが可能となり、実際の粘弾性材料の特性に即した精度の高いシミュレーションが可能になる。

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

図２には、シミュレーション方法の処理手順の一例が示される。本実施形態では、先ずフィラーが配合されたゴム材料モデルが設定される（ステップＳ１）。図３には、微視構造としてのゴム材料モデル２の一例が視覚化して示されている。

該ゴム材料モデル２は、解析しようとするゴム材料（実在するか否かを問わない）の微小領域が、有限個の小さな要素２ａ、２ｂ、２ｃ…に置き換えられたものである。各要素２ａ、２ｂ、２ｃ…は、数値解析が可能に定義される。数値解析が可能とは、例えば有限要素法、有限体積法、差分法又は境界要素法といった数値解析法により、各要素ないし系全体についての変形計算が可能なことを意味する。具体的には、各要素２ａ、２ｂ、２ｃ…について、座標系における節点座標値、要素形状、材料特性などが定義される。各要素２ａ、２ｂ、２ｃ…には、例えば２次元平面としての三角形ないし四辺形の要素、３次元要素としては、例えば４ないし６面体の要素が好ましく用いられる。これにより、ゴム材料モデル２は、前記コンピュータ装置１にて取り扱い可能な数値データを構成しうる。

この実施形態のゴム材料モデル２は、後述する変形シミュレーションにおいて平面ひずみ状態で解析が行われる。したがってＺ方向にはひずみを持たない。ゴム材料モデル２の一辺の長さは、例えば縦横それぞれ３００nm×３００nmとしている。

この実施形態のゴム材料モデル２は、ゴムマトリックス部分がモデル化されたマトリックスモデル３（最も濃い部分）と、このマトリックスモデル３の中に分散して配されかつフィラー（充填材）がモデル化されたフィラーモデル４（白い部分）と、前記マトリックスモデル３と前記フィラーモデル４との間に介在しかつ両者の間の界面を形成する界面モデル５（やや黒い部分）とからなるものが例示される。

前記マトリックスモデル３は、ゴム材料モデル２の主要部を構成し、かつ、例えば三角形ないし四辺形の要素を用いて表現されている。マトリックスモデル３を構成する各要素には、材料特性として、下記式（３）の応力とひずみとの関係が定義される。したがって、マトリックスモデル３は各要素について、それぞれ応力とひずみとの関係が得られる。

ここで、式（３）の導出過程について触れる。ゴム材料は、変形による体積変化が小さいため、それを無視することができる。このため、非圧縮性の粘弾性材料では、非圧縮条件により密度を一定とすることができ、kirchhoff 応力Ｓ_ijは、式（４）で表すことができる。なお、Ｅ_ijはGreen ひずみの成分、ｐは静水圧（境界条件により定まる）、Ｘ_i は、応力もひずみも０である状態Ｃ0 における任意の物体点Ｐの位置、ｘ_j は変形した状態Ｃにおける物体点Ｐの位置である。

Cauchy応力の成分σ_ijと、kirchhoff 応力の成分Ｓ_mnとは下記式（５）の関係を有する。従って、式（４）から式（６）を得ることができる。

ここで、体積が一定の場合、物体の体積変化率を示す式（５）のＪは１とみなせる。また、式（２）のひずみエネルギー関数は、左Cauchy-Green変形テンソルＡ_ijのひずみの１次の不変量Ｉ₁ だけの関数であるため、式（６）から式（７）が得られる。

また、下記の関係式（８ないし１０）を用いることにより、式（７）の速度形式表示は式（１１）になる。

そして、体積一定条件下で、式（１１）に示すCauchy応力のJaumann 速度をkirchhoff 応力のJaumann 速度に置き換えても本質的な差異はないことと、変形速度テンソルＤをひずみ速度テンソルに置き換えることとにより、非圧縮性の粘弾性材料の構成式の速度形式表示を前記式（３）で表すことができる。

ところが、単にこの式（３）だけを用いた場合では、Aruudaらの非特許文献１と同様、ゴム材料の特徴の一つでもあるエネルギーロスを考慮することができない。即ち図１７の状態Ｑから荷重を除荷するシミュレーションを行うと、該曲線と実質的に同じ経路を通ってひずみの回復が行われる。これでは、ヒステリシスループが形成されずエネルギーロスの計算を行うことができない。

発明者らは、Aruudaらのモデルを前提として種々の実験を試みた。先ずゴム材料は、上述のように複雑に絡み合った長い分子鎖ｃが伸びることによって数百％にも達し得る大ひずみが許容されると考えられる。そこで、発明者らは、荷重の負荷における変形過程において、粘弾性材料の分子鎖ｃの互いに絡み合った部分がほどけたり（接合点ｂの数の減少）、或いは荷重の除荷によってさらなる絡み合い（接合点ｂの数の増加）が生じ得るとの仮定を立てた。換言すれば、この仮定は、前記式（１）における１本の分子鎖当たりの平均セグメント数Ｎは、負荷ないし除荷変形時において可変のパラメータであることを意味している。発明者らは、この仮定を検証したところ、粘弾性材料のヒステリシスを表現可能であることが分かった。

このように、本実施形態のゴム材料モデル２は、前記式（３）において、１本の分子鎖当たりの平均セグメント数Ｎが負荷変形時と除荷変形時とで異なるパラメータとして定義されている。これにより、マトリックスモデル３の各要素は、エネルギーロスを計算上取り込むことができる。ここで、負荷変形時とは、微小時間の間でモデルのひずみが増大する変形であり、逆に除荷変形時とは、ひずみが減少する変形を意味する。

図４（Ａ）には、粘弾性材料の巨視的な３次元の網目構造体ｈを、また図４（Ｂ）には、この網目構造体ｈを構成する８鎖モデルｇの１つを拡大して示した。この例の網目構造体ｈは、幅方向、高さ方向及び奥行き方向にそれぞれ８鎖モデルｇがｋ個結合したものとする。

いま、網目構造体ｈに含まれる接合点ｂの総数を「からみ数」として符号ｍで表すと、からみ数ｍは、式（１２）で表すことができる。同様に、網目構造体ｈに含まれる分子鎖ｃの数（即ち、マトリックスモデル３の単位体積中に含まれる分子鎖の数）をｎとすると、このｎは、式（１３）で表すことができる。
ｍ＝（ｋ＋１）³ ＋ｋ³ …式（１２）
ｎ＝８ｋ³ …式（１３）

ここで、ｋは十分に大きい数とすると、ｋの３次項以外を省略して上式はそれぞれ次のような式（１４）、（１５）で表すことができ、さらにこれらの関係から、からみ数ｍは、ｎを用いて、式（１６）で表すことができる。
ｍ＝２ｋ³ …式（１４）
ｎ＝８ｋ³ …式（１５）
ｍ＝ｎ／４ …式（１６）

さらに、マトリックスモデル３は、変形してもその中に含まれる分子鎖のセグメントの総数Ｎ_A は変化しないと仮定すると、式（１７）〜（１８）が成り立つ。
Ｎ_A ＝ｎ・Ｎ …式（１７）
Ｎ＝Ｎ_A ／ｎ＝Ｎ_A ／４ｍ …式（１８）

発明者らは、シミュレーションにおいてマトリックスモデル３のエネルギーロスを表現するために、荷重の負荷過程及びこれに続くひずみの回復過程において、粘弾性材料の分子鎖ｃのからみ数ｍが変化するとの着想を試みた。即ち、図５（Ａ）に示すように、例えば一つの接合点ｂで接合されている分子鎖ｃ１ないしｃ４に矢印方向の引張応力が作用すると、各分子鎖ｃ１ないしｃ４は伸び、接合点ｂは大きなひずみを受けて破損するものと考えた。

この結果、図５（Ｂ）に示すように、これまで２本であった分子鎖ｃ１及びｃ２は、あたかも１本の長い分子鎖ｃ５のように振る舞うものと考えられる。分子鎖ｃ３及びｃ４についても同様である。そして、このような現象は、ゴム材料の負荷変形が進むにつれて逐次発生していくものと考えられる。

以上のような接合点ｂの破損は、マトリックスモデル３におけるからみ数ｍの減少に他ならない。マトリックスモデル３自体は材料の出入りがないため、その中に含まれる分子鎖ｃのセグメントの総数Ｎ_A が変化しないと仮定すると、マトリックスモデル３の負荷変形が進むにつれて、１本の分子鎖ｃに含まれる平均セグメント数Ｎが増加することは式（１８）から明らかである。つまり、前記平均セグメント数Ｎは、一定ではなく、負荷変形時と除荷変形時とで異なる可変のパラメータとすることが適切と考えられる。これによって、マトリックスモデル３の変形時におけるエネルギーロスが再現できる。これについては、さらに詳細を後で述べる。

このような様子をコンピュータ上でより好ましく再現するためには、前記平均セグメント数Ｎは、負荷変形時では、そのひずみに関するパラメータに応じて増大させることが特に好適である。前記ひすみに関するパラメータとしては、特に制限されるものではないが、例えば、ひずみ量、ひずみ速度又はひずみの１次の不変量Ｉ₁ などが挙げられる。本実施形態では、前記平均セグメント数Ｎが、下記式（１９）に示すように、当該マトリックスモデル３の各要素個々において、それぞれひずみの１次の不変量Ｉ₁ を平方根であるパラメータλc の関数となるものを示す。

なお式（１９）は、種々の実験によって定めた一例であり、上記ＡないしＥは、いずれも定数であって、ゴム試験片の単純な１軸引張試験などの実測結果からように定めることができる。この例では、上記定数を次のように設定している。
Ａ＝＋２．９４９３
Ｂ＝−５．８０２９
Ｃ＝＋５．５２２０
Ｄ＝−１．３５８２
Ｅ＝＋０．１３２５

図６には、マトリックスモデル３の各要素の負荷変形時における平均セグメント数Ｎとパラメータλc との関係を示す。ひずみに関するパラメータであるλc が大きくなると、平均セグメント数Ｎは滑らかに増大する。この例では、パラメータλc の上限は２．５である。後述するゴム材料モデル２の変形シミュレーションでは、マトリックスモデル３の各要素について、負荷変形時においてはパラメータλc が常時計算される。計算されたλc は、式（１９）に代入される。これにより、当該要素の当該ひずみ状態における平均セグメント数Ｎが計算できる。

他方、マトリックスモデル３の除荷変形時には、平均セグメント数Ｎは一定値としている。平均セグメント数Ｎを決定する方法としては、例えば次の方法がある。解析対象となる粘弾性材料において、一つの応力−ひずみ曲線がある場合、先ず除荷時の曲線に合うように前記ｎ、Ｎを定める（つまり、Ｎ_A ＝ｎ・Ｎが決まる。）。次に、負荷時、除荷時とも分子鎖の総セグメント数Ｎ_A は同一であるため、負荷時の曲線に整合するよう、各ひずみにおける平均セグメント数Ｎを求める（このＮは変化させる。）。そして、決定された負荷時の平均セグメントＮに一致するよう、式（１９）のパラメータＡないしＥを決定する。本実施形態では、Ｎ＝６．６を使用し、かつ、負荷終了時のＮが除荷時のＮとなるように設定している。

前記フィラーモデル４は、本実施形態ではフィラーとしてのカーボンブラックをモデル化したものが示されている。但しフィラーは、カーボンブラックに限定されるものではなく、例えばシリカ等でも良い。この例では、フィラーモデル４の物理的形状は、実際のゴム中に充填されたカーボンブラックの形状を電子顕微鏡にて撮像した形状に基づいて定められている。カーボンブラック６の形状（最小単位）は、具体的には、図７に略示するように、炭素原子からなる直径数１０nm程度の球状の一次粒子７が不規則に３次元的に結合した葡萄の房状構造をなしている。

またカーボンブラックは、ゴム等の粘弾性材料に比べると数百倍の硬さ（縦弾性係数）を持つため、フィラーモデル４は、本実施形態では粘弾性体ではなく弾性体として取り扱われる。したがって、フィラーモデル４は、材料特性として縦弾性係数が与えられ、変形計算上、応力とひずみとが比例する。フィラーモデル４の量は、ゴムへの配合量を考慮して適宜定められる。

またフィラーモデル４の周りには、前記界面モデル５が設定される。物理的な構造として、フィラーとゴムマトリックスとの界面にこのような特殊な層が形成されていることは確認されていない。しかしながら、上で述べたように、フィラーとマトリックスゴムとの界面では、様々な現象が生じており、とりわけ界面における滑りないし摩擦現象は、比較的大きなエネルギーロスを生じさせる。

本実施形態の界面モデル５は、フィラーモデル４を連続して取り囲み、かつ、薄い厚さで設定されるとともに、マトリックスモデル３とは異なる粘弾性特性が定義される。これは、ゴム材料モデル２において、界面での大きなエネルギーロスの再現を可能とするのに役立つ。なお界面モデル５は必ずしも連続してフィラーモデル４を取り囲む形態に限定されるものではないが、その大部分を囲むことが望ましい。また界面モデル５の厚さｔも特に限定はされないが、大きすぎるとエネルギーロスが過大に表現されてしまう傾向があり、逆に小さすぎても界面でのエネルギーロスの再現が困難となりやすい。種々の実験の結果、前記厚さｔは、例えば１〜２０nm、より好ましくは５〜１０nmであることが望ましい。

また界面モデル５は、粘弾性材料として取り扱われる。従って、マトリックスモデル３と同様に前記式（３）に基づいて応力とひずみとの関係が定義される。本実施形態の界面モデル５は、前記マトリックスモデル３よりも軟い粘弾性特性が定義される。具体的には、応力−ひずみ（又は伸び）曲線における弾性成分の傾き（弾性体の縦弾性係数に相当）をマトリックスモデル３のそれよりも小としている。さらに本実施形態の界面モデル５は、前記マトリックスモデル３よりもヒステリシスロス（ひずみの１サイクル当たりに生じるエネルギーロス）が大きい粘弾性特性が定義される。

次に本実施形態のシミュレーション方法では、設定されたゴム材料モデル２を用いて変形シミュレーションが行われる（ステップＳ３）。変形シミュレーションの具体的な処理手順は、図８に示される。変形シミュレーションでは、先ずデータがコンピュータ装置１に入力される（ステップＳ３１）。入力されるデータには、各要素に定義された節点の位置や材料特性といった情報が含まれる。

コンピュータ装置１では、入力されたデータに基づいて各要素の剛性マトリックスを作成し（ステップＳ３２）、しかる後、全体構造の剛性マトリックスを組み立てる（ステップＳ３３）。全体構造の剛性マトリックスには、既知節点の変位、節点力が導入され（ステップＳ３４）、剛性方程式の解析が行われる。そして、未知節点変位が決定され（ステップＳ３５）、各要素のひずみ、応力、主応力といった物理量を計算し、出力する（ステップＳ３６ないし３７）。ステップＳ３８では、計算を終了させるか否かの判定がなされ、否定的である場合には、ステップＳ３２以降を繰り返す。

シミュレーションは、例えば有限要素法を用いたエンジニアリング系の解析アプリケーションソフトウエア（例えば米国リバモア・ソフトウェア・テクノロジー社で開発・改良されたＬＳ−ＤＹＮＡ等）を用いて行うことができる。

本シミュレーションは、均質化法（漸近展開均質化法）に基づいて行われる態様を例示する。均質化法は、図９に示すように、図３に示した微視構造（ユニットセル）を周期的に持っているゴム材料全体Ｍを表現するｘ_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.
具体的には前記式（３）に加え下記式（２０）の巨視的平衡方程式及び式（２１）の特性変位関数（Ｙ−periodic）が用いられる。

また、本実施形態では、式（３）の定数等を次のように設定した。
Ｃ^R ＝０．２６８
Ｎ＝６．６
Ｔ＝２９６
ｋ_B ＝１．３８０６６×１０^-29
ｎ＝Ｃ^R ／Ｋ_B ／Ｔ＝６．５５８×１０²⁵
Ｎ_A ＝ｎ・Ｎ＝４．３２８×１０²⁶
フィラーモデルの体積含有率 ３０％
フィラーモデルの縦弾性係数Ｅ：１００ＭＰａ
フィラーモデルのポアソン比ν：０．３

また変形シミュレーションは、巨視的モデルＭとして２mm×２mmの矩形状の解析領域を設定し、この巨視的モデルＭに一様な一軸引張変形を発生させるため、図９のｘ２方向に平均ひずみ速度１．０×１０^-5／ｓを加え、ひずみを０．０１きざみで零から２．５まで漸増させるとともに、各ひずみ毎にゴム材料モデル２に作用する応力を計算している。またひずみが２．５に達した後は、逆に０．０１きざみで前記と同じひずみ速度でひずみを零まで漸減させた。各ひずみ状態において、それぞれ平均セグメント数Ｎが計算され、この値は式（３）へ代入され逐次計算が行われる。なおゴム材料モデルは、厚さ方向（図３のＺ軸方向）に変化しないようにArrudaらの３次元８鎖モデルが用いられている。また、マトリックスモデル３及び界面モデル５の平均セグメント数Ｎは次のように設定した。

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

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

前記変形計算が行われると、その結果から必要な物理量を取得することができる（ステップＳ４）。図１０には、マトリックスモデル３、フィラーモデル４及び界面モデル５について、それぞれ単独で変形シミュレーションを行った結果を示している。フィラーモデル４が最も高い弾性を示し、エネルギーロスは生じていない。界面モデル５は、同じ応力でもマトリックスモデル３よりも大きなひずみが生じる粘弾性特性が表れており、かつ、エネルギーロス（曲線のループ面積）も大きくなっている。

また図１１の曲線Ｌａは、上記各モデルを図３の微視構造（セルユニット）として組み入れたゴム材料モデル全体における真応力とひずみとの関係を示している。負荷変形時には、非線形の第１の曲線Ｌａ１が得られた。また除荷変形時では、第１の曲線Ｌａ１とは異なる第２の曲線Ｌａ２が得られた。第２の曲線Ｌａ２は、第１の曲線Ｌａ１よりも軟化しており（低弾性）、ヒステリシスループが再現できた。このループの閉面積を計算することにより、変形１サイクルでのエネルギーロスを計算することができる。

ちなみに、分子鎖１本当たりの平均セグメント数Ｎを一定として Arruda らと同様の変形シミュレーションを行った場合、負荷変形時及び除荷変形時とも、同じ曲線を通り、エネルギーロスを表現し得ない結果となる。また上記実施形態では、平均セグメント数Ｎがパラメータλc の関数である態様を示したが、λc に代えてひずみ速度やひずみ量を用いることができる。この場合、式（１９）と同様に、負荷変形時と除荷変形時とで係数を違えることにより、エネルギーロスを表現することができる。

各曲線Ｌａ１、Ｌａ２の形状は、式（１９）において、係数ＡないしＥを適宜設定し分子鎖１本当たりの平均セグメント数Ｎを定める関数を変えることにより種々設定しうる。したがって、例えば解析しようとするゴム材料に応じて、係数ＡないしＥを種々設定することにより、材料毎にエネルギーロスなどを詳細に検討することができる。これは、粘弾性材料を主要部として用いるタイヤ、ゴルフボールの性能改善などに大いに役立つ。また、シミュレーション対象の材料の応力−ひずみ曲線が分かっている場合、それに応じて前記係数ＡないしＥ等を調整することにより、実際の粘弾性材料と同じ応力ーひずみの関係を表す粘弾性材料のシミュレーションを簡単に行うことができる。

また、図１１には、界面モデル５にマトリックスモデル３と同一の材料特性を与えたモデル（つまり、界面を考慮していないモデル）で同一の実験を行った結果が、曲線Ｌｂとして示されている。曲線Ｌｂは、ゴム材料モデル２において、ひずみが大きく生じる界面の軟質領域がないため、応力ーひずみ曲線の傾きが曲線Ｌａよりも立ち上がっている。また、界面モデル５のように界面で大きなエネルギーロスが生じにくいため、モデル全体のエネルギーロスも小さい。

さらに図１１には、界面モデル５に、曲線Ｌａとは逆の材料特性を与えてシミュレーションを行った結果が、曲線Ｌｃとして示されている。即ち、曲線Ｌｃを示すゴム材料モデル２は、図１２に示すように、界面モデル５の方がマトリックスモデル３よりも高弾性を示すように粘弾性特性を定義したものである。この結果、ゴム材料モデル２の系全体の粘弾性特性が曲線Ｌｂよりも非常に高い勾配を示し、高弾性化されていることが確認できる。なお、この実施形態では、各モデルの平均セグメント数Ｎは次のように設定している。

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

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

図１３には、前記ゴム材料モデル２（曲線Ｌａのもの）の単軸引張の変形シミュレーションにおける一つの微視構造（ユニットセル）について、変形〜回復状態の進行過程を間欠的に視覚化して示す。図１３（Ａ）〜（Ｅ）は、負荷変形時を、同（Ｆ）〜（Ｊ）は除荷変形をそれぞれ示している。また色彩が白く変化しているところは、ひずみの大きい領域を示している。この結果では、フィラーモデル４の界面及びフィラーモデル４、４間に大きなひずみが集中していることが正確に再現されている。特に、フィラーモデル４、４間の距離が小さい部分に大きなひずみが生じていることが分かる。

図１４には、最大ひずみ時における各要素の引張方向のひずみ分布を示し、（Ａ）は、界面を考慮していない図１１の曲線Ｌｂを、（Ｂ）は界面を考慮した同曲線Ｌａのものである。色彩が薄くなるところほど大きなひずみが生じていることを示す。界面を考慮していない（Ａ）のものでは、ひずみの発生が全体的に広範囲にかつ穏やかに生じているのに対して、界面を考慮した（Ｂ）のものでは、フィラーモデル４の周辺に大きなひずみが集中していることが良く再現できている。

また図１５には、各要素毎にひずみの１サイクルのエネルギーロスを計算し、その大小を色彩にて示す。色彩が薄いところほど大きなエネルギーロスが生じていることを示す。またモデルの形状は、図１４の最大変形形状で表しており、（Ａ）は、図１１の曲線Ｌｂを、（Ｂ）は同曲線Ｌａのものである。界面を考慮した（Ｂ）のものでは、フィラーモデル４の界面に集中して大きなエネルギーロスが生じていることが確認できる。

このようなシミュレーションは、例えばゴム材へのフィラーの配合割合、フィラーの種類を種々違えたゴム材について、エネルギーロスを含めた粘弾性特性を精度良く評価することができる。

本実施形態で用いたコンピュータ装置の一例を示す斜視図である。 本実施形態の処理手順を示すフローチャートである。 ゴム材料モデル（微視構造）の一実施形態を示す線図である。 （Ａ）は粘弾性材料の網目構造体を示す斜視図、（Ｂ）は８鎖モデルの一例を示す斜視図である。 （Ａ）、（Ｂ）は分子鎖の接合点の破断を説明する線図である。 パラメータλc と分子鎖１本当たりの平均セグメント数Ｎとの関係を示すグラフである。 カーボンブラックの形状を示す線図である。 変形シミュレーションの手順を示すフローチャートである。 均質化法を説明する微視構造と全体構造との関係を示す。 マトリックスモデル、フィラーモデル及び界面モデルの個々の応力ーひずみ曲線を示す。 ゴム材料モデル（全体構造）のシミュレーション結果として真応力−ひずみとの関係を示すグラフである。 マトリックスモデル及び界面モデルの他の例を示す応力ーひずみ曲線を示す。 （Ａ）〜（Ｊ）は、ゴム材料モデルの変形過程を視覚化して示す線図である。 （Ａ）は界面層未考慮のモデル、（Ｂ）は、界面層を考慮したモデルの引張ひずみの分布図である。 （Ａ）は界面層未考慮のモデル、（Ｂ）は、界面層を考慮したモデルのエネルギーロスの分布図である。 （Ａ）は粘弾性材料、（Ｂ）はその分子鎖１構造を説明する線図、（Ｃ）は１本の分子鎖の拡大図、（Ｄ）はセグメントの拡大図である。 Arrudaらによるシミュレーション結果として、応力ーひずみ曲線を示す。

符号の説明

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

Claims (6)

1. フィラーが配合されたゴム材料の変形をシミュレーションするゴム材料のシミュレーション方法であって、
   前記ゴム材料を数値解析が可能な要素でモデル化したゴム材料モデルを設定するステップと、
   前記ゴム材料モデルに条件を設定して変形計算を行うステップと、
   前記変形計算から必要な物理量を取得するステップとを含むとともに、
   前記ゴム材料モデルは、ゴムマトリックスをモデル化したマトリックスモデルと、フィラーをモデル化したフィラーモデルと、前記マトリックスモデルと前記フィラーモデルとの間の界面を形成する界面モデルとを含み、
   前記界面モデルは、前記フィラーモデルを連続して取り囲みかつ厚さを有するとともに前記マトリックスモデルとは異なる粘弾性特性が定義されることを特徴とするゴム材料のシミュレーション方法。
2. 前記界面モデルは、前記マトリックスモデルよりも軟い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法。
3. 前記界面モデルは、前記マトリックスモデルよりも硬い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法。
4. 前記界面モデルは、前記マトリックスモデルよりもヒステリシスロスが大きい粘弾性特性が定義されることを特徴とする請求項１乃至３のいずれかに記載のゴム材料のシミュレーション方法。
5. 前記界面モデルの厚さは１〜２０nmである請求項１記載のゴム材料のシミュレーション方法。
6. 前記フィラーモデルは、ゴムよりも硬い弾性率が定義される請求項１記載のゴム材料のシミュレーション方法。