(57)【特許請求の範囲】
【請求項１】
弾性要素と粘弾性要素とが並列配置された粘弾性材料構成則に基いて、粘弾性材料の特性を解析する解析装置であって、
節点を有する有限数の要素に分割された粘弾性材料モデルに所定の入力条件を設定して前記節点の変位量を計算する第１計算部と、
前記変位量を用いて前記節点におけるひずみ速度を計算する第２計算部と、
前記ひずみ速度を底とする冪乗の値及び温度−時間換算則のシフトファクターを底とする冪乗の値に比例する値を前記粘弾性要素の緩和時間として計算する第３計算部と、
前記緩和時間を用いて前記節点における応力を計算する第４計算部と、
を備えることを特徴とする解析装置。
【請求項２】
請求項１に記載の解析装置において、
前記シフトファクターを計算するシフトファクター計算部をさらに備え、
前記シフトファクター計算部は、
解析対象である前記粘弾性材料を用いて複数の環境温度にて実施された調和振動試験の試験結果から、前記複数の環境温度のそれぞれにおける前記粘弾性材料の弾性率を算出し、
前記複数の環境温度の中の基準温度での基準弾性率を用いて前記複数の環境温度のそれぞれの前記弾性率を標準化し、
当該標準化された弾性率の対数値と温度との関係を表す関数ｆ（Ｔ）を求め、
前記粘弾性材料を用いた調和振動試験における弾性率と振幅との関係から得られた材料定数ｍｉと、前記関数ｆ（Ｔ）と、以下の（Ａ）式とに基いてシフトファクターα（Ｔ）を算出する
ように構成された、
解析装置。
【数１】
0006988599.tif 000046
【請求項３】
請求項１に記載の解析装置において、
前記シフトファクターを算出するシフトファクター計算部をさらに備え、
前記シフトファクター計算部は、
解析対象である前記粘弾性材料を用いて複数の環境温度にて実施された応力緩和試験の試験結果から、前記複数の環境温度のそれぞれについて応力緩和グラフを作成し、
前記複数の環境温度の中の基準温度の最大応力を用いて前記複数の環境温度のそれぞれの応力緩和グラフを標準化して、前記複数の環境温度のそれぞれにおける緩和時間を算出し、
前記基準温度の緩和時間を用いて、前記複数の環境温度のそれぞれにおける緩和時間の対数値を標準化し、当該標準化された緩和時間の対数値と温度との関係を表す関数ｆ（Ｔ）を求め、
前記粘弾性材料を用いた調和振動試験における弾性率と振幅との関係から得られた材料定数ｍｉと、前記関数ｆ（Ｔ）と、以下の（Ｂ）式とに基いてシフトファクターα（Ｔ）を算出する
ように構成された、
解析装置。
【数２】
0006988599.tif 000047
【請求項４】
請求項１乃至請求項３のいずれか一項に記載の解析装置において、
前記シフトファクターを底とする冪乗の値の冪数は、前記ひずみ速度を底とする冪乗の値の冪数に１を足した値である、
解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、有限要素法により粘弾性材料の特性を解析する解析装置に関する。
【背景技術】
【０００２】
従来から、粘弾性材料（例えば、ゴム）の特性を解析するためのモデルとして、弾性要素と粘弾性要素を並列に配置した粘弾性材料構成則が知られている。従来の装置の一つ（以下、「従来装置」と称呼する。）は、粘弾性要素の減衰特性を示す緩和時間を定数とすることにより、粘弾性材料の応力−ひずみ特性を解析する（例えば、特許文献１）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特許第６０４８３５８号
【発明の概要】
【０００４】
粘弾性材料の挙動は環境温度によって変化する場合がある。しかしながら、従来装置においては、緩和時間が、環境温度に依存しない関数として定義されている。従って、従来装置は、粘弾性材料の環境温度に対する依存性を考慮して、粘弾性材料の特性（例えば、応力−ひずみ特性）を解析することができない。
【０００５】
本発明は、上記課題を解決するためになされた。即ち、本発明の目的の一つは、粘弾性材料の環境温度に対する依存性を再現することが可能な解析装置を提供することである。
【０００６】
本発明の解析装置（以下、「本発明装置」と称呼される場合がある。）は、弾性要素と粘弾性要素とが並列配置された粘弾性材料構成則に基いて、粘弾性材料の特性を解析する装置である。
当該解析装置は、
節点を有する有限数の要素に分割された粘弾性材料モデルに所定の入力条件を設定して前記節点の変位量を計算する第１計算部（１０４１）と、
前記変位量を用いて前記節点におけるひずみ速度を計算する第２計算部（１０４２）と、
前記ひずみ速度を底とする冪乗の値及び温度−時間換算則のシフトファクター（α（Ｔ））を底とする冪乗の値に比例する値を前記粘弾性要素の緩和時間（τｉ）として計算する第３計算部（１０４３）と、
前記緩和時間を用いて前記節点における応力を計算する第４計算部（１０４４）と、
を備える。
【０００７】
係る構成を有する本発明装置では、粘弾性要素の減衰特性を表す緩和時間が、温度−時間換算則のシフトファクターの冪関数を用いて算出される。従って、本発明装置は、粘弾性材料の環境温度に対する依存性を再現することができる。これにより、環境温度が変化した際の粘弾性材料の特性（例えば、応力−ひずみ特性）の予測精度を高めることができる。
【０００８】
本発明装置の一の態様は、前記シフトファクターを算出するシフトファクター計算部（１０３０）をさらに備える。
前記シフトファクター計算部は、
解析対象である前記粘弾性材料を用いて複数の環境温度にて実施された調和振動試験の試験結果から、前記複数の環境温度のそれぞれにおける前記粘弾性材料の弾性率を算出し、
前記複数の環境温度の中の基準温度（Ｔ０）での基準弾性率（Ｇ（Ｔ０））を用いて前記複数の環境温度のそれぞれの前記弾性率を標準化し、
当該標準化された弾性率の対数値と温度との関係を表す関数ｆ（Ｔ）を求め、
前記粘弾性材料を用いた調和振動試験における弾性率と振幅との関係から得られた材料定数ｍｉと、前記関数ｆ（Ｔ）と、以下の（Ａ）式とに基いてシフトファクターα（Ｔ）を算出する
ように構成されている。
【０００９】
【数１】
0006988599.tif 000002
【００１０】
本態様によれば、複数の環境温度にて実施された調和振動試験の試験結果から、温度−時間換算則のシフトファクターを算出することができる。
【００１１】
本発明装置の一の態様は、前記シフトファクターを算出するシフトファクター計算部（１０３０）をさらに備える。
前記シフトファクター計算部は、
解析対象である前記粘弾性材料を用いて複数の環境温度にて実施された応力緩和試験の試験結果から、前記複数の環境温度のそれぞれについて応力緩和グラフを作成し、
前記複数の環境温度の中の基準温度（Ｔ０）の最大応力（σ０）を用いて前記複数の環境温度のそれぞれの応力緩和グラフを標準化して、前記複数の環境温度のそれぞれにおける緩和時間を算出し、
前記基準温度の緩和時間を用いて、前記複数の環境温度のそれぞれにおける緩和時間の対数値を標準化し、当該標準化された緩和時間の対数値と温度との関係を表す関数ｆ（Ｔ）を求め、
前記粘弾性材料を用いた調和振動試験における弾性率と振幅との関係から得られた材料定数ｍｉと、前記関数ｆ（Ｔ）と、以下の（Ｂ）式とに基いてシフトファクターα（Ｔ）を算出する
ように構成されている。
【００１２】
【数２】
0006988599.tif 000003
【００１３】
本態様によれば、複数の環境温度にて実施された応力緩和試験の試験結果から、温度−時間換算則のシフトファクターを算出することができる。
【００１４】
本発明装置の一の態様において、前記シフトファクターを底とする冪乗の値の冪数は、前記ひずみ速度を底とする冪乗の値の冪数に１を足した値である。
【００１５】
上記説明においては、本発明の理解を助けるために、後述する実施形態に対応する発明の構成に対し、その実施形態で用いた名称及び／又は符号を括弧書きで添えている。しかしながら、本発明の各構成要素は、前記名称及び／又は符号によって規定される実施形態に限定されるものではない。
【図面の簡単な説明】
【００１６】
【図１】粘弾性材料構成則を模式的に示す図である。
【図２】調和振動試験における試験条件を示す図である。
【図３】材料定数ｍｉの算出方法を模式的に示すグラフである。
【図４】応力緩和試験における応力緩和曲線を模式的に示す図である。
【図５】入力ひずみを変化させた場合における応力緩和曲線を示す図である。
【図６】ひずみと応力の静特性曲線を示す図である。
【図７】一定ひずみ速度試験におけるひずみと応力の関係曲線及び静特性曲線を示す図である。
【図８】本発明の実施形態に係る温度−時間換算則のシフトファクターα（Ｔ）を算出する流れを示す図である。
【図９】本発明の実施形態に係る解析装置のハードウェア構成を示すブロック図である。
【図１０】本発明の実施形態に係る解析装置の機能構成を示すブロック図である。
【図１１】本発明の実施形態に係るＣＰＵが実行する「解析ルーチン」を示したフローチャートである。
【図１２】（ａ）は、ゴム材料に対して調和振動試験を行った際の温度と弾性率との関係を表すグラフであり、（ｂ）は、本発明の実施形態に係る解析装置によって温度と弾性率との関係を算出した結果を表すグラフである。
【図１３】（ａ）は、ゴム材料に対して調和振動試験を行った際の「温度」と「応力とひずみの位相差」との関係を表すグラフであり、（ｂ）は、本発明の実施形態に係る解析装置によって「温度」と「応力とひずみの位相差」との関係を算出した結果を表すグラフである。
【図１４】本発明の変形例に係る温度−時間換算則のシフトファクターを算出する流れを示す図である。
【発明を実施するための形態】
【００１７】
以下、添付図面を参照して本発明の実施形態について説明する。なお、添付図面は本発明の原理に則った具体的な実施形態を示しているが、これらは本発明を理解するための例あり、本発明を限定的に解釈するために用いられるべきでない。
【００１８】
＜粘弾性材料構成則を用いた解析の概要＞
図１を用いて、本実施形態において用いられる粘弾性材料構成則を説明する。図１に示す粘弾性材料モデルでは、弾性率Ｇ０の弾性要素と、複数の粘弾性要素とが並列に接続される。複数の粘弾性要素のそれぞれは、弾性率Ｇｉ（ｉ＝１〜Ｎ、Ｎは自然数）の弾性要素と、当該弾性要素に直列に接続された粘性係数ηｉの粘性要素とを有している。弾性要素と粘弾性要素とを組み合わせた本モデルは、タイヤ及びゴムブッシュ等のゴム部品の特性を表すためのモデルとして用いられる。
【００１９】
ここで、並列接続される弾性要素及び粘弾性要素における剛性割合をγｉ（ｉ＝０〜Ｎ）、粘弾性要素における緩和時間をτｉ＝ηｉ／Ｇｉとする。この粘弾性材料モデルにおける時刻ｔでの応力Ｓは、下記（１）式及び（２）式により表される。
【００２０】
【数３】
0006988599.tif 000004
【００２１】
（１）式において、Ｓは、第２Ｐｉｏｌａ−Ｋｉｒｃｈｈｏｆｆ応力を示し、上付き記号を有するＳ〇は、粘性力成分を除去した弾性成分のみの応力を示す。Ｊは、粘弾性材料の体積変化率を示す。体積変化率Ｊは、ある物質点における変形前と変形後の位置の線形変換の関係を示す変形勾配テンソルＦのデターミナント（ｄｅｔ）を用いて、Ｊ＝ｄｅｔ［Ｆ］により表される。演算子ＤＥＶは、右Ｃａｕｃｈｙ−ＧｒｅｅｎテンソルＣ＝ＦＴ・Ｆを用いて、下記（３）式で表される。なお、（３）式における［・］は、演算子ＤＥＶの演算対象となる変数を表す。
【００２２】
【数４】
0006988599.tif 000005
【００２３】
（１）式におけるＱｉは、それぞれの粘弾性要素における粘性力を示す。Ｑｉは（２）式に示す発展方程式により表される。
【００２４】
【数５】
0006988599.tif 000006
は、超弾性体におけるひずみポテンシャルエネルギーの偏差成分を表す。
【００２５】
【数６】
0006988599.tif 000007
は、体積成分を除去した修正右Ｃａｕｃｈｙ−Ｇｒｅｅｎテンソルであり、下記（４）式で表される。
【００２６】
【数７】
0006988599.tif 000008
【００２７】
上記（１）式及び（２）式により表される第２Ｐｉｏｌａ−Ｋｉｒｃｈｈｏｆｆ応力Ｓは、下記（５）式に示す積分形で表すことができる。
【００２８】
【数８】
0006988599.tif 000009
【００２９】
【数９】
0006988599.tif 000010
は、ひずみポテンシャルエネルギーの体積成分である。更に、ｇ（ｔ）は、緩和関数であり、上記（６）式により表される。
【００３０】
ここで、時間ｔの関数として表される（５）式を下記（７）式に変形することで、時刻ｔｎ＋１における第２Ｐｉｏｌａ−Ｋｉｒｃｈｈｏｆｆ応力Ｓｎ＋１を得ることができる。なお、計算ステップｎにおける関数を（・）ｎ、計算ステップｎ＋１における関数を（・）ｎ＋１と表記している。
【００３１】
【数１０】
0006988599.tif 000011
【００３２】
（７）式におけるＨ（ｉ）ｎ＋１は、中点定理を用いて［ｔｎ，ｔｎ＋１］の時間間隔で積分した近似解として得られる中間関数であり、下記（８）式で表される。
【００３３】
【数１１】
0006988599.tif 000012
【００３４】
【数１２】
0006988599.tif 000013
は、それぞれ、下記（９）式及び（１０）式で定義される。
【００３５】
【数１３】
0006988599.tif 000014
【００３６】
なお、Ｋｉｒｃｈｈｏｆｆ弾性応力
【数１４】
0006988599.tif 000015
を、下記（１１）式によって定義する。
【００３７】
【数１５】
0006988599.tif 000016
【００３８】
これにより、（９）式は、下記（１２）式により表すこともできる。
【００３９】
【数１６】
0006988599.tif 000017
【００４０】
更に、Ｋｉｒｃｈｈｏｆｆ応力テンソルτｎ＋１は、第２Ｐｉｏｌａ−Ｋｉｒｃｈｈｏｆｆ応力Ｓｎ＋１を用いて下記（１３）式により表すことができる。
【００４１】
【数１７】
0006988599.tif 000018
【００４２】
従って、（７）式及び（１３）式より、Ｋｉｒｃｈｈｏｆｆ応力テンソルを下記（１４）式により表すことができる。
【００４３】
【数１８】
0006988599.tif 000019
【００４４】
なお、（１４）式における演算子ｄｅｖは、下記（１５）式により定義される。（１５）式における（・）は、演算子ＤＥＶの演算対象となる変数を表す。
【００４５】
【数１９】
0006988599.tif 000020
【００４６】
更に、（１４）式における緩和関数ｇ＊は、下記（１６）式により定義される。
【００４７】
【数２０】
0006988599.tif 000021
【００４８】
以上より、中間関数として、Ｈ（ｉ）ｎ＋１、及び、
【数２１】
0006988599.tif 000022
を計算ステップ毎に保持することにより、Ｋｉｒｃｈｈｏｆｆ応力テンソルの値を得ることができる。
【００４９】
＜具体的な解析処理内容＞
従来装置の粘弾性モデルでは、緩和時間τｉが、環境温度に依存しない関数として定義されている。従って、環境温度が変化した際の粘弾性材料の特性（例えば、応力−ひずみ特性）の予測精度が高くないという課題があった。
【００５０】
そこで、本願発明者は、緩和時間τｉを、温度−時間換算則のシフトファクターの冪関数として定義することを見出した。これにより、粘弾性材料の環境温度に対する依存性を再現し、粘弾性材料の特性（例えば、応力−ひずみ特性）の予測精度を高めることができる。
【００５１】
本実施形態では、下記（１７）式に示す関係式を満たす緩和時間τｉを用いる。
【００５２】
【数２２】
0006988599.tif 000023
【００５３】
【数２３】
0006988599.tif 000024
は、偏差成分を除去したＧｒｅｅｎ−Ｌａｇｒａｎｇｅひずみテンソルであり、修正右Ｃａｕｃｈｙ−Ｇｒｅｅｎテンソル
【数２４】
0006988599.tif 000025
を用いて下記（１８）式で表される。
【００５４】
【数２５】
0006988599.tif 000026
【００５５】
【数２６】
0006988599.tif 000027
は、ひずみ速度であり、ひずみテンソルの時間微分値を表す。
【００５６】
【数２７】
0006988599.tif 000028
は、ひずみ速度の大きさを表し、３次元のひずみテンソルを用いる場合には、下記（１９）式により表される。
【００５７】
【数２８】
0006988599.tif 000029
【００５８】
α（Ｔ）は、本実施形態で新たに導入する値であり、温度−時間換算則のシフトファクターである。Ｔは環境温度である。Ａｉ及びｍｉは、それぞれ、材料定数である。従って、（１７）式に示すように、緩和時間τｉは、ひずみ速度を底とする冪乗の値及び温度−時間換算則のシフトファクターを底とする冪乗の値に比例する値として定義される。更に、温度−時間換算則のシフトファクターを底とする冪乗の値の冪数は、ひずみ速度を底とする冪乗の値の冪数に１を足した値である。
【００５９】
なお、シフトファクターα（Ｔ）の算出方法は後述する。更に、材料定数Ａｉ及びｍｉの算出方法についても後述する。
【００６０】
次に、（１７）式に示す緩和時間τｉの導入に際し、数値解析上必要となる式の変形について述べる。（１４）式を用いてＫｉｒｃｈｈｏｆｆ応力を算出するためには、緩和時間τｉを含む式を変形する必要がある。具体的には、（８）式を下記（２０）式に変形する。
【００６１】
【数２９】
0006988599.tif 000030
【００６２】
更に、（１６）式を下記（２２）式に変形する。
【００６３】
【数３０】
0006988599.tif 000031
【００６４】
ここで、緩和時間τｉの関数は、下記（２３）式及び（２４）式により定義される。
【００６５】
【数３１】
0006988599.tif 000032
【００６６】
なお、（２０）式の右辺第２項及び（２１）式の右辺第２項、並びに、（２２）式では、緩和時間τｉとして（２４）式が用いられる。これは、中点定理を用いて［ｔｎ，ｔｎ＋１］の時間間隔で積分した近似解を用いるためである。
【００６７】
次に、本実施形態に示す粘弾性材料モデルにおける材料定数ｍｉ及びＡｉの算出方法について述べる。
【００６８】
図２は、調和振動試験における試験条件を示す図である。本図は、横軸を時間ｔとして、予ひずみＥｐｒｅ、振幅ε、周波数ωとした場合における粘弾性材料の試験体に加振されるひずみＥ＝Ｅｐｒｅ＋εｓｉｎ（ωｔ）のグラフを示す。このようなひずみＥを加振した場合における試験体の動的弾性率Ｇを測定する。そして、下記（２５）式の関係を用いて材料定数ｍｉを導出することができる。
【００６９】
【数３２】
0006988599.tif 000033
【００７０】
ここで、（２）式に（１７）式を代入することにより、（２）式を下記（２６）式の発展方程式に変形する。（２５）式は（２６）式を解くことにより得られる。
【００７１】
【数３３】
0006988599.tif 000034
【００７２】
図３は、材料定数ｍｉの算出方法を模式的に示すグラフである。本図は、調和振動試験において加振する振幅εの対数値ｌｏｇ（ε）を横軸として、予ひずみＥｐｒｅ及び周波数ωを変化させたときの試験結果値である動的弾性率Ｇの対数値ｌｏｇ（Ｇ）をプロットしたものである。グラフ上の近似線Ｌ１は、予ひずみＥｐｒｅ１と周波数ω１に対して、振幅の値ε１〜ε３を変化させたときのプロット値に対応する傾きｌ１の近似線である。同様に、近似線Ｌ２は、予ひずみＥｐｒｅ２と周波数ω２に対して、振幅の値ε１〜ε３を変化させたときのプロット値に対応する傾きｌ２の近似線である。近似線Ｌ３は、予ひずみＥｐｒｅ３と周波数ω３に対して、振幅の値ε１〜ε３を変化させたときのプロット値に対応する傾きｌ３の近似線である。（２５）式の両辺に対数をとることにより、下記（２７）式を得ることができる。従って、傾きｌｉから材料定数ｍｉを得ることができる。
【００７３】
【数３４】
0006988599.tif 000035
【００７４】
（２７）式により、材料定数ｍｉは、以下の関係式（ｍｉ＝−１−ｌｉ）から得ることができる。なお、式（２７）におけるβは、図３に示す近似線の切片である。これにより、調和振動試験における加振周波数ωｉに対応する材料定数ｍｉを得ることができる。なお、図３に示す振幅εの値は例示である。振幅εとして、２つ値が用いられてもよいし、４以上の値が用いられてもよい。
【００７５】
なお、ｍｉの値が加振周波数ωｉによってそれほど変化しない場合には、それぞれの周波数に対応するｍｉを得る代わりに、それぞれのｍｉの平均値として得られる共通の冪数ｍが用いられてもよい。共通の冪数ｍを用いることにより、周波数ごとに異なるｍｉを用いる場合と比べて、計算負荷を減らすことができる。
【００７６】
材料定数Ａｉは、上記の方法で得られた材料定数ｍｉを用いて、下記（２８）式の関係から得ることができる。（２８）式における周波数ωｉは、調和振動試験における加振周波数である。この周波数を変化させることにより周波数成分ωｉが異なるそれぞれの粘弾性要素に対応する材料定数Ａｉを得ることができる。
【００７７】
【数３５】
0006988599.tif 000036
【００７８】
次に、弾性要素の剛性割合γ０の算出方法について述べる。図４は、応力緩和試験における応力緩和曲線を模式的に示す図である。本図は、入力ひずみＥｐｒｅを一定にした場合における応力Ｑの時間変化を測定した試験結果を示している。この試験結果に示す時刻ｔ＝０の瞬間応力Ｑ０と、時刻ｔ＝∞における緩和応力Ｑ∞とから、弾性要素の剛性割合を（γ０＝Ｑ∞／Ｑ０）により得ることができる。
【００７９】
次に、弾性要素の弾性率Ｇ０の算出方法について述べる。本実施形態では、超弾性の材料モデルであるＹｅｏｈ材料モデルに定められる超弾性係数Ｃ１０、Ｃ２０及びＣ３０を定めることにより弾性要素の弾性率Ｇ０を決定する。これらの超弾性係数は、図４に示した応力緩和試験において、複数の入力ひずみＥｐｒｅに対して緩和応力Ｑ∞を得ることにより導出することができる。
【００８０】
図５は、入力ひずみＥｋ（ｋ＝１〜Ｎ、Ｎは自然数）を変化させた場合における応力緩和曲線を示す図である。図６は、ひずみＥと応力Ｑの静特性曲線ＱＲを示す図である。図５に示すように、値の異なる複数の負荷ひずみＥ１〜ＥＮに対して応力緩和試験を実施し、各負荷ひずみＥｋに対応する緩和応力ＱＲ（Ｅ１）〜ＱＲ（ＥＮ）を得る。この関係式をひずみ−応力曲線としてグラフ化することにより、図６に示す静特性曲線ＱＲが得られる。この静特性曲線を近似することのできる超弾性係数Ｃ１０、Ｃ２０、Ｃ３０を公知の技術を用いて定めることにより、弾性要素の弾性率Ｇ０を決定することができる。
【００８１】
次に、粘弾性要素の剛性割合γｉの算出方法について述べる。図７は、一定ひずみ速度試験におけるひずみＥと応力Ｑの関係曲線Ｑｉ及び静特性曲線ＱＲを示している。静特性曲線ＱＲは、図６に示す静特性曲線と同じである。関係曲線Ｑ１〜ＱＮは、一定ひずみ速度試験においてひずみ速度ＶｉをＶ１〜ＶＮに変化させた場合に対応するひずみＥと応力Ｑの試験結果を示す。ここで、ｉ番目の関係曲線Ｑｉは、ｉ番目の粘弾性要素の動特性に対応し、加振周波数ωｉから算出した材料定数Ａｉの粘弾性要素に対応する。このとき、調和振動試験における加振周波数ωｉと振幅εとから、ひずみ速度Ｖｉは以下の関係式（Ｖｉ＝ωｉε）を満たす。
【００８２】
各粘弾性要素の剛性割合γｉは、図６に示すひずみＥと応力Ｑの関係曲線Ｑｉにより囲まれる面積Ｚｉから得ることができる。静特性曲線ＱＲにより囲まれる面積をＺ０とし、ｉ−１番目からｉ番目にひずみ速度が増加することによる面積の増加分をＺｉとすると、剛性割合γｉは以下の関係式（γｉ＝γ０×Ｚｉ／Ｚ０）により算出することができる。
【００８３】
以上の方法により、材料定数ｍｉ及びＡｉ、剛性割合γ０及びγｉ、並びに、弾性率Ｇ０を同定することができる。
【００８４】
次に、温度−時間換算則のシフトファクターα（Ｔ）の算出方法について説明する。図８（ａ）に示すように、複数の環境温度（Ｔ＝Ｔ０、Ｔ１、Ｔ２、…、ＴＮ）にて粘弾性材料の試験体に対して調和振動試験を実施する。このとき、振幅εは一定であり、且つ、予ひずみも一定である。
【００８５】
ここで、調和振動試験における弾性率Ｇ（ｔ）とシフトファクターα（Ｔ）との関係について、下記（２９）式が成り立つ。
【００８６】
【数３６】
0006988599.tif 000037
【００８７】
従って、まず、調和振動試験の試験結果から弾性率Ｇ（Ｔ）を算出する。ここで、複数の環境温度（Ｔ０、Ｔ１、Ｔ２、…、ＴＮ）のうち、Ｔ０を基準温度として用いる。基準温度Ｔ０の弾性率Ｇ（Ｔ０）を用いて、複数の環境温度での弾性率のそれぞれを標準化する。次に、図８（ｂ）に示すように、温度Ｔを横軸とし、「標準化された弾性率の対数値ｌｏｇ（Ｇ（Ｔ）／Ｇ（０））」をプロットした散布図を作成する。そして、図８（ｂ）に示すように、最小二乗法により回帰関数ｆ（Ｔ）を求める。回帰関数ｆ（Ｔ）は、回帰直線又は回帰曲線の何れであってもよい。ここで、求められた関数ｆ（Ｔ）とシフトファクターα（Ｔ）との関係について、下記（３０）式が成り立つ。
【００８８】
【数３７】
0006988599.tif 000038
【００８９】
従って、材料定数ｍｉを用いてシフトファクターα（Ｔ）を算出することができる。
【００９０】
＜解析装置のハードウェア構成＞
次に、有限要素法により粘弾性材料の特性を解析する解析装置（以下、「本実施装置」と称呼する場合がある。）ついて説明する。図９に示すように、解析装置は、情報処理装置９１０と、入力装置９２０と、出力装置９３０とを備える。
【００９１】
情報処理装置９１０は、ＣＰＵ９１１と、ＲＡＭ９１２と、ＲＯＭ９１３と、ハードディスク（ＨＤＤ）９１４と、Ｉ／Ｏインタフェース９１５とを備える。ＲＯＭ９１３は、ＣＰＵ９１１が実行するインストラクション（プログラム、ルーチン）を格納している。ＣＰＵ９１１は、当該インストラクションを実行することにより以降で説明する各種機能を実現するようになっている。
【００９２】
情報処理装置９１０は、Ｉ／Ｏインタフェース９１５を介して入力装置９２０及び出力装置９３０に接続されている。入力装置９２０は、ユーザからの各種要求を受け付ける装置であり、キーボード及びマウスを含む。出力装置９３０は、情報処理装置９１０が行った処理結果を出力するディスプレイを含む。
【００９３】
＜解析装置の機能構成＞
ＣＰＵ９１１は、ＲＯＭ９１３に格納されているインストラクションを読み出して実行することにより、図１０に示す「構築部１０１０、材料定数計算部１０２０、シフトファクター（ＳＦ）計算部１０３０及び応力計算部１０４０」の各機能を実現するようになっている。なお、図１０に示す「試験結果記憶部１０５０、モデル記憶部１０６０、パラメータ記憶部１０７０及び計算結果記憶部１０８０」は、上記のＲＡＭ９１２及び／又はＨＤＤ９１４によって実現される。
【００９４】
構築部１０１０は、解析対象とする粘弾性材料の粘弾性材料モデルを構築する。粘弾性材料モデルでは、解析対象となる粘弾性材料が要素（Ｅｌｅｍｅｎｔ）及び節点（Ｇｒｉｄ）の集合体として表現される。例えば、粘弾性材料を２次元モデルによって表現する場合、各要素として、３節点を有する三角形状の要素又は４節点を有する四角形状の要素が用いられる。粘弾性材料を３次元モデルによって表現する場合、各要素として、四面体形状の要素又は六面体形状の要素が用いられてもよい。粘弾性材料モデルでは、各要素について、座標系における節点座標値、要素形状及び材料特性等が定義される。構築部１０１０は、例えば、プリプロセッサと呼ばれる公知の技術に基づいた汎用ソフトウェアを用いることにより実現することができる。構築部１０１０は、節点を有する有限数の要素に分割された粘弾性材料モデルをモデル記憶部１０６０に格納する。
【００９５】
試験結果記憶部１０５０は、材料定数ｍｉ及びＡｉ、剛性割合γ０及びγｉ、弾性要素における超弾性係数Ｃ１０、Ｃ２０及びＣ３０、並びに、シフトファクターα（Ｔ）を算出するのに必要な試験結果（調和振動試験及び応力緩和試験等）のデータを記憶している。
【００９６】
材料定数計算部１０２０は、試験結果記憶部１０５０から試験結果のデータを取得して、粘弾性材料モデルの材料定数ｍｉ及びＡｉ、剛性割合γ０及びγｉ、並びに、弾性要素における超弾性係数Ｃ１０、Ｃ２０及びＣ３０（以下、これらをまとめて「パラメータ」と称呼する場合がある。）を上述したように算出する。材料定数計算部１０２０は、当該算出したパラメータをパラメータ記憶部１０７０に格納する。なお、パラメータ記憶部１０７０は、複数種類の粘弾性材料の解析を可能とするために、粘弾性材料ごとにパラメータを記憶してもよい。
【００９７】
ＳＦ計算部１０３０は、試験結果記憶部１０５０から、複数の環境温度にて実施された調和振動試験の試験結果データを取得する。ＳＦ計算部１０３０は、複数の環境温度のそれぞれにおける弾性率Ｇ（Ｔ）を算出する。次に、ＳＦ計算部１０３０は、複数の環境温度の中の基準温度Ｔ０の基準弾性率Ｇ（Ｔ０）を用いて、複数の環境温度のそれぞれの弾性率Ｇ（Ｔ）を標準化する。ＳＦ計算部１０３０は、図８（ｂ）に示すように、「標準化された弾性率の対数値ｌｏｇ（Ｇ（Ｔ）／Ｇ（０））」と温度Ｔとの関係を表す関数ｆ（Ｔ）を求める。ＳＦ計算部１０３０は、材料定数ｍｉと、関数ｆ（Ｔ）と、（３０）式の関係とに基いてシフトファクターα（Ｔ）を算出する。ＳＦ計算部１０３０は、算出したシフトファクターα（Ｔ）をパラメータ記憶部１０７０に格納する。なお、パラメータ記憶部１０７０は、複数種類の粘弾性材料の解析を可能とするために、粘弾性材料ごとにシフトファクターα（Ｔ）を記億してもよい。
【００９８】
応力計算部１０４０は、粘弾性材料モデルを用いて、各要素の節点における変位量、ひずみ量及び応力を算出する。応力計算部１０４０は、第１計算部１０４１と、第２計算部１０４２と、第３計算部１０４３と、第４計算部１０４４とを有する。
【００９９】
第１計算部１０４１は、粘弾性材料モデルをモデル記憶部１０６０から取得する。第１計算部１０４１は、粘弾性材料モデルに境界条件を設定し、計算ステップ毎の入力条件から各要素における節点の変位量Ｕを計算する。境界条件とは、粘弾性材料の挙動をシミュレートする上で粘弾性材料モデルに与える各種条件である。第１計算部１０４１は、各要素における剛性方程式を解くための剛性マトリックスを作成する。そして、第１計算部１０４１は、粘弾性材料モデルの全体構造を表す全体剛性マトリックスを作成する。第１計算部１０４１は、全体剛性マトリックスに対して、入力条件（例えば、既知節点の変位量及び節点力）を導入して解析処理を実行することにより、未知節点の変位量Ｕを計算する。なお、第１計算部１０４１は、例えば、ソルバーと呼ばれる公知の技術に基づいた汎用ソフトウェアを用いることにより実現することができる。
【０１００】
第２計算部１０４２は、第１計算部１０４１により得られた変位量Ｕを入力値として受け取る。第２計算部１０４２は、変位量Ｕを用いて各要素の節点におけるひずみ速度を算出する。なお、第２計算部１０４２は、計算ステップｎにおいて得られたひずみ量を用いることにより、次の計算ステップであるステップｎ＋１におけるひずみ速度を算出する。
【０１０１】
第２計算部１０４２は、計算ステップｎにおける節点の全変位量ψｎと、第１計算部１０４１により得られた変位量Ｕから、計算ステップｎ＋１における全変位量をψｎ＋１＝ψｎ＋Ｕにより算出する。ここで、下記の（３１）式乃至（３６）式の関係が成り立つ。
【０１０２】
【数３８】
0006988599.tif 000039
【０１０３】
第２計算部１０４２は、全変位量ψｎ＋１、及び、（３２）式乃至（３６）式の関係から、計算ステップｎ＋１における変形勾配テンソルＦｎ＋１、体積変化率（ヤコビアン）Ｊｎ＋１、右Ｃａｕｃｈｙ−ＧｒｅｅｎテンソルＣｎ＋１、偏差成分を除去した修正変形勾配テンソル
【数３９】
0006988599.tif 000040
及び、修正右Ｃａｕｃｈｙ−Ｇｒｅｅｎテンソル
【数４０】
0006988599.tif 000041
を算出する。
【０１０４】
なお、（３２）式におけるＤは、変形勾配テンソルＦを得るための微分演算子である。
【０１０５】
第２計算部１０４２は、（３６）式により得られた修正右Ｃａｕｃｈｙ−Ｇｒｅｅｎテンソルと、計算ステップｎにおけるひずみ量から、下記（３７）式及び（３８）式により計算ステップｎ＋１におけるひずみ量とひずみ速度とを算出する。なお、Δｔｎは、各計算ステップに対応して離散化された時間ステップである。
【０１０６】
【数４１】
0006988599.tif 000042
【０１０７】
第３計算部１０４３は、第２計算部１０４２により得られたひずみ速度を入力値として受け取る。更に、第３計算部１０４３は、材料定数ｍｉ及びＡｉをパラメータ記憶部１０７０から取得する。更に、第３計算部１０４３は、環境温度Ｔに対応するシフトファクターα（Ｔ）をパラメータ記憶部１０７０から取得する。第３計算部１０４３は、材料定数ｍｉ及びＡｉ、シフトファクターα（Ｔ）、並びに、ひずみ速度から、応力の計算に必要となる緩和時間τｉを算出する。具体的には、第３計算部１０４３は、（２３）式及び（２４）式により緩和時間τｉを計算する。
【０１０８】
第４計算部１０４４は、第３計算部１０４３により得られた緩和時間τｉを入力値として受け取る。第４計算部１０４４は、緩和時間τｉを用いて節点におけるＫｉｒｃｈｈｏｆｆ応力を求める。第４計算部１０４４は、まず、（１１）式を用いて計算ステップｎ＋１におけるＫｉｒｃｈｈｏｆｆ弾性応力を求める。次に、第４計算部１０４４は、（２０）式及び（２１）式を用いて中間関数Ｈ（ｉ）ｎ＋１を求めるとともに、（１２）式を用いて
【数４２】
0006988599.tif 000043
を求める。
【０１０９】
次に、第４計算部１０４４は、（１４）式及び（２２）式を用いてＫｉｒｃｈｈｏｆｆ応力を算出する。第４計算部１０４４は、各計算ステップで得られた中間関数及びＫｉｒｃｈｈｏｆｆ応力の値を計算結果記憶部１０８０に記憶する。応力計算部１０４０は、上述した処理を所定の計算ステップ数だけ繰り返し実行する。最終的に、第４計算部１０４４は、計算結果記憶部１０８０に記憶されている計算結果から、粘弾性材料モデルにおける応力−ひずみ曲線及び応力緩和曲線を作成する。第４計算部１０４４は、応力−ひずみ曲線及び応力緩和曲線を出力装置９３０に表示させる。
【０１１０】
＜実際の作動＞
ＣＰＵは、所定時間が経過する毎に、図１１により示した「解析ルーチン」を実行するようになっている。所定のタイミングになると、ＣＰＵは、ステップ１１００から図１１のルーチンを開始して以下のステップ１１０１乃至ステップ１１０８の処理を順に行い、その後、ステップ１１０９に進む。
【０１１１】
ステップ１１０１：ＣＰＵは、解析対象とする粘弾性材料の粘弾性材料モデルを構築する。
ステップ１１０２：ＣＰＵは、粘弾性材料モデルの材料定数ｍｉ及びＡｉ、剛性割合γ０及びγｉ、並びに、弾性要素における超弾性係数Ｃ１０、Ｃ２０及びＣ３０を上述したように算出する。
ステップ１１０３：ＣＰＵは、粘弾性材料モデルのシフトファクターα（Ｔ）を上述したように算出する。
【０１１２】
ステップ１１０４：ＣＰＵは、計算ステップ数ｎを「１」に設定する。
ステップ１１０５：ＣＰＵは、上述したように、粘弾性材料モデルの各要素の節点における変位量Ｕを算出する。
ステップ１１０６：ＣＰＵは、上述したように、変位量Ｕを用いて各要素の節点におけるひずみ速度を算出する。
ステップ１１０７：ＣＰＵは、上述したように、材料定数ｍｉ及びＡｉ、シフトファクターα（Ｔ）、並びに、ひずみ速度から緩和時間τｉを算出する。
ステップ１１０８：ＣＰＵは、上述したように、緩和時間τｉを用いて節点におけるＫｉｒｃｈｈｏｆｆ応力を算出する。
【０１１３】
ＣＰＵは、ステップ１１０９に進むと、所定の終了条件が成立しているか否かを判定する。終了条件は、上述のステップ１１０５乃至ステップ１１０８の処理が予め決められた所定の計算ステップ数Ｎｔｈだけ実行されたとき（即ち、ｎ＝Ｎｔｈであるとき）に成立する。
【０１１４】
終了条件が成立していない場合、ＣＰＵは、そのステップ１１０９にて「Ｎｏ」と判定してステップ１１１０に進み、計算ステップ数ｎをインクリメントする。その後、ＣＰＵは、ステップ１１０５に戻る。
【０１１５】
終了条件が成立している場合、ＣＰＵは、そのステップ１１０９にて「Ｙｅｓ」と判定してステップ１１１１に進む。ステップ１１１１にて、ＣＰＵは、計算結果から、粘弾性材料モデルにおける応力−ひずみ曲線及び応力緩和曲線を作成して、これらを出力装置９３０に表示させる。その後、ＣＰＵは、ステップ１１９５に進んで本ルーチンを一旦終了する。
【０１１６】
＜実験結果＞
図１２（ａ）は、ゴム材料に対して調和振動試験を行った際の温度と弾性率との関係を表すグラフである。なお、調和振動試験の条件は以下の通りである。
振幅：０．２μｍ
周波数：１Ｈｚ
予ひずみ：０．２５Ｎ
【０１１７】
図１２（ｂ）は、本実施装置において温度と弾性率の関係を算出した結果を表すグラフである。本図に示すように、本実施装置は、環境温度の変化に対するゴム材料の弾性率の挙動を再現することができる。
【０１１８】
図１３（ａ）は、上記と同じ条件で調和振動試験を行った際の「温度」と「応力とひずみの位相差ｔａｎδ」との関係を表すグラフである。図１３（ｂ）は、本実施装置において「温度」と「応力とひずみの位相差ｔａｎδ」との関係を算出した結果を表すグラフである。本実施装置は、環境温度の変化に対する「応力とひずみの位相差」の挙動も再現することができる。
【０１１９】
本実施装置では、粘弾性要素の減衰特性を表す緩和時間τｉが、温度−時間換算則のシフトファクターの冪関数として定義される。従って、本実施装置は、粘弾性材料の環境温度に対する依存性を再現することができる。これにより、粘弾性材料の特性（例えば、応力−ひずみ特性）の予測精度を高めることができる。
【０１２０】
なお、本発明は上記実施形態に限定されることはなく、本発明の範囲内において種々の変形例を採用することができる。
【０１２１】
温度−時間換算則のシフトファクターα（Ｔ）の算出方法は、上記の例に限定されない。ＳＦ計算部１０３０は、以下の流れでシフトファクターα（Ｔ）求めてもよい。例えば、複数の環境温度（Ｔ＝Ｔ０、Ｔ１、Ｔ２、…、ＴＮ）にて粘弾性材料の試験体に対して応力緩和試験を実施する。このとき、予ひずみは一定とし、且つ、ひずみ速度も一定とする。
【０１２２】
図１４（ａ）に示すように、ＳＦ計算部１０３０は、応力緩和試験の試験結果データから、複数の環境温度のそれぞれについて応力緩和グラフを作成する。ここで、複数の環境温度のうち、Ｔ０を基準温度として用いる。次に、ＳＦ計算部１０３０は、基準温度Ｔ０の最大応力σ０を用いて、複数の環境温度のそれぞれの応力緩和グラフを標準化する。そして、図１４（ｂ）に示すように、ＳＦ計算部１０３０は、複数の環境温度のそれぞれにおける緩和時間（τ０、τ１、τ２、…）を算出する。
【０１２３】
次に、図１４（ｃ）に示すように、ＳＦ計算部１０３０は、緩和時間の対数値と温度との関係を表すグラフを作成する。次に、ＳＦ計算部１０３０は、基準温度Ｔ０の緩和時間τ０を用いて、複数の環境温度のそれぞれにおける緩和時間の対数値を標準化する。そして、図１４（ｄ）に示すように、ＳＦ計算部１０３０は、温度Ｔを横軸とし、「標準化された緩和時間の対数値ｌｏｇ（τ（Ｔ）／τ０）」をプロットした散布図を作成する。ＳＦ計算部１０３０は、最小二乗法により回帰関数ｆ（Ｔ）を求める。回帰関数ｆ（Ｔ）は、回帰直線又は回帰曲線の何れであってもよい。求められた関数ｆ（Ｔ）とシフトファクターα（Ｔ）との関係について、（３０）式が成り立つ。従って、ＳＦ計算部１０３０は、材料定数ｍｉと、関数ｆ（Ｔ）と、（３０）式の関係とに基いてシフトファクターα（Ｔ）を算出する。ＳＦ計算部１０３０は、算出したシフトファクターα（Ｔ）をパラメータ記憶部１０７０に格納する。
【０１２４】
なお、一定ひずみ速度試験の結果から材料定数Ａｉを算出してもよい。（２８）式を変形することにより、下記（３９）式が得られる。
【０１２５】
【数４３】
0006988599.tif 000044
【０１２６】
ここで、分子のωｉεは調和振動試験における最大ひずみ速度であることから、一定ひずみ速度試験におけるひずみ速度Ｖｎｏｍに対応付けることができる。また、分母の振幅εは調和振動試験における最大ひずみ量であることから、一定ひずみ速度試験における測定ひずみＥ＊に対応付けることができる。従って、一定ひずみ速度試験におけるひずみ速度Ｖｎｏｍと測定ひずみＥ＊を用いて、（３９）式は、下記（４０）式のように書き換えることができる。
【０１２７】
【数４４】
0006988599.tif 000045
【０１２８】
従って、式（４０）を用いて一定ひずみ速度試験における試験結果より材料定数Ａｉを算出することができる。
【符号の説明】
【０１２９】
９１０…情報処理装置、９２０…入力装置、９３０…出力装置、９１１…ＣＰＵ、９１２…ＲＡＭ、９１３…ＲＯＭ、９１４…ＨＤＤ、９１５…Ｉ／Ｏインタフェース、１０１０…構築部、１０２０…材料定数計算部、１０３０…シフトファクター（ＳＦ）計算部、１０４０…応力計算部、１０４１…第１計算部、１０４２…第２計算部、１０４３…第３計算部、１０４４…第４計算部、１０５０…試験結果記憶部、１０６０…モデル記憶部、１０７０…パラメータ記憶部、１０８０…計算結果記憶部。
【図１】
0006988599.tif 000048
【図２】
0006988599.tif 000049
【図３】
0006988599.tif 000050
【図４】
0006988599.tif 000051
【図５】
0006988599.tif 000052
【図６】
0006988599.tif 000053
【図７】
0006988599.tif 000054
【図８】
0006988599.tif 000055
【図９】
0006988599.tif 000056
【図１０】
0006988599.tif 000057
【図１１】
0006988599.tif 000058
【図１２】
0006988599.tif 000059
【図１３】
0006988599.tif 000060
【図１４】
0006988599.tif 000061