特許 6048358 解析装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6048358

(24)【登録日】2016年12月2日

(45)【発行日】2016年12月21日

(54)【発明の名称】解析装置

(51)【国際特許分類】

   G01L 1/00 20060101AFI20161212BHJP

   G01N 3/00 20060101ALI20161212BHJP

   G06F 17/50 20060101ALI20161212BHJP

【ＦＩ】

   G01L1/00 M

   G01N3/00 K

   G06F17/50 612H

【請求項の数】4

【全頁数】19

(21)【出願番号】特願2013-211369(P2013-211369)

(22)【出願日】2013年10月8日

(65)【公開番号】特開2015-75383(P2015-75383A)

(43)【公開日】2015年4月20日

【審査請求日】2015年12月17日

(73)【特許権者】

【識別番号】000003207

【氏名又は名称】トヨタ自動車株式会社

(74)【代理人】

【識別番号】100105924

【弁理士】

【氏名又は名称】森下 賢樹

(74)【代理人】

【識別番号】100109047

【弁理士】

【氏名又は名称】村田 雄祐

(74)【代理人】

【識別番号】100109081

【弁理士】

【氏名又は名称】三木 友由

(72)【発明者】

【氏名】ファン ビン ロン

【審査官】 森 雅之

(56)【参考文献】

【文献】 特許第４８５２６２６（ＪＰ，Ｂ２）

【文献】 特許第４５９４０４３（ＪＰ，Ｂ２）

【文献】 特許第３４６６５８３（ＪＰ，Ｂ２）

【文献】 特許第３４６６５８４（ＪＰ，Ｂ２）

【文献】 特許第３４６６５８５（ＪＰ，Ｂ２）

【文献】 特開２０１４−１００４７（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｌ １／

Ｇ０１Ｎ ３／

Ｇ０６Ｆ １７／５０

(57)【特許請求の範囲】

【請求項1】

弾性要素と粘弾性要素とが並列配置された非線形粘弾性材料構成則に基づいて、粘弾性材料の応力−ひずみ特性を解析する装置であって、
節点を境界とする有限数の要素に分割された粘弾性材料モデルに条件を設定して前記節点の変位量を計算する第１計算部と、
前記変位量を用いて前記節点におけるひずみ速度を計算する第２計算部と、
前記ひずみ速度を底とする冪乗の値に比例する値を前記粘弾性要素の緩和時間として計算する第３計算部と、
前記ひずみ速度より計算された緩和時間を用いて前記節点における応力を計算する第４計算部と、
を備えることを特徴とする解析装置。

【請求項2】

前記緩和時間を得るために解析対象とする粘弾性材料に応じて定められる定数として、前記ひずみ速度を底とする冪数ｍｉと、前記冪乗の値に掛け合わされる比例定数１／Ａｉと、を記憶する記憶部をさらに備え、
前記第３計算部は、前記記憶部に記憶される冪数ｍｉ及び比例定数１／Ａｉを読み出して前記緩和時間を計算することを特徴とする請求項１に記載の解析装置。

【請求項3】

前記記憶部は、
解析対象とする粘弾性材料を用いた調和振動試験における絶対弾性率Ｇと振幅εの測定結果を入力として以下の（１）式により得られる冪数ｍｉを記憶し、
前記調和振動試験における加振周波数ωと振幅εの測定結果および前記冪数ｍｉの値を入力として以下の（２）式により得られる比例定数１／Ａｉを記憶することを特徴とする請求項２に記載の解析装置。

【数41】

【請求項4】

弾性要素と粘弾性要素とが並列配置された非線形粘弾性材料構成則に基づいて、粘弾性材料の応力−ひずみ特性を得るためのプログラムであって、
コンピュータに、
節点を境界とする有限数の要素に分割された粘弾性材料モデルに条件を設定して前記節点の変位量を計算する機能と、
前記変位量を用いて前記節点におけるひずみ速度を計算する機能と、
前記ひずみ速度を底とする冪乗の値に比例する値を前記粘弾性要素の緩和時間として計算する機能と、
前記ひずみ速度より計算された緩和時間を用いて前記節点における応力を計算する機能と、
を実現させるためのプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、解析装置に関し、特に、有限要素法により粘弾性材料の特性を解析する装置に関する。

