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

Description

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

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

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

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

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

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

しかしながら、ゴム材料の重要な物性値である前記貯蔵弾性率Ｅ’、損失弾性率Ｅ”及び／又は損失正接ｔａｎδを計算するためには、応力の振幅σ0と位相差δとが分かれば良い。従って、このようなシミュレーションを短時間で終わらせるためには、応力の前記近似波形の全て求めるのは得策ではなく、応力のピーク近傍領域のみの必要最低限の計算結果を得ることがより能率的である。

本発明は、以上のような問題点に鑑み案出なされたもので、１計算ステップ当たりの時間増分Δｔを変化させることを基本として、貯蔵弾性率Ｅ’、損失弾性率Ｅ”及び／又は損失正接ｔａｎδといった物性値を短時間でかつ精度良く計算しうるゴム材料のシミュレーション方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、ゴム材料に正弦波状に変化する歪を与えたときの応力の変化をコンピュータを用いて計算するゴム材料のシミュレーション方法であって、前記コンピュータに、ゴム材料を有限個の要素で分割したゴム材料モデルを入力するステップと、前記コンピュータが、入力された前記ゴム材料モデルの変形条件として、正弦波状に変化する歪を、１計算ステップあたりの時間増分の間隔で離散的に付与する歪付与ステップと、前記ゴム材料モデルと前記歪の値とに基づいて、各歪状態での前記ゴム材料モデルに生じる応力を計算する応力計算ステップとを行うとともに、前記歪付与ステップは、前記時間増分が変化するとともに、前記時間増分は、前記応力のピーク近傍で小さく、かつ、前記ピーク近傍以外で大きいことを特徴とする。

また請求項２記載の発明は、前記時間増分は、前記歪のピークから０．３８ラジアンの領域で小さく、それ以外の領域で大きい請求項１に記載のゴム材料のシミュレーション方法である。

また請求項３記載の発明は、前記コンピュータは、応力の振幅と、応力と前記歪との位相差を得るステップと、これらの値を用いて貯蔵弾性率Ｅ’、損失弾性率Ｅ”及び／又は損失正接ｔａｎδを計算するステップとをさらに含む請求項１又は２に記載のゴム材料のシミュレーション方法である。

本発明のゴム材料のシミュレーション方法によれば、コンピュータは、ゴム材料モデルの変形条件として、正弦波状に変化する歪を、１計算ステップ当たりの時間増分の間隔で離散的に付与する歪付与ステップと、前記ゴム材料モデルと前記歪の値とに基づいて、前記ゴム材料モデルに生じる応力を計算する応力計算ステップとを行う。そして、前記歪付与ステップでは、時間増分が変化することを特徴としている。

従って、例えば、時間増分を、応力のピーク近傍において小さくすることにより、計算精度が要求されるピーク近傍での応力の正弦波状の計算結果を多数得ることができる。他方、ピーク近傍以外では、時間増分を前記ピーク近傍よりも、より大きく設定することにより、計算回数を減らすことができる。従って、より少ない計算量で精度良くゴム材料モデルの応力の応答等を調べることができる。

本実施形態で用いたコンピュータ装置の一例を示す斜視図である。 本実施形態の処理手順を示すフローチャートである。 （ａ）はゴム材料モデル（微視構造）の一実施形態を示す線図、（ｂ）はそのＸ部拡大図である。 歪と応力との関係を示すグラフである。 歪付与ステップで与えられる歪と時間増分Δｔとの関係を示すグラフである。 他の実施形態の歪付与ステップで与えられる歪と時間増分Δｔとの関係を示すグラフである。 変形シミュレーションの処理手順を示すフローチャートである。 均質化法を説明する概念図である。 （ａ）は本実施形態の歪付与ステップを示す歪と時間との関係を示すグラフ、（ｂ）は、それに基づいて計算された応力と時間との関係を示すグラフである。 （ａ）は従来の歪付与ステップを示す歪と時間との関係を示すグラフ、（ｂ）は、その条件に基づいて計算された応力と時間との関係を示すグラフである。 ゴム材料の正弦波状の歪と、これに対応した応力との関係を示すグラフである。

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

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

