(57)【特許請求の範囲】
【請求項１】
コンピュータが実行する方法であって、
予め設定された解析条件下において予め設定された分子モデルデータを用いた分子動力学計算を行い、各解析時点における応力データを算出するステップと、
前記応力データに基づき各解析時点毎の緩和弾性率を算出するステップと、
前記緩和弾性率を補間関数で近似した近似値を算出するステップと、
前記緩和弾性率の近似値を微分処理して緩和スペクトルを算出するステップと、
を含む緩和スペクトルを算出する方法。
【請求項２】
前記補間関数は、Ｂスプライン関数である、請求項１に記載の方法。
【請求項３】
前記緩和弾性率の近似値から緩和スペクトルを算出する演算は、式（２）に表されるＳｃｈｗａｒｚｌ−Ｓｔａｖｅｒｍａｎによる近似式を用いており、次数ｋは２次以上である、請求項１又は２に記載の方法。
【数１】
0006579904.tif 000012
Ｇは緩和弾性率、τは緩和時間である。
【請求項４】
請求項１〜３のいずれかに記載の方法であって、
前記緩和スペクトルに基づき、貯蔵弾性率、損失弾性率及び損失正接のうち少なくとも１つの物性値を算出するステップを更に含む、物性値を算出する方法。
【請求項５】
予め設定された解析条件下において予め設定された分子モデルデータを用いた分子動力学計算の実行結果に基づき、各解析時点における応力データを算出する応力算出部と、
前記応力データに基づき各解析時点毎の緩和弾性率を算出する緩和弾性率算出部と、
前記緩和弾性率を補間関数で近似した近似値を算出する近似部と、
前記緩和弾性率の近似値を微分処理して緩和スペクトルを算出する緩和スペクトル算出部と、
を備える緩和スペクトルを算出する装置。
【請求項６】
前記補間関数は、Ｂスプラインである、請求項５に記載の装置。
【請求項７】
前記緩和弾性率の近似値から緩和スペクトルを算出する演算は、式（２）に表されるＳｃｈｗａｒｚｌ−Ｓｔａｖｅｒｍａｎによる近似式を用いており、次数ｋは２次以上である、請求項５又は６に記載の装置。
【数２】
0006579904.tif 000013
Ｇは緩和弾性率、τは緩和時間である。
【請求項８】
請求項５〜７に記載の装置であって、
前記緩和スペクトルに基づき、貯蔵弾性率、損失弾性率及び損失正接のうち少なくとも１つの物性値を算出する物性値算出部を更に備える、物性値を算出する装置。
【請求項９】
請求項１〜４のいずれかに記載の方法をコンピュータに実行させるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本開示は、分子動力学計算により緩和スペクトルを算出する方法、物性値（貯蔵弾性率、損失弾性率及び損失正接など）を算出する方法、装置及びプログラムに関する。
【背景技術】
【０００２】
貯蔵弾性率、損失弾性率及び損失正接等の物性値の計算は、次のように行われる。
（１）分子動力学計算を実行し、応力データを得る。
（２）グリーン・久保公式に基づき応力データから緩和弾性率を算出する。
（３）Ｓｃｈｗａｒｚｌ−Ｓｔａｖｅｒｍａｎによる近似式を用い、緩和弾性率に基づき緩和スペクトルを算出する。
（４）緩和スペクトルをＳｔｉｅｌｔｊｅｓ変換し、貯蔵弾性率、損失弾性率及び損失正接を算出する。
【０００３】
例えば、特許文献１には、高分子モデルを用いて分子動力学計算を行い、その結果からグリーン・久保公式に基づき緩和弾性率を算出することが開示されている。緩和弾性率Ｇ（ｔ）は、式（１）を用いて算出できる。
【数１】
0006579904.tif 000002
σ（ｔ）は時刻ｔのせん断応力テンソル、Ｖは体積、Ｔは温度、ＫＢはボルツマン定数である
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１４−２０３２６２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
緩和スペクトルＨ（τ）を算出するためのＳｃｈｗａｒｚｌ−Ｓｔａｖｅｒｍａｎによる近似式は、式（２）に表される。式（２）の次数ｋを無限大にする計算は現実的に不可能であるので、式（３）に示す１次の式が利用しやすい。すなわち、緩和スペクトルの算出には、式（３）に示すように微分演算が含まれる。τは緩和時間である。
【数２】
0006579904.tif 000003
【数３】
0006579904.tif 000004
【０００６】
緩和スペクトルＨ（τ）は、貯蔵弾性率、損失弾性率及び損失正接等の物性値を算出する基礎となるデータであるので、正しく計算されていることが望まれる。しかし、緩和スペクトルＨ（τ）がどの程度正しいかを見積もることは難しい。そこで、式（４）を用いて緩和スペクトルＨ（τ）をラプラス変換すれば、緩和弾性率Ｇ（ｔ）に戻すことができることを利用する。式（４）を用いて復元した緩和弾性率Ｇ（ｔ）復元と元の緩和弾性率Ｇ（ｔ）が同じになれば、緩和スペクトルＨ（τ）が正しく算出されていることになる。図５に示すように、Ｇ（ｔ）とＧ（ｔ）復元との残差の最小二乗和を算出して評価することが考えられる。
【数４】
0006579904.tif 000005
【０００７】
しかしながら、図５に示すように、短時間領域と長時間領域において両者が大きくことなっていることが分かる。すなわち、緩和スペクトルの算出精度が十分であるとはいえない。
【０００８】
緩和スペクトルＨ（τ）を精度良く求めるためのアプローチの一つとして、Ｓｃｈｗａｒｚｌ−Ｓｔａｖｅｒｍａｎによる近似式（２）の次数を上げることが考えられる。しかしながら、図６Ａに示す１次から図６Ｂに示す２次になれば改良されるが、図６Ｃに示す３次になれば、逆に悪化してしまう場合があることが判明した。図６Ａに示す例では、式（２）の１次による算出結果の残差二乗和が0.0455である。図６Ｂに示す例では、式（２）の２次による算出結果の残差二乗和が0.0247である。図６Ｃに示す例では、式（２）の３次による算出結果の残差二乗和が0.829であった。３次で悪化している。
【０００９】
上記した緩和スペクトルの算出精度が十分でないという認識、精度を向上させる課題、その解決手段に言及している文献は見つかっていない。
【００１０】
本開示は、上記新たな知見に基づく課題に着目してなされたものであって、その目的は、算出精度を向上させた緩和スペクトル及びそれを用いた物性値を算出する方法、装置及びプログラムを提供することである。
【課題を解決するための手段】
【００１１】
本開示は、上記目的を達成するために、次のような手段を講じている。
【００１２】
すなわち、本開示の緩和スペクトルを算出する方法は、
予め設定された解析条件下において予め設定された分子モデルデータを用いた分子動力学計算を行い、各解析時点における応力データを算出するステップと、
前記応力データに基づき各解析時点毎の緩和弾性率を算出するステップと、
前記緩和弾性率を補間関数で近似した近似値を算出するステップと、
前記緩和弾性率の近似値を微分処理して緩和スペクトルを算出するステップと、
を含む。
【００１３】
本開示の緩和スペクトルを算出する装置は、
予め設定された解析条件下において予め設定された分子モデルデータを用いた分子動力学計算の実行結果に基づき、各解析時点における応力データを算出する応力算出部と、
前記応力データに基づき各解析時点毎の緩和弾性率を算出する緩和弾性率算出部と、
前記緩和弾性率を補間関数で近似した近似値を算出する近似部と、
前記緩和弾性率の近似値を微分処理して緩和スペクトルを算出する緩和スペクトル算出部と、
を備える。
【００１４】
このように、離散データである緩和弾性率を補間関数で近似するので、緩和弾性率のバラツキが平滑化され、微分処理による精度の悪化を抑制し、精度を向上させることが可能となる。
【図面の簡単な説明】
【００１５】
【図１】本開示の緩和スペクトル及び物性値を算出する装置を示すブロック図。
【図２】本開示の装置で実行される算出処理ルーチンを示すフローチャート。
【図３】緩和弾性率及びそれを補間関数で近似した近似値を示す図。
【図４Ａ】Ｂスプライン関数で近似した緩和弾性率Ｇ（ｔ）を１次の式（３）で緩和スペクトルＨ（τ）に変換し、式（４）のラプラス変換で緩和弾性率Ｇ（ｔ）に戻したものを示す図。
【図４Ｂ】Ｂスプライン関数で近似した緩和弾性率Ｇ（ｔ）を２次の式（５）で緩和スペクトルＨ（τ）に変換し、式（４）のラプラス変換で緩和弾性率Ｇ（ｔ）に戻したものを示す図。
【図４Ｃ】Ｂスプライン関数で近似した緩和弾性率Ｇ（ｔ）を３次の式（６）で緩和スペクトルＨ（τ）に変換し、式（４）のラプラス変換で緩和弾性率Ｇ（ｔ）に戻したものを示す図。
【図５】緩和弾性率Ｇ（ｔ）をそのまま１次の式（３）で緩和スペクトルＨ（τ）に変換し、式（４）のラプラス変換で緩和弾性率Ｇ（ｔ）に戻したものを示す図。
【図６Ａ】緩和弾性率Ｇ（ｔ）をそのまま１次の式（３）で緩和スペクトルＨ（τ）に変換し、式（４）のラプラス変換で緩和弾性率Ｇ（ｔ）に戻したものを示す図。
【図６Ｂ】緩和弾性率Ｇ（ｔ）をそのまま２次の式（５）で緩和スペクトルＨ（τ）に変換し、式（４）のラプラス変換で緩和弾性率Ｇ（ｔ）に戻したものを示す図。
【図６Ｃ】緩和弾性率Ｇ（ｔ）をそのまま３次の式（６）で緩和スペクトルＨ（τ）に変換し、式（４）のラプラス変換で緩和弾性率Ｇ（ｔ）に戻したものを示す図。
【発明を実施するための形態】
【００１６】
以下、本開示の一実施形態を、図面を参照して説明する。
【００１７】
［緩和スペクトル及び物性値を算出する装置］
本実施形態の装置は、予め設定された分子モデルデータを用いた分子動力学計算を行い、緩和スペクトルを算出する装置である。また、装置は、算出した緩和スペクトルを用いて物性値（貯蔵弾性率、損失弾性率及び損失正接）を算出する。勿論、緩和スペクトルのみを算出するように構成してもよいし、貯蔵弾性率、損失弾性率及び損失正接のうち少なくとも１つを算出するように構成してもよい。
【００１８】
図１に示すように、装置１は、設定部１０と、分子動力学計算実行部１１と、応力算出部１２と、緩和弾性率算出部１３と、近似部１４と、緩和スペクトル算出部１５と、物性値算出部１６と、を有する。これら各部１０〜１６は、ＣＰＵ、メモリ、各種インターフェイス等を備えたパソコン等の情報処理装置において予め記憶されている図３に示す物性値算出処理ルーチンをＣＰＵが実行することによりソフトウェア及びハードウェアが協働して実現される。
【００１９】
図１に示す設定部１０は、キーボードやマウス等の既知の操作部を介してユーザからの操作を受け付け、解析対象となる分子モデル（ゴムモデルなど）に関する情報の設定、温度・体積などの分子動力学計算に用いる各種解析条件の設定を実行し、これらをメモリに記憶する。
【００２０】
図１に示す分子動力学計算実行部１１は、予め設定された分子モデルデータを用いて予め設定された解析条件下において分子動力学計算を行い、各粒子の座標及び速度の時系列データを算出する。一例として、ランジュバン動力学のアルゴリズムによる定温における分子動力学計算などが挙げられる。
【００２１】
図１に示す応力算出部１２は、分子動力学計算実行部１１の算出結果（時系列データ）に基づき、各解析時間における応力データを算出する。具体的には、各解析時点における応力テンソルを算出する。
【００２２】
図１に示す緩和弾性率算出部１３は、応力データに基づき各解析時点毎の緩和弾性率を算出する。具体的には、時系列データから、せん断応力テンソルの自己相関関数を計算し、式（１）に示すグリーン・久保公式に代入して緩和弾性率Ｇ（ｔ）を算出する。図３は、鎖長１００の高分子１００本で構成される高分子メルトモデルを用いて算出した緩和弾性率Ｇ（ｔ）をバツ印で示す。同図は縦軸及び横軸共にログスケールに設定している。同図に示すように、解析時点毎に緩和弾性率の値が存在する離散的な時系列データである。
【数５】
0006579904.tif 000006
σ（ｔ）は時刻ｔのせん断応力テンソル、Ｖは体積、Ｔは温度、ＫＢはボルツマン定数である。
【００２３】
図１に示す近似部１４は、緩和弾性率Ｇ（ｔ）を補間関数で近似した近似値を算出する。補間関数としては、Ｂスプライン、ベジェ、非一様有理Ｂスプライン（Non-Uniform Rational Basis Spline）、スプライン、ラグランジュ補間が挙げられる。本実施形態では、３次のＢスプラインを用いているが、これに限定されず、２次のＢスプラインも使用できる。この中では最もＢスプラインが取り扱いしやすい。ベジェは多数の制御点を持つ場合には接続条件が煩雑になるからである。スプラインは、滑らかな曲線が得られるものの、弾性緩和率の値をそのまま制御点にすれば計算結果の振動（バラツキ）をダイレクトに拾ってしまうからである。ラグランジュ補間は制御点が多くなると振動性が出てしまうからである。
【００２４】
図３は、緩和弾性率Ｇ（ｔ）をバツ印で示すと共に、緩和弾性率Ｇ（ｔ）を３次のＢスプライン関数で近似し、得られる近似値（Ｂスプライン関数の値）を実線で描いた図である。バツ印で示す離散的データである緩和弾性率Ｇ（ｔ）が平滑化され、バラツキが低減されていることが分かる。
【００２５】
図１に示す緩和スペクトル算出部１５は、近似値を微分処理して緩和スペクトルＨ（τ）を算出する。具体的には、Ｓｃｈｗａｒｚｌ−Ｓｔａｖｅｒｍａｎによる近似式（２）を用い、緩和弾性率Ｇ（ｔ）の近似値（Ｂスプライン関数の値）を微分処理して緩和スペクトルＨ（τ）を算出する。式（３）は１次、式（５）は２次、式（６）は３次である。
【数６】
0006579904.tif 000007
【数７】
0006579904.tif 000008
【数８】
0006579904.tif 000009
【００２６】
図４Ａ〜Ｃは、Ｂスプライン関数で近似した緩和弾性率Ｇ（ｔ）を、１次の式（３）、２次の式（５）、３次の式（６）で緩和スペクトルＨ（τ）に変換し、それらを式（４）のラプラス変換で緩和弾性率Ｇ（ｔ）復元に戻したものを示す。同４Ａに示す１次の式（２）による算出結果の残差二乗和が0.0324である。図４Ｂに示す２次の式（２）による算出結果の残差二乗和が0.0173である。図４Ｃに示す３次の式（２）による算出結果の残差二乗和が0.0133であった。次数が上がるほど残差が小さくなり、精度が向上している。図６Ａ〜Ｃに比べても残差が小さくなり、精度が向上している。精度が向上した原因について、発明者は、次のように考えている。従来では、緩和弾性率データにはバラツキが含まれており、緩和弾性率Ｇ（ｔ）から緩和スペクトルＨ（τ）を算出する際に行う微分処理によって、算出精度が悪化していた。本開示の方法では、緩和弾性率Ｇ（ｔ）を補間関数で近似しているので、緩和弾性率Ｇ（ｔ）の振動成分（バラツキ）が平滑化され、悪化の原因が低減されたためと考えている。
【００２７】
したがって、式（２）に表されるＳｃｈｗａｒｚｌ−Ｓｔａｖｅｒｍａｎによる近似式を用いており、次数は１以上であればよい。より一層精度を向上させるのであれば２次以上がよい。１次から２次への向上代が最も大きく、２次から３次への向上代は小さくなる。よって、演算コストに対する精度の向上代を考慮すれば、２次が好ましい。
【００２８】
図１に示す物性値算出部１６は、緩和スペクトルＨ（τ）に基づき、貯蔵弾性率Ｇ’（ω）、損失弾性率Ｇ’’（ω）及び損失正接ｔａｎδのうち少なくとも１つの物性値を算出する。具体的には、貯蔵弾性率Ｇ’（ω）を式（７）で、損失弾性率Ｇ’’（ω）を式（８）で算出する。損失正接ｔａｎδは式（９）で算出する。ωは周波数であり、δは歪と応力の位相差を示す。なお、物性値算出部１６は実装に応じて省略可能である。
【数９】
0006579904.tif 000010
【００２９】
［緩和スペクトル及び物性値を算出する装置］
上記装置１を用い、緩和スペクトルを算出する方法を説明する。また、緩和スペクトルを用いて物性値を算出する方法を、図２を参照して説明する。
【００３０】
まず、ステップＳＴ１において、初期設定部１０は、解析対象となる分子モデルデータ、温度・体積などの分子動力学計算に用いる各種解析条件の設定を実行し、これらをメモリに記憶する。
【００３１】
ステップＳＴ２において、図１に示す分子動力学計算実行部１１は、メモリに記憶されている分子モデルデータを用いて予め設定された解析条件下において分子動力学計算を行い、各粒子の座標及び速度の時系列データを算出する。
【００３２】
次のステップＳＴ３において、応力算出部１２は、時系列データに基づき各解析時間における応力データを算出する。次のステップＳＴ４において、緩和弾性率算出部１３は、応力データに基づき緩和弾性率Ｇ（ｔ）を算出する。
【００３３】
次のステップＳＴ５において、近似部１４は、緩和弾性率Ｇ（ｔ）を補間関数で近似した近似値を算出する。
【００３４】
次のステップＳＴ６において、緩和スペクトル算出部１５は、緩和弾性率Ｇ（ｔ）の近似値を微分処理して緩和スペクトルＨ（τ）を算出する。
【００３５】
次のステップＳＴ７において、物性値算出部１６は、緩和スペクトルＨ（τ）に基づき、貯蔵弾性率Ｇ’（ω）、損失弾性率Ｇ’’（ω）及び損失正接ｔａｎδのうち少なくとも１つの物性値を算出する。
【００３６】
以上のように、本実施形態の緩和スペクトルを算出する方法は、
予め設定された解析条件下において予め設定された分子モデルデータを用いた分子動力学計算を行い、各解析時点における応力データを算出するステップ（ＳＴ２〜ＳＴ３）と、
応力データに基づき各解析時点毎の緩和弾性率Ｇ（ｔ）を算出するステップ（ＳＴ４）と、
緩和弾性率Ｇ（ｔ）を補間関数で近似した近似値を算出するステップ（ＳＴ５）と、
緩和弾性率Ｇ（ｔ）の近似値を微分処理して緩和スペクトルＨ（τ）を算出するステップ（ＳＴ６）と、
を含む。
【００３７】
本実施形態の緩和スペクトルを算出する装置１は、
予め設定された解析条件下において予め設定された分子モデルデータを用いた分子動力学計算の実行結果に基づき、各解析時点における応力データを算出する応力算出部１２と、
応力データに基づき各解析時点毎の緩和弾性率Ｇ（ｔ）を算出する緩和弾性率算出部１３と、
緩和弾性率Ｇ（ｔ）を補間関数で近似した近似値を算出する近似部１４と、
緩和弾性率Ｇ（ｔ）の近似値を微分処理して緩和スペクトルＨ（τ）を算出する緩和スペクトル算出部１５と、
を備える。
【００３８】
このように、離散データである緩和弾性率Ｇ（ｔ）を補間関数で近似するので、緩和弾性率Ｇ（ｔ）のバラツキが平滑化され、微分処理による精度の悪化を抑制し、精度を向上させることが可能となる。
【００３９】
本実施形態において、補間関数はＢスプライン関数である。
この構成によれば、緩和弾性率Ｇ（ｔ）の各値をそのまま制御点に使用できるので、実装が容易である。また、滑らかな曲線が得られ、且つ、他の補間関数に比べてバラツキを抑制できるので、他の補間関数に比べて、精度を向上させることが可能となる。
【００４０】
本実施形態において、緩和弾性率Ｇ（ｔ）の近似値から緩和スペクトルＨ（τ）を算出する演算は、式（２）に表されるＳｃｈｗａｒｚｌ−Ｓｔａｖｅｒｍａｎによる近似式を用いており、次数ｋは２次以上が好ましい。
【数１０】
0006579904.tif 000011
Ｇは緩和弾性率である。
この構成によれば、１次よりも精度を向上させることが可能となる。
【００４１】
本実施形態の方法において、緩和スペクトルＨ（τ）に基づき、貯蔵弾性率Ｇ’（ω）、損失弾性率Ｇ’’（ω）及び損失正接ｔａｎδのうち少なくとも１つの物性値を算出するステップ（ＳＴ７）を更に含む。
本実施形態の装置において、緩和スペクトルＨ（τ）に基づき、貯蔵弾性率Ｇ’（ω）、損失弾性率Ｇ’’（ω）及び損失正接ｔａｎδのうち少なくとも１つの物性値を算出する物性値算出部１６を更に備える。
この構成によれば、従来に比べて精度のよい緩和スペクトルＨ（τ）を用いるので、物性値の算出精度を向上させることが可能となる。
【００４２】
本実施形態に係るコンピュータプログラムは、上記方法を構成する各ステップをコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。言い換えると、上記方法を使用しているとも言える。
【００４３】
以上、本開示の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本開示の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。
【００４４】
例えば、図１に示す各部１０〜１６は、所定プログラムをコンピュータのＣＰＵで実行することで実現しているが、各部を専用メモリや専用回路で構成してもよい。
【００４５】
上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本開示の趣旨を逸脱しない範囲で種々変形が可能である。
【符号の説明】
【００４６】
１２…応力算出部
１３…緩和弾性率算出部
１４…近似部
１５…緩和スペクトル算出部
１６…物性値算出部
【図１】
0006579904.tif 000014
【図２】
0006579904.tif 000015
【図３】
0006579904.tif 000016
【図４Ａ】
0006579904.tif 000017
【図４Ｂ】
0006579904.tif 000018
【図４Ｃ】
0006579904.tif 000019
【図５】
0006579904.tif 000020
【図６Ａ】
0006579904.tif 000021
【図６Ｂ】
0006579904.tif 000022
【図６Ｃ】
0006579904.tif 000023