【背景技術】

【0002】

ゴム材料などの粘弾性材料の特性を解析するためのモデルとして、弾性要素と粘弾性要素を並列配置した粘弾性材料構成則が用いられる。例えば、特許文献１や非特許文献１では、粘弾性要素の減衰特性を示す緩和時間を定数とすることにより、ゴム弾性および粘弾性を加味した応力とひずみの相関関係を解析する方法が示されている。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特開２０１０−２５６２９３号公報

【非特許文献】

【0004】

【非特許文献1】J.C. Simo and T.J.R Hughes著「Computational Inelasticity」 Interdisciplinary Applied Mathematics volume 7, Springer Verlag(1998)

【発明の概要】

【発明が解決しようとする課題】

【0005】

一般に、ゴム材料に対する調和振動試験では加振する振幅の大きさによって応力とひずみの相関関係が変化することとなるため、振幅依存性を再現することのできる解析モデルを用いることが望ましい。

【0006】

本発明は、こうした状況に鑑みなされたものであり、振幅依存性に対する再現性を高めた粘弾性材料の解析技術を提供することを目的とする。

【課題を解決するための手段】

【0007】

上記課題を解決するために、本発明のある態様の解析装置は、弾性要素と粘弾性要素とが並列配置された非線形粘弾性材料構成則に基づいて、粘弾性材料の応力−ひずみ特性を解析する装置であって、節点を境界とする有限数の要素に分割された粘弾性材料モデルに条件を設定して節点の変位量を計算する第１計算部と、変位量を用いて節点におけるひずみ速度を計算する第２計算部と、ひずみ速度を底とする冪乗の値に比例する値を粘弾性要素の緩和時間として計算する第３計算部と、ひずみ速度より計算された緩和時間を用いて節点における応力を計算する第４計算部と、を備える。

【0008】

この態様によると、解析装置は、ひずみ速度の冪関数として表される緩和時間を用いて粘弾性材料の応力−ひずみ特性を解析することができる。これにより、粘弾性材料に加振される振幅量に応じた応力−ひずみ特性の再現性を高めることができる。

【発明の効果】

【0009】

本発明によれば、粘弾性要素の緩和時間として、ひずみ速度を底とする冪乗の値を用いることにより、振幅依存性に対する再現性を高めた解析技術を提供することができる。

【図面の簡単な説明】

【0010】

【図1】本発明の実施の形態に係る粘弾性材料構成則を模式的に示す図である。

【図2】調和振動試験における試験条件を示す図である。

【図3】冪数の算出方法を模式的に示すグラフである。

【図4】応力緩和試験における応力緩和曲線を模式的に示す図である。

【図5】入力ひずみを変化させた場合における応力緩和曲線を示す図である。

【図6】ひずみと応力の静特性曲線を示す図である。

【図7】一定ひずみ速度試験におけるひずみと応力の関係曲線および静特性曲線を示す図である。

【図8】実施の形態に係る解析装置の機能構成を示すブロック図である。

【図9】解析装置により得られる応力−ひずみ曲線の一例を示すグラフである。

【図10】解析装置により得られる応力緩和曲線の一例を示すグラフである。

【図11】解析装置により得られる周波数特性の一例を示すグラフである。

【図12】解析装置の動作を示すフローチャートである。

【発明を実施するための形態】

【0011】

以下、図面を参照して本発明の実施の形態について説明する。

【0012】

図１は、本発明の実施の形態に係る粘弾性材料構成則を模式的に示す図である。本図に示す粘弾性材料モデルでは、弾性率Ｇ０の弾性要素と、弾性率Ｇｉ（ｉ＝１〜Ｎ、Ｎは自然数）の弾性要素と粘性係数ηｉの粘性要素とが直列接続された複数の粘弾性要素とが並列接続される。弾性要素と粘弾性要素とを組み合わせた本モデルは、タイヤやゴムブッシュなどのゴム部品の動特性を表すためのモデルとして用いられる。

【0013】

