JP6604555B2 - 粘弾性材料のシミュレーション方法、構造体のシミュレーション方法、粘弾性材料のシミュレーション装置、及び、プログラム - Google Patents

Description

本発明は、粘弾性材料の力学変形を再現する粘弾性材料のシミュレーション方法、構造体のシミュレーション方法、粘弾性材料のシミュレーション装置、及び、プログラムに関する。

タイヤの転がり抵抗は、粘弾性特性を示すゴムの力学特性の影響を大きく受ける。このゴムの粘弾性特性は、転がり抵抗の他に、タイヤの騒音、振動、操縦安定性、及び耐久性などのタイヤ様々な特性にも影響を与える。上記ゴムの粘弾性特性は、ゴムの温度、ゴムの歪み速度、歪み振幅などによって変化する。このため、優れた騒音や振動特性を有する低燃費タイヤの開発を、コンピュータシミュレーションを利用して行なう場合、転がり抵抗のシミュレーション計算と、振動や騒音のシミュレーション計算を行なう。しかし、騒音や振動特性に影響を与える歪みは微小（歪みが１％未満）であり、転がり抵抗に影響を与える歪みは大きい（歪みが１０％以上）。このため、歪みが微小（歪みが１％未満）な条件から歪みが大きな条件（歪みが１０％以上）までの広い歪みの範囲で、ゴムの力学変形の粘弾性特性を、統一的に再現することが好ましい。

粘弾性材料の振幅依存性に対する再現性を高めた技術が知られている（特許文献１）。
当該技術では、弾性要素と粘弾性要素とが並列配置された非線形粘弾性材料構成則に基づいて、粘弾性材料の応力−歪み特性を解析するとき、節点を境界とする有限数の要素に分割された有限要素モデルに条件を設定して節点の変位量を計算し、この変位量を用いて節点における歪み速度を計算し、歪み速度を底とする冪乗の値に比例する値を粘弾性要素の緩和時間として計算し、歪み速度より計算された緩和時間を用いて節点における応力を計算する。

また、ゴム製品の粘弾性特性の歪依存性を各々表す非線形弾性方程式及び非線形粘弾性方程式を用いて、ゴム製品の粘弾性応答性能を、有限要素法を用いて予測する粘弾性応答性能予測方法も知られている（特許文献２）。
当該方法では、非線形弾性方程式は、貯蔵弾性率を表す方程式であり、貯蔵弾性率Ｇ’は、γをせん断歪、Ａ、Ｃ、ｎを、フィラーネットワーク構造に依存する物理定数とすると、Ｇ’＝Ａ・（γ^２＋Ｃ^２）^−ｎ／２で表される。

特開２０１５−７５３８３号公報 特許第４２９９７３３号公報

上記方法では、ゴムの変形状態に応じて粘弾性特性を定めることができるが、歪みが１％未満から１０％以上の広い歪みの範囲の変形状態におけるゴムの粘弾性特性を統一的に再現して計算することは難しい。

そこで、本発明は、ゴム等の粘弾性材料の力学変形をコンピュータで再現する際に、粘弾性材料の様々な大きさの力学変形であっても、力学変形の大きさに応じて変化する粘弾性特性を統一的に再現することができる粘弾性モデルを用いた粘弾性材料のシミュレーション方法、構造体のシミュレーション方法、粘弾性材料のシミュレーション装置、及び、プログラムを提供することを目的とする。

本発明の一態様は、粘弾性材料の力学変形をコンピュータで再現するシミュレーション方法であって、
粘弾性材料の力学変形を非線形の弾性特性で表した第１ネットワーク要素と、前記力学変形の粘弾性特性を、粘弾性パラメータを用いて表したネットワーク要素であって、前記力学変形の大きさに応じて前記粘弾性特性が変化する第２ネットワーク要素と、を並列接続した粘弾性モデルをコンピュータに作成させるステップと、
前記粘弾性モデルを用いてコンピュータに、前記粘弾性材料の前記力学変形を再現するための演算を実行させるステップと、を含み、
前記第２ネットワーク要素は、前記粘弾性パラメータの値を前記力学変形に応じて設定するために用いるスイッチング関数を含む、ことを特徴とする。

前記第２ネットワーク要素は、大きさの異なる少なくとも２つの力学変形を再現した小変形ネットワーク要素及び大変形ネットワーク要素を含み、
前記小変形ネットワーク要素と前記大変形ネットワーク要素の前記粘弾性パラメータの値は、前記スイッチング関数を用いて前記粘弾性材料の前記力学変形に応じて設定される、ことが好ましい。

前記小変形ネットワーク要素の前記粘弾性パラメータの値は、予め設定された第１の値に前記スイッチング関数の値を乗算した値であり、
前記大変形ネットワーク要素の前記粘弾性パラメータの値は、予め設定された第２の値に、前記スイッチング関数の値を１から差し引いた値を乗算した値である。

前記小変形ネットワーク要素及び前記大変形ネットワーク要素のそれぞれは、前記力学変形を粘弾性特性で表した、前記粘弾性パラメータの値が互いに異なるサブネットワーク要素を少なくとも２つ以上含む、ことが好ましい。

前記スイッチング関数は、前記力学変形における、前記粘弾性材料の歪みの情報あるいは前記粘弾性材料の変形勾配の情報に応じて前記スイッチング関数の値を定める、ことが好ましい。

前記第２ネットワーク要素の粘弾性特性は、前記粘弾性材料に作用する粘性歪みの時間緩和を表す粘性歪みの発展式を用いて表され、前記粘性歪みの発展式は、前記粘弾性材料の歪み速度の関数、前記粘弾性材料の変形速度の関数、あるいは前記粘弾性材料の全歪みと前記粘弾性材料の粘性歪みの関数である、ことが好ましい。

前記粘性歪みの発展式は、前記粘弾性材料に与えられる歪みと前記粘性歪みの差分に対する分数階微分の演算を含む、ことが好ましい。

前記粘性歪みの発展式は、分数階微分を用いて下記式で表される、ことが好ましい。
（上記式において、ｅは前記粘弾性材料の力学変形における歪みを表し、ｅ_ｉｎ,ｉは前記第２ネットワーク要素のうちｉ番目のサブネットワーク要素における前記粘性歪みを表し、Ｚ_ｉ、ｄ_ｉ、λ_ｉ、α_ｉは、ｉ番目の前記サブネットワーク要素における粘弾性パラメータであり、ｑ_ｉはｉ番目の前記サブネットワーク要素における時間ｔに依存した内部変数を表し、α_ｉは０より大きく１より小さい実数であり、Ｄは微分演算子を表す）

前記粘性歪みの発展式は、増分形式の粘性歪みの発展式である、ことが好ましい。

前記スイッチング関数を用いて値を設定する粘弾性パラメータは、前記粘性歪みの時間緩和を表す緩和時間パラメータを含む、ことが好ましい。

前記スイッチング関数は、前記力学変形における、前記粘弾性材料の歪みの情報あるいは前記粘弾性材料の変形勾配の情報に応じて前記スイッチング関数の値を定める関数であって、前記歪みの情報あるいは前記変形勾配の情報と、基準値との差が大きくなる程、前記スイッチング関数の値が小さくなる関数である、ことが好ましい。

また、本発明の他の一態様は、ゴム部材を含む構造体のシミュレーション方法であって、前記構造体のシミュレーションを行うとき、前記粘弾性材料のシミュレーション方法を前記ゴム部材に適用する、ことを特徴とする構造体のシミュレーション方法である。
前記構造体は、例えば、空気入りタイヤである。

