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

Description

本発明は、マトリックスゴムとフィラーとを含むゴム材料モデルから損失正接ｔａｎδ等の物性値を短時間で精度良く予測・計算するのに役立つシミュレーション方法に関する。

従来、マトリックスゴムとフィラーとを含むゴム材料の物性値を把握するために、粘弾性試験機を用いて歪と応力との関係が測定される。粘弾性試験機では、図１２に示されるように、ゴム材料に、変形条件として一定の周波数及び振幅で変化する歪ａを与え、そのときに発生する応力ｂが測定される。一例を挙げると、ゴム材料（試験片）には、例えば１０％の初期伸長引張歪が与えられるとともに、この状態から片振幅１％かつ１０Hzの周期で正弦波のひずみを負荷したときの応力が測定される。

測定される応力は、通常、歪ａとある位相差δを持って正弦波状に変化する。また、この応力の振幅σ0と、前記位相差δとを用いて、例えば、以下の重要な物性値を調べることができる。
貯蔵弾性率Ｅ’＝｛σ０ ／ε０｝・ｃｏｓδ
損失弾性率Ｅ”＝｛σ０ ／ε０｝・ｓｉｎδ
損失正接ｔａｎδ＝Ｅ”／Ｅ’
なお、ε０ は歪ａの片振幅である。

また、近年では、マトリックスゴムとフィラーとを含むゴム材料の設計・開発のために、コンピュータを用いたシミュレーション方法が、例えば下記特許文献１乃至２のように種々提案されている。このようなコンピュータシミュレーションでは、ゴムの分子鎖及びフィラー等を計算に織り込むことができ、実際のゴム材料を試作することなく上述の物性値の計算が可能になってきた。従って、上述のようなシミュレーションを利用することにより、例えばカーボンやシリカ等のフィラーの割合を異ならせた種々のゴム材料について、その物性値等がどのように変化するか等を、実際にゴム材料を試作することなく予測することができる。

特開２００６−１７５９３７号公報 特開２００９−２５９０４３号公報

ところで、上述のシミュレーションでは、ゴム材料モデルの変形条件として、図１２に示したような正弦波状に変化する歪ａを連続的に与えることができない。このため、シミュレーションでは、図１１（ａ）に示されるように、歪が、一定の時間増分Δｔの間隔で離散的にかつ正弦波状に付与される。そして、この離散的に与えられる各歪に応じた応力が収束計算等によりそれぞれ求められる。応力の応答は、例えば図１１（ｂ）に示される。また、応力の変化を調べるために、離散的に得られた応力の各計算点に正弦波でカーブフィットを施して近似波形を求め、この近似波形に基づいて、例えば応力の振幅σ0等が求められていた。

しかしながら、上述のようなシミュレーション方法では、応力の波形を精度良く得るためには、図１１（ｂ）に示されるように１周期中に応力の計算数を比較的密に行う必要があり、多くの計算時間が必要であった。また、応力と歪との位相差が安定するためには、複数の周期の計算が必要であるため、やはり多くの計算時間が必要であった。

本発明は、以上のような問題点に鑑み案出なされたもので、各要素の歪変化量に基づくエネルギーロスに着目し、これを利用して、損失正接ｔａｎδを短時間でかつ精度良く計算しうるゴム材料のシミュレーション方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、マトリックスゴムとフィラーとを含むゴム材料の損失正接ｔａｎδをコンピュータを用いて計算するゴム材料のシミュレーション方法であって、前記コンピュータに、前記マトリックスゴム及び前記フィラーをそれぞれ有限個の要素で分割した得たゴム材料モデルを入力するステップと、前記コンピュータが、入力された前記ゴム材料モデルの変形条件として、大小２種類の歪を付与するステップと、前記ゴム材料モデルの大小２種類の各歪状態の変形シミュレーションを行い、各要素について、それぞれの状態での垂直歪及びせん断歪を計算し、前記計算された歪を用いて各要素のエネルギーロスを計算するステップと、前記各要素のエネルギーロスの総和、前記大小２種類の歪の差から得られる歪振幅、及びゴム材料モデルの弾性率に基づいてゴム材料モデルの損失正接ｔａｎδを計算するステップとを含むことを特徴とする。

また請求項２記載の発明は、前記コンピュータは、前記各要素の単位体積当たりのエネルギーロスを下記式（１）で計算する請求項１記載のゴム材料のシミュレーション方法である。
Ｗ＝Σπ・Ｅ・（εp ／２）² ・ｔａｎδ …（１）
ここで、Σは、垂直歪及びせん断歪の各成分の全ての和を示し、Ｅは各要素の貯蔵弾性率、εpは大小２種類の歪が作用したときの要素の歪変化量、ｔａｎδは各要素の損失正接である。