ここで、並列接続される弾性要素および粘弾性要素における剛性割合をγｉ（ｉ＝０〜Ｎ）、粘弾性要素における緩和時間をτｉ＝ηｉ／Ｇｉとすると、この粘弾性材料モデルにおける時刻ｔでの応力Ｓは、下記の式（１）、（２）により表される。

【0014】

【数1】

【0015】

式（１）において、Ｓは、第２Ｐｉｏｌａ−Ｋｉｒｃｈｈｏｆｆ応力を示し、上付き記号を有するＳ○は、粘性力成分を除去した弾性成分のみの応力であることを示す。Ｊは、粘弾性材料の体積変化率を示す。体積変化率Ｊは、ある物質点における変形前と変形後の位置の線形変換の関係を示す変形勾配テンソルＦのデターミナント（ｄｅｔ）を用いて、Ｊ＝ｄｅｔ［Ｆ］により表される。演算子ＤＥＶは、右Ｃａｕｃｈｙ−ＧｒｅｅｎテンソルＣ＝ＦＴ・Ｆを用いて、下記式（３）で表される。なお、式（３）における［・］は、演算子ＤＥＶの演算対象となる変数を表す。

【0016】

【数2】

【0017】

また、Ｑｉは、それぞれの粘弾性要素における粘性力を示し、式（１）におけるＱｉは式（２）に示す発展方程式により表される。式（２）において、

【0018】

【数3】

【0019】

は、超弾性体におけるひずみポテンシャルエネルギーの偏差成分を表す。また、

【0020】

【数4】

【0021】

は、体積成分を除去した修正右Ｃａｕｃｈｙ−Ｇｒｅｅｎテンソルであり、式（４）で表される。

【0022】

【数5】

【0023】

上記の式（１）、（２）により表される第２Ｐｉｏｌａ−Ｋｉｒｃｈｈｏｆｆ応力Ｓは、積分因子ｅｘｐ（ｔ／τｉ）を用いた畳み込み積分形式を積分することで、下記の式（５）に示す積分形で表すことができる。

【0024】

【数6】

【0025】

なお、式（５）における

【0026】

【数7】

【0027】

は、ひずみポテンシャルエネルギーの体積成分である。また、ｇ（ｔ）は、緩和関数であり、上記の式（６）により表される。

【0028】

ここで、時間ｔの関数として表される式（５）を下記の式（７）に変形することで、時刻ｔｎ＋１における第２Ｐｉｏｌａ−Ｋｉｒｃｈｈｏｆｆ応力Ｓｎ＋１を得ることができる。なお、計算ステップｎにおける関数を（・）ｎ、計算ステップｎ＋１における関数を（・）ｎ＋１と表記している。

【0029】

【数8】

【0030】

式（７）における

【0031】

【数9】

【0032】

は、中点定理を用いて［ｔｎ，ｔｎ＋１］の時間間隔で積分した近似解として得られる中間関数であり、下記式（８）で表される。

【0033】

【数10】

【0034】

ここで、

【0035】

【数11】

【0036】

は、下記式（９）、（１０）で定義される。

【0037】

【数12】

【0038】

なお、式（９）は、Ｋｉｒｃｈｈｏｆｆ弾性応力

【0039】

【数13】

【0040】

を、下記式（１１）

【0041】

【数14】

【0042】

により定義することにより、下記式（１２）により表すこともできる。

【0043】

【数15】

【0044】

また、Ｋｉｒｃｈｈｏｆｆ応力テンソル

【0045】

【数16】

【0046】

は、第２Ｐｉｏｌａ−Ｋｉｒｃｈｈｏｆｆ応力Ｓｎ＋１を用いて下記の式（１３）により表すことができる。

【0047】

【数17】

【0048】

したがって、式（７）、（１３）より、Ｋｉｒｃｈｈｏｆｆ応力テンソルを下記の式（１４）により表すことができる。

【0049】

【数18】

【0050】

なお、式（１４）における演算子ｄｅｖは、下記の式（１５）により定義される。なお、式（１５）における（・）は、演算子ＤＥＶの演算対象となる変数を表す。

【0051】

【数19】

【0052】

また、式（１４）における緩和関数ｇ^＊は、下記の式（１６）により定義される。

【0053】

【数20】

【0054】