図３（ａ）及びそのＸ部拡大図である（ｂ）に視覚化して示されるように、ゴム材料モデル２は、例えばシリカ等のフィラーがマトリックスゴム中に分散されたフィラー入りのゴム材料がモデル化されている。また、図３（ａ）のモデルは、繰り返し最小単位の微視構造としてのゴム材料モデル２の一例が視覚化して示されている。該ゴム材料モデル２は、解析しようとするフィラー配合ゴム材料の微小領域が、有限個の小さな要素（メッシュ）２ａ、２ｂ、２ｃ…に置き換えられたものである。各要素２ａ、２ｂ、２ｃ…は、数値解析が可能に定義される。

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

この実施形態のゴム材料モデル２は、後述する変形シミュレーションにおいて平面ひずみ状態の解析、さらに詳しくはｙ軸方向の引張変形シミュレーションが行われる。つまり、本実施形態では、ｚ方向にはひずみを持たない２次元のシミュレーションが行われる。この実施形態において、微視構造としてのゴム材料モデル２は、例えば３００nm×３００nmの正方形である。

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

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

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

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

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

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

また、マトリックスモデル３、シリカモデル４及び界面モデル５には、それぞれ、これまでの実験結果等に基づいて、物性についてのパラメータ（密度、弾性率等）が入力される。

次に、ゴム材料モデル２を変形させるための変形条件が設定される（ステップＳ２）。本実施形態では、図３のｙ方向に正弦波状の歪を加えてゴム材料モデル２に歪が正弦波状に変化する動的な引張変形を与える条件が定義される。また、正弦波状に変化する歪は、１計算ステップあたりの時間増分Δｔの間隔で離散的に付与されるとともに、該時間増分Δｔが変化することを特徴としている。

例えば、上述の時間増分Δｔは、正弦波状に変化するものとして得られる応力のピーク近傍において小さくすることが良い。これにより、応力のピーク近傍において、より多くの計算結果を得ることができ、カーブフィット等の近似曲線の精度を高め、ひいてはこれらから応力の振幅σ_０などをより正確に求めることができる。他方、応力のピーク近傍以外では、時間増分Δｔは、前記ピーク近傍よりも、より大きく設定することが良い。これにより、本来あまり必要ではない領域での計算の回数を減らすことができる。従って、より少ない計算量で精度良くゴム材料モデルの応力の応答を調べることができる。

応力のピーク近傍において前記時間増分Δｔを小さくするためには、歪と、これによって生じるであろう応力との大凡の位相差δを予め把握しておく必要がある。例えば、シミュレーションとして、任意の測定条件において、ｔａｎδが０．１〜０．４０の範囲のゴム材料を対象とする場合、ｔａｎ^-1δは最大でも０．３８ラジアンである。従って、図４に示されるように、ゴム材料に歪を正弦波状に変化させて与えると、それに対応して生じる応力（収束した応力）は、０よりも大きくかつ０．３８ラジアン以下の位相差δが生じると予測できる。

つまり、応力のピークＰ２も、歪のピークＰ１から０．３８ラジアンの領域Ａ１で生じる。従って、応力のピーク付近を含んでその近傍領域での計算頻度を高めるためには、歪の前記時間増分Δｔを、歪のピークＰ１から領域Ａ１の範囲で小さくする一方、それ以外の領域で時間増分Δｔを大きく設定することによって、計算精度と計算時間とをバランス良く両立することができる。

図５には、このような方法に基づいて、変化する時間増分Δｔで正弦波状に変化する歪を離散的に付与する具体例が示されている。図５の実施形態では、歪のピークＰ１から領域Ａの範囲の時間増分Δｔ１だけが小さく設定される。それ以外の領域の時間増分Δｔ２は、時間増分Δｔ１よりも大きくかつ一定で設定されている。なお、時間増分Δｔ１、Δｔ２の値自体については、測定条件の周波数等に応じて適宜設定できるが、前記時間増分Δｔ１については、好ましくは０．１〜０．３ラジアン程度が望ましい。