また請求項３記載の発明は、前記コンピュータは、前記ゴム材料モデルの損失正接ｔａｎδを下記式（２）で計算する請求項１又は２記載のゴム材料のシミュレーション方法である。
ｔａｎδ＝Ｗａ／｛π・Ｅａ・εa²｝ …（２）
ここで、Ｗａは各要素の単位体積当たりのエネルギーロスの総和、Ｅａはゴム材料モデルの弾性率、εaは大小２種類の歪の差から得られる歪振幅である。

本発明では、マトリックスゴムとフィラーとを含むゴム材料の損失正接ｔａｎδをコンピュータを用いて計算するゴム材料のシミュレーション方法において、コンピュータに、マトリックスゴム及びフィラーをそれぞれ有限個の要素で分割した得たゴム材料モデルを入力するステップと、前記コンピュータが、入力された前記ゴム材料モデルの変形条件として、大小２種類の歪を付与するステップと、前記ゴム材料モデルの大小２種類の各歪状態において、それぞれ垂直歪とせん断歪の各成分を計算し、前記計算された歪を用いて各要素の単位体積当たりのエネルギーロスを計算するステップと、前記各要素のエネルギーロスの総和、前記大小２種類の歪の差から得られる歪振幅、及びゴム材料モデルの弾性率とに基づいてゴム材料モデルの損失正接ｔａｎδを計算するステップとを含むことを特徴とする。従って、本発明によれば、大小２つの歪の状態での各要素のエネルギーロスから、ゴム材料モデル全体のエネルギーロスを計算することができる。従って、より少ない計算量でゴム材料モデルの損失正接ｔａｎδを計算することができる。

本実施形態で用いられるコンピュータ装置の一例を示す斜視図である。 本実施形態の処理手順を示すフローチャートである。 （ａ）はゴム材料モデル（微視構造）の一実施形態を示す線図、（ｂ）はそのＸ部拡大図である。 周期的な歪を説明するグラフである。 本実施形態で与えられる歪を示すグラフである。 変形シミュレーションの処理手順を示すフローチャートである。 均質化法を説明する概念図である。 （ａ）乃至（ｃ）は、要素の垂直歪及びせん断歪を説明する線図である。 ｔａｎδの計算手順を示すフローチャートである。 本発明と従来手法との計算結果を比較するグラフである。 （ａ）は従来のシミュレーションにおける歪と時間との関係を示すグラフ、（ｂ）は、その条件で計算された応力と時間との関係を示すグラフである。 ゴム材料に与えられる正弦波状の歪と、これに対応する応力との関係を示すグラフである。

以下、本発明の実施の一形態が図面に基づき説明する。
本発明では、マトリックスゴムとフィラーとを含むゴム材料の損失正接ｔａｎδをコンピュータを用いて計算するゴム材料のシミュレーション方法である。図１には、このようなシミュレーション方法を実施するためのコンピュータ装置１が示されている。該コンピュータ装置１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含んで構成されている。本体１ａの内部には、ＣＰＵ、ＲＯＭ、作業用メモリー及び磁気ディスク等の大容量記憶装置が設けられる。また、本体１ａには、ＣＤ−ＲＯＭやフレキシブルディスクのドライブ装置１ａ１、１ａ２が設けられる。そして、前記大容量記憶装置には後述する本発明のシミュレーション方法を実行するための処理手順（プログラム）が記憶されている。

図２には、本実施形態のシミュレーション方法の処理手順の一例が示される。本実施形態では、先ず、ゴム材料モデルが設定される（ステップＳ１）。本実施形態のゴム材料モデルは、例えばタイヤ用のゴム材料を模したものであって、このようなゴム材料は、例えば初期歪１０％、周波数１０Hz、歪片振幅１％、温度３０℃の測定条件で、損失正接ｔａｎδが０．１〜０．４程度の値を示すものが一般的である。

図３（ａ）及びそのＸ部拡大図である（ｂ）に視覚化して示されるように、ゴム材料モデル２は、例えばカーボンやシリカといったフィラーがマトリックスゴム中に分散して配されたフィラー入りのゴム材料がモデル化される。また、図３（ａ）のモデルは、繰り返し最小単位の微視構造としてのゴム材料モデル２の一例が視覚化して示されており、この実施形態において、微視構造としてのゴム材料モデル２は、例えば３００nm×３００nmの正方形である。