以上より、中間関数として

【0055】

【数21】

【0056】

を計算ステップ毎に保持することにより、Ｋｉｒｃｈｈｏｆｆ応力テンソルの値を数値解析により得ることができる。なお、上述の式変形の詳細については、例えば、非特許文献１に記載される。

【0057】

本実施の形態では、図１に示す粘弾性モデルにおいて、粘弾性要素の減衰特性を表す緩和時間τｉをひずみ速度に依存した変数とすることにより、振幅依存性の再現性を高めた粘弾性材料の解析方法を示す。一方、非特許文献１などに示される粘弾性モデル（Ｓｉｍｏモデル）では、緩和時間τｉを定数として計算することにより、ゴム材料における様々な動特性を表すことができるものの、調和振動試験などで見られる振幅依存性の予測精度が高くないという課題があった。本発明者は、振幅依存性の精度を高める手法として、緩和時間τｉをひずみ速度の冪関数とすることが有効であることを見出した。本実施の形態では、下記式（１７）に示す関係式を満たす緩和時間τｉを用いることにより、振幅依存性の再現性を高めた解析手法を提案する。

【0058】

【数22】

【0059】

式（１７）において、

【0060】

【数23】

【0061】

は、偏差成分を除去したＧｒｅｅｎ−Ｌａｇｒａｎｇｅひずみテンソルであり、修正右Ｃａｕｃｈｙ−Ｇｒｅｅｎテンソル

【0062】

【数24】

【0063】

を用いて下記の式（１８）で表される。

【0064】

【数25】

【0065】

また、

【0066】

【数26】

【0067】

は、ひずみ速度であり、ひずみテンソルの時間微分値を表す。また、

【0068】

【数27】

【0069】

は、ひずみ速度の大きさを表し、３次元のひずみテンソルを用いる場合には、下記式（１９）により表される。

【0070】

【数28】

【0071】

なお、冪数ｍｉおよび比例定数１／Ａｉは、本実施の形態で導入する値であり、調和振動試験の試験結果を用いて後述する関係式により導出される。

【0072】

つづいて、式（１７）に示す緩和時間τｉの導入に際し、数値解析上必要となる式の変形について述べる。上述の式（１４）を用いてＫｉｒｃｈｈｏｆｆ応力を算出するためには、緩和時間τｉを含む式を変形する必要がある。具体的には、式（８）を下記の式（１８）に、式（９）を下記の式（１９）に変形すればよい。

【0073】

【数29】

【0074】

ここで、緩和時間τｉの関数は、下記の式（２０）、（２１）により定義される。

【数30】

【0075】

なお、式（１８）の右辺第２項および式（１９）における緩和時間τｉとして式（２１）を用いるのは、中点定理を用いて［ｔｎ，ｔｎ＋１］の時間間隔で積分した近似解を用いたためである。

【0076】

次に、本実施の形態に示す粘弾性材料モデルにおける材料定数の同定方法について述べる。まず、冪数ｍｉおよび比例定数１／Ａｉの同定方法について示した後に、並列配置される弾性要素の剛性割合γ０および弾性率Ｇ０の同定方法を示し、最後に粘弾性要素の剛性割合γｉの同定方法を示す。

【0077】

図２は、調和振動試験における試験条件を示す図である。本図は、横軸を時間ｔとして、予ひずみＥｐｒｅ、振幅ε、周波数ωのとした場合における粘弾性材料の試験体に加振されるひずみＥ＝Ｅｐｒｅ＋εｓｉｎ（ωｔ）のグラフを示す。このようなひずみＥを加振した場合における試験体の動的弾性率Ｇを測定することにより、下記式（２２）の関係を用いて冪数ｍｉを導出することができる。

【0078】

【数31】

【0079】

なお、式（２２）は、図２に示す調和振動試験の入力に対して、式（２）に式（１７）を代入した発展方程式である下記式（２３）を解くことにより得られる。

【0080】

【数32】

【0081】