図６には、時間増分Δｔの変化について、さらに他の実施形態が示される。この実施形態では、前記領域Ａ１を含む歪のピークＰ２の前後の近傍領域Ｂの時間増分Δｔ１が、それ以外の領域の時間増分Δｔ２よりも小さく設定されている。この近傍領域Ｂは、歪の１周期の１／４程度の範囲とするのが好適である。

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

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

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

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

前記変形計算が行われると、その結果から必要な物理量を取得することができる（ステップＳ４）。物理量としては、計算によって得られた応力の振幅σ_０と、応力と前記歪との位相差δとを用いて、貯蔵弾性率Ｅ’、損失弾性率Ｅ”及び／又は損失正接ｔａｎδなどを計算するステップが含まれるのが望ましい。

図９には、このようなシミュレーション方法の具体例が示される。図９（ａ）は、ゴム材料モデルに与えられる変形条件としての歪について、その値と時間との関係が示されている。この例の各条件は次の通りである。
ゴム材料モデルの初期伸長：１０％
周波数：１０Hz
ひずみの片振幅：１％
時間増分Δｔ１：０．００２７５sec
時間増分Δｔ２：０．００５５sec
変形１周期中の計算回数：１４ステップ

図９（ｂ）には、このような変形条件によって計算された応力と時間との関係が示されている（結果は、位相差δが収束した変形７周期後の状態が示されている。）。図９（ｂ）から明らかなように、時間増分は、応力のピーク近傍で小さく、かつ、ピーク近傍以外で大きくなっている。図９（ｂ）の計算結果に正弦波のカーブフィットを施し、応力振幅σ0及び位相差δを計算し、それらに基づいて貯蔵弾性率Ｅ’、損失弾性率Ｅ”及び損失正接ｔａｎδを計算した。結果は、下記の通りである。
Ｅ’＝４．７１（ＭＰａ）
Ｅ”＝０．８４（ＭＰａ）
ｔａｎδ＝０．１８

他方、図１０には、従来のシミュレーション方法の結果を示し、（ａ）は変形条件としての歪、（ｂ）はそれに基づいて計算された応力（結果が収束した７周期後の結果を示している）を示す。時間増分Δｔは０．００２５secで一定とし、歪の１周期中の計算回数は４０ステップである。

図１０（ｂ）から、貯蔵弾性率Ｅ’、損失弾性率Ｅ”及び損失正接ｔａｎδを計算したところ、図９（ｂ）から計算した値と一致することが確認できた。一方、図１０のシミュレーション方法の計算時間を１００とすると、図９のシミュレーション方法の計算時間は約３５でなり、大幅な短縮が確認できた。

このように、本発明のシミュレーション方法によれば、より少ない計算量で精度良くゴム材料モデルの応力の応答を調べることができる。

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

Claims (3)

1. ゴム材料に正弦波状に変化する歪を与えたときの応力の変化をコンピュータを用いて計算するゴム材料のシミュレーション方法であって、
   前記コンピュータに、ゴム材料を有限個の要素で分割したゴム材料モデルを入力するステップと、
   前記コンピュータが、入力された前記ゴム材料モデルの変形条件として、正弦波状に変化する歪を、１計算ステップあたりの時間増分の間隔で離散的に付与する歪付与ステップと、
   前記ゴム材料モデルと前記歪の値とに基づいて、各歪状態での前記ゴム材料モデルに生じる応力を計算する応力計算ステップとを行うとともに、
   前記歪付与ステップは、前記時間増分が変化するとともに、前記時間増分は、前記応力のピーク近傍で小さく、かつ、前記ピーク近傍以外で大きいことを特徴とするゴム材料のシミュレーション方法。
2. 前記時間増分は、前記歪のピークから０．３８ラジアンの領域で小さく、それ以外の領域で大きい請求項１に記載のゴム材料のシミュレーション方法。
3. 前記コンピュータは、応力の振幅と、応力と前記歪との位相差を得るステップと、
   これらの値を用いて貯蔵弾性率Ｅ’、損失弾性率Ｅ”及び／又は損失正接ｔａｎδを計算するステップとをさらに含む請求項１又は２に記載のゴム材料のシミュレーション方法。
   

Images