本発明の他の一態様は、粘弾性材料の力学変形を再現するシミュレーション装置であって、
前記力学変形を非線形の弾性特性で表した第１ネットワーク要素と、前記力学変形の粘弾性特性を、粘弾性パラメータを用いて表したネットワーク要素であって、前記力学変形の大きさに応じて前記粘弾性特性が変化する第２ネットワーク要素と、を並列接続したモデルであって、前記第２ネットワーク要素は、前記粘弾性パラメータの値を、前記力学変形に応じて設定するスイッチング関数を含む、粘弾性モデルを作成するモデル作成部と、
前記粘弾性モデルを用いて、前記粘弾性材料の前記力学変形を再現するための演算を行う演算部と、を含む、ことを特徴とする。

さらに、本発明のさらに他の一態様は、粘弾性材料の力学変形をコンピュータに再現させるプログラムであって、
前記力学変形を非線形の弾性特性で表した第１ネットワーク要素と、前記力学変形の粘弾性特性を、粘弾性パラメータを用いて表したネットワーク要素であって、前記力学変形の大きさに応じて前記粘弾性特性が変化する第２ネットワーク要素と、を並列接続した粘弾性モデルをコンピュータに設定させる手順と、
前記粘弾性モデルを用いてコンピュータに、前記粘弾性材料の前記力学変形を再現するための演算を実行させる手順と、を含み、
前記第２ネットワーク要素は、前記粘弾性パラメータの値を、前記力学変形に応じて設定するスイッチング関数を含む、ことを特徴とする。

上述の粘弾性材料のシミュレーション方法、構造体のシミュレーション方法、粘弾性材料のシミュレーション装置、及び、プログラムによれば、粘弾性材料の様々な大きさの力学変形に応じて変化する粘弾性特性を統一的に再現することができる。

本実施形態の粘弾性材料のシミュレーション方法を実施するシミュレーション装置の構成の一例を示す図である。 本実施形態の粘弾性モデルの一例を説明する図である。 （ａ），（ｂ）は、本実施形態の粘弾性モデルの好ましい形態の一例を示す図である。 本実施形態の粘弾性モデルで用いるスイッチング関数ｈの一例を示す図である。 本実施形態における粘弾性パラメータの値を算出する方法を説明する図である。 本実施形態の粘弾性材料のシミュレーション方法のフローの一例を説明する図である。 本実施形態のスイッチング関数ｈを含んだ実施例の粘弾性モデルを用いた結果の一例を示す図である。 本実施形態のスイッチング関数ｈを含んだ実施例の粘弾性モデルを用いた結果の一例を示す図である。 （ａ）は、粘弾性試験の一例を説明する図であり、（ｂ）は、実施例の粘弾性モデルを用いた結果の一例を示す図である。 （ａ）は、粘弾性試験の一例を説明する図であり、（ｂ）は、実施例の粘弾性モデルを用いた結果の一例を示す図である。 （ａ）、（ｂ）は、実施例の粘弾性モデルを用いた結果の一例を示す図である。 スイッチング関数ｈを含まない比較例１の粘弾性モデルを用いた結果の一例を示す図である。 スイッチング関数ｈを含まない比較例１の粘弾性モデルを用いた結果の一例を示す図である。 （ａ）〜（ｄ）は、スイッチング関数ｈを含まない比較例１の粘弾性モデルを用いた結果の一例を示す図である。 スイッチング関数ｈを含まない比較例２の粘弾性モデルを用いた結果の一例を示す図である。 スイッチング関数ｈを含まない比較例２の粘弾性モデルを用いた結果の一例を示す図である。 （ａ）〜（ｄ）は、スイッチング関数ｈを含まない比較例２の粘弾性モデルを用いた結果の一例を示す図である。

以下、本発明の粘弾性材料のシミュレーション方法、構造体のシミュレーション方法、粘弾性材料のシミュレーション装置、及びプログラムについて添付の図面を参照しながら詳細に説明する。図１は、本実施形態の粘弾性材料のシミュレーション方法及び構造体のシミュレーション方法を実施する粘弾性材料のシミュレーション装置１０の構成の一例を示す図である。

シミュレーション装置１０は、ＣＰＵ１２、ＲＡＭ及びＲＯＭを含む記憶部１４、及び入出力部１６を含むコンピュータにより構成される。入出力部１６は、粘弾性測定試験機２６、マウスやキーボード等の入力操作系２８、及びディスプレイ・プリンタ３０に接続されている。シミュレーション装置１０は、記憶部１４に記憶されたプログラムを読みだして起動することにより、シミュレーション装置１０は後述する機能を有するソフトウェアモジュールを形成する。すなわち、プログラムを起動することによって、シミュレーション装置１０には、パラメータ設定部１８、モデル作成部２０、演算部２２、及び制御管理部２４がソフトウェアモジュールとして設けられる。パラメータ設定部１８、モデル作成部２０、演算部２２、及び制御管理部２４の機能は、本実施形態の粘弾性材料のシミュレーションを説明した後、詳細に説明する。

（本実施形態の粘弾性材料のシミュレーション）
本実施形態の粘弾性材料のシミュレーションでは、歪みが１％未満から１０％以上の広い範囲で変形する粘弾性試験の実験データに一致するように、粘弾性モデルの粘弾性パラメータの値を設定し、この値を用いて、粘弾性材料の粘弾性特性を反映した有限要素法等におけるシミュレーション用のマトリクスを作成して粘弾性材料のシミュレーション及び構造体のシミュレーションを行う。
以降、有限要素法を用いたシミュレーション方法を例に説明するが、有限要素法に限らない。例えば応力と歪みの関係を用いてシミュレーションを行うメッシュフリー法や差分法等のシミュレーションに適用することができる。

図２は、本実施形態の粘弾性材料の粘弾性特性を再現する粘弾性モデルの一例を説明する図である。図２に示すように、粘弾性モデル４０は、両端１，２の間で粘弾性特性を示すモデルである。粘弾性モデル４０は、力学変形を非線形の弾性特性で表した第１ネットワーク要素ＥＡと、力学変形の粘弾性特性を、粘弾性パラメータを用いて表したネットワーク要素であって、力学変形の大きさに応じて粘弾性特性が変化する第２ネットワーク要素ＥＢ，ＥＣと、を並列接続したモデルである。
ここで、第１ネットワーク要素ＥＡの弾性特性は、一定の値を有する材料定数で表されてもよいし、歪みエネルギー密度関数で表したムーニー−リブリンモデルやアルダ−ボイスモデル等の公知の弾性ポテンシャルを用いて超弾性体の特性として表されてもよい。すなわち、弾性ポテンシャルを用いる場合、弾性ポテンシャルによって弾性特性を表す材料定数が定式化される。

一方、第２ネットワーク要素ＥＢは、粘弾性材料の微小変形を再現する微小変形ネットワーク要素である。微小変形の歪みは、例えば、歪み振幅が１％未満の歪みであることが好ましい。第２ネットワーク要素ＥＣは、粘弾性材料の大変形を再現する大変形ネットワーク要素である。大変形の歪みは、歪み振幅が１０％以上の歪みであることが好ましい。本実施形態では、第１ネットワーク要素ＥＡと、第２ネットワーク要素ＥＢと、第２ネットワーク要素ＥＣとが並列接続されている。ここで、第２ネットワーク要素ＥＢ、ＥＣのそれぞれは、粘弾性材料の粘弾性特性に大きな影響を与える粘性歪みを備える。この粘性歪みは、後述する粘性歪みの発展方程式を用いて、歪みと応力の関係を定めることができる。粘性歪みの発展方程式は、粘性歪みを定めるための複数の粘弾性パラメータを含み、この粘弾性パラメータに適切な値を与えることにより、実際の粘弾性材料の粘弾性特性を再現することができる。このように、第２ネットワーク要素は、大きさの異なる少なくとも２つの力学変形を再現した小変形ネットワーク要素（第２ネットワーク要素ＥＢ）及び大変形ネットワーク要素（第２ネットワーク要素ＥＣ）を含む。小変形ネットワーク要素と大変形ネットワーク要素の粘弾性パラメータの値は、後述するスイッチング関数を用いて力学変形に応じて設定されることが好ましい。