図３は、冪数ｍｉの算出方法を模式的に示すグラフである。本図は、調和振動試験において加振する振幅εの対数値ｌｏｇ（ε）を横軸として、予ひずみＥｐｒｅおよび周波数ωを変化させたときの試験結果値である動的弾性率Ｇの対数値ｌｏｇ（Ｇ）をプロットしたものである。グラフ上の近似線Ｌ１は、予ひずみＥｐｒｅ１と周波数ω１に対して、振幅の値ε１〜ε３を変化させたときのプロット値に対応する傾きｌ１の近似線である。同様に、近似線Ｌ２、Ｌ３は、予ひずみＥｐｒｅ２、Ｅｐｒｅ３と周波数ω２、ω３に対して、振幅の値ε１〜ε３を変化させたときのプロット値に対応する傾きｌ２、ｌ３の近似線である。式（２２）の両辺に対数をとることにより、下記式（２４）を得ることができることから、この傾きｌｉから冪数ｍｉを得ることができる。

【0082】

【数33】

【0083】

式（２４）より、冪数はｍｉ＝−１−ｌｉの関係式により得ることができる。なお、式（２４）におけるαは、図３に示す近似線の切片を示す。これにより、調和振動試験における加振周波数ωｉに対応する冪数ｍｉを得ることができる。なお、図３に示す振幅εの値は例示であり、冪数ｍｉを得るために４以上の振幅εの値を用いてもよいし、振幅εの値を２としてよい。

【0084】

なお、冪数ｍｉの値が加振周波数ωｉによってそれほど変化しない場合には、それぞれの周波数に対応するｍｉを得る代わりに、それぞれのｍｉの平均値として得られる共通の冪数ｍを用いることとしてもよい。共通の冪数ｍを用いることにより、周波数ごとに異なる冪数ｍｉを用いる場合と比べて、計算負荷を減らすことができる。

【0085】

また、比例定数１／Ａｉは、上記方法で得られた冪数ｍｉを用いて、下記式（２５）の関係から得ることができる。上述の冪数ｍｉにおける導出方法と同様、式（２５）における周波数ωｉは、調和振動試験における加振周波数であり、この周波数を変化させることにより周波数成分ωｉが異なるそれぞれの粘弾性要素に対応する比例定数１／Ａｉを得ることができる。

【0086】

【数34】

【0087】

つづいて、弾性要素の剛性割合をγ０の算出方法について述べる。図４は、応力緩和試験における応力緩和曲線を模式的に示す図である。本図は、入力ひずみＥｐｒｅを一定にした場合における応力Ｑの時間変化を測定した試験結果を示している。この試験結果に示す時刻ｔ＝０の瞬間応力Ｑ０と、時刻ｔ＝∞における緩和応力Ｑ∞とから、弾性要素の剛性割合をγ０＝Ｑ∞／Ｑ０により得ることができる。

【0088】

つづいて、弾性要素の弾性率Ｇ０の算出方法について述べる。本実施の形態では、超弾性の材料モデルであるＹｅｏｈ材料モデルに定められる超弾性係数Ｃ１０、Ｃ２０、Ｃ３０を定めることにより弾性要素の弾性率を決定する。これらの超弾性係数は、図４に示した応力緩和試験において、複数の入力ひずみＥｐｒｅに対して緩和応力Ｑ∞を得ることにより導出することができる。

【0089】

図５は、入力ひずみＥｋを変化させた場合における応力緩和曲線を示す図である。図６は、ひずみＥと応力Ｑの静特性曲線ＱＲを示す図である。図５に示すように、値の異なる複数の負荷ひずみＥ１〜ＥＮに対して応力緩和試験を実施し、各負荷ひずみＥｋに対応する緩和応力ＱＲ（Ｅ１）〜ＱＲ（ＥＮ）を得る。この関係式をひずみ−応力曲線としてグラフ化することにより、図６に示す静特性曲線ＱＲが得られる。この静特性曲線を近似することのできる超弾性係数Ｃ１０、Ｃ２０、Ｃ３０を公知の技術を用いて定めることにより、弾性要素の弾性率を決定することができる。

【0090】