該ゴム材料モデル２は、解析しようとするフィラー配合のゴム材料の微小領域が、有限個の小さな要素（メッシュ）２ａ、２ｂ、２ｃ…に置き換えられたものである。各要素２ａ、２ｂ、２ｃ…は、数値解析が可能に定義される。

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

本実施形態のゴム材料モデル２は、後述する変形シミュレーションにおいて平面ひずみ状態の解析、より詳しくはｙ軸方向の引張変形シミュレーションが行われる。つまり、本実施形態では、ｚ方向にはひずみを持たない２次元のシミュレーションが行われる。

また、本実施形態のゴム材料モデル２は、ゴムマトリックス部分がモデル化されたマトリックスモデル３と、このマトリックスモデル３の中に配置されかつ前記フィラーがモデル化されたフィラーモデル４と、各フィラーモデル４の周りを環状に取り囲みかつ前記マトリックスモデル３よりも硬い物性が定義された界面モデル５（図３（ｂ）に示す）とを含む。

前記マトリックスモデル３は、ゴム材料モデル２の主要部を構成し、かつ、例えば三角形ないし四辺形の複数個の要素を用いて表現されている。変形計算を行うために、マトリックスモデル３を構成する各要素には、例えば、その物性として応力と伸びとの関係を表す関数が定義される。本実施形態のゴム材料のシミュレーション方法では、ゴム弾性応答を表現するために、前記マトリックスモデル３及び界面モデル５のゴム部分は、いずれも分子鎖網目理論に基づいて計算が行われる。なお、分子鎖網目理論については、本件出願人らが提案している特開２０１０−２０５１６５号公報や特開２００９−２１６６１２号公報に記載されており、これらの文献に記載の方法で本シミュレーションの計算過程に導入されている。従って、分子鎖網目理論については、ここでの詳細な説明は省略する。

前記フィラーモデル４は、シリカからなるフィラーが四辺形の複数個の要素を用いてモデル化されたもので、全体として円形に形成されている。三次元モデルの場合、フィラーモデル４は、球形にモデル化されるのが望ましい。また、シリカは、直径約１０〜３００nm程度であり、ゴムに比べて非常に硬い粒子からなる。このため、フィラーモデル４には、このような解析対象となるシリカの物性とほぼ等しい物性が設定される。即ち、本実施形態において、フィラーモデル４は、粘弾性体ではなく弾性体として取り扱われる。また、フィラーモデル４の粒子の個数は、例えば、解析対象のゴム材料のシリカ配合量などに基づいて適宜決定される。

前記界面モデル５は、シリカからなるフィラーとマトリックスゴムとを化学的に結合させるシランカップリング剤の働きをシミュレーションに取り込むためにモデル化されたものである。本実施形態の界面モデル５は、フィラーモデル４の周りを小さい厚さｔで環状に連続して取り囲むように設定されている。従って、界面モデル５の内周面はフィラーモデル４の外周面に接触しかつ固着されている。本実施形態では、界面モデル５の内周面とフィラーモデル４の外周面とは、互いに剥離しない条件が設定されるが、必要に応じて、予め定めた値以上の応力が生じたときに、フィラーモデル４と界面モデル５との境界を分離させるような条件が設定されても良い。なお、界面モデル５の外周面は、マトリックスゴムモデル３に接触しかつ固着されている。

前記界面モデル５の厚さｔは、特に限定されるものではないが、種々の実験結果等に鑑み、フィラーモデル４の直径の１０〜３０％程度、より好ましくは１５〜２５％程度に設定されるのが実際のゴム材料と整合する点で望ましい。

また、前記界面モデル５にも、応力と伸びとの関係が定義される。実際のシリカの界面結合材の物性に鑑み、この界面モデル５はマトリックスゴムよりも硬い物性、即ち、マトリックスモデル３よりも伸び難く定義される。ただし、界面モデル５は、フィラーモデル４よりは軟らかいのは言うまでもない。

また、マトリックスモデル３、フィラーモデル４及び界面モデル５には、それぞれ、これまでの実験結果等に基づいて、物性についてのパラメータ（例えば、密度、弾性率、損失正接ｔａｎδ等）が入力される。

次に、上記マトリックスモデル３、フィラーモデル４及び界面モデル５の複合体であるゴム材料モデル２を変形させるための変形条件が設定される（ステップＳ２）。従来、ゴム材料からｔａｎδを実測する場合、ゴム材料モデル２に図４に示されるような周期的な歪が与えられる。また、従来、シミュレーションでゴム材料モデル２のｔａｎδを計算する場合、図１１（ａ）に示したように、擬似的に、周期的な歪がゴム材料モデルに離散的に与えられる。本発明では、前記ゴム材料モデル２からｔａｎδを計算するためのゴム材料モデルの変形条件として、大小２種類の歪が与えられる。