本実施形態の好ましい形態として、第２ネットワーク要素ＥＢ、ＥＣのそれぞれは、力学変形を粘弾性特性で表した、粘弾性パラメータの値が互いに異なるサブネットワーク要素を少なくとも２つ以上含む。図３（ａ），（ｂ）は、本実施形態の好ましい形態の一例を示す図である。第２ネットワーク要素ＥＢは、図３（ａ）に示すように、粘弾性パラメータの値が互いに異なるサブネットワーク要素ＥＢ１〜ＥＢｊ（ｊは２以上の自然数）で構成される。第２ネットワーク要素ＥＣは、図３（ｂ）に示すように、粘弾性パラメータの値が互いに異なるサブネットワーク要素ＥＣ（ｊ＋１）〜ＥＣｋ（ｋはｊより大きい２以上の自然数）で構成される。

このような粘弾性モデルにおける両端１，２の間に働く応力σ_１２は、両端１，２間の歪みをｅ_１２、サブネットワーク要素の粘性歪みをｅ_in,ｉ（ｉは、サブネットワーク要素を識別する番号で、ｉ＝１〜ｋの自然数）として、下記式（１）で表される。下記式（１）では、理解のために単純化してテンソルの代わりにスカラーで表している。式（１）の第１項のＡは第１ネットワーク要素ＥＡの非線形の弾性特性を表す係数であり、第２項のＢ_ｉは第２ネットワーク要素ＥＢのサブネットワーク要素ＥＢｉの粘弾特性を表す係数であり、第３項のＣ_ｉは第２ネットワーク要素ＥＣのサブネットワーク要素ＥＣｉの粘弾特性を表す係数である。係数Ａ、Ｂ_ｉ、Ｃ_ｉは、実験データに合わせて値が設定される材料定数である。

ここで、粘性歪みｅ_in,ｉは、下記式（２）で表される粘性歪みの発展方程式に従って算出されることが好ましい。式（２）では、歪みｅ_１２をｅと記している。

ここで、式（２）中のＺ_ｉ、ｄ_ｉ、λ_ｉは、ｉ番目のサブネットワーク要素の粘弾性パラメータであり、サブネットワーク要素毎に値が設定される。この値は、シミュレーション装置１０のパラメータ設定部１８において、予め粘弾性測定試験機２６で得られた実験データに一致するように、この実験データから抽出され、記憶部１４に記憶されている。なお、ｑ_ｉはｉ番目のサブネットワーク要素における時間ｔに依存した内部変数を表し、ｑ_ｉ（０）＝０である。式（２）の粘性歪みの発展方程式から、例えば、後退オイラー法を用いて粘性歪みｅ_ｉｎ，iを時間増分形式で求めることができる。ここで、歪みｅは、例えば、Piola歪み（＝（Ｃ^−１−Ｉ）／２：Ｃは右コーシー−グリーン変形テンソル、Ｉは単位テンソルである）あるいはグリーン−ラグランジェ歪み（＝（Ｃ−Ｉ）／２）が用いられる。粘弾性パラメータＺ_ｉは、サブネットワーク要素の緩和時間を示すパラメータである。また、Ｈ_ｉ（ｔ）＝１の場合、線形の粘弾性特性を表す。また、ｄｅ（ｔ）／ｄｔは、歪み速度である。式（２）では、歪み速度の内積の平方根（||Ａ||＝√（Ａ：Ａ）、Ａは歪み速度を表す）を用いている。また、ｄｅ（ｔ）／ｄｔの代わりに、変形速度テンソルＤの内積の平方根を用いることもできる。変形速度テンソルＤは、速度勾配テンソルＬをｄｖ／ｄｘで表すと、Ｌの対称成分、すなわち、Ｄ＝（Ｌ^Ｔ＋Ｌ）／２で表すことができる。このように、本実施形態では、第２ネットワーク要素ＥＢ，ＥＣの粘弾性特性は、粘性歪みの時間変動を表す粘性歪みの発展式を用いて表され、粘性歪みの発展式は、粘弾性材料の歪み速度ｄｅ（ｔ）／ｄｔあるいは変形速度（変形速度テンソルＤ）の関数である、ことが好ましい。さらに、粘性歪みの発展式は、粘弾性材料の全歪みと粘弾性材料の粘性歪みの関数であることも好ましい。粘弾性材料の全歪みは上記式（２）のｅに相当し、粘弾性材料の粘性歪みは上記式（２）のｅ_in,ｉに相当する。したがって、上記式（２）からわかるように、粘性歪みの発展式は、粘弾性材料の全歪みと粘弾性材料の粘性歪みの差分の関数でもある。粘弾性材料の全歪みと粘弾性材料の粘性歪みの差分は、弾性歪みを表し、弾性歪みと応力は、材料定数（ヤング率や横せん断係数）や歪みエネルギー関数を用いて一意に対応付けることができるため、粘性歪みの発展式は、応力の関数でもある。

さらに、好ましい形態として、粘性歪みの発展式は、式（２）の右辺の（ｅ−ｅ_ｉｎ，ｉ）に代えて粘弾性材料に与えられる歪みと粘性歪みの差分（ｅ−ｅ_ｉｎ，ｉ）に対する分数階微分の演算を含む式である。例えば、下記式（３）に示すように、粘性歪みの発展式を定めることができる。式（３）において、Ｄ^１−αｉ（Ｄは微分演算子）は、分数階微分の演算子である。ここで、α_ｉは、０より大きく１より小さい実数である。α_ｉは、分数階微分の分数階を定めるパラメータであり、α_ｉ＝１の場合、０階微分になり、式（２）と一致する。

このような、歪みｅと粘性歪みｅ_ｉｎ，ｉとの差分の分数階微分の演算結果は、歪みｅと粘性歪みｅ_ｉｎ，ｉとの差分の位相情報を含み、粘性歪みｅ_ｉｎ，ｉの挙動を含んだ粘弾性特性を、実際の実験データとよく一致させることができる。分数階微分については、“分数階微分モデル粘弾性の基礎理論”、清水信行他、日本機械学会D&D Conference 2000 ＣＤ−ＲＯＭ論文集に詳細に説明されている。

式（３）で表される式を増分形式で記載すると、下記式（４）で表すことができる。式（４）では、現在の時間ｔ_n＋１における粘性歪みをｅ_ｉｎ，ｉ ^ｔn＋１と表している。式（４）に示されるように、式（４）の右辺は全て既知であり、発散する部分がないので、左辺の現在（時間ｔ_n＋１）の粘性歪みｅ_ｉｎ，ｉは、安定した値として算出され得る。しかも、粘性歪みの各時間ステップの算出では、以前に算出した計算結果を用いることができるので、効率的に粘性歪みｅ_ｉｎ，ｉを算出することができる。すなわち、式（４）では、両端１，２間に歪みｅが時間によって変化するように与えられたとき、過去（時間ｔ_n以前の時間）の歪みｅと算出された過去（時間ｔ_n以前の時間）の粘性歪みｅ_ｉｎ，ｉから、時間ｔ_n＋１における粘性歪みｅ_ｉｎ，ｉを算出することができる。このように、本実施形態では、粘性歪みの発展式を、増分形式の粘性歪みの発展式として表して、粘性歪みを算出することが、短時間で効率よくシミュレーションを行う点から好ましい。