最後に、粘弾性要素の剛性割合γｉの算出方法について述べる。図７は、一定ひずみ速度試験におけるひずみＥと応力Ｑの関係曲線Ｑｉおよび静特性曲線ＱＲを示す図である。静特性曲線ＱＲは、図６に示す静特性曲線と同じである。関係曲線Ｑ１〜ＱＮは、一定ひずみ速度試験においてひずみ速度ＶｉをＶ１〜ＶＮに変化させた場合に対応するひずみＥと応力Ｑの試験結果を示す。ここで、ｉ番目の関係曲線Ｑｉは、ｉ番目の粘弾性要素の動特性に対応し、加振周波数ωｉから算出した比例定数１／Ａｉの粘弾性要素に対応する。このような対応関係とするためには、調和振動試験における加振周波数ωｉと振幅εを用いて、ひずみ速度Ｖｉ＝ωｉεの関係を満たす値を定めればよい。

【0091】

このとき、各粘弾性要素の剛性割合γｉは、図６に示すひずみＥと応力Ｑの関係曲線Ｑｉにより囲まれる面積Ｚｉから得ることができる。静特性曲線ＱＲにより囲まれる面積をＺ０とし、ｉ−１番目からｉ番目にひずみ速度が増加することによる面積の増加分をＺｉとすると、剛性割合はγｉ＝γ０×Ｚｉ／Ｚ０により同定される。

【0092】

以上の方法により、本実施の形態の粘弾性材料モデルにおける材料定数ｍｉ、Ａｉ、γ０、Ｇ０、γｉを同定することができる。つまり、本実施の形態においては、材料試験の結果から材料定数を一意に決めることができる。

【0093】

なお、その他の材料定数の同定方法として、最適化計算による同定方法が挙げられるが、最適化計算による同定方法では計算が収束せず解が得られない、得られた最適解が局所解であるために再現精度が低下するなどの課題があった。一方、本実施の形態では、最適化計算を必要としないため、解が得られないという課題を解消することができる。また、材料試験の結果から材料定数を同定するため、解析の再現精度を高めることができる。

【0094】

次に、上述の式（２３）に示す発展方程式を利用した有限要素解析法について述べる。図８は、実施の形態に係る解析装置１００の機能構成を示すブロック図である。解析装置１００は、構築部１０と、計算部２０と、記憶部３０と、出力部４０を備える。

【0095】

本明細書のブロック図において示される各ブロックは、ハードウェア的には、コンピュータのＣＰＵをはじめとする素子や機械装置で実現でき、ソフトウェア的にはコンピュータプログラム等によって実現されるが、ここでは、それらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックはハードウェア、ソフトウェアの組み合わせによっていろいろなかたちで実現できることは、当業者に理解されるところである。

【0096】

構築部１０は、解析しようとする粘弾性材料を有限個の小さな要素に置き換えて、計算部２０の計算対象となる粘弾性材料モデルを構築する。粘弾性材料モデルを構成する各要素は数値解析が可能に定義され、具体的には、各要素について、座標系における節点座標値、要素形状、材料特性などが定義される。２次元モデルを対象とする場合には、各要素として三角形状を有する３節点要素や、四角形状を有する４節点要素を用いればよい。また、３次元モデルを対象とする場合には、四面体形状を有する４節点要素や、六面体形状を有する６節点要素などを用いればよい。構築部１０は、例えば、プリプロセッサと呼ばれる公知の技術に基づいた汎用ソフトウェアを用いることにより実現することができる。

【0097】

記憶部３０は、計算部２０による処理に必要な入力値を記憶するとともに、計算部２０より計算された各計算ステップにおける出力値を記憶する。記憶部３０は、入力値として粘弾性材料モデルの材料定数である冪数ｍｉ、比例定数１／Ａｉ、剛性割合γ０、γｉ、弾性要素における超弾性係数Ｃ１０、Ｃ２０、Ｃ３０を少なくとも記憶する。これらの材料定数は、解析対象とする粘弾性材料の特性に応じて、上述の材料定数同定方法により設定される値である。記憶部３０には、例えば、複数種類の粘弾性材料の解析を可能とするために複数種類の材料定数が記憶されており、解析対象とする粘弾性材料の種類に応じて、対応する材料定数が計算部２０により読み出される。

【0098】