設定される上記大小２種類の歪の値は、特に限定されるものではないが、物性値として本来得たいｔａｎδの測定条件に基づいて決定することが望ましい。例えば、図５に仮想線で示されるように、本来得たいｔａｎδの測定条件が、初期歪１０％、振幅歪±１％の振幅歪であった場合、シミュレーションにおいては、前記振幅歪の最大ピーク歪（εmax＝１１％）と、最小ピーク歪（εmin＝９％）との値が、上記大小２種類の歪の値として設定されるのが好ましい。これにより、上記測定条件により近い変形条件でシミュレーションを行うことが可能になる。

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

コンピュータ装置１では、入力されたデータに基づいて各要素の剛性マトリックスを作成し（ステップＳ３２）、しかる後、全体構造の剛性マトリックスを組み立てる（ステップＳ３３）。全体構造の剛性マトリックスには、既知節点の変位、節点力が導入され（ステップＳ３４）、剛性方程式の解析が行われる。そして、未知節点変位が決定され（ステップＳ３５）、先に定められた変形条件（各ひずみ）に基づいた応力を計算する（ステップＳ３６）。この際、コンピュータ装置１は、先に述べたように、ゴム材料モデル２に、変形条件として２種類の歪の一つを与え、計算を行う。そして、コンピュータ装置１は、それらの結果を記憶装置等に出力する（ステップ３７）。

ステップＳ３８では、計算を終了させるか否かの判定がなされる。設定された２種類の歪について全て計算が終えていると判断された場合、図１のステップＳ４に戻る。他方、設定された２種類の歪について全て計算が終えていないと判断された場合には、新たなあ歪の条件でステップＳ３２以降が繰り返される。このようなシミュレーション（変形計算）は、例えば有限要素法を用いたエンジニアリング系の解析アプリケーションソフトウエア（例えば米国リバモア・ソフトウェア・テクノロジー社で開発・改良されたＬＳ−ＤＹＮＡ等）を用いて行うことができる。

また、本シミュレーションは、均質化法（漸近展開均質化法）に基づいて行われる。均質化法は、図６に示されるように、図３に示した微視構造（均質化法では「ユニットセル」とも呼ばれる）を周期的に持っているゴム材料全体Ｍを表現するｘI と、前記微視構造を表現するｙI との独立した２変数が用いられる。微視的スケールと巨視的スケールという異なる尺度の場におけるそれぞれ独立した変数を漸近展開することにより、図３に示した微視構造のモデル構造を反映させたゴム材料全体の平均的な力学応答を近似的に求めることができる。

前記変形計算が行われると、コンピュータ装置１は、その結果から必要な物理量を取得する（メモリー等に読む込み）ことができる（ステップＳ４）。前記物理量としては、前記大小２種類の各歪条件における各要素の垂直歪及びせん断歪が含まれる。図８には、二次元のソリッド要素からなる要素２ａを例に挙げ、（ａ）にはその垂直歪が示される。該垂直歪は、ｘ軸方向の垂直歪ε1とｙ軸方向の垂直歪ε2とを含む。三次元の場合、垂直歪には、さらにｚ軸方向の歪が含まれる。また、図８（ｂ）及び（ｃ）には、せん断歪が示され、これには、ｘ軸方向のせん断歪ε3とｙ軸方向のせん断歪ε4とが含まれる。三次元の場合、せん断歪には、さらにｚ軸方向の歪が含まれる。

次に、ゴム材料モデル２のｔａｎδが計算される（ステップＳ５）。ｔａｎδの計算の一例は、図９に示されている。本実施形態において、コンピュータ装置１は、各要素２ａ、２ｂ、２ｃ…の単位体積当たりのエネルギーロスＷｉ（ｉ＝１〜要素数）を下記式（１）で計算する（ステップＳ５１）。
Ｗｉ＝Σπ・Ｅ・（εp ／２）²・ｔａｎδi …（１）
式（１）は、歪の変化に着目して実験によって得られた公知の理論式であり、式中、”Σ”は要素の垂直歪及びせん断歪の各成分の全ての和、”Ｅ”は要素の貯蔵弾性率、εpは大小２種類の歪（εmax、εmin）が作用したときの要素の歪の変化量、ｔａｎδiは要素の損失正接をそれぞれ表している。