このような微小変形を再現する第２ネットワーク要素ＥＢの粘弾性パラメータの値と大変形を再現する第２ネットワーク要素ＥＣの粘弾性パラメータの値を、両端１，２の間の力学変形に応じて設定するために用いるスイッチング関数を、第２ネットワーク要素ＥＢ，ＥＣは含んでいる。このとき、第２ネットワーク要素ＥＢの粘弾性パラメータの値は、予め設定された第１の値、すなわち微小変形を代表する歪み（例えば０．２％歪み）におけるパラメータの値にスイッチング関数の値を乗算した値であり、第２ネットワーク要素ＥＣの粘弾性パラメータの値は、予め設定された第２の値、すなわち大変形を代表する歪み（例えば１５％歪み）におけるパラメータの値に、スイッチング関数の値を１から差し引いた値を乗算した値であることが好ましい。これによって微小変形と大変形の間の変形（中程度の変形）における粘弾性特性を再現することができる。
なお、第２ネットワーク要素ＥＢ，ＥＣがそれぞれサブネットワーク要素を含む場合、各サブネットワーク要素の粘弾性パラメータの値は、スイッチング関数を用いて両端１，２間の力学変形に応じて設定することが好ましい。

上記実施形態では、式（１）に示すように応力と歪みの関係を前提として粘性歪みを定めて第２ネットワーク要素ＥＢ，ＥＣにおける粘弾性特性を定めるが、マイクロ応力とアフィン変形によるマイクロ歪みの関係を利用して粘性歪みを定めて第２ネットワーク要素ＥＢ，ＥＣにおける粘弾性特性を定めることもできる。この場合の粘性歪みの発展方程式は、式（２）と同じ内容で表すことができる。このとき、上記式（２）で用いる歪みｅには、アフィン変形によるマイクロ歪みε_i（＝ｌｎ（λｆ_i）：λｆはアフィン変形におけるマイクロストレッチ（＝｜Ｆ・r_i｜：Ｆは変形勾配テンソル、r_iは離散化方向ベクトル））が用いられる。また、式（３）と同じように、分数階微分の演算を含む式で表すこともできる。アフィン変形によるマイクロ歪みについては、“C. Miehe and S. Goktepe, A micro-macro approach to rubber-like materials-Part II: The micro-sphere model of finite rubber viscoelasticity, Journal of the Mechanics and Physics of Solids, 53 (2005) 2231-2258”に詳細に説明されている。

図４は、スイッチング関数ｈの一例を示す図である。図４に示すスイッチング関数ｈは一例である。図４に示すスイッチング関数ｈは、一定の値を定める領域と、力学変形の量Ｉに従がって単調減少する領域とを有するが、単調増加する領域等があってもよく、極大点や極小点を有する領域があってもよい。図４に示すスイッチング関数ｈは、両端１，２間の力学変形の量Ｉに応じて値が変化する関数である。スイッチング関数ｈは、歪みが０．２％未満の力学変形から歪みが０．２％の力学変形までの微小変形時のスイッチング関数ｈの値を１．０とし、１５％歪みに相当する力学変形以上の大変形時のスイッチング関数ｈの値を０とし、その間の範囲を、ｌｏｇ（Ｉ−３）が大きくなるに連れてスイッチング関数ｈの値が単調減少するように設定されている。ここで、複数の粘弾性パラメータの少なくとも１つに関して、スイッチング関数ｈが１．０のとき、０．２％歪みに相当する力学変形における粘弾性パラメータの値を用い、スイッチング関数ｈが０．０のとき、１５％歪みに相当する力学変形における粘弾性パラメータの値を用い、その間の範囲では、ｌｏｇ（Ｉ−３）の大きさに応じて定まるスイッチング関数ｈの値を用いて、粘弾性パラメータの値を設定する。具体的には、０．２％の歪みに相当する力学変形における粘弾性パラメータの値と１５％の歪みに相当する力学変形における粘弾性パラメータの値との間で、スイッチング関数ｈの値を用いて算出した値を、粘弾性パラメータの値として設定する。図５は、粘弾性パラメータの値を算出する方法を説明する図である。図５に示す例では、粘弾性パラメータとしてＺ_ｉを用いている。図５に示すように、ｈをスイッチング関数の値として、ｈ及び（１−ｈ）を重み付け係数とし、具体的には、ｈ×（０．２％歪みに相当する力学変形における粘弾性パラメータＺ_ｉの値）を第２ネットワーク要素ＥＢの粘弾性パラメータとし、（１−ｈ）×（１５％歪みに相当する力学変形における粘弾性パラメータＺ_ｉの値）を第２ネットワーク要素ＥＣの粘弾性パラメータとし、それぞれのネットワーク要素ごとに応力を計算し、それらの応力を足し合わせる。また、周期的に繰返し負荷を与えて粘弾性材料を変形させる場合、スイッチング関数ｈの値は負荷時に経験した最大変形における値を固定値として用いてもよい。すなわち、除荷時にもスイッチング関数ｈの値は、最大変形時における値を用いてもよい。あるいは、周期的に繰返し負荷が与えられる場合、負荷方向が反転した（負荷時から除荷時に変化、あるいは除荷時から負荷時に変化した）際に、基準とする変形状態を変更し、例えば、反転した時点の変形状態を基準状態とし、力学変形の量の、現在の変形状態と基準状態との差に応じて値が変化するスイッチング関数ｈを作成してもよい。すなわち、負荷・除荷時におけるスイッチング関数ｈの値は、力学変形の量の、基準とする変形と現在の変形との間の最大差を用いて定まるものでもよい。ここで、基準とする変形と、現在の変形との間の差とは、例えば無変形状態の力学変形の量（例えば、歪みあるいは変形勾配テンソルの不変量）と、現在の力学変形の量との差、あるいは、周期的に繰返し負荷が与えられる場合における、負荷方向が反転した（負荷時から除荷時に変化、あるいは除荷時から負荷時に変化）した時点の力学変形の量（例えば、歪みあるいは変形勾配テンソルの不変量）と、現在の力学変形の量との差を表す。

図４中のグラフの横軸を定めるｌｏｇ（Ｉ−３）における力学変形の量Ｉは、右コーシー−グリーン変形テンソルの第１不変量であることが好ましい。右コーシー−グリーン変形テンソルは、変形勾配テンソルをＦで表したときＦ^Ｔ・Ｆで表されるため、右コーシー−グリーン変形テンソルの第１不変量は、力学変形に伴う主伸張比をλ_１、λ_２、λ_３としたとき、λ_１ ^２＋λ_２ ^２＋λ_３ ^２と表すことができる。すなわち、力学変形に伴う伸張の大きさに応じて、スイッチング関数ｈの値を変化させることが好ましい。
このようなスイッチング関数ｈによって値が設定される粘弾性パラメータは、Ｚ_ｉ、ｄ_ｉ、λ_ｉ等であるが、特に、粘性歪みの時間緩和を表す緩和時間パラメータＺ_ｉを含むことが好ましい。この場合、緩和時間パラメータＺ_ｉには、第２ネットワーク要素ＥＢの緩和時間パラメータＺ_ｉの値には、微小変形を代表する歪み（例えば０．５％歪み）における緩和時間パラメータＺ_ｉの値にスイッチング関数ｈの値を乗算した値を用い、第２ネットワーク要素ＥＣの緩和時間パラメータＺ_ｉの値には、大変形を代表する歪み（例えば１５％歪み）における緩和時間パラメータＺ_ｉの値に、（１−スイッチング関数ｈ）の値を乗算した値を用いることが好ましい。