計算部２０は、構築部１０により構築された粘弾性材料モデルを用いて、各要素の節点における変位量、ひずみ量、応力などの値を計算ステップ毎に算出する。計算部２０は、所定の計算終了条件を満たすまで計算処理を繰り返し実行し、得られた値を記憶部３０に記憶させる。計算部２０は、第１計算部２１と、第２計算部２２と、第３計算部２３と、第４計算部２４を有する。

【0099】

第１計算部２１は、有限数の要素に分割された粘弾性材料モデルに境界条件を設定し、計算ステップ毎の入力条件から各要素における節点の変位量Ｕを計算する。第１計算部２１は、各要素における剛性方程式を解くための剛性マトリックスを作成した後、粘弾性材料モデルの全体構造を表す全体構成マトリックスを作成する。全体剛性マトリックスに対して、入力条件となる既知節点の変位量および節点力を導入して解析処理を施すことにより、未知節点の変位量Ｕが計算される。なお、第１計算部２１は、例えば、ソルバーと呼ばれる公知の技術に基づいた汎用ソフトウェアを用いることにより実現することができる。

【0100】

第２計算部２２は、第１計算部２１により得られた変位量Ｕを入力値として各要素の節点におけるひずみ速度を計算する。なお、第２計算部２２は、計算ステップｎにおいて得られたひずみ量を用いることにより、次の計算ステップであるステップｎ＋１におけるひずみ速度を算出する。

【0101】

第２計算部２２は、計算ステップｎにおける節点の全変位量ψｎと、第１計算部２１により得られた変位量Ｕから、計算ステップｎ＋１における全変位量をψｎ＋１＝ψｎ＋Ｕにより算出する。この全変位量ψｎ＋１を用いて、下記の式（２６）−（３１）の関係により、計算ステップｎ＋１における変形勾配テンソルＦｎ＋１、体積変化率（ヤコビアン）Ｊｎ＋１、右Ｃａｕｃｈｙ−ＧｒｅｅｎテンソルＣｎ＋１、偏差成分を除去した修正変形勾配テンソル

【0102】

【数35】

【0103】

及び修正右Ｃａｕｃｈｙ−Ｇｒｅｅｎテンソル

【0104】

【数36】

【0105】

を計算する。

【0106】

【数37】

【0107】

なお、式（２７）におけるＤは、変形勾配テンソルＦを得るための微分演算子である。

【0108】

式（３１）により得られた修正右Ｃａｕｃｈｙ−Ｇｒｅｅｎテンソルと、計算ステップｎにおけるひずみ量から、下記の式（３２）、（３３）により計算ステップｎ＋１におけるひずみ量とひずみ速度を得ることができる。なお、式（１８）におけるΔｔｎは、各計算ステップに対応して離散化された時間ステップである。

【0109】

【数38】

【0110】

第３計算部２３は、第２計算部２２により得られたひずみ速度から、応力の計算に必要となる緩和時間τｉを算出する。第３計算部２３は、記憶部３０に記憶される冪数ｍｉおよび比例定数１／Ａｉを読み出し、上述の式（２０）、（２１）により緩和時間τｉを計算する。

【0111】

第４計算部２４は、ひずみ速度より計算された緩和時間τｉを用いて節点におけるＫｉｒｃｈｈｏｆｆ応力を求める。第４計算部２４は、まず、式（１１）を用いて計算ステップｎ＋１におけるＫｉｒｃｈｈｏｆｆ弾性応力を求める。次に、式（１２）、（１８）を用いて中間変数を求める。最後に、式（１４）、（１９）を用いてＫｉｒｃｈｈｏｆｆ応力を得る。なお、第４計算部２４は、各計算ステップで得られた中間変数および応力の値を記憶部３０に記憶する。

【0112】

出力部４０は、計算部２０により得られた計算結果を出力する。出力部４０は、計算結果として、粘弾性材料モデルにおける応力−ひずみ曲線や応力緩和曲線などを得るための数値データを出力する。

【0113】

