JP2005146146A

JP2005146146A - ゴム材料のシミュレーション方法

Info

Publication number: JP2005146146A
Application number: JP2003386992A
Authority: JP
Inventors: Masato Naito; 正登内藤; Masaki Shiraishi; 正貴白石
Original assignee: Sumitomo Rubber Industries Ltd
Current assignee: Sumitomo Rubber Industries Ltd
Priority date: 2003-11-17
Filing date: 2003-11-17
Publication date: 2005-06-09
Anticipated expiration: 2023-11-17
Also published as: JP3660932B2

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 Mechani cs 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としている。

この実施形態のゴム材料モデル２は、ゴムマトリックス部分がモデル化されたマトリックスモデル３と、距離を隔てて向き合う少なくとも２つのフィラー粒子をモデル化したフィラーモデル４Ａ、４Ｂ（総称するとき、単に「フィラーモデル４」と言うことがある。）と、前記２つのフィラーモデル４Ａ、４Ｂ間に設けられ、かつ、該フィラーモデル４Ａ、４Ｂ間の距離に基づくフィラー間引力が定義されたフィラー間モデル５とからなるものが例示される。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

フィラーが配合されたゴム材料の変形シミュレーションにおいてフィラー間の引力を考慮する場合、カーボンブラックの一次粒子又はその凝集体である二次粒子のいずれをモデル化しても良い。この実施形態では、フィラーモデル４は、ラグビーボール状に定義された２つが示される。各フィラーモデル４Ａ、４Ｂは同一の構成で定義され、距離ｄ（ただし、ｄ≠０）を隔てて向き合っている。

例えば下記に記すneedleman らの実験を参考にすれば、２物体間の距離が大きすぎてもまた小さすぎても、それらの間の引力が作用しなくなることが報告されている。
A Continuum Model for Void Nucleation by Inclusion Debonding / A. needleman (Journa l of Applied Mechanics SEP. 1987 Vol.54 P.525-531)
このシミュレーションにおいて、フィラーモデル４Ａとフィラーモデル４Ｂとの間の初期の距離ｄの値は最も小さい所で数ナノメータとし、より具体的には１ないし３ナノメータ程度の範囲で選択している。ただし、初期の距離ｄの値は、フィラーモデル４の位置やシミュレーションの条件などに応じて種々異なる。

また図９には、上記needleman らが提案した剥離を開始した界面ａ、ｂ間に働く相互引力に基づく引力Ｔｎと、距離ｄとの関係が示されている。なお横軸のδは、界面ａ、ｂ間の距離ｄを増大させたときに引力が零となるときの限界の長さ（以下、この長さを「特徴距離」と呼ぶ）である。図から明らかなように、界面ａ、ｂ間の引力Ｔｎは、距離ｄが増大するにつれて徐々に増大してピークσmax に至り、さらに距離ｄが増大することによって滑らかに減少し、かつ特徴距離で零に至る放物状曲線を描く距離ｄの関数である。needleman らは、引力Ｔｎの関数を式（２０）のように定めている（二次元形）。法線方向ｎと接線方向ｓとは、図１０のように定める。

ゴム材料のシミュレーションの場合、フィラーモデル４Ａ、４Ｂ間に初期状態でマトリックスモデルが介在していること、及びこの初期状態ではフィラー間引力は零であるとの前提を置く。このため、式（２０）の関数は、本実施形態にそのまま適用することはできない。この実験結果をフィラーモデル４Ａ、４Ｂ間に適用するために、次の修正が加えられる。

先ず、初期のフィラーモデル間の距離を導入するため、図９のグラフを横軸方向に−ｄ（初期のフィラーモデル間の距離）、縦軸方向に−Ｔｄ（引力）だけ移動させる。この場合、式（２０）は、式（２１）のように修正できる。本シミュレーションでは、フィラーモデル間の距離ｄに応じたフィラー間引力を式（２１）から求め得る。なおこの例では二次元形を示すが、３次元形としても良い。

図１１には、フィラーモデル間の初期の距離ｄ０を種々違えた場合の引力Ｔｎの一例を示している。前記特徴距離δは、例えば初期の距離ｄ０＝０の場合には約１１nm、また初期の距離ｄ０＝２．０nmの場合には約４．８nmになる。このように、初期の距離ｄ０、特徴長さδは、条件に応じて適宜設定される。