このように、本実施形態のスイッチング関数ｈは、力学変形の量として、歪みの情報として扱われる力学変形における右コーシー−グリーン（Cauchy-Green）変形テンソルの不変量に基づいてスイッチング関数ｈの値を定めるが、歪みの情報である左Cauchy-Green変形テンソルの不変量、Green-Lagrange歪み、対数歪み、アフィン変形によるマイクロ歪みに基づいてスイッチング関数ｈの値を定めることも好ましい。また、粘弾性材料の変形勾配の情報に応じてスイッチング関数ｈの値を定めてもよい。この場合、スイッチング関数ｈは、歪みの情報（例えば、右コーシー−グリーン変形テンソルの不変量、左コーシー−グリーン変形テンソルの不変量、Green-Lagrange歪み、対数歪み、アフィン変形によるマイクロ歪みの値）あるいは変形勾配の情報とある基準値との差が大きくなる程、スイッチング関数ｈの値が小さくなる関数であることが好ましい。
上述したように、式（２）で用いる歪みｅにアフィン変形によるマイクロ歪みを用いる場合、右コーシー−グリーン変形テンソルの不変量の代わりに、アフィン変形によるマイクロ歪みが用いられる。この場合、スイッチング関数ｈは、２つの力学変形（微小変形及び大変形）におけるアフィン変形によるマイクロ歪みの間の範囲で、アフィン変形によるマイクロ歪みが大きくなり、ある基準値からの差が大きくなる程、スイッチング関数ｈの値が小さくなる関数であることが好ましい。

このように、本実施形態の粘弾性モデル４０の第２ネットワーク要素ＥＢ，ＥＣは、粘弾性パラメータの値を、力学変形に応じて設定するスイッチング関数ｈを含むので、ゴム等の粘弾性材料の力学変形をコンピュータで計算する際に、粘弾性材料の様々な大きさの力学変形であっても、力学変形の大きさに応じて変化する粘弾性特性を統一的に再現することができる。

このようなネットワーク要素を含んだモデルを、以下説明する粘弾性材料のシミュレーションに好適に適用することができる。

（シミュレーション装置の説明）
図１に示すパラメータ設定部１８は、粘弾性測定試験機２６から得られる各種の粘弾性材料の実験データ（大変形の実験データと微小変形の実験データ）を用いて、上述した図２，３に示すネットワーク要素を含む粘弾性モデル４０を用いて、粘弾性モデル４０で規定されている粘弾性パラメータである係数Ａ、Ｂ_ｉ、Ｃ_ｉ、Ｚ_ｉ、ｄ_ｉ、λ_ｉ、α_ｉの値を抽出する。この値の抽出に関しては、具体的には、粘弾性パラメータである係数Ａ、Ｂ_ｉ、Ｃ_ｉ、Ｚ_ｉ、ｄ_ｉ、λ_ｉ、α_ｉの値を種々変更しながら、応力を算出し、算出結果が粘弾性試験の応力の実験データに一致するように、係数Ａ、Ｂ_ｉ、Ｃ_ｉ、Ｚ_ｉ、ｄ_ｉ、λ_ｉ、α_ｉの値を探索することにより、値を抽出する。このようにして抽出された値は、記憶部１４に記憶される。

粘弾性材料の力学変形のシミュレーションを行う場合、モデル作成部２０は、ディスプレイに表示されたモデル作成画面に従がってオペレータがマウスやキーボード等の入力操作系２８から情報を入力することにより、粘弾性材料のシミュレーションモデルを作成する。モデル作成部２０は、例えば、有限要素法に基づいた構造体のシミュレーションモデルを作成する。構造体のシミュレーションモデルには、粘弾性材料のシミュレーションモデルが含まれている。すなわち、本実施形態では、力学変形のシミュレーションを行う粘弾性材料のシミュレーションモデルを含み、金属部材や有機繊維材等の弾性材料のシミュレーションモデルを含んだ構造体のシミュレーションモデルを作成する。粘弾性材料としては、ゴム部材が挙げられる。また、構造体として、空気入りタイヤが挙げられる。勿論、構造体は、粘弾性材料のみで構成された構造体であってもよい。シミュレーションモデルの作成では、具体的には、複数の節点と、節点で定義される要素を定め、各要素に材料定数を定める。ここで、シミュレーションモデルには、第１ネットワーク要素ＥＡ及び第２ネットワーク要素ＥＢ，ＥＣを含む粘弾性モデル４０が設けられる。粘弾性モデル４０の作成では、粘弾性材料の粘弾性特性が再現されるように、パラメータ設定部１８は、サブネットワーク要素の数と、各サブネットワーク要素に対応した係数Ｂ_ｉ、Ｃ_ｉ、Ｚ_ｉ、ｄ_ｉ、λ_ｉ、α_ｉの値を記憶部１４から読み出して設定する。

演算部２２は、静解析及び動解析のシミュレーションを行うことができるような公知のソルバー、例えばＡＢＡＱＵＳが用いられる。演算部２２は、まず、シミュレーションモデルを用いて静解析を行う。静解析では、時間変化をしない変形を構造体のシミュレーションモデルに与えて時間変動のない歪み及び応力を算出する解析である。したがって、静解析では、演算部２２は、粘性歪み-ｅ_ｉｎ，ｉを＝０として演算する。この静解析は、動解析を開始するときの粘弾性パラメータの値を設定するために行なわれる。例えば、静解析では、空気入りタイヤの場合、空気入りタイヤとリムで囲まれたタイヤ空洞領域の雰囲気圧力が所定の空気圧になるようにタイヤ空洞領域に空気を充填すること、及び空気を充填した空気入りタイヤを設定した荷重で路面に接地させること、を構造体のシミュレーションモデルで再現する。

演算部２２は、静解析のシミュレーションで得られたシミュレーションモデルの各要素の力学変形から、例えば右コーシー−グリーン変形テンソルの第１不変量を求め、この不変量を力学変形の量Ｉとして、各要素のスイッチング関数ｈの値を求める。求めた値に基づいて、各要素の粘弾性パラメータの値を定めて、動解析のシミュレーションを開始する。動解析のシミュレーションでは、演算部２２は、時間ステップ毎に、変形を与えたシミュレーションモデルの歪み及び応力を算出する。動解析のシミュレーション中、力学変形のレベルは変化するので、力学変形の量Ｉも変化する。このため、動解析における時間ステップ毎に、力学変形の量Ｉを求め、求めた結果に応じてスイッチング関数ｈの値を求めることができる。これにより、動解析中、粘弾性パラメータの値も変化させることができる。動解析のシミュレーションで時間ステップ毎に歪み及び応力を算出するとき、式（２）に従がって、あるいは式（３）、（４）に従がって粘性歪みｅ_ｉｎ，ｉが算出され、応力が式（１）を用いて算出される。こうして、演算部２２は、シミュレーションモデルに与えた変形に対する歪みと応力を、粘弾性材料の力学変形を表す物理量として算出することにより、粘弾性材料の力学変形を再現する。例えば、動解析では、空気入りタイヤの場合、接地した空気入りタイヤを路面上で転動させることを再現するシミュレーションを行なって転がり抵抗（タイヤ回転軸の前後方向に作用する力）を算出する。このシミュレーションでは、結果にゴム部材の大変形時の特性が大きく影響する。また、空気入りタイヤの振動や騒音のシミュレーションを行う場合、振動や騒音の発生を再現するシミュレーションを行って、固有振動数や振動や騒音の周波数スペクトル波形、振動レベルあるいは音圧レベルを算出する。このシミュレーションでは、結果にゴム部材の微小変形時の特性が大きく影響する。

制御管理部２４は、パラメータ設定部１８、モデル作成部２０、演算部２２の動作を管理し制御する部分である。また、制御管理部２４は、演算部２２で演算された粘弾性材料の再現結果、すなわち歪みと応力の結果をディスプレイ・プリンタ３０に出力するように制御する。
このように本実施形態では、粘弾性材料としてゴム部材を含む構造体のシミュレーション方法を行うことができる。この場合、上述した粘弾性材料のシミュレーション方法を、構造体のゴム部材に適用する。弾性部材については、粘性歪み-ｅ_ｉｎ，ｉ＝０であるので、動解析では、時間に対して一定の材料定数が用いられる。