図９は、解析装置１００により得られる応力−ひずみ曲線の一例を示すグラフであり、一定ひずみ速度試験の試験結果に対応する。なお、本図では、出力部４０からの出力結果である「計算値」とともに、対応する一定ひずみ速度試験で得られた「試験値」を示している。計算値Ｃ１〜Ｃ４は、それぞれ異なるひずみ速度に対応する計算結果であり、試験値Ｔ１〜Ｔ４は、計算の条件と同じひずみ速度にて一定ひずみ速度試験を行った場合の試験結果である。本実施の形態の解析装置を用いることにより、図示されるように、一定ひずみ速度試験の結果を精度良く再現する計算結果を得ることができる。

【0114】

図１０は、解析装置１００により得られる応力緩和曲線の一例を示すグラフである。なお、本図においても、出力部４０からの出力結果である「計算値」を示す実線とともに、応力緩和試験で得られた「試験値」を破線で示している。それぞれの応力緩和曲線は、異なる瞬間応力を加えた場合における計算値および試験値を示している。本実施の形態の解析装置を用いることにより、図示されるように、応力緩和試験の結果を精度良く再現する計算結果を得ることができる。

【0115】

図１１は、解析装置１００により得られる周波数特性の一例を示すグラフである。本図は、調和振動試験において加振周波数ωと振幅εを変化させた場合における絶対弾性率の値を求めた結果を示す。本図では、調和振動試験における振幅εの大きさを予ひずみに対して０．２％とした場合と、２．０％とした場合の計算結果を示している。本図に示すように、加振周波数の変化による絶対弾性率の変化を表すとともに、振幅の増大に伴う絶対弾性率の減少を表すことができた。一般に、ゴム材料などの粘弾性材料では、振幅が大きくなると絶対弾性率が減少することが知られていることから、本実施の形態における解析装置１００によれば、粘弾性材料における振幅依存性を再現できる。

【0116】

以上の構成による解析装置１００の動作を説明する。
図１２は、監視装置２の動作を示すフローチャートである。まず、解析対象とする粘弾性材料モデルを構築し（Ｓ１０）、解析対象の粘弾性材料モデルに対応した材料定数を設定する（Ｓ１２）。次に、入力条件に基づいて粘弾性材料モデルを構成する各要素の節点変位量を計算ステップ毎に算出する（Ｓ１４）。得られた設定変位量からひずみ速度を算出するとともに（Ｓ１６）、ひずみ速度から緩和時間を算出し（Ｓ１８）、ひずみ速度から得られる緩和時間を用いて応力を算出する（Ｓ２０）。計算の終了条件を満たしていない場合には、Ｓ１４〜Ｓ２０のステップを繰り返して計算ステップ毎の応力を算出する（Ｓ２２のＮ）。終了条件を満たした場合には（Ｓ２２のＹ）、得られた計算結果を出力する（Ｓ２４）。

【0117】

以上、本発明を実施形態にもとづいて説明した。本発明は上記実施形態に限定されず、種々の設計変更が可能であり、様々な変形例が可能であること、またそうした変形例も本発明の範囲にあることは、当業者に理解されるところである。

【0118】

上述の実施の形態では、調和振動試験の結果から式（２５）の関係式を用いて比例定数１／Ａｉを同定する方法を述べた。変形例においては、一定ひずみ速度試験の結果から比例定数１／Ａｉを同定することとしてもよい。以下、変形例における比例定数１／Ａｉの同定方法について述べる。

【0119】

上述の式（２５）の左辺を変形することにより、下記式（３４）が得られる。

【0120】

【数39】

【0121】

ここで、分母のωｉεは調和振動試験における最大ひずみ速度であることから、一定ひずみ速度試験におけるひずみ速度Ｖｎｏｍに対応付けることができる。また、分子の振幅εは調和振動試験における最大ひずみ量であることから、一定ひずみ速度試験における測定ひずみＥ＊に対応付けることができる。したがって、一定ひずみ速度試験におけるひずみ速度Ｖｎｏｍと測定ひずみＥ＊を用いて、式（３４）は、下記式（３５）のように書き換えることができる。

【0122】

【数40】

【0123】

したがって、変形例においては、式（３５）を用いて一定ひずみ速度試験における試験結果より比例定数１／Ａｉを同定することができる。

【符号の説明】

【0124】

１０…構築部、２０…計算部、２１…第１計算部、２２…第２計算部、２３…第３計算部、２４…第４計算部、３０…記憶部、１００…解析装置。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】