また図１２には、フィラー間引力Ｔｎの例として、初期の距離ｄ０＝１．０nmの場合を取り出して示す。フィラー間引力Ｔｎを示す前記式（２１）は、フィラーモデル間の負荷変形（フィラーモデル４Ａ、４Ｂ間の距離を増大させる変形）時においては、距離ｄが増大するにつれて、ｄ＝１．０nmの点Ｐ１から徐々に増大し、かつピークσmax に至る漸増領域Ｌａと、さらなる距離ｄの増大によって滑らかに減少し、かつ、予め定めた特徴距離δの点Ｐ２で零に至る漸減領域Ｌｂとを含む放物状曲線を描く関数である。距離ｄが特徴距離δよりも大きくなるとフィラーモデル４Ａ、４Ｂ間にはフィラー間引力は作用しない。

また前記引力Ｔｎは、負荷変形時の後に続く除荷変形時においては、除荷開始時のフィラーモデル４Ａ、４Ｂ間の距離ｄのときの値から初期の状態まで線形に減少する。例えば漸減領域Ｌｂ上における任意の点Ｐ３から除荷が開始される場合、図の鎖線の矢印のように、フィラー間引力Ｔｎは、距離ｄに応じて、当該点Ｐ３と点Ｐ１とを結ぶ直線上の値をとりながら線形に減少する。なお、負荷変形がピークσmax を超えていない場合には、負荷変形と同じルートを通って初期の状態まで復元できる。この場合、可逆変形となる。これにより、フィラー間引力は、ピークσmax を超えた場合、荷重の負荷及び除荷の１サイクルで閉ループを描く非可逆変形をマトリックスモデル３に生じさせる。この閉ループの面積は熱として消費されるエネルギーロスを表す。このような除荷時の引力Ｔｎは、下記式（２２）により求めることができる。

また本発明では、このようなフィラー間引力を生じさせるために、マトリックスモデル３とは別のフィラー間モデル５がゴム材料モデル２の中に定義される。図４には、やや太い実線でこのフィラー間モデル５が拡大表示される。フィラー間モデル５は、２つのフィラーモデル４Ａ、４Ｂの間に設けられる。各フィラー間モデル５は、この例ではいずれも四辺形である。四辺形をなす各フィラー間モデル５において、２つの節点は、一方のフィラーモデル４Ａの外表面の節点と共有され、残り２つの節点は、他方のフィラーモデル４Ａの外表面の節点と共有されている。計算効率を向上させるために、このフィラー間モデル５を含めて、ゴム材料モデル２の微視構造（ユニットセル）は、そのｘ軸及びｙ軸の各中心線に対して対称にモデリングされている。

またフィラー間モデル５は、マトリックスモデル３のようなそれ自身の剛性を持っていない。換言すれば、式（２１）、同（２２）を用いてフィラー間引力を計算するために設けられたフィラーモデル間の微細領域である。フィラー間モデル５は、フィラーモデル４Ａ、４Ｂの間の距離の方向（ｙ軸方向）には分割されていない。従って、フィラー間モデル５の形状は、変形の過程度、その都度、４つの節点（換言すればフィラーモデルの相対位置）よって単純に決定される。また式（２１）、同（２２）によってフィラーモデルの単位面積当たりのフィラー間引力（相互作用力）が各フィラー間モデル５において計算されるが、それは図１０に示した４つの節点に作用する力へと置き換えられる。これにより、フィラー間引力をマトリックスモデル３に作用させることができる。

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

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

これらの処理を行うための有限要素方程式は式（２３）の通りである。

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

本シミュレーションは、均質化法（漸近展開均質化法）に基づいて行われる態様を例示する。均質化法は、図１４に示されるように、図３に示したゴム材料モデル２の微視構造（ユニットセル）を周期的に持っているゴム材料全体Ｍを表現するｘ_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＝ｎ・Ｎ＝４．３２８×１０²⁶
フィラーモデルの縦弾性係数Ｅ：１００ＭＰａ
フィラーモデルのポアソン比：０．３
特徴長さδ（δｎ＝δｓ）＝２２．７５nm
ピークσmax ＝０．００１ＭＰａ
法線方向とせん断方向との強さ比γ＝１．０

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

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

