(57)【特許請求の範囲】
【請求項１】
粘弾性材料の変形をシミュレーションする粘弾性材料のシミュレーション方法であって、
粘弾性材料を数値解析が可能な要素でモデル化した粘弾性材料モデルを設定するステップと、
前記粘弾性材料モデルに条件を設定して変形計算を行うステップと、
前記変形計算から必要な物理量を取得するステップとを含むとともに、
前記粘弾性材料モデルは、ゴムマトリックス部分がモデル化されたマトリックスモデルと、このマトリックスモデルの中に分散して配されかつフィラーがモデル化されたフィラーモデルと、前記マトリックスモデルと前記フィラーモデルとの間に介在しかつ両者の間の界面を形成する界面モデルとを含み、しかも前記界面モデルは、前記フィラーモデルを連続して取り囲みかつ厚さを有するとともに前記マトリックスモデルとは異なる粘弾性特性が定義されるとともに、
前記マトリックスモデル及び前記界面モデルは、下記式（１）で示される応力とひずみとの関係が定義され、かつ、該式（１）における１本の分子鎖当たりの平均セグメント数Ｎが、負荷変形時と除荷変形時とで異なることによりエネルギーロスを表現可能なことを特徴とする粘弾性材料のシミュレーション方法。
【数１】
0003668239.tif 000012
【請求項２】
前記平均セグメント数Ｎは、負荷変形時ではひずみに関するパラメータに応じて増大し、除荷変形時では一定値をなすことを特徴とする請求項１記載の粘弾性材料のシミュレーション方法。
【請求項３】
前記ひずみに関するパラメータは、ひずみの１次の不変量Ｉ1 であることを特徴とする請求項２記載の粘弾性材料のシミュレーション方法。
【請求項４】
前記ひずみに関するパラメータは、ひずみ量又はひずみ速度であることを特徴とする請求項２記載の粘弾性材料のシミュレーション方法。
【請求項５】
前記界面モデルは、前記マトリックスモデルよりも軟い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法。
【請求項６】
前記界面モデルは、前記マトリックスモデルよりも硬い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法。
【請求項７】
前記界面モデルは、前記マトリックスモデルよりもヒステリシスロスが大きい粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法。
【請求項８】
前記界面モデルの厚さは１〜２０nmである請求項１記載のゴム材料のシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ゴムなどの粘弾性材料の変形状態を精度良く解析するのに役立つ粘弾性材料のシミュレーション方法に関する。
【背景技術】
【０００２】
ゴムをはじめとする粘弾性材料は、タイヤ、スポーツ用品、その他各種の工業製品に使用される。粘弾性材料は、負荷を受けると大きく変形し、負荷を完全に取り除くと元の状態へと復元しうる。また粘弾性材料は、変形に際して応力とひずみとが比例せず、かつ、荷重の負荷時と除荷時とでは応力−ひずみ曲線の経路を異にする非線形性を有している。試作の手間とコストとを減じるために、粘弾性材料の変形過程などを予めコンピュータを用いてシミュレーションすることが行われている。従来の粘弾性材料のシミュレーション方法としては、下記の特許文献１や非特許文献１等が知られている。
【０００３】
【特許文献１】特開２００２−３６５２０５号公報
【非特許文献１】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を生じるようなひずみエネルギー関数Ｗが存在する物質として定義される。換言すれば、ひずみエネルギー関数は、粘弾性材料が変形したときに蓄えられたポテンシャルエネルギの存在を仮定的に示す。従って、ひずみエネルギー関数Ｗの微分勾配から応力とひずみ関係を得ることができる。
【数２】
0003668239.tif 000002
【００１１】
Aruudaらは、非ガウス鎖理論により、粘弾性材料の変形が大きくなるとエントロピー変化が急激に大きくなる（分子鎖が伸びきって配向する）と考え、式（３）の粘弾性材料（ゴム弾性体）のひずみエネルギ関数Ｗを示した。そして、このひずみエネルギー関数Ｗを上記式に代入することにより、粘弾性材料の応力とひずみとの関係を取り出すことができる。
【数３】
0003668239.tif 000003
【発明の開示】
【発明が解決しようとする課題】
【００１２】
Aruudaらの応力とひずみとの関係を用いて、粘弾性材料の変形シミュレーションを行うことにより、例えば図２０に示すように、１軸引張変形シミュレーションにおいて、非線形な応力とひずみとの関係を得ることができる。この結果は、荷重の負荷変形時における実測値と良い相関を示す。
【００１３】
しかしながら、Aruudaらの非特許文献１は、粘弾性材料の特徴の一つでもあるエネルギーロスが考慮されていない。このため、図２０の状態Ｑから荷重を除荷するシミュレーションを行うと、該曲線と実質的に同じ経路を通ってひずみの回復が行われる。これでは、ヒステリシスループが形成されずエネルギーロスの計算を行うことができない。
【００１４】
発明者らは、Aruudaらのモデルを前提として種々の実験を試みた。先ず粘弾性材料は、上述のように複雑に絡み合った長い分子鎖ｃが伸びることによって数百％にも達し得る大ひずみが許容されると考えられる。そこで、発明者らは、荷重の負荷における変形過程において、粘弾性材料の分子鎖ｃの互いに絡み合った部分がほどけたり（接合点ｂの数の減少）、或いは荷重の除荷によってさらなる絡み合い（接合点ｂの数の増加）が生じ得るとの仮定を立てた。換言すれば、この仮定は、前記式（１）における１本の分子鎖当たりの平均セグメント数Ｎは、負荷ないし除荷変形時において可変のパラメータであることを意味している。発明者らは、この仮定を検証したところ、粘弾性材料のヒステリシスを表現可能であることが分かった。
【００１５】
以上のように、本発明は、エネルギーロスを発現させることが可能な特性を粘弾性材料モデル自体に組み入れることを基本として、より精度良く粘弾性材料を変形シミュレーションしうる粘弾性材料のシミュレーション方法を提供することを目的としている。
【課題を解決するための手段】
【００１６】
本発明のうち請求項１記載の発明は、粘弾性材料の変形をシミュレーションする粘弾性材料のシミュレーション方法であって、粘弾性材料を数値解析が可能な要素でモデル化した粘弾性材料モデルを設定するステップと、前記粘弾性材料モデルに条件を設定して変形計算を行うステップと、前記変形計算から必要な物理量を取得するステップとを含むとともに、前記粘弾性材料モデルは、ゴムマトリックス部分がモデル化されたマトリックスモデルと、このマトリックスモデルの中に分散して配されかつフィラーがモデル化されたフィラーモデルと、前記マトリックスモデルと前記フィラーモデルとの間に介在しかつ両者の間の界面を形成する界面モデルとを含み、しかも前記界面モデルは、前記フィラーモデルを連続して取り囲みかつ厚さを有するとともに前記マトリックスモデルとは異なる粘弾性特性が定義されるとともに、前記マトリックスモデル及び前記界面モデルは、下記式（１）で示される応力とひずみとの関係が定義され、かつ、該式（１）における１本の分子鎖当たりの平均セグメント数Ｎが、負荷変形時と除荷変形時とで異なることによりエネルギーロスを表現可能なことを特徴とする粘弾性材料のシミュレーション方法である。
【数４】
0003668239.tif 000004
【００１７】
また請求項２記載の発明は、前記平均セグメント数Ｎは、負荷変形時ではひずみに関するパラメータに応じて増大し、除荷変形時では一定値をなすことを特徴とする請求項１記載の粘弾性材料のシミュレーション方法である。
【００１８】
また請求項３記載の発明は、前記ひずみに関するパラメータは、ひずみの１次の不変量Ｉ1 であることを特徴とする請求項２記載の粘弾性材料のシミュレーション方法である。
【００１９】
また請求項４記載の発明は、前記ひずみに関するパラメータは、ひずみ量又はひずみ速度であることを特徴とする請求項２記載の粘弾性材料のシミュレーション方法である。また請求項５記載の発明は、前記界面モデルは、前記マトリックスモデルよりも軟い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法であり、請求項６記載の発明は、前記界面モデルは、前記マトリックスモデルよりも硬い粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法であり、請求項７記載の発明は、前記界面モデルは、前記マトリックスモデルよりもヒステリシスロスが大きい粘弾性特性が定義されることを特徴とする請求項１記載のゴム材料のシミュレーション方法であり、請求項８記載の発明は、前記界面モデルの厚さは１〜２０nmである請求項１記載のゴム材料のシミュレーション方法である。
【発明の効果】
【００２０】
本発明の粘弾性材料のシミュレーション方法では、粘弾性材料のエネルギーロスが表現できる。従って、実際の粘弾性材料の特性に即した精度の高いシミュレーションが可能である。またエネルギーロスは、式（１）における１本の分子鎖当たりの平均セグメント数Ｎを、負荷変形時と除荷変形時とで異なるパラメータとして定めることによって表現しうる結果、特許文献１のような数多くの実験データの収集工程を不要とし、より一層、解析コストと手間を減じることができる。
【発明を実施するための最良の形態】
【００２１】
以下、本発明の実施の一形態を図面に基づき説明する。
図１には、本発明のシミュレーション方法を実施するためのコンピュータ装置１が示されている。このコンピュータ装置１は、本体１ａと、入力手段としてのキーボード１ｂ、マウス１ｃと、出力手段としてのディスプレイ装置１ｄとから構成されている。本体１ａには、図示していないが、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリー、磁気ディスクなどの大容量記憶装置、ＣＤ−ＲＯＭやフレキシブルディスクのドライブ１ａ１、１ａ２が設けられている。そして、前記大容量記憶装置には後述する本発明のシミュレーション方法を実行するための処理手順（プログラム）が記憶されている。コンピュータ装置１にはＥＷＳなどが好適である。
【００２２】
図２には、シミュレーション方法の処理手順の一例が示される。本実施形態では、先ず粘弾性材料モデルが設定される（ステップＳ１）。図３には、粘弾性材料モデル２の一例が視覚化して示されている。
【００２３】
該粘弾性材料モデル２は、解析しようとする粘弾性材料（実在するか否かを問わない）が、有限個の小さな要素２ａ、２ｂ、２ｃ…に置き換えられたものである。各要素２ａ、２ｂ、２ｃ…は、数値解析が可能に定義される。数値解析が可能とは、例えば有限要素法、有限体積法、差分法又は境界要素法といった数値解析法により、各要素ないし系全体についての変形計算が可能なことを意味する。具体的には、各要素２ａ、２ｂ、２ｃ…について、座標系における節点座標値、要素形状、材料特性などが定義される。各要素２ａ、２ｂ、２ｃ…には、例えば２次元平面としての三角形ないし四角形の要素、３次元要素としては、例えば４ないし６面体の要素が好ましく用いられる。これにより、粘弾性材料モデル２は、前記コンピュータ装置１にて取り扱い可能な数値データを構成しうる。
【００２４】
この実施形態の粘弾性材料モデル２は、後述する変形シミュレーションにおいて平面ひずみ状態で解析が行われる。したがってＺ方向にはひずみを持たない。粘弾性材料モデル２の一辺の長さは、例えば縦横それぞれ２mmである。該粘弾性材料モデル２は、四辺形をなす９つの平面要素２ａないし２ｉにより定義される。
【００２５】
前記各要素２ａないし２ｉには、材料特性として、前記式（１）の応力とひずみとの関係が定義される。したがって、各要素の２ａないし２ｉについて、それぞれ応力とひずみとの関係が得られる。式（１）の導出過程について簡単に触れる。粘弾性材料は、変形による体積変化が小さいため、それを無視することができる。このため、非圧縮性の粘弾性材料では、非圧縮条件により密度を一定とすることができ、kirchhoff 応力Ｓijは、式（４）で表すことができる。なお、ＥijはGreen ひずみの成分、ｐは静水圧（境界条件により定まる）、Ｘi は、応力もひずみも０である状態Ｃ0 における任意の物体点Ｐの位置、ｘj は変形した状態Ｃにおける物体点Ｐの位置である。
【数５】
0003668239.tif 000005
【００２６】
Cauchy応力の成分σijと、kirchhoff 応力の成分Ｓmnとは下記式（５）の関係を有    するから、式（４）から式（６）を得ることができる。
【数６】
0003668239.tif 000006
【００２７】
ここで、体積が一定の場合、物体の体積変化率を示す式（５）のＪは１とみなせる。また、式（３）のひずみエネルギー関数は、左Cauchy-Green変形テンソルＡijのひずみの  １次の不変量Ｉ1 だけの関数であるため、式（６）から式（７）が得られる。
【数７】
0003668239.tif 000007
【００２８】
また、下記の関係式（８ないし１０）を用いることにより、式（７）の速度形式表示は式（１１）になる。
【数８】
0003668239.tif 000008
【００２９】
そして、体積一定条件下で、式（１１）に示すCauchy応力のJaumann 速度をkirchhoff 応力のJaumann 速度に置き換えても本質的な差異はないことと、変形速度テンソルＤをひずみ速度テンソルに置き換えることとにより、非圧縮性の粘弾性材料の構成式の速度形式表示を前記式（１）で表すことができる。
【００３０】
本実施形態の粘弾性材料モデル２は、前記式（１）において、１本の分子鎖当たりの平均セグメント数Ｎが負荷変形時と除荷変形時とで異なるパラメータとして定義される。これにより、粘弾性材料モデル２の各要素２ａないし２ｉは、エネルギーロスを表現することができる。この点がArrudaらの提案とは基本的に異なる。ここで、負荷変形時とは、微小時間の間で粘弾性材料モデル２のひずみが増大する変形であり、逆に除荷変形時とは、ひずみが減少する変形である。
【００３１】
図４（Ａ）には、粘弾性材料の巨視的な３次元の網目構造体ｈ、また図４（Ｂ）には、この網目構造体ｈを構成する８鎖モデルｇの１つを拡大して示す。この例の網目構造体ｈは、幅方向、高さ方向及び奥行き方向にそれぞれ８鎖モデルｇがｋ個結合したものとする。
【００３２】
いま、網目構造体ｈに含まれる接合点ｂの総数を「からみ数」として符号ｍで表すと、からみ数ｍは、式（１２）で表すことができる。同様に、網目構造体ｈに含まれる分子鎖ｃの数（即ち、粘弾性材料モデル２の単位体積中に含まれる分子鎖の数）をｎとすると、このｎは、式（１３）で表すことができる。
ｍ＝（ｋ＋１）3 ＋ｋ3       …式（１２）
ｎ＝８ｋ3                   …式（１３）
【００３３】
ここで、ｋは十分に大きい数とすると、ｋの３次項以外を省略して上式はそれぞれ次のような式（１４）、（１５）で表すことができ、さらにこれらの関係から、からみ数ｍは、ｎを用いて、式（１６）で表すことができる。
ｍ＝２ｋ3   …式（１４）
ｎ＝８ｋ3   …式（１５）
ｍ＝ｎ／４  …式（１６）
【００３４】
さらに、粘弾性材料モデルは、変形してもその中に含まれる分子鎖のセグメントの総数ＮA は変化しないと仮定すると、式（１７）〜（１８）が成り立つ。
ＮA ＝ｎ・Ｎ              …式（１７）
Ｎ＝ＮA ／ｎ＝ＮA ／４ｍ  …式（１８）
【００３５】
発明者らは、シミュレーションにおいて粘弾性材料モデル２のエネルギーロスを表現するために、荷重の負荷過程及びこれに続くひずみの回復過程において、粘弾性材料の分子鎖ｃのからみ数ｍが変化するとの着想を試みた。即ち、図５（Ａ）に示すように、例えば一つの接合点ｂで接合されている分子鎖ｃ１ないしｃ４に矢印方向の引張応力が作用すると、各分子鎖ｃ１ないしｃ４は伸び、接合点ｂは大きなひずみを受けて破損するものと考えた。
【００３６】
この結果、図５（Ｂ）に示すように、これまで２本であった分子鎖ｃ１及びｃ２は、あたかも１本の長い分子鎖ｃ５のように振る舞うものと考えられる。分子鎖ｃ３及びｃ４についても同様である。そして、このような現象は、粘弾性材料モデル２の負荷変形が進むにつれて逐次発生していくものと考えられる。
【００３７】
以上のような接合点ｂの破損は、粘弾性材料モデル２におけるからみ数ｍの減少に他ならない。粘弾性材料モデル自体は材料の出入りがないため、その中に含まれる分子鎖ｃのセグメントの総数ＮA が変化しないと仮定すると、粘弾性材料モデル２の負荷変形が進むにつれて、１本の分子鎖ｃに含まれる平均セグメント数Ｎが増加することは式（１８）から明らかとなる。つまり、前記平均セグメント数Ｎは、一定ではなく、負荷変形時と除荷変形時とで異なる可変のパラメータとすることが適切と考えられる。これによって、粘弾性材料モデルの変形時におけるエネルギーロスが再現できる。これについては、さらに詳細を後で述べる。
【００３８】
このような様子をコンピュータ上でより好ましく再現するためには、前記平均セグメント数Ｎは、負荷変形時では、そのひずみに関するパラメータに応じて増大させることが特に好適である。前記ひすみに関するパラメータとしては、特に制限されるものではないが、例えば、ひずみ量、ひずみ速度又はひずみの１次の不変量Ｉ1 などが挙げられる。本実施形態では、前記平均セグメント数Ｎが、下記式（１９）に示すように、当該粘弾性材料モデル２の各要素２ａないし２ｉ個々において、それぞれ１次のひずみの不変量Ｉ1 を平方根であるパラメータλc の関数となるものを示す。
【００３９】
【数９】
0003668239.tif 000009
【００４０】
なお式（１９）は、種々の実験によって定めた一例であり、上記ＡないしＥは、いずれも定数であって、ゴム試験片の単純な１軸引張試験などの実測結果からように定めることができる。この例では、上記定数を次のように設定している。
Ａ＝＋２．９４９３
Ｂ＝−５．８０２９
Ｃ＝＋５．５２２０
Ｄ＝−１．３５８２
Ｅ＝＋０．１３２５
【００４１】
図６には、粘弾性材料モデル２の各要素２ａないし２ｉの負荷変形時における平均セグメント数Ｎとパラメータλc との関係を示す。ひずみに関するパラメータであるλc が大きくなると、平均セグメント数Ｎは滑らかに増大する。この例では、パラメータλc の上限は２．５である。後述する粘弾性材料モデル２の変形シミュレーションでは、粘弾性材料モデル２の各要素２ａないし２ｉについて、負荷変形時においてはパラメータλc が常時計算される。計算されたλc は、式（１９）に代入される。これにより、当該要素の当該ひずみ状態における平均セグメント数Ｎが計算できる。
【００４２】
他方、粘弾性材料モデル２の除荷変形時には、平均セグメント数Ｎは一定値としている。平均セグメント数Ｎを決定する方法としては、例えば次の方法がある。解析対象となる粘弾性材料において、一つの応力−ひずみ曲線がある場合、先ず除荷時の曲線に合うように前記ｎ、Ｎを定める（つまり、ＮA ＝ｎ・Ｎが決まる。）。次に、負荷時、除荷時とも分子鎖の総セグメント数ＮA は同一であるため、負荷時の曲線に整合するよう、各ひずみにおける平均セグメント数Ｎを求める（このＮは変化させる。）。そして、決定された負荷時の平均セグメントＮに一致するよう、式（１９）のパラメータＡないしＥを決定する。本実施形態では、Ｎ＝６．６を使用し、かつ、負荷終了時のＮが除荷時のＮとなるように設定している。
【００４３】
次に本実施形態のシミュレーション方法では、設定された粘弾性材料モデル２を用いて変形シミュレーションが行われる（ステップＳ３）。変形シミュレーションの具体的な処理手順は、図７に示される。変形シミュレーションでは、先ずデータがコンピュータ装置１に入力される（ステップＳ３１）。入力されるデータには、各要素に定義された節点の位置や材料特性といった情報が含まれる。
【００４４】
コンピュータ装置１では、入力されたデータに基づいて各要素の剛性マトリックスを作成し（ステップＳ３２）、しかる後、全体構造の剛性マトリックスを組み立てる（ステップＳ３３）。全体構造の剛性マトリックスには、既知節点の変位、節点力が導入され（ステップＳ３４）、剛性方程式の解析が行われる。そして、未知節点変位が決定される（ステップＳ３５）。そして、各要素のひずみ、応力、主応力といった物理量を計算し、同出力する（ステップＳ３６ないし３７）。ステップＳ３８では、計算を終了させるか否かの判定がなされ、否定的である場合には、ステップＳ３２以降を繰り返す。
【００４５】
本実施形態では、前記粘弾性材料モデル２に１軸の引張試験を行い、引張のひずみを０．０１きざみで０〜２．０まで漸増させるとともに、各ひずみ毎に粘弾性材料モデル２に作用する応力を計算する。またひずみが２．０に達した後は、逆に０．０１きざみでひずみを０まで漸減させる。各ひずみ状態において、それぞれ平均セグメント数Ｎが計算され、この値は式（１）のひずみエネルギー関数へと代入され逐次計算が行われる。なお粘弾性材料モデルは、厚さ方向（図３のＺ軸方向）に変化しないようにArrudaらの３次元８鎖モデルが用いられている。
【００４６】
シミュレーションは、コンピュータ装置１を用いて行われる。上述のように粘弾性材料モデル２の材料特性と、変形条件とが決定できれば、ステップＳ３２〜Ｓ３７の計算自体は、例えば有限要素法を用いたエンジニアリング系の解析アプリケーションソフトウエア（例えば米国リバモア・ソフトウェア・テクノロジー社で開発・改良されたＬＳ−ＤＹＮＡ等）を用いて行うことができる。
【００４７】
また、本実施形態では、式（１）の定数等を次のように設定した。
ＣR ＝０．２６８                                                                Ｎ＝６．６                                                                      Ｔ＝２９６
ｋB ＝１．３８０６６×１０-29                                                   ｎ＝ＣR ／ＫB ／Ｔ＝６．５５８×１０25                                          ＮA ＝ｎ・Ｎ＝４．３２８×１０26
【００４８】
前記変形計算が行われると、その結果から必要な物理量を取得することができる（ステップＳ４）。図８には、粘弾性材料モデルの物理量として、真応力とひずみとの関係を示す。粘弾性材料モデル２の負荷変形時には、非線形の第１の曲線Ｌ１が得られた。また除荷変形時では、第１の曲線Ｌ１とは異なる第２の曲線Ｌ２が得られた。第２の曲線Ｌ２は、第１の曲線Ｌ１よりも低弾性を示し、ヒステリシスループが再現できた。このループの閉面積を計算することにより、変形１サイクルでのエネルギーロスを計算しうる。
【００４９】
各曲線Ｌ１、Ｌ２の形状は、式（１９）において、係数ＡないしＥを適宜設定し分子鎖１本当たりの平均セグメント数Ｎを表す関数を変えることにより種々設定しうる。したがって、例えば解析しようとするゴム材料に応じて、係数ＡないしＥを種々設定することにより、材料毎にエネルギーロスなどを詳細に検討することができる。これは、粘弾性材料を主要部として用いるタイヤ、ゴルフボールの性能改善などに大いに役立つ。また、シミュレーション対象の材料の応力−ひずみ曲線が分かっている場合、それに応じて前記係数ＡないしＥ等を調整することにより、実際の粘弾性材料と同じ応力ーひずみの関係を表す粘弾性材料のシミュレーションを簡単に行うことができる。
【００５０】
なお本発明の比較例として、分子鎖１本当たりの平均セグメント数Ｎを一定としてArrudaらと同様の変形シミュレーションを行った。この例では平均セグメント数Ｎの値は、次のように設定した。
Ｎ＝６．６
【００５１】
図９には、比較例の結果として、粘弾性材料モデル２の真応力とひずみとの関係を示す。この場合、粘弾性材料モデル２は、負荷変形時及び除荷変形時とも、同じ曲線Ｌ３を通る。したがって、粘弾性材料の特徴であるエネルギーロスを表現し得ない。これは、Ｎの値とは関係がなく、Ｎが負荷変形時と除荷変形時とで同じであることに基づいている。
【００５２】
また上記実施形態では、平均セグメント数Ｎが１次のひずみの不変量Ｉ1 の平方根であるλ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）が用いられる。
【００６１】
【数１０】
0003668239.tif 000010
【数１１】
0003668239.tif 000011
【００６２】
また、本実施形態では、式（１）の定数等を次のように設定した。
ＣR ＝０．２６８
Ｎ＝６．６
Ｔ＝２９６
ｋB ＝１．３８０６６×１０-29 
ｎ＝ＣR ／ＫB ／Ｔ＝６．５５８×１０25
ＮA ＝ｎ・Ｎ＝４．３２８×１０26
フィラーモデルの体積含有率  ３０％
フィラーモデルの縦弾性係数Ｅ：１００ＭＰａ
フィラーモデルのポアソン比ν：０．３
【００６３】
また変形シミュレーションは、巨視的モデルＭとして２mm×２mmの矩形状の解析領域を設定し、この巨視的モデルＭに一様な一軸引張変形を発生させるため、図１２のｘ２方向に平均ひずみ速度１．０×１０-5／ｓを加え、ひずみを０．０１きざみで零から２．５まで漸増させるとともに、各ひずみ毎にゴム材料モデル２に作用する応力を計算している。またひずみが２．５に達した後は、逆に０．０１きざみで前記と同じひずみ速度でひずみを零まで漸減させた。各ひずみ状態において、それぞれ平均セグメント数Ｎが計算され、この値は式（１）へ代入され逐次計算が行われる。なおゴム材料モデルは、厚さ方向に変化しないようにArrudaらの３次元８鎖モデルが用いられている。また、マトリックスモデル３及び界面モデル５の平均セグメント数Ｎは次のように設定した。
【００６４】
＜マトリックスモデル＞
・負荷変形時の平均セグメント数Ｎ
Ｎ＝-3.2368+20.6175 λc-21.8168 λc2+10.8227λc3-1.9003 λc4
・除荷変形時の平均セグメント数Ｎ（一定）
Ｎ＝６．６
・分子鎖のセグメントの総数ＮA （一定）
ＮA ＝４．３２８１×１０26
【００６５】
＜界面モデル＞
・負荷変形時の平均セグメント数Ｎ
Ｎ＝-5.9286+20.6175 λc-21.8168 λc2+10.8227λc3-1.9003 λc4
・除荷変形時の平均セグメント数Ｎ（一定）
Ｎ＝３．９１
・分子鎖のセグメントの総数ＮA （一定）
ＮA ＝３．２０３×１０25
【００６６】
前記変形計算が行われると、その結果から必要な物理量を取得することができる（ステップＳ４）。図１３には、マトリックスモデル３、フィラーモデル４及び界面モデル５について、それぞれ単独で変形シミュレーションを行った結果を示している。フィラーモデル４が最も高い弾性を示し、エネルギーロスは生じていない。界面モデル５は、同じ応力でもマトリックスモデル３よりも大きなひずみが生じる粘弾性特性が表れており、かつ、エネルギーロス（曲線のループ面積）も大きくなっている。
【００６７】
また図１４の曲線Ｌａは、上記各モデルを図１０の微視構造（セルユニット）として組み入れたゴム材料モデル全体Ｍにおける真応力とひずみとの関係を示している。負荷変形時には、非線形の第１の曲線Ｌａ１が得られた。また除荷変形時では、第１の曲線Ｌａ１とは異なる第２の曲線Ｌａ２が得られた。第２の曲線Ｌａ２は、第１の曲線Ｌａ１よりも軟化しており（低弾性）、ヒステリシスループが再現できた。このループの閉面積を計算することにより、変形１サイクルでのエネルギーロスを計算することができる。
【００６８】
ちなみに、分子鎖１本当たりの平均セグメント数Ｎを一定として Arruda らと同様の変形シミュレーションを行った場合、負荷変形時及び除荷変形時とも、同じ曲線を通り、エネルギーロスを表現し得ない結果となる。また上記実施形態では、平均セグメント数Ｎがパラメータλc の関数である態様を示したが、λc に代えてひずみ速度やひずみ量を用いることができる。この場合、式（１９）と同様に、負荷変形時と除荷変形時とで係数を違えることにより、エネルギーロスを表現することができる。
【００６９】
各曲線Ｌａ１、Ｌａ２の形状は、式（１９）において、係数ＡないしＥを適宜設定し分子鎖１本当たりの平均セグメント数Ｎを定める関数を変えることにより種々設定しうるのは前記同様である。したがって、解析しようとするゴム材料に応じて、係数ＡないしＥを種々設定することにより、材料毎により詳細にエネルギーロスなどを詳細に検討することができる。
【００７０】
また、曲線Ｌｂは、図１０の粘弾性材料モデル２において、界面モデル５にマトリックスモデル３と同一の材料特性を与えたモデル（つまり、界面を考慮していないモデル）で同一の実験を行った結果を示している。界面の軟質領域がないため、応力ーひずみ曲線がＬａよりも立ち上がっている。また、界面モデル５の大きなエネルギーロスがないため、モデル全体のエネルギーロスも小さいことが分かる。
【００７１】
さらに曲線Ｌｃは、界面モデル５に、前記とは逆の材料特性を与えてシミュレーションを行った結果を示している。即ち、曲線Ｌｃを得た粘弾性材料モデル２は、図１５に示すように、界面モデル５の方がマトリックスモデル３よりも高弾性を示すように粘弾性特性を定義したものである。なお、各モデルの平均セグメント数Ｎは次のように設定している。
【００７２】
＜マトリックスモデル＞
・負荷変形時の平均セグメント数Ｎ
Ｎ＝-3.2368+20.6175 λc-21.8168 λc2+10.8227λc3-1.9003 λc4
・除荷変形時の平均セグメント数Ｎ（一定）
Ｎ＝６．６
・分子鎖のセグメントの総数ＮA （一定）
ＮA ＝４．３２８１×１０26
【００７３】
＜界面モデル＞
・負荷変形時の平均セグメント数Ｎ
Ｎ＝-5.4800+20.6175 λc-21.8168 λc2+10.8227λc3-1.9003 λc4
・除荷変形時の平均セグメント数Ｎ（一定）
Ｎ＝４．３６
・分子鎖のセグメントの総数ＮA （一定）
ＮA ＝８．５６９×１０26
【００７４】
この結果、図１４に示したように、粘弾性材料モデル２の系全体の粘弾性特性が曲線Ｌｂよりも高い勾配、高弾性化されていることが確認できる。
【００７５】
図１６には、前記粘弾性材料モデル２（曲線Ｌａのもの）の一軸引張の変形シミュレーションにおける変形〜回復状態の進行過程を間欠的に視覚化して示す。図１６（Ａ）〜（Ｅ）は、負荷変形時を、同（Ｆ）〜（Ｊ）は除荷変形をそれぞれ示している。また色彩が白く変化しているところは、ひずみの大きい領域を示している。この結果では、フィラーモデル４の界面及びフィラーモデル４、４間に大きなひずみが集中していることが正確に再現されている。
【００７６】
図１７には、最大ひずみ時における各要素の引張方向のひずみ分布を示し、（Ａ）は、図１３の曲線Ｌｂを、（Ｂ）は同曲線Ｌａのものである。色彩が薄くなるところほど大きなひずみが生じていることを示す。界面を考慮していない（Ａ）のものでは、ひずみの発生が全体的に広範囲にかつ穏やかに生じているのに対して、界面を考慮した（Ｂ）のものでは、フィラーモデル４の周辺にひずみが集中していることが良く再現できている。
【００７７】
また図１８には、各要素毎にひずみの１サイクルのエネルギーロスを計算し、その大小を色彩にて示す。色彩が薄いところほど大きなエネルギーロスが生じていることを示す。また粘弾性材料モデル２は、図１７の最大変形形状で表しており、（Ａ）は、図１４の曲線Ｌｂを、（Ｂ）は同曲線Ｌａのものである。界面を考慮した（Ｂ）のものでは、フィラーモデル４の界面に集中して大きなエネルギーロスが生じていることが確認できる。
【図面の簡単な説明】
【００７８】
【図１】本実施形態で用いたコンピュータ装置の一例を示す斜視図である。
【図２】本実施形態の処理手順を示すフローチャートである。
【図３】粘弾性材料モデルの一実施形態を示す線図である。
【図４】（Ａ）は粘弾性材料の網目構造体を示す斜視図、（Ｂ）は８鎖モデルの一例を示す斜視図である。
【図５】（Ａ）、（Ｂ）は分子鎖の接合点の破断を説明する線図である。
【図６】パラメータλc と分子鎖１本当たりの平均セグメント数Ｎとの関係を示すグラフである。
【図７】変形シミュレーションの手順を示すフローチャートである。
【図８】粘弾性材料モデルのシミュレーション結果として真応力−ひずみとの関係を示すグラフである。
【図９】比較例の粘弾性材料モデルのシミュレーション結果として真応力−ひずみとの関係を示すグラフである。
【図１０】粘弾性材料モデルの他の実施形態を示す線図である。
【図１１】カーボンブラックの形状を示す線図である。
【図１２】均質化法を説明する微視構造と全体構造との関係を示す。
【図１３】マトリックスモデル、フィラーモデル及び界面モデルの応力ーひずみ曲線を示す。
【図１４】他の実施形態として、粘弾性材料モデルの変形シミュレーション結果から得られた真応力ーひずみ曲線である。
【図１５】マトリックスモデル及び界面モデルの応力ーひずみ曲線を示す。
【図１６】（Ａ）〜（Ｊ）は、粘弾性材料モデルの変形過程を視覚化して示す線図である。
【図１７】（Ａ）は界面層未考慮のモデル、（Ｂ）は、界面層を考慮したモデルの引張ひずみの分布図である。
【図１８】（Ａ）は界面層未考慮のモデル、（Ｂ）は、界面層を考慮したモデルのエネルギーロスの分布図である。
【図１９】（Ａ）は粘弾性材料、（Ｂ）はその分子鎖１構造を説明する線図、（Ｃ）は１本の分子鎖の拡大図、（Ｄ）はセグメントの拡大図である。
【図２０】Arrudaらによるシミュレーション結果として、応力ーひずみ曲線を示す。
【符号の説明】
【００７９】
１  コンピュータ装置
２  粘弾性材料モデル
３  マトリックスモデル
４  フィラーモデル
５  界面モデル
【図１】
0003668239.tif 000013
【図２】
0003668239.tif 000014
【図３】
0003668239.tif 000015
【図４】
0003668239.tif 000016
【図５】
0003668239.tif 000017
【図６】
0003668239.tif 000018
【図７】
0003668239.tif 000019
【図８】
0003668239.tif 000020
【図９】
0003668239.tif 000021
【図１１】
0003668239.tif 000022
【図１２】
0003668239.tif 000023
【図１３】
0003668239.tif 000024
【図１４】
0003668239.tif 000025
【図１５】
0003668239.tif 000026
【図１９】
0003668239.tif 000027
【図２０】
0003668239.tif 000028
【図１０】
0003668239.tif 000029
【図１６】
0003668239.tif 000030
【図１７】
0003668239.tif 000031
【図１８】
0003668239.tif 000032