ゴム材料に正弦波状に変化する歪を与えたときの応力の変化をコンピュータを用いて計算するゴム材料のシミュレーション方法であって、前記コンピュータに、ゴム材料を有限個の要素で分割したゴム材料モデルを入力するステップと、前記コンピュータが、入力された前記ゴム材料モデルの変形条件として、正弦波状に変化する歪を、1計算ステップあたりの時間増分の間隔で離散的に付与する歪付与ステップと、前記ゴム材料モデルと前記歪の値とに基づいて、各歪状態での前記ゴム材料モデルに生じる応力を計算する応力計算ステップとを行うとともに、前記歪付与ステップは、前記時間増分が変化することを特徴とするゴム材料のシミュレーション方法。前記時間増分は、前記応力のピーク近傍で小さく、かつ、前記ピーク近傍以外で大きい請求項1記載のゴム材料のシミュレーション方法。前記時間増分は、前記歪のピークから0.38ラジアンの領域で小さく、それ以外の領域で大きい請求項1又は2に記載のゴム材料のシミュレーション方法。前記コンピュータは、応力の振幅と、応力と前記歪との位相差を得るステップと、これらの値を用いて貯蔵弾性率E’、損失弾性率E及び/又は損失正接tanδを計算するステップとをさらに含む請求項1乃至3のいずれかに記載のゴム材料のシミュレーション方法。ところで、上述のシミュレーションでは、ゴム材料モデルの変形条件として、図11に示したような正弦波状に変化する歪aを連続的に与えることができない。このため、シミュレーションでは、図10aに示されるように、歪が、一定の時間増分Δtの間隔で離散的にかつ正弦波状に付与される。そして、この離散的に与えられる各歪に応じた応力が収束計算等によりそれぞれ求められる。応力の応答は、例えば図10bに示される。また、応力の変化を調べるために、離散的に得られた応力の各計算点に正弦波でカーブフィットを施して近似波形を求め、この近似波形に基づいて、例えば応力の振幅σ0等が求められていた。しかしながら、ゴム材料の重要な物性値である前記貯蔵弾性率E’、損失弾性率E及び/又は損失正接tanδを計算するためには、応力の振幅σ0と位相差δとが分かれば良い。従って、このようなシミュレーションを短時間で終わらせるためには、応力の前記近似波形の全て求めるのは得策ではなく、応力のピーク近傍領域のみの必要最低限の計算結果を得ることがより能率的である。本発明は、以上のような問題点に鑑み案出なされたもので、1計算ステップ当たりの時間増分Δtを変化させることを基本として、貯蔵弾性率E’、損失弾性率E及び/又は損失正接tanδといった物性値を短時間でかつ精度良く計算しうるゴム材料のシミュレーション方法を提供することを主たる目的としている。本発明のうち請求項1記載の発明は、ゴム材料に正弦波状に変化する歪を与えたときの応力の変化をコンピュータを用いて計算するゴム材料のシミュレーション方法であって、前記コンピュータに、ゴム材料を有限個の要素で分割したゴム材料モデルを入力するステップと、前記コンピュータが、入力された前記ゴム材料モデルの変形条件として、正弦波状に変化する歪を、1計算ステップあたりの時間増分の間隔で離散的に付与する歪付与ステップと、前記ゴム材料モデルと前記歪の値とに基づいて、各歪状態での前記ゴム材料モデルに生じる応力を計算する応力計算ステップとを行うとともに、前記歪付与ステップは、前記時間増分が変化することを特徴とする。また請求項2記載の発明は、前記時間増分は、前記応力のピーク近傍で小さく、かつ、前記ピーク近傍以外で大きい請求項1記載のゴム材料のシミュレーション方法である。また請求項3記載の発明は、前記時間増分は、前記歪のピークから0.38ラジアンの領域で小さく、それ以外の領域で大きい請求項1又は2に記載のゴム材料のシミュレーション方法である。また請求項4記載の発明は、前記コンピュータは、応力の振幅と、応力と前記歪との位相差を得るステップと、これらの値を用いて貯蔵弾性率E’、損失弾性率E及び/又は損失正接tanδを計算するステップとをさらに含む請求項1乃至3のいずれかに記載のゴム材料のシミュレーション方法である。本発明のゴム材料のシミュレーション方法によれば、コンピュータは、ゴム材料モデルの変形条件として、正弦波状に変化する歪を、1計算ステップ当たりの時間増分の間隔で離散的に付与する歪付与ステップと、前記ゴム材料モデルと前記歪の値とに基づいて、前記ゴム材料モデルに生じる応力を計算する応力計算ステップとを行う。そして、前記歪付与ステップでは、時間増分が変化することを特徴としている。従って、例えば、時間増分を、応力のピーク近傍において小さくすることにより、計算精度が要求されるピーク近傍での応力の正弦波状の計算結果を多数得ることができる。他方、ピーク近傍以外では、時間増分を前記ピーク近傍よりも、より大きく設定することにより、計算回数を減らすことができる。従って、より少ない計算量で精度良くゴム材料モデルの応力の応答等を調べることができる。以下、本発明の実施の一形態が図面に基づき説明する。図1には、本発明のシミュレーション方法を実施するためのコンピュータ装置1が示されている。このコンピュータ装置1は、本体1a、キーボード1b、マウス1c及びディスプレイ装置1dを含んで構成されている。本体1aの内部には、CPU、ROM、作業用メモリー及び磁気ディスク等の大容量記憶装置が設けられる。また、本体1aには、CD-ROMやフレキシブルディスクのドライブ装置1a1、1a2が設けられる。そして、前記大容量記憶装置には後述する本発明のシミュレーション方法を実行するための処理手順プログラムが記憶されている。図2には、本実施形態のシミュレーション方法の処理手順の一例が示される。本実施形態では、先ず、ゴム材料モデルが設定されるステップS1。本実施形態のゴム材料は、タイヤ用のゴム材料である。このようなゴム材料は、例えば初期歪10、周波数10Hz、歪片振幅1、温度30°Cの測定条件で、損失正接tanδが0.1~0.4程度の値を示すものが一般的である。図3a及びそのX部拡大図であるbに視覚化して示されるように、ゴム材料モデル2は、例えばシリカ等のフィラーがマトリックスゴム中に分散されたフィラー入りのゴム材料がモデル化されている。また、図3aのモデルは、繰り返し最小単位の微視構造としてのゴム材料モデル2の一例が視覚化して示されている。該ゴム材料モデル2は、解析しようとするフィラー配合ゴム材料の微小領域が、有限個の小さな要素メッシュ2a、2b、2c...に置き換えられたものである。各要素2a、2b、2c...は、数値解析が可能に定義される。前記数値解析が可能とは、例えば有限要素法、有限体積法、差分法又は境界要素法といった数値解析法により、各要素ないし系全体についての変形計算が可能なことを意味する。具体的には、各要素2a、2b、2c...について、座標系における節点座標値、要素形状、材料特性などが定義される。各要素2a、2b、2c...には、例えば2次元平面としての三角形ないし四辺形の要素、3次元要素としては、例えば4ないし6面体の要素が好ましく用いられる。これにより、ゴム材料モデル2は、前記コンピュータ装置1にて取り扱い可能な数値データを構成する。この実施形態のゴム材料モデル2は、後述する変形シミュレーションにおいて平面ひずみ状態の解析、さらに詳しくはy軸方向の引張変形シミュレーションが行われる。つまり、本実施形態では、z方向にはひずみを持たない2次元のシミュレーションが行われる。この実施形態において、微視構造としてのゴム材料モデル2は、例えば300nm×300nmの正方形である。また、本実施形態のゴム材料モデル2は、ゴムマトリックス部分がモデル化されたマトリックスモデル3と、このマトリックスモデル3の中に配置されかつシリカがモデル化されたシリカモデル4と、各シリカモデル4の周りを環状に取り囲みかつ前記マトリックスモデル3よりも硬い物性が定義された界面モデル5図3bに示すとを含む。前記マトリックスモデル3は、ゴム材料モデル2の主要部を構成し、かつ、例えば三角形ないし四辺形の複数個の要素を用いて表現されている。変形計算を行うために、マトリックスモデル3を構成する各要素には、その物性として応力と伸びとの関係を表す関数が定義される。本実施形態のゴム材料のシミュレーション方法では、ゴム弾性応答を表現するために、前記マトリックスモデル3及び界面モデル5のゴム部分は、いずれも分子鎖網目理論に基づいて計算が行われる。なお、分子鎖網目理論については、本件出願人らが提案している特開2010-205165号公報や特開2009-216612号公報に記載されており、これらの文献に記載の方法で本シミュレーションの計算過程に導入されている。従って、分子鎖網目理論については、ここでの詳細な説明は省略する。前記シリカモデル4は、シリカを四辺形の複数個の要素を用いてモデル化したもので、全体として円形に形成されている。三次元モデルの場合、シリカモデル4は、球形にモデル化されるのが望ましい。また、シリカは、直径約10~300nm程度であり、ゴムに比べて非常に硬い粒子からなる。シリカモデル4には、このような解析対象となるシリカの物性とほぼ等しい物性が設定される。即ち、本実施形態において、シリカモデル4は、粘弾性体ではなく弾性体として取り扱われる。また、シリカモデル4の粒子の個数は、例えば、解析対象のゴム材料のシリカ配合量に基づいて適宜決定される。前記界面モデル5は、シリカとマトリックスゴムとを化学的に結合させるシランカップリング剤の働きをシミュレーションに取り込むためにモデル化したものである。本実施形態の界面モデル5は、シリカモデル4の周りを小さい厚さtで環状に連続して取り囲むように設定されている。従って、界面モデル5の内周面はシリカモデル4の外周面に接触して固着されている。本実施形態では、界面モデル5の内周面とシリカモデル4の外周面とは、互いに剥離しない条件が設定されるが、必要に応じて、予め定めた値以上の応力が生じたときに、シリカモデル4と界面モデル5との境界を分離させるような条件が設定されても良い。なお、界面モデル5の外周面は、マトリックスゴムモデル3に接触している固着されている。前記界面モデル5の厚さtは、特に限定されるものではないが、種々の実験結果などに鑑み、シリカモデル4の直径の10~30程度、より好ましくは15~25程度に設定されるのが実際のゴム材料と整合する点で望ましい。前記界面モデル5にも、応力と伸びとの関係が定義される。実際の界面結合材の物性に鑑み、この界面モデル5はマトリックスゴムよりも硬い物性、即ち、マトリックスモデル3よりも伸び難く定義される。ただし、界面モデル5は、シリカモデル4よりは軟らかいのは言うまでもない。また、マトリックスモデル3、シリカモデル4及び界面モデル5には、それぞれ、これまでの実験結果等に基づいて、物性についてのパラメータ密度、弾性率等が入力される。次に、ゴム材料モデル2を変形させるための変形条件が設定されるステップS2。本実施形態では、図3のy方向に正弦波状の歪を加えてゴム材料モデル2に歪が正弦波状に変化する動的な引張変形を与える条件が定義される。また、正弦波状に変化する歪は、1計算ステップあたりの時間増分Δtの間隔で離散的に付与されるとともに、該時間増分Δtが変化することを特徴としている。例えば、上述の時間増分Δtは、正弦波状に変化するものとして得られる応力のピーク近傍において小さくすることが良い。これにより、応力のピーク近傍において、より多くの計算結果を得ることができ、カーブフィット等の近似曲線の精度を高め、ひいてはこれらから応力の振幅σ0などをより正確に求めることができる。他方、応力のピーク近傍以外では、時間増分Δtは、前記ピーク近傍よりも、より大きく設定することが良い。これにより、本来あまり必要ではない領域での計算の回数を減らすことができる。従って、より少ない計算量で精度良くゴム材料モデルの応力の応答を調べることができる。応力のピーク近傍において前記時間増分Δtを小さくするためには、歪と、これによって生じるであろう応力との大凡の位相差δを予め把握しておく必要がある。例えば、シミュレーションとして、任意の測定条件において、tanδが0.1~0.40の範囲のゴム材料を対象とする場合、tan1δは最大でも0.38ラジアンである。従って、図4に示されるように、ゴム材料に歪を正弦波状に変化させて与えると、それに対応して生じる応力収束した応力は、0よりも大きくかつ0.38ラジアン以下の位相差δが生じると予測できる。つまり、応力のピークP2も、歪のピークP1から0.38ラジアンの領域A1で生じる。従って、応力のピーク付近を含んでその近傍領域での計算頻度を高めるためには、歪の前記時間増分Δtを、歪のピークP1から領域A1の範囲で小さくする一方、それ以外の領域で時間増分Δtを大きく設定することによって、計算精度と計算時間とをバランス良く両立することができる。図5には、このような方法に基づいて、変化する時間増分Δtで正弦波状に変化する歪を離散的に付与する具体例が示されている。図5の実施形態では、歪のピークP1から領域Aの範囲の時間増分Δt1だけが小さく設定される。それ以外の領域の時間増分Δt2は、時間増分Δt1よりも大きくかつ一定で設定されている。なお、時間増分Δt1、Δt2の値自体については、測定条件の周波数等に応じて適宜設定できるが、前記時間増分Δt1については、好ましくは0.1~0.3ラジアン程度が望ましい。図6には、時間増分Δtの変化について、さらに他の実施形態が示される。この実施形態では、前記領域A1を含む歪のピークP1の前後の近傍領域Bの時間増分Δt1が、それ以外の領域の時間増分Δt2よりも小さく設定されている。この近傍領域Bは、歪の1周期の1/4程度の範囲とするのが好適である。次に本実施形態のシミュレーション方法では、上述のように設定されたゴム材料モデル2を用いて変形シミュレーションが行われるステップS3。コンピュータ装置1が行う変形シミュレーションの具体的な処理手順は、図7に示される。変形シミュレーションでは、先ずゴム材料モデル2の各種のデータがコンピュータ装置1に入力されるステップS31。入力されるデータには、各要素に定義された節点の位置や材料特性といった情報が含まれる。コンピュータ装置1では、入力されたデータに基づいて各要素の剛性マトリックスを作成しステップS32、しかる後、全体構造の剛性マトリックスを組み立てるステップS33。全体構造の剛性マトリックスには、既知節点の変位、節点力が導入されステップS34、剛性方程式の解析が行われる。そして、未知節点変位が決定されステップS35、先に定められた変形条件各ひずみに基づいた応力を計算するステップS36。この際、コンピュータ装置1は、先に述べたように、ゴム材料モデル2に、正弦波状に変化する歪を、離散的かつ時間増分Δtを変化させて与える歪付与ステップ。そして、コンピュータ装置1は、その結果を記憶装置等に出力するステップ37。ステップS38では、計算を終了させるか否かの判定がなされ、否定的である場合には、ステップS32以降が繰り返される。このようなシミュレーション変形計算は、例えば有限要素法を用いたエンジニアリング系の解析アプリケーションソフトウエア例えば米国リバモアソフトウェアテクノロジー社で開発改良されたLS-DYNA等を用いて行うことができる。また、本シミュレーションは、均質化法漸近展開均質化法に基づいて行われる。均質化法は、図8に示されるように、図3に示した微視構造均質化法ではユニットセルとも呼ばれるを周期的に持っているゴム材料全体Mを表現するxIと、前記微視構造を表現するyIとの独立した2変数が用いられる。微視的スケールと巨視的スケールという異なる尺度の場におけるそれぞれ独立した変数を漸近展開することにより、図3に示した微視構造のモデル構造を反映させたゴム材料全体の平均的な力学応答を近似的に求めることができる。前記変形計算が行われると、その結果から必要な物理量を取得することができるステップS4。物理量としては、計算によって得られた応力の振幅σ0と、応力と前記歪との位相差δとを用いて、貯蔵弾性率E’、損失弾性率E及び/又は損失正接tanδなどを計算するステップが含まれるのが望ましい。図9には、このようなシミュレーション方法の具体例が示される。図9aは、ゴム材料モデルに与えられる変形条件としての歪について、その値と時間との関係が示されている。この例の各条件は次の通りである。ゴム材料モデルの初期伸長:10周波数:10Hzひずみの片振幅:1時間増分Δt1:0.00275sec時間増分Δt2:0.0055sec変形1周期中の計算回数:14ステップ図9bには、このような変形条件によって計算された応力と時間との関係が示されている結果は、位相差δが収束した変形7周期後の状態が示されている。。図9bから明らかなように、時間増分は、応力のピーク近傍で小さく、かつ、ピーク近傍以外で大きくなっている。図9bの計算結果に正弦波のカーブフィットを施し、応力振幅σ0及び位相差δを計算し、それらに基づいて貯蔵弾性率E’、損失弾性率E及び損失正接tanδを計算した。結果は、下記の通りである。E’=4.71MPaE=0.84MPatanδ=0.18他方、図10には、従来のシミュレーション方法の結果を示し、aは変形条件としての歪、bはそれに基づいて計算された応力結果が収束した7周期後の結果を示しているを示す。時間増分Δtは0.0025secで一定とし、歪の1周期中の計算回数は40ステップである。図10bから、貯蔵弾性率E’、損失弾性率E及び損失正接tanδを計算したところ、図9bから計算した値と一致することが確認できた。一方、図10のシミュレーション方法の計算時間を100とすると、図9のシミュレーション方法の計算時間は約35でなり、大幅な短縮が確認できた。このように、本発明のシミュレーション方法によれば、より少ない計算量で精度良くゴム材料モデルの応力の応答を調べることができる。