（粘弾性材料のシミュレーション方法）
図６は、本実施形態の粘弾性材料のシミュレーション方法のフローの一例を説明する図である。図６に示すフローの例は、有限要素法を用いて解析を行なう有限要素モデルを作成する例である。
モデル作成部２０は、オペレータの入力指示に応じて、粘弾性材料を含んだ構造体の有限要素モデルを作成する（ステップＳ１０）。このとき、制御管理部２４は、オペレータの入力指示によって、粘弾性試験を再現するために、静解析及び動解析のシミュレーションの内容を設定する。また、パラメータ設定部１８は、粘弾性材料の種類に応じてシミュレーションモデルに用いる粘弾性モデル４０の各粘弾性パラメータの値を設定する。

この後、演算部２２は、静解析のシミュレーションを行う（ステップＳ１２）。この静解析のシミュレーションにより、所定の変形が与えられた静的な（時間によって変化しない）力学変形を再現する。すなわち、静解析のシミュレーションでは、動解析のシミュレーションを行うときの粘弾性材料の係数Ｂ_ｉ、Ｃ_ｉ、Ｚ_ｉ、ｄ_ｉ、λ_ｉである粘弾性パラメータの値をスイッチング関数ｈの値に基づいて設定するために、粘弾性試験の開始時の粘弾性材料の力学変形を再現する。パラメータ設定部１８は、粘弾性パラメータである係数Ｂ_ｉ、Ｃ_ｉ、Ｚ_ｉ、ｄ_ｉ、λ_ｉの値を、力学変形に応じて値が設定されるスイッチング関数ｈを用いて設定する（ステップＳ１４）。

この後、演算部２２は、動解析のシミュレーションを行う（ステップＳ１６）。例えば、静解析で与えられた歪みを初期状態として、粘弾性材料のシミュレーションモデルを時間的に変化する変形を与えて（振動させて）動解析のシミュレーションを行う。このとき、動解析のシミュレーションの各時間ステップにおいて、スイッチング関数ｈの値は、力学変形のレベルに応じて変化するが、このスイッチング関数ｈの値が変化する度に、係数Ｂ_ｉ、Ｃ_ｉ、Ｚ_ｉ、ｄ_ｉ、λ_ｉの値を設定し直すことができる。こうして、演算部２２は、シミュレーションモデルを用いて動解析のシミュレーションを行なうことにより、粘弾性材料を含む構造体の力学変形を再現することができる。
なお、本実施形態のフローでは、粘弾性試験の開始時の力学変形を再現するために予め静解析を行なうが、必ずしも静解析のシミュレーションは行わなくてもよい。例えば、試験開始時の粘弾性材料が力学変形を受けていないような粘弾性試験を再現する場合、シミュレーションモデルに予め変形を与えない状態から、時間とともに変形が大きくなるような動解析のシミュレーションを行えばよく、静解析のシミュレーションを行わなくてもよい。この場合、動解析のシミュレーションの各時間ステップにおいて、スイッチング関数ｈの値を用いて、スイッチング関数ｈの値が変化する度に、係数Ｂ_ｉ、Ｃ_ｉ、Ｚ_ｉ、ｄ_ｉ、λ_ｉの値を設定し直せばよい。
また、試験開始時の粘弾性材料が力学変形を予め受けているような粘弾性試験を再現する場合でも、静解析を行わず、動解析のシミュレーションの各時間ステップにおいて、スイッチング関数ｈの値を用いて、スイッチング関数ｈの値が変化する度に、係数Ｂ_ｉ、Ｃ_ｉ、Ｚ_ｉ、ｄ_ｉ、λ_ｉの値を設定し直してもよい。

この後、制御管理部２４は、算出されたシミュレーションの結果、すなわち応力の算出結果を歪みの算出結果とともに、粘弾性材料の力学変形を表す物理量としてディスプレイ・プリンタ３０に出力させる（ステップＳ１８）。
このようなステップＳ１０〜１８は、シミュレーション装置１０を構成するコンピュータにより実質的に行なわれる。

（実施例、比較例）
本実施形態は、第２ネットワーク要素ＥＢ，ＥＣと、スイッチング関数ｈを用いることにより、大変形及び微小変形の幅広い範囲の粘弾性特性の実験データを統一的に再現することができる粘弾性材料のシミュレーションモデルを開示する。この点を、実施例、比較例を用いて説明する。
図７〜１１は、本実施形態のスイッチング関数ｈを含んだ粘弾性モデル（実施例）の効果の一例を示す図である。粘弾性モデル（実施例）は、第１ネットワーク要素ＥＡと第２ネットワーク要素ＥＢ，ＥＣとを並列接続したモデルであって、第２ネットワーク要素ＥＢ，ＥＣは、スイッチング関数ｈを用いて粘弾性パラメータの値を、力学変形に応じて設定することが可能に構成されている。本実施例では、第２ネットワーク要素ＥＢのサブネットワーク要素の総数を２とし、第２ネットワーク要素ＥＣのサブネットワーク要素の総数を１とした。また、微小変形を再現する第２ネットワーク要素ＥＢのサブネットワーク要素には、式（３）に示す分数階微分の演算子を含む粘性歪み発展方程式を用い、大変形を再現する第２ネットワーク要素ＥＣのサブネットワーク要素には、式（２）に示す粘性歪み発展方程式を用いた。この粘弾性モデルを用いた有限要素モデルであるシミュレーションモデルを作成した。シミュレーションモデルは、粘弾性試験に用いる粘弾性材料のサンプルを再現した有限要素モデルである。
具体的には、図７〜１１は、微小変形の粘弾性試験の実験データに一致するように、第１ネットワーク要素ＥＡと第２ネットワーク要素ＥＢを含む粘弾性モデルの粘弾性パラメータの値を設定し、さらに、大変形の粘弾性試験の実験データに一致するように、第１ネットワーク要素ＥＡと第２ネットワーク要素ＥＣを含む粘弾性モデルの粘弾性パラメータの値を設定した後、図４に示すスイッチング関数ｈを設定し、シミュレーションモデルに微小変形および大変形を与えたとき、微小変形および大変形の粘弾性試験の実験データに略一致した応力を算出することができることを示している。粘弾性材料として、カーボンブラックを２５％（体積分率）含むスチレン・ブタジエンゴムのサンプルを用いた。

より具体的には、図７では、温度２５℃、歪み振幅０．１％（平均歪み１０％）の条件下の、微小変形の粘弾性試験における貯蔵弾性率Ｅ’と損失弾性率Ｅ’’の角周波数”ω’ｒａｄ／ｓｅｃ”依存性を示している。図７に示すように、貯蔵弾性率Ｅ’と損失弾性率Ｅ’’の角周波数”ω’ｒａｄ／ｓｅｃ”依存性に関して、シミュレーションモデルの算出結果“●Ｅ’（実施例）”及び“■Ｅ’ ’（実施例）”は、実験データ“−−−Ｅ’(実験)”及び“―Ｅ’’(実験)”に略一致している。図８では、周波数０．３１２５Ｈｚにおける微小変形の粘弾性試験における貯蔵弾性率Ｅ’の歪み振幅”strain amplitude”依存性に関して、シミュレーションモデルの算出結果“○DMA（実施例）”は、実験データ“●DMA(実験)”に略一致している。