次に、コンピュータ装置１は、各要素のエネルギーロスＷｉを計算し終えると、これらの総和Ｗａを計算する（ステップＳ５２）。

次に、コンピュータ装置１は、ゴム材料モデル２の損失正接ｔａｎδを下記式（２）で計算する（ステップＳ５３）。
ｔａｎδ＝Ｗａ／｛π・Ｅａ・εa²｝ …（２）
式（２）は、上記式（１）の理論式を、ゴム材料モデルの系全体に適用したもので、Ｗａは前記各要素の単位体積当たりのエネルギーロスＷｉの総和、Ｅａはゴム材料モデルの弾性率、εaは大小２種類の歪の差から得られる歪振幅｛（εmax−εmin）／２｝である。なお、ゴム材料モデルの弾性率Ｅａは、その応力σと歪εとの関係から、下式によって求める。
Ｅａ＝（σmax−σmin）／（εmax−εmin）
上記σは、全要素の応力と要素の体積から求めた平均の応力で、σmax−σminは、その最大と最小の差を意味している。

このように、本発明によれば、大小２つの歪の状態での各要素のエネルギーロスから、ゴム材料モデル全体のエネルギーロスを計算し、さらにそこからｔａｎδを計算することができる。従って、より少ない計算量でゴム材料モデルのｔａｎδを計算することができる。

図１０には、図３に示したゴム材料モデル２を使用して、本発明のシミュレーション方法で求めたｔａｎδと従来手法で求めたｔａｎδとを比較したグラフが示される。これらのシミュレーションの各条件は次の通りである。
ゴム材料モデルの初期伸長：１０％
周波数：１０Hz
ひずみの片振幅：１％

従来手法では、ゴム材料モデルの変形条件として、図１１（ａ）で示したように、短い時間増分Δｔ１で徐々に変化する歪が複数回離散的に与えられた。なお、従来手法では、計算される応力が、歪に対する位相差が安定するまで繰り返し計算された。本実施例では、この計算に７周期を要した。

他方、本発明では、ゴム材料モデルの変形条件として、図５に示したように、前記振幅歪の最大ピーク歪（εmax＝１１％）と、最小ピーク歪（εmin＝９％）との２種類だけが与えられた。

図１０から明らかなように、本発明と従来手法とは、ｔａｎδに関して、正の相関を有していることが確認できる。また、ｔａｎδの計算に要した時間は、従来手法の役１／１５であり、大幅な短縮が確認できた。従って、本発明によれば、より少ない計算時間で精度良くゴム材料モデルのｔａｎδを計算・予測することができる。

１ コンピュータ装置
２ ゴム材料モデル

Claims (3)

1. マトリックスゴムとフィラーとを含むゴム材料の損失正接ｔａｎδをコンピュータを用いて計算するゴム材料のシミュレーション方法であって、
   前記コンピュータに、前記マトリックスゴム及び前記フィラーをそれぞれ有限個の要素で分割した得たゴム材料モデルを入力するステップと、
   前記コンピュータが、入力された前記ゴム材料モデルの変形条件として、大小２種類の歪を付与するステップと、
   前記ゴム材料モデルの大小２種類の各歪状態の変形シミュレーションを行い、各要素について、それぞれの状態での垂直歪及びせん断歪を計算し、前記計算された歪を用いて各要素のエネルギーロスを計算するステップと、
   前記各要素のエネルギーロスの総和、前記大小２種類の歪の差から得られる歪振幅、及びゴム材料モデルの弾性率に基づいてゴム材料モデルの損失正接ｔａｎδを計算するステップとを含むことを特徴とするゴム材料のシミュレーション方法。
2. 前記コンピュータは、前記各要素の単位体積当たりのエネルギーロスを下記式（１）で計算する請求項１記載のゴム材料のシミュレーション方法。
   Ｗ＝Σπ・Ｅ・（εp ／２）² ・ｔａｎδ …（１）
   （ここで、Σは、垂直歪及びせん断歪の各成分の全ての和を示し、Ｅは各要素の貯蔵弾性率、εpは大小２種類の歪が作用したときの要素の歪変化量、ｔａｎδは各要素の損失正接である。）
3. 前記コンピュータは、前記ゴム材料モデルの損失正接ｔａｎδを下記式（２）で計算する請求項１又は２記載のゴム材料のシミュレーション方法。
   ｔａｎδ＝Ｗａ／｛π・Ｅａ・εa²｝ …（２）
   （ここで、Ｗａは各要素の単位体積当たりのエネルギーロスの総和、Ｅａはゴム材料モデルの弾性率、εaは大小２種類の歪の差から得られる歪振幅である。）

Images