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

しかしながら、Aruudaらの非特許文献１は、粘弾性材料の特徴の一つでもあるエネルギーロスが考慮されていない。このため、図２０の状態Ｑから荷重を除荷するシミュレーションを行うと、該曲線と実質的に同じ経路を通ってひずみの回復が行われる。これでは、ヒステリシスループが形成されずエネルギーロスの計算を行うことができない。

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

以上のように、本発明は、エネルギーロスを発現させることが可能な特性を粘弾性材料モデル自体に組み入れることを基本として、より精度良く粘弾性材料を変形シミュレーションしうる粘弾性材料のシミュレーション方法を提供することを目的としている。

本発明のうち請求項１記載の発明は、粘弾性材料の変形をシミュレーションする粘弾性材料のシミュレーション方法であって、粘弾性材料を数値解析が可能な要素でモデル化した粘弾性材料モデルを設定するステップと、前記粘弾性材料モデルに条件を設定して変形計算を行うステップと、前記変形計算から必要な物理量を取得するステップとを含むとともに、前記粘弾性材料モデルは、ゴムマトリックス部分がモデル化されたマトリックスモデルと、このマトリックスモデルの中に分散して配されかつフィラーがモデル化されたフィラーモデルと、前記マトリックスモデルと前記フィラーモデルとの間に介在しかつ両者の間の界面を形成する界面モデルとを含み、しかも前記界面モデルは、前記フィラーモデルを連続して取り囲みかつ厚さを有するとともに前記マトリックスモデルとは異なる粘弾性特性が定義されるとともに、前記マトリックスモデル及び前記界面モデルは、下記式（１）で示される応力とひずみとの関係が定義され、かつ、該式（１）における１本の分子鎖当たりの平均セグメント数Ｎが、負荷変形時と除荷変形時とで異なることによりエネルギーロスを表現可能なことを特徴とする粘弾性材料のシミュレーション方法である。

また請求項２記載の発明は、前記平均セグメント数Ｎは、負荷変形時ではひずみに関するパラメータに応じて増大し、除荷変形時では一定値をなすことを特徴とする請求項１記載の粘弾性材料のシミュレーション方法である。

また請求項３記載の発明は、前記ひずみに関するパラメータは、ひずみの１次の不変量Ｉ₁ であることを特徴とする請求項２記載の粘弾性材料のシミュレーション方法である。

また請求項４記載の発明は、前記ひずみに関するパラメータは、ひずみ量又はひずみ速度であることを特徴とする請求項２記載の粘弾性材料のシミュレーション方法である。また請求項５記載の発明は、前記界面モデルは、前記マトリックスモデルよりも軟い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法であり、請求項６記載の発明は、前記界面モデルは、前記マトリックスモデルよりも硬い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法であり、請求項７記載の発明は、前記界面モデルは、前記マトリックスモデルよりもヒステリシスロスが大きい粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法であり、請求項８記載の発明は、前記界面モデルの厚さは１〜２０nmである請求項１記載のゴム材料のシミュレーション方法である。

本発明の粘弾性材料のシミュレーション方法では、粘弾性材料のエネルギーロスが表現できる。従って、実際の粘弾性材料の特性に即した精度の高いシミュレーションが可能である。またエネルギーロスは、式（１）における１本の分子鎖当たりの平均セグメント数Ｎを、負荷変形時と除荷変形時とで異なるパラメータとして定めることによって表現しうる結果、特許文献１のような数多くの実験データの収集工程を不要とし、より一層、解析コストと手間を減じることができる。

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

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

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

この実施形態の粘弾性材料モデル２は、後述する変形シミュレーションにおいて平面ひずみ状態で解析が行われる。したがってＺ方向にはひずみを持たない。粘弾性材料モデル２の一辺の長さは、例えば縦横それぞれ２mmである。該粘弾性材料モデル２は、四辺形をなす９つの平面要素２ａないし２ｉにより定義される。