図９（ａ）、（ｂ）は、大変形の粘弾性試験の内容とその結果を示す。図９（ａ）に示すように、歪み速度１．３×１０^−２（／ｓｅｃ）の条件下、歪み振幅４０％（１０サイクル）、歪み振幅２０％（１０サイクル）の歪みを、粘弾性材料のＪＩＳ−１号ダンベル試験片（JIS K 6251）に与えた。このときの応力”stress”に関して、図８（ｂ）に示すように、本実施形態の算出結果“−−−Uni_multi_cycle（実施例）”は、実験データ“―Uni_multi_cycle(実験）”に略一致している。

図１０（ａ）は、大変形の粘弾性試験として粘弾性材料に与えた公称歪み“Nominal strain”の時間変化を示す図である。図１０（ｂ）は、このときの応力“stress”の変化を示す図である。応力“stress”に関して、図１０（ｂ）に示すように、シミュレーションモデルの算出結果“−−−Relax（実施例）”は、実験データ“―Relax(実験）”に略一致している。

図１１（ａ）は、粘弾性材料に０から５０％の伸張歪みを繰り返し与える、大変形の粘弾性試験である繰り返し伸張試験（歪み速度２．８×１０^−２（／ｓ））の結果を示す図である。図１１（ａ）に示すように、応力“stress”に関して、シミュレーションモデルの算出結果“−−−Uni_1_cycle_High（実施例）”は、実験データ“―Uni_1_cycle_High(実験）”に略一致している。

図１１（ｂ）は、粘弾性材料に０から５０％の伸張歪みを繰り返し与える、大変形の粘弾性試験である繰返し伸張試験（歪み速度５．７×１０^−５（／ｓ））の結果を示す図である。図１１（ｂ）に示すように、応力“stress”に関して、シミュレーションモデルの算出結果“−−−Uni_1_cycle_low（実施例）”は、実験データ“―Uni_1_cycle_low(実験）”に略一致している。

このように、微小変形と大変形のそれぞれで粘弾性試験の実験データに一致するように、第１ネットワーク要素ＥＡと第２ネットワーク要素ＥＢ,ＥＣを含む粘弾性モデル４０を用いた粘弾性材料のシミュレーションモデルにおいて、第２ネットワーク要素ＥＢ,ＥＣの粘弾性パラメータの値を設定し、このモデルに微小から大変形を与えたとき、このモデルは、微小から大変形の粘弾性試験の実験データに略一致した応力を算出することができる。

図１２〜図１４は、比較例１の粘弾性モデルを用いた結果の一例を示す図である。比較例１の粘弾性モデルは、本実施形態と同様に、第１ネットワーク要素ＥＡ、第２ネットワーク要素ＥＢ，ＥＣを含むが、スイッチング関数ｈを含まない粘弾性モデルである。そのため、第２ネットワーク要素ＥＢと第２ネットワーク要素ＥＣは区別されておらず同一である。比較例１においては、第２ネットワーク要素の総数を１７とした。この粘弾性モデルを含んだ有限要素モデルであるシミュレーションモデルを作成した。シミュレーションモデルは、粘弾性試験に用いる粘弾性材料のサンプルを再現した有限要素モデルである。

図１２〜１４では、微小変形の粘弾性特性の実験データに一致するように、粘弾性パラメータの値を設定した粘弾性モデル（比較例１）を用いて大変形の粘弾性特性の応力の結果を算出した。

図１２は、図７と同様の微小変形の粘弾性試験を行なったときの貯蔵弾性率Ｅ’と損失弾性率Ｅ’’の角周波数”ω’ｒａｄ／ｓｅｃ”依存性を示す図である。図１２に示すように、シミュレーションモデルの算出結果“●Ｅ’(比較例１)”及び“■Ｅ’ ’(比較例１)”は、実験データ“−−−Ｅ’(実験)”及び“―Ｅ’’(実験)”に略一致している。図１３は、図８に示す例と同様に、周波数０．３１２５Ｈｚにおける貯蔵弾性率Ｅ’の歪み振幅”strain amplitude”依存性の結果を示す図である。この場合、シミュレーションモデルの算出結果“○DMA(比較例１)”は、実験データ“●DMA(実験)”と一致していない。

さらに、図１４（ａ）〜（ｄ）は、図９（ｂ）、図１０（ｂ）、及び図１１（ａ），（ｂ）の結果に対応する大変形の粘弾性試験における応力の結果を示す。いずれの場合も、比較例１の応力の結果は、大変形の粘弾性試験の実験データの結果と一致しない。このように、スイッチング関数ｈを用いない比較例１では、微小変形の粘弾性試験の実験データに一致するように粘弾性パラメータの値を定めると、大変形の粘弾性試験の実験データに一致しないことがわかる。

図１５〜１７は、比較例２の粘弾性モデルを用いた結果の一例を示す図である。比較例２の粘弾性モデルは、比較例１と同様に、第１ネットワーク要素ＥＡ、第２ネットワーク要素ＥＢ，ＥＣを含むが、スイッチング関数ｈを含まないモデルである。そのため、第２ネットワーク要素ＥＢと第２ネットワーク要素ＥＣは区別されておらず同一である。比較例２においては、第２ネットワーク要素の総数を２とした。この粘弾性モデルを含んだ有限要素モデルであるシミュレーションモデルを作成した。シミュレーションモデルは、粘弾性試験に用いる粘弾性材料のサンプルを再現した有限要素モデルである。

図１５〜１７は、比較例２の粘弾性モデルを用いたときの結果を、実験データとともに示す。この場合、大変形の粘弾性試験の実験データに一致するように、粘弾性パラメータの値を設定した粘弾性モデルを用いて微小変形の粘弾性試験における応力の結果を算出した。

図１７（ａ）〜（ｄ）は、図９（ｂ）、図１０（ｂ）、及び図１１（ａ），（ｂ）の結果に対応する大変形の粘弾性試験の応力の結果を示す。いずれの場合も、比較例２の結果は、大変形の粘弾性試験の実験データの結果と略一致している。図１５、図１６は、比較例２における、図７、図８に対応する微小変形の粘弾性試験における貯蔵弾性率Ｅ’あるいは損失弾性率Ｅ’’の結果を示す。いずれの場合も、比較例２の貯蔵弾性率Ｅ’あるいは損失弾性率Ｅ’’の結果は、実施例に比べて粘弾性試験の実験データと一致しない。このように、スイッチング関数ｈを用いない比較例２のモデルでは、大変形の粘弾性試験の実験データに一致するように、粘弾性パラメータの値を定めると、微小変形の粘弾性試験の実験データに一致しないことがわかる。

したがって、比較例１及び比較例２では、大変形の粘弾性試験の実験データ及び微小変形の粘弾性試験の実験データのいずれも満足するような粘弾性パラメータの値を得ることができない。これに対して、スイッチング関数ｈを用いることで、変形の粘弾性試験の実験データ及び微小変形の粘弾性試験の実験データのいずれも満足するような粘弾性パラメータの値を得ることができる。すなわち、本実施形態のスイッチング関数ｈを含む粘弾性モデル４０の効果は明らかである。
したがって、本実施形態の粘弾性材料のシミュレーション方法を用いれば、粘弾性材料を含む構造体において粘弾性材料が大変形の力学変形をするシミュレーションと、粘弾性材料が微小変形の力学変形をするシミュレーションを、１つの構造体のシミュレーションモデルを用いて統一して行なうことができる。例えば、空気入りタイヤの場合、ゴム部材の大変形時の特性が大きく影響する転がり抵抗の評価と、ゴム部材の微小変形時の特性が大きく影響する振動レベルや音圧レベルの評価とを、１つの粘弾性材料のシミュレーションモデルを含む空気入りタイヤのシミュレーションモデルを用いて統一して行うことができる。

本実施形態は、粘弾性材料のシミュレーション方法及びシミュレーション装置を含むが、さらに、コンピュータに実行させるプログラムも含む。このプログラムは、記録媒体に記録されたものでもよいし、インターネット等を介して提供されるものでもよい。
すなわち、粘弾性材料の力学変形をコンピュータに計算させるプログラムは、
力学変形を非線形の弾性特性で表した第１ネットワーク要素ＥＡと、力学変形の粘弾性特性を、粘弾性パラメータを用いて表したネットワーク要素であって、力学変形の大きさに応じて粘弾性特性が変化する第２ネットワーク要素ＥＢ，ＥＣと、を並列接続した粘弾性モデルをコンピュータに設定させる手順と、
前記粘弾性モデルを用いてコンピュータに、前記粘弾性材料の前記力学変形を再現するための演算を実行させる手順と、を含む。
このとき、第２ネットワーク要素ＥＢ、ＥＣは、粘弾性パラメータの値を、力学変形に応じて設定するスイッチング関数ｈを含む。
また、本実施形態は、上記手順を含むコンピュータが実行可能なプログラムを記録した、非一時的（non-transitory）なコンピュータの読み取り可能な記録媒体を提供することができる。

以上、本発明の粘弾性材料のシミュレーション方法、構造体のシミュレーション方法、粘弾性材料のシミュレーション装置、及び、プログラムについて詳細に説明したが、本発明は上記実施形態に限定されず、本発明の主旨を逸脱しない範囲において、種々の改良や変更をしてもよいのはもちろんである。

１０ シミュレーション装置
１２ ＣＰＵ
１４ 記憶部
１６ 入出力部
１８ パラメータ設定部
２０ モデル作成部
２２ 演算部
２４ 制御管理部
２６ 粘弾性測定試験機
２８ 入力操作系
３０ ディスプレイ・プリンタ
４０ 粘弾性モデル

Claims (15)

1. 粘弾性材料の力学変形をコンピュータで再現するシミュレーション方法であって、
   粘弾性材料の力学変形を非線形の弾性特性で表した第１ネットワーク要素と、前記力学変形の粘弾性特性を、粘弾性パラメータを用いて表したネットワーク要素であって、前記力学変形の大きさに応じて前記粘弾性特性が変化する第２ネットワーク要素と、を並列接続した粘弾性モデルをコンピュータに作成させるステップと、
   前記粘弾性モデルを用いてコンピュータに、前記粘弾性材料の前記力学変形を再現するための演算を実行させるステップと、を含み、
   前記第２ネットワーク要素は、前記粘弾性パラメータの値を前記力学変形に応じて設定するために用いるスイッチング関数を含む、ことを特徴とする粘弾性材料のシミュレーション方法。
2. 前記第２ネットワーク要素は、大きさの異なる少なくとも２つの力学変形を再現した小変形ネットワーク要素及び大変形ネットワーク要素を含み、
   前記小変形ネットワーク要素と前記大変形ネットワーク要素の前記粘弾性パラメータの値は、前記スイッチング関数を用いて前記粘弾性材料の前記力学変形に応じて設定される、請求項１に記載の粘弾性材料のシミュレーション方法。
3. 前記小変形ネットワーク要素の前記粘弾性パラメータの値は、予め設定された第１の値に前記スイッチング関数の値を乗算した値であり、
   前記大変形ネットワーク要素の前記粘弾性パラメータの値は、予め設定された第２の値に、前記スイッチング関数の値を１から差し引いた値を乗算した値である、請求項２に記載の粘弾性材料のシミュレーション方法。
4. 前記小変形ネットワーク要素及び前記大変形ネットワーク要素のそれぞれは、前記力学変形を粘弾性特性で表した、前記粘弾性パラメータの値が互いに異なるサブネットワーク要素を少なくとも２つ以上含む、請求項３に記載の粘弾性材料のシミュレーション方法。
5. 前記スイッチング関数は、前記力学変形における、前記粘弾性材料の歪みの情報あるいは前記粘弾性材料の変形勾配の情報に応じて前記スイッチング関数の値を定める、請求項１〜４のいずれか１項に記載の粘弾性材料のシミュレーション方法。
6. 前記第２ネットワーク要素の粘弾性特性は、前記粘弾性材料に作用する粘性歪みの時間緩和を表す粘性歪みの発展式を用いて表され、前記粘性歪みの発展式は、前記粘弾性材料の歪み速度の関数、前記粘弾性材料の変形速度の関数、あるいは前記粘弾性材料の全歪みと前記粘弾性材料の粘性歪みの関数である、請求項１〜５のいずれか１項に記載の粘弾性材料のシミュレーション方法。
7. 前記粘性歪みの発展式は、前記粘弾性材料に与えられる歪みと前記粘性歪みの差分に対する分数階微分の演算を含む、請求項６に記載の粘弾性材料のシミュレーション方法。
8. 前記粘性歪みの発展式は、分数階微分を用いて下記式で表される、請求項６または７に記載の粘弾性材料のシミュレーション方法。
   （上記式において、ｅは前記粘弾性材料の力学変形における歪みを表し、ｅ_ｉｎ,ｉは前記第２ネットワーク要素のうちｉ番目のサブネットワーク要素における前記粘性歪みを表し、Ｚ_ｉ、ｄ_ｉ、λ_ｉ、α_ｉは、ｉ番目の前記サブネットワーク要素における粘弾性パラメータであり、ｑ_ｉはｉ番目の前記サブネットワーク要素における時間ｔに依存した内部変数を表し、α_ｉは０より大きく１より小さい実数であり、Ｄは微分演算子を表す）
9. 前記粘性歪みの発展式は、増分形式の粘性歪みの発展式である、請求項６〜８のいずれか１項に記載の粘弾性材料のシミュレーション方法。
10. 前記スイッチング関数を用いて値を設定する粘弾性パラメータは、前記粘性歪みの時間緩和を表す緩和時間パラメータを含む、請求項６〜９のいずれか１項に記載の粘弾性材料のシミュレーション方法。
11. 前記スイッチング関数は、前記力学変形における、前記粘弾性材料の歪みの情報あるいは前記粘弾性材料の変形勾配の情報に応じて前記スイッチング関数の値を定める関数であって、前記歪みの情報あるいは前記変形勾配の情報と、ある基準値との差が大きくなる程、前記スイッチング関数の値が小さくなる関数である、請求項６〜１０のいずれか１項に記載の粘弾性材料のシミュレーション方法。
12. ゴム部材を含む構造体のシミュレーション方法であって、
    前記構造体のシミュレーションを行うとき、請求項１〜１１のいずれか１項に記載の粘弾性材料のシミュレーション方法を前記ゴム部材に適用する、ことを特徴とする構造体のシミュレーション方法。
13. 前記構造体は、空気入りタイヤである、請求項１２に記載の構造体のシミュレーション方法。
14. 粘弾性材料の力学変形を再現するシミュレーション装置であって、
    前記力学変形を非線形の弾性特性で表した第１ネットワーク要素と、前記力学変形の粘弾性特性を、粘弾性パラメータを用いて表したネットワーク要素であって、前記力学変形の大きさに応じて前記粘弾性特性が変化する第２ネットワーク要素と、を並列接続したモデルであって、前記第２ネットワーク要素は、前記粘弾性パラメータの値を、前記力学変形に応じて設定するスイッチング関数を含む、粘弾性モデルを作成するモデル作成部と、
    前記粘弾性モデルを用いて、前記粘弾性材料の前記力学変形を再現するための演算を行う演算部と、を含む、ことを特徴とする粘弾性材料のシミュレーション装置。
15. 粘弾性材料の力学変形をコンピュータに再現させるプログラムであって、
    前記力学変形を非線形の弾性特性で表した第１ネットワーク要素と、前記力学変形の粘弾性特性を、粘弾性パラメータを用いて表したネットワーク要素であって、前記力学変形の大きさに応じて前記粘弾性特性が変化する第２ネットワーク要素と、を並列接続した粘弾性モデルをコンピュータに設定させる手順と、
    前記粘弾性モデルを用いてコンピュータに、前記粘弾性材料の前記力学変形を再現するための演算を実行させる手順と、を含み、
    前記第２ネットワーク要素は、前記粘弾性パラメータの値を、前記力学変形に応じて設定するスイッチング関数を含む、ことを特徴とするプログラム。