前記変形計算が行われると、その結果から必要な物理量を取得することができる（ステップＳ４）。図１５には、マトリックスモデル３とフィラーモデル４とについて、それぞれ単独で変形シミュレーションを行った結果を示している。

また図１６には、図３の微視構造（セルユニット）として組み入れたゴム材料モデル全体における真応力とひずみとの関係を示している。実線のグラフは、フィラー間モデル５を定義していないため、フィラー間引力の相互作用が考慮されていない結果であり、破線のグラフは、フィラー間モデル５を定義した本実施形態の結果を示す。フィラー間引力を考慮した本実施形態のシミュレーションでは、フィラー間引力を無視したものに比べて、やや硬化していることが分かる。また、フィラー間引力を考慮したものでは、エネルギーロスを表す閉ループ面積が、フィラー間引力を無視したものに対して０．９２９％アップしており、僅かではあるがフィラー間引力に伴うエネルギーロス分大きくなっていることが確認できた。この値はパラメータの設定によって変えることができる。また、この実施形態では、セルユニットの中に２つのフィラーモデル４Ａ、４Ｂを設定したが、例えば解析しようとするゴム材料のフィラーの配合割合などに応じてフィラーモデル４の密度を変えることにより、より正確なエネルギーロスの評価ないし比較を行うことができる。

また図１７には、各フィラー間モデル５に生じた引力と距離ｄとの関係が示される。要素ア、イ、ウ、エ及びオは、それぞれ図４の符号ア、イ、ウ、エ及びオが示す要素に対応している。要素ア、イ及びウについては、距離ｄの変化が小さかったためピークσmax を超えていない。このため、引力Ｔｎは、負荷及び除荷変形において同一の経路を通ったものと考えられる。また要素オについては、初期の距離ｄが当初から大きいため、フィラー間引力が生じなかったと考えられる。これに対して要素エは、引力Ｔｎがピークを超えて復元したため、負荷と除荷とで異なる経路を通ったと考えられる。

このようなシミュレーションは、例えばゴム材へのフィラーの配合割合、フィラーの種類を種々違えたゴム材について、エネルギーロスを含めた粘弾性特性をより詳しく評価することができる。これは、とりわけ粘弾性材料を主要部として用いるタイヤ、ゴルフボールの工業製品に性能改善に大いに役立つ。

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

符号の説明

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

Claims

フィラーが配合されたゴム材料の変形をシミュレーションするゴム材料のシミュレーション方法であって、
前記ゴム材料を数値解析が可能な要素でモデル化したゴム材料モデルを設定するステップと、
前記ゴム材料モデルに条件を設定して変形計算を行うステップと、
前記変形計算から必要な物理量を取得するステップとを含むとともに、
前記ゴム材料モデルは、ゴムマトリックスをモデル化したマトリックスモデルと、
距離を隔てて向き合う少なくとも２つのフィラー粒子をモデル化したフィラーモデルと、
前記２つのフィラーモデル間に設けられ、かつ、該フィラーモデル間の距離に基づくフィラー間
引力が定義されたフィラー間モデルとを含むことを特徴とするゴム材料のシミュレーション方法。
前記フィラーモデルは、カーボンブラックの一次粒子又はその凝集体である二次粒子をモデル化したカーボンブラックモデルであることを特徴とする請求項１記載のゴム材料のシミュレーション方法。
前記フィラー間引力は、フィラーモデル間の領域の負荷変形時において、フィラーモデルの間の距離が増大するにつれて徐々に増大しかつピークに至る漸増領域と、
さらなる距離の増大によって滑らかに減少し、かつ、予め定めた特徴距離で零に至る漸減領域とを含む放物状曲線を描く関数に基づき定義されることを特徴とする請求項１又は２に記載のゴム材料のシミュレーション方法。
前記引力は、負荷変形時の後に続く除荷変形時において、除荷開始時から初期の状態まで線形に減少することにより、荷重の負荷及び除荷の１サイクルでループを描くことを特徴とする請求項３に記載のゴム材料のシミュレーション方法。