前記各要素２ａないし２ｉには、材料特性として、前記式（１）の応力とひずみとの関係が定義される。したがって、各要素の２ａないし２ｉについて、それぞれ応力とひずみとの関係が得られる。式（１）の導出過程について簡単に触れる。粘弾性材料は、変形による体積変化が小さいため、それを無視することができる。このため、非圧縮性の粘弾性材料では、非圧縮条件により密度を一定とすることができ、kirchhoff 応力Ｓ_ijは、式（４）で表すことができる。なお、Ｅ_ijはGreen ひずみの成分、ｐは静水圧（境界条件により定まる）、Ｘ_i は、応力もひずみも０である状態Ｃ₀ における任意の物体点Ｐの位置、ｘ_j は変形した状態Ｃにおける物体点Ｐの位置である。

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

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

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

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

本実施形態の粘弾性材料モデル２は、前記式（１）において、１本の分子鎖当たりの平均セグメント数Ｎが負荷変形時と除荷変形時とで異なるパラメータとして定義される。これにより、粘弾性材料モデル２の各要素２ａないし２ｉは、エネルギーロスを表現することができる。この点がArrudaらの提案とは基本的に異なる。ここで、負荷変形時とは、微小時間の間で粘弾性材料モデル２のひずみが増大する変形であり、逆に除荷変形時とは、ひずみが減少する変形である。

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

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

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

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

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

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

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

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

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

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

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

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

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

本実施形態では、前記粘弾性材料モデル２に１軸の引張試験を行い、引張のひずみを０．０１きざみで０〜２．０まで漸増させるとともに、各ひずみ毎に粘弾性材料モデル２に作用する応力を計算する。またひずみが２．０に達した後は、逆に０．０１きざみでひずみを０まで漸減させる。各ひずみ状態において、それぞれ平均セグメント数Ｎが計算され、この値は式（１）のひずみエネルギー関数へと代入され逐次計算が行われる。なお粘弾性材料モデルは、厚さ方向（図３のＺ軸方向）に変化しないようにArrudaらの３次元８鎖モデルが用いられている。

シミュレーションは、コンピュータ装置１を用いて行われる。上述のように粘弾性材料モデル２の材料特性と、変形条件とが決定できれば、ステップＳ３２〜Ｓ３７の計算自体は、例えば有限要素法を用いたエンジニアリング系の解析アプリケーションソフトウエア（例えば米国リバモア・ソフトウェア・テクノロジー社で開発・改良されたＬＳ−ＤＹＮＡ等）を用いて行うことができる。

また、本実施形態では、式（１）の定数等を次のように設定した。
Ｃ^R ＝０．２６８ Ｎ＝６．６ Ｔ＝２９６
ｋ_B ＝１．３８０６６×１０^-29 ｎ＝Ｃ^R ／Ｋ_B ／Ｔ＝６．５５８×１０²⁵ Ｎ_A ＝ｎ・Ｎ＝４．３２８×１０²⁶

前記変形計算が行われると、その結果から必要な物理量を取得することができる（ステップＳ４）。図８には、粘弾性材料モデルの物理量として、真応力とひずみとの関係を示す。粘弾性材料モデル２の負荷変形時には、非線形の第１の曲線Ｌ１が得られた。また除荷変形時では、第１の曲線Ｌ１とは異なる第２の曲線Ｌ２が得られた。第２の曲線Ｌ２は、第１の曲線Ｌ１よりも低弾性を示し、ヒステリシスループが再現できた。このループの閉面積を計算することにより、変形１サイクルでのエネルギーロスを計算しうる。

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

なお本発明の比較例として、分子鎖１本当たりの平均セグメント数Ｎを一定としてArrudaらと同様の変形シミュレーションを行った。この例では平均セグメント数Ｎの値は、次のように設定した。
Ｎ＝６．６

図９には、比較例の結果として、粘弾性材料モデル２の真応力とひずみとの関係を示す。この場合、粘弾性材料モデル２は、負荷変形時及び除荷変形時とも、同じ曲線Ｌ３を通る。したがって、粘弾性材料の特徴であるエネルギーロスを表現し得ない。これは、Ｎの値とは関係がなく、Ｎが負荷変形時と除荷変形時とで同じであることに基づいている。

また上記実施形態では、平均セグメント数Ｎが１次のひずみの不変量Ｉ₁ の平方根であるλc の関数である態様を示したが、λc に代えてひずみ速度やひずみ量を用いることができる。この場合、式（１９）と同様に、負荷変形時と除荷変形時とで係数を違えることにより、エネルギーロスを表現することができる。

図１０には、本発明の他の実施形態を示す。
この実施形態では、粘弾性材料モデル２が、ゴムマトリックス中にフィラー（充填材）を配合したゴムをモデル化したものが例示される。即ち、粘弾性材料モデル２は、ゴムマトリックス部分がモデル化されたマトリックスモデル３（最も濃い部分）と、このマトリックスモデル３の中に分散して配されかつフィラー（充填材）がモデル化されたフィラーモデル４（白い部分）と、前記マトリックスモデル３と前記フィラーモデル４との間に介在しかつ両者の間の界面を形成する界面モデル５（やや黒い部分）とを含むものが例示される。図１０の粘弾性材料モデル２の縦、横の長さは、３００ｎｍ×３００ｎｍとしている。

前記マトリックスモデル３は、粘弾性材料モデル２の主要部を構成し、かつ、例えば三角形ないし四辺形の要素を用いて表現されている。マトリックスモデル３を構成する各要素には、材料特性として、前述した式（１）の応力とひずみとの関係が定義される。式（１）中の定数などは条件に応じ前記実施形態とは異ならせ得るのは言うまでもない。

前記フィラーモデル４は、本実施形態ではフィラーとしてのカーボンブラックをモデル化したものが示されている。但し、フィラーは、カーボンブラックに限定されるものではなく、例えばシリカ等でも良い。この例では、フィラーモデル４の物理的形状は、実際のゴム中に充填されたカーボンブラックの形状を電子顕微鏡にて撮像した形状に基づいて定められている。カーボンブラック６の形状（最小単位）は、具体的には、図１１に略示するように、炭素原子からなる直径数１０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 (8)

1. 粘弾性材料の変形をシミュレーションする粘弾性材料のシミュレーション方法であって、
   粘弾性材料を数値解析が可能な要素でモデル化した粘弾性材料モデルを設定するステップと、
   前記粘弾性材料モデルに条件を設定して変形計算を行うステップと、
   前記変形計算から必要な物理量を取得するステップとを含むとともに、
   前記粘弾性材料モデルは、ゴムマトリックス部分がモデル化されたマトリックスモデルと、このマトリックスモデルの中に分散して配されかつフィラーがモデル化されたフィラーモデルと、前記マトリックスモデルと前記フィラーモデルとの間に介在しかつ両者の間の界面を形成する界面モデルとを含み、しかも前記界面モデルは、前記フィラーモデルを連続して取り囲みかつ厚さを有するとともに前記マトリックスモデルとは異なる粘弾性特性が定義されるとともに、
   前記マトリックスモデル及び前記界面モデルは、下記式（１）で示される応力とひずみとの関係が定義され、かつ、該式（１）における１本の分子鎖当たりの平均セグメント数Ｎが、負荷変形時と除荷変形時とで異なることによりエネルギーロスを表現可能なことを特徴とする粘弾性材料のシミュレーション方法。
   
2. 前記平均セグメント数Ｎは、負荷変形時ではひずみに関するパラメータに応じて増大し、除荷変形時では一定値をなすことを特徴とする請求項１記載の粘弾性材料のシミュレーション方法。
3. 前記ひずみに関するパラメータは、ひずみの１次の不変量Ｉ₁ であることを特徴とする請求項２記載の粘弾性材料のシミュレーション方法。
4. 前記ひずみに関するパラメータは、ひずみ量又はひずみ速度であることを特徴とする請求項２記載の粘弾性材料のシミュレーション方法。
5. 前記界面モデルは、前記マトリックスモデルよりも軟い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法。
6. 前記界面モデルは、前記マトリックスモデルよりも硬い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法。
7. 前記界面モデルは、前記マトリックスモデルよりもヒステリシスロスが大きい粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法。
8. 前記界面モデルの厚さは１〜２０nmである請求項１記載のゴム材料のシミュレーション方法。