JP6385113B2 - 高分子モデルの緩和弾性率を解析する装置、方法及びコンピュータプログラム。 - Google Patents

Description

本発明は、粘弾性（複素弾性率）を算出するために有用な、高分子モデルの緩和弾性率を解析する装置、方法及びコンピュータプログラムに関する。

高分子モデルなどの粘弾性材料の複素弾性率（貯蔵弾性率、損失弾性率、損失正接）は、緩和弾性率を表す式をフーリエ変換することで得られることが知られている。緩和弾性率は、グリーン・久保公式によりシステム全体の応力に基づいて算出される。緩和弾性率を表現するための級数近似として、例えば特許文献１のように、プロニー級数近似を行うことが知られている。特許文献１では、各級数の係数及び緩和時間の決定方法について開示しているようである。

さらに、進んだ解析方法として、本発明の発明者らは、非特許文献１において、システムをいくつかのグループに分割しておき、緩和弾性率を級数近似したときの各係数に対する各グループの寄与を算出する可能に、固有値問題を解くことを述べている。

特開２００３−３４５７８０号公報

分子動力学による高分子粘弾性の解析法に関する研究，日野 理、高橋 宏幸，日本ゴム協会２０１３年年次大会講演要旨集，８７頁

高分子モデルの緩和弾性率を解析するにあたり、上記特許文献１の方法では、システム全体の緩和弾性率を表現する各級数の緩和時間τ_ｐと係数を最小二乗法などによりフィティングしているため、フーリエ変換により複素弾性率を適切に得ることができる。しかし、緩和弾性率を級数近似したときの各係数に対する、システムを構成する各グループの寄与を算出することはできない。

一方、非特許文献１の方法では、固有値問題を解くことで得られた解を用いれば、緩和弾性率を級数近似したときの各係数に対する、システムを構成する各グループの寄与を算出する可能である。しかし、級数近似した結果をフーリエ変換して粘弾性（複素弾性率）を算出し、算出結果をプロットした場合に、図６Ａに示すように、本来では有るはずのない凹凸部がプロット結果に表れてしまう。このような凹凸部は、固有値問題に対応可能に近似する課程において人為的に生じたと考えられる。事実と異なり、このような間違ったデータに基づき設計を進めると不具合が生じるおそれがある。

本発明は、このような課題に着目してなされたものであって、その目的は、緩和弾性率を級数近似したときの各係数に対する、システムを構成する各グループの寄与を算出可能であり、且つ、フーリエ変換による適切な複素弾性率を算出可能に、高分子モデルの緩和弾性率を解析する装置、方法及びコンピュータプログラムを提供することである。

本発明は、上記目的を達成するために、次のような手段を講じている。

すなわち、本発明の高分子モデルの緩和弾性率を解析する装置は、高分子モデルをＮ個のグループに分割するグループ設定部と、予め設定された解析条件を用いて平衡状態における前記高分子モデルの分子動力学計算を行い、前記グループ単位での応力の時間発展データを算出し、前記グループ単位の緩和弾性率Ｇ_ｊｉ（ｔ）［ｔは時間、ｉ，ｊはグループに付けた番号を表す］を算出する緩和弾性率算出部と、時刻ｔにおけるグループ毎の応力Ｇ_ｊｉ（ｔ）を要素とするＮ×Ｎ個の行列Ａ（ｔ）を生成する行列生成部と、終了時刻（ｔ_０＋ｔ）における行列Ａ（ｔ_０＋ｔ）を、初期時刻ｔ_０における行列Ａ（ｔ_０）と指数関数と固有ベクトルｇ_ｐを含むプロニー級数形式で表した固有値方程式を解き、行列Ａ（ｔ）を近似する固有値問題解決部と、を備え、前記指数関数として、ｅｘｐ｛−［（ｔ_０＋ｔ）^β＋ｔ_０ ^β］／τ_ｐ｝で表されるＫＷＷ（Kohlausch-Williams-Watts）関数を用いることを特徴とする。
ただし、τ_ｐは緩和時間を表す。

本発明の高分子モデルの緩和弾性率を解析する方法は、コンピュータが実行する方法であって、高分子モデルをＮ個のグループに分割するステップと、予め設定された解析条件を用いて平衡状態における前記高分子モデルの分子動力学計算を行い、前記グループ単位での応力の時間発展データを算出し、前記グループ単位の緩和弾性率Ｇ_ｊｉ（ｔ）［ｔは時間、ｉ，ｊはグループに付けた番号を表す］を算出するステップと、時刻ｔにおけるグループ毎の応力Ｇ_ｊｉ（ｔ）を要素とするＮ×Ｎ個の行列Ａ（ｔ）を生成するステップと、終了時刻（ｔ_０＋ｔ）における行列Ａ（ｔ_０＋ｔ）を、初期時刻ｔ_０における行列Ａ（ｔ_０）と指数関数と固有ベクトルｇ_ｐを含むプロニー級数形式で表した固有値方程式を解き、行列Ａ（ｔ）を近似するステップと、を含み、前記指数関数として、ｅｘｐ｛−［（ｔ_０＋ｔ）^β＋ｔ_０ ^β］／τ_ｐ｝で表されるＫＷＷ（Kohlausch-Williams-Watts）関数を用いることを特徴とする。
ただし、τ_ｐは緩和時間を表す。

このように、周波数空間において、従来の指数関数よりも広いＫＷＷ関数を用いるので、得られた近似式をフーリエ変換して複素弾性率を算出した場合に、従来の指数関数で生じていたデータの欠けが生じることなく、複素弾性率を適切に得ることが可能となる。したがって、複素弾性率を適切に得ることができる点と、級数近似したときの各成分がどういう緩和時間に寄与しているかを表現できる点とを両立できる。

具体的な前記固有値方程式は、次の式（６）で表される。
Ａ（ｔ_０＋ｔ）ｇ_ｐ＝ｅｘｐ｛−［（ｔ_０＋ｔ）^β＋ｔ_０ ^β］／τ_ｐ｝Ａ（ｔ_０）ｇ_ｐ …（６）

ＫＷＷ関数のパラメータβを適切な値に設定するためには、前記ＫＷＷ関数のパラメータβは、次の式（９）を満足する値に設定されることが好ましい。
Ｇ（ｔ）＝Ｇ（ｔ_０）ｅｘｐ｛−［（ｔ＋ｔ_０）^β＋ｔ_０ ^β］／２ｔ_０｝ …（９）

本発明は、上記方法を構成するステップをコンピュータに実行させるプログラムとして特定可能である。

本発明の高分子モデルの緩和弾性率を解析する装置を示すブロック図。 装置が実行する解析処理ルーチンを示すフローチャート。 高分子モデルの一例を示す図。 図３Ａに示すモデルを用いた分子動力学計算により得られた緩和弾性率を示す図。 応力緩和を表現する指数緩和に関する説明図。 装置が実行する固有値解析処理ルーチンを示すフローチャート。 従来法の級数近似をフーリエ変換して算出した複素弾性率を示す図。 本発明の級数近似をフーリエ変換して算出した複素弾性率を示す図。

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

［高分子モデルの緩和弾性率を解析する装置］
本実施形態の装置は、粘弾性（複素弾性率）を算出するために有用な、高分子モデルの緩和弾性率を解析する装置である。具体的には、装置は、分子動力学計算を行い、時間と共に変化する応力を算出し、当該応力に基づいて緩和弾性率を算出する。得られる緩和弾性率は離散的なデータであるため、プロニー級数近似及び固有値問題を解き、時間的に連続的なデータを得る。得られた結果をフーリエ変換すれば、粘弾性、すなわち複素弾性率（貯蔵弾性率、損失弾性率、損失正接）を得ることが可能となる。

図１に示すように、装置１は、グループ設定部１０と、緩和弾性率算出部１１と、行列生成部１２と、固有値問題解決部１３と、β決定部１４と、を有する。これら各部１０〜１４は、ＣＰＵ、メモリ、各種インターフェイス等を備えたパソコン等の情報処理装置において予め記憶されている図２に示す解析処理ルーチンをＣＰＵが実行することによりソフトウェア及びハードウェアが協働して実現される。

図１に示すグループ設定部１０は、メモリに記憶されている高分子モデルをＮ個のグループに分割する設定を行う（図３Ａ参照）。図３Ａの例は、５００個の粒子で構成される複数の分子鎖を、Ｎ＝１０個のグループに分割した図である。なお、メモリには、キーボードやマウス等の既知の操作部を介してユーザからの操作を受け付け、高分子モデルに関する情報の設定、緩和応力を計算するための分子動力学計算に用いる各種解析条件が設定されている。

図１に示す緩和弾性率算出部１１は、予め設定された解析条件を用いて平衡状態における高分子モデルの分子動力学計算を行い、グループ単位での応力σ_ｉ（ｔ）の時間発展データを算出し、グループ単位の緩和弾性率Ｇ_ｊｉ（ｔ）［ｔは時間、ｉ，ｊはグループに付けた番号を表す］を算出する。

緩和弾性率Ｇ（ｔ）は、グリーン・久保公式に基づき次の式（１）で算出可能である。
Ｖは体積、Ｔは温度、σはシステム（系）全体の応力、Ｓ_全時間は全解析時間における単位時間の数を表す。
ここで、システム全体の応力σ（ｔ）は、次の式（２）のように、グループ毎の部分応力σ_ｉ（ｔ）として表せる。
上記２式は、次の式（３），（４）のように相関行列として表現できる。

図３Ｂは、図３Ａに示すモデルを用いて分子動力学計算により算出した緩和弾性率を示す。横軸が時間［τ］を示し、縦軸が弾性率［ε／σ^３］を示す。

図１に示す行列生成部１２は、時刻ｔにおけるグループ毎の応力Ｇ_ｊｉ（ｔ）を要素とするＮ×Ｎ個の行列Ａ（ｔ）を生成する。行列Ａ（ｔ）は、次の式（５）で表される。

図１に示す固有値問題解決部１３は、終了時刻（ｔ_０＋ｔ）における行列Ａ（ｔ_０＋ｔ）と固有ベクトルｇ _ｐ の積を、初期時刻ｔ_０における行列Ａ（ｔ_０）と指数関数で表される固有値λ _ｐ と固有ベクトルｇ_ｐ との積で表した固有値方程式を解き、行列Ａ（ｔ）を近似する。上記固有値方程式は、次の式（６）〜（８）で表される。

ここで、λ_ｐは固有値、ｇ_ｐは固有ベクトル、τ_ｐは緩和時間である。

ｅｘｐ｛−［（ｔ_０＋ｔ）^β −ｔ_０ ^β］／τ_ｐ｝は、拡張された指数関数（Stretched exponential function）と呼ばれるＫＷＷ（Kohlausch-Williams-Watts）関数である。ＫＷＷ関数は、通常の指数関数に比べて減衰が遅い関数であり、フーリエ変換した周波数領域において広いスペクトルを持つ特徴を有する。

ＫＷＷ関数のパラメータβは１よりも小さい値であるが、初期時刻ｔ_０と終了時刻（ｔ_０＋ｔ）に応じた適切な値を採用するためには、次の式（９）を満足する値に設定されるのが好ましい。図１に示すβ決定部１４は、次の式（９）を数値的に解き、同式（９）を満足するβを決定する。
パラメータβは小さいほど減衰が遅くなり、周波数領域でスペクトルが広くなる。本実施形態の例では、０．５〜０．７の値が好ましかった。なお、本実施形態では、β決定部１４を設けてβを動的に決定しているが、パラメータβを予め決定しておき、メモリに設定しておいてもよい。

具体的に、固有値問題解決部１３は、式（６）に示す固有値問題を解くにあたり、長時間部分から順次緩和解析を行い、全領域の応力緩和を、個々の緩和解析結果の和で表現する（プロニー級数表示）。すなわち、図４に示すように、緩和時間（全体）は、緩和解析（１）〜（３）の総和で再現される。固有値問題の解法は、非特許文献１と同じであるので、詳細は述べないが、簡単に説明すれば、図５のフローに基づく。すなわち、ステップＳＴ１１において、長時間部分（終端緩和付近）に初期時刻ｔ_０及び終了時刻（ｔ_０＋ｔ）を設定し、固有値問題解決部１３が、ステップＳＴ１２において式（６）の固有値問題を解き、最大Ｎ組の固有値λ_ｐと固有ベクトルｇ_ｐを得る。次のステップＳＴ１３において、式（１０）でＡ（ｔ）を近似した結果をＡ’（ｔ）とする。
ここで、緩和弾性率は、次の式（１３）で近似される。
なお、グループｉのｇ_ｐへの寄与は、式（１３）の第ｉ行の要素の和で定義できる。

次のステップＳＴ１４において、全時間を近似できたか否かを判定する。近似できていない場合には、ステップＳＴ１５において、Ａ（ｔ）−Ａ’（ｔ）を新たなＡ（ｔ）に更新し、ステップＳＴ１１の処理に戻る。

［高分子モデルの緩和弾性率を解析する方法］
上記装置１を用いて緩和弾性率を解析する方法について、図２を用いて説明する。
まず、ステップＳＴ１において、グループ設定部１０が、高分子モデルをＮ個のグループに分割する。次のステップＳＴ２において、緩和弾性率算出部１１が、予め設定された解析条件を用いて平衡状態における高分子モデルの分子動力学計算を行い、グループ単位での応力の時間発展データを算出し、グループ単位の緩和弾性率を算出する。次のステップＳＴ３において、行列生成部１２が、時刻ｔにおけるグループ毎の応力Ｇ_ｊｉ（ｔ）を要素とするＮ×Ｎ個の行列Ａ（ｔ）を生成する。次のステップＳＴ４において、固有値問題解決部１３が、終了時刻（ｔ_０＋ｔ）における行列Ａ（ｔ_０＋ｔ）を、初期時刻ｔ_０における行列Ａ（ｔ_０）と指数関数と固有ベクトルｇ_ｐを含むプロニー級数形式で表した固有値方程式を解き、行列Ａ（ｔ）を近似する。

本発明の効果を示すために、従来（非特許文献１）との比較を示す。従来では、次の式（１４）を用いて固有値問題を解いていた。
従来の方法で得られた級数近似の結果をフーリエ変換すると、図６Ａに示すように、損失正接（ｔａｎδ）において、本来では有るはずのない凹凸部（図中にて円で囲んで示す）がプロット結果に表れてしまう。これは、応力緩和を指数関数の足し算で表現するにあたり、関数の緩和時間軸における分布に粗密があると周波数空間において強度がない部分ができやすく、その結果、スペクトルが無い部分ができてしまい、従来方法では凹凸部が生じてしまうと考えられる。従来のように一般的な指数関数を用いた緩和解析法では、上記凹凸部のような欠けが生じやすい、緩和弾性率のような実時間の値では問題がないが、フーリエ変換で得られる複素弾性率のような周波数空間では問題となる。

一方、本発明のようにＫＷＷ関数を用いれば、図６Ｂに示すように、上記凹凸部は出現しない。すなわち、本発明は、特許文献１に記載の利点（複素弾性率を適切に得ることができる点）と、非特許文献１に記載の利点（級数近似したときの各成分がどういう緩和時間に寄与しているかを表現できる点）とを両立する解析方法であることが分かる。

以上のように、本実施形態の高分子モデルの緩和弾性率を解析する装置は、高分子モデルをＮ個のグループに分割するグループ設定部１０と、予め設定された解析条件を用いて平衡状態における高分子モデルの分子動力学計算を行い、グループ単位での応力の時間発展データを算出し、グループ単位の緩和弾性率Ｇ_ｊｉ（ｔ）［ｔは時間、ｉ，ｊはグループに付けた番号を表す］を算出する緩和弾性率算出部１１と、時刻ｔにおけるグループ毎の応力Ｇ_ｊｉ（ｔ）を要素とするＮ×Ｎ個の行列Ａ（ｔ）を生成する行列生成部１２と、終了時刻（ｔ_０＋ｔ）における行列Ａ（ｔ_０＋ｔ）を、初期時刻ｔ_０における行列Ａ（ｔ_０）と指数関数と固有ベクトルｇ_ｐを含むプロニー級数形式で表した固有値方程式を解き、行列Ａ（ｔ）を近似する固有値問題解決部１３と、を有する。前記指数関数として、ｅｘｐ｛−［（ｔ_０＋ｔ）^β＋ｔ_０ ^β］／τ_ｐ｝で表されるＫＷＷ（Kohlausch-Williams-Watts）関数を用いる。

本実施形態の高分子モデルの緩和弾性率を解析する方法は、コンピュータが実行する方法であって、高分子モデルをＮ個のグループに分割するステップ（ＳＴ１）と、予め設定された解析条件を用いて平衡状態における高分子モデルの分子動力学計算を行い、グループ単位での応力の時間発展データを算出し、グループ単位の緩和弾性率Ｇ_ｊｉ（ｔ）［ｔは時間、ｉ，ｊはグループに付けた番号を表す］を算出するステップ（ＳＴ２）と、時刻ｔにおけるグループ毎の応力Ｇ_ｊｉ（ｔ）を要素とするＮ×Ｎ個の行列Ａ（ｔ）を生成するステップ（ＳＴ３）と、終了時刻（ｔ_０＋ｔ）における行列Ａ（ｔ_０＋ｔ）を、初期時刻ｔ_０における行列Ａ（ｔ_０）と指数関数と固有ベクトルｇ_ｐを含むプロニー級数形式で表した固有値方程式を解き、行列Ａ（ｔ）を近似するステップ（ＳＴ４）と、を含む。前記指数関数として、ｅｘｐ｛−［（ｔ_０＋ｔ）^β＋ｔ_０ ^β］／τ_ｐ｝で表されるＫＷＷ（Kohlausch-Williams-Watts）関数を用いる。

このように、周波数空間において、従来の指数関数よりも広いＫＷＷ関数を用いるので、得られた近似式をフーリエ変換して複素弾性率を算出した場合に、従来の指数関数で生じていたデータの欠けが生じることなく、複素弾性率を適切に得ることが可能となる。したがって、複素弾性率を適切に得ることができる点と、級数近似したときの各成分がどういう緩和時間に寄与しているかを表現できる点とを両立できる。

本実施形態では、具体的に、前記固有値方程式は、次の式（６）で表される。
Ａ（ｔ_０＋ｔ）ｇ_ｐ＝ｅｘｐ｛−［（ｔ_０＋ｔ）^β＋ｔ_０ ^β］／τ_ｐ｝Ａ（ｔ_０）ｇ_ｐ …（６）

本実施形態では、前記ＫＷＷ関数のパラメータβは、次の式（９）を満足する値に設定している。
Ｇ（ｔ）＝Ｇ（ｔ_０）ｅｘｐ｛−［（ｔ＋ｔ_０）^β＋ｔ_０ ^β］／２ｔ_０｝ …（９）
これによれば、ＫＷＷ関数のパラメータβを適切な値に設定可能となる。

本実施形態に係るコンピュータプログラムは、上記方法を構成する各ステップをコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。言い換えると、上記方法を使用しているとも言える。

以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。

上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。

１０…グループ設定部
１１…緩和弾性率算出部
１２…行列生成部
１３…固有値問題解決部
１４…β決定部

Claims (5)

1. 高分子モデルをＮ個のグループに分割するグループ設定部と、
   予め設定された解析条件を用いて平衡状態における前記高分子モデルの分子動力学計算を行い、前記グループ単位での応力の時間発展データを算出し、前記グループ単位の緩和弾性率Ｇ_ｊｉ（ｔ）［ｔは時間、ｉ，ｊはグループに付けた番号を表す］を算出する緩和弾性率算出部と、
   時刻ｔにおけるグループ毎の緩和弾性率Ｇ_ｊｉ（ｔ）を要素とするＮ×Ｎ個の行列Ａ（ｔ）を生成する行列生成部と、
   終了時刻（ｔ_０＋ｔ）における行列Ａ（ｔ_０＋ｔ）と固有ベクトルｇ _ｐ の積を、初期時刻ｔ_０における行列Ａ（ｔ_０）と指数関数で表される固有値λ _ｐ と固有ベクトルｇ_ｐ との積で表した固有値方程式を解いて最大Ｎ組の固有値λ _ｐ 及び固有ベクトルｇ _ｐ を取得し、取得した固有値λ _ｐ 及び固有ベクトルｇ _ｐ を用いて行列Ａ（ｔ）を近似する固有値問題解決部と、を備え、
   前記指数関数として、ＫＷＷ（Kohlausch-Williams-Watts）関数を用いることを特徴とする高分子モデルの緩和弾性率を解析する装置。
2. 前記固有値方程式は、次の式（６）で表される請求項１に記載の装置。
   Ａ（ｔ_０＋ｔ）ｇ_ｐ＝ｅｘｐ｛−［（ｔ_０＋ｔ）^β −ｔ_０ ^β］／τ_ｐ｝Ａ（ｔ_０）ｇ_ｐ …（６）
   ただし、τ _ｐ は緩和時間を表す。
3. コンピュータが実行する方法であって、
   高分子モデルをＮ個のグループに分割するステップと、
   予め設定された解析条件を用いて平衡状態における前記高分子モデルの分子動力学計算を行い、前記グループ単位での応力の時間発展データを算出し、前記グループ単位の緩和弾性率Ｇ_ｊｉ（ｔ）［ｔは時間、ｉ，ｊはグループに付けた番号を表す］を算出するステップと、
   時刻ｔにおけるグループ毎の緩和弾性率Ｇ_ｊｉ（ｔ）を要素とするＮ×Ｎ個の行列Ａ（ｔ）を生成するステップと、
   終了時刻（ｔ_０＋ｔ）における行列Ａ（ｔ_０＋ｔ）と固有ベクトルｇ _ｐ の積を、初期時刻ｔ_０における行列Ａ（ｔ_０）と指数関数で表される固有値λ _ｐ と固有ベクトルｇ_ｐ との積で表した固有値方程式を解いて最大Ｎ組の固有値λ _ｐ 及び固有ベクトルｇ _ｐ を取得し、取得した固有値λ _ｐ 及び固有ベクトルｇ _ｐ を用いて行列Ａ（ｔ）を近似するステップと、を含み、
   前記指数関数として、ＫＷＷ（Kohlausch-Williams-Watts）関数を用いることを特徴とする高分子モデルの緩和弾性率を解析する方法。
4. 前記固有値方程式は、次の式（６）で表される請求項３に記載の方法。
   Ａ（ｔ_０＋ｔ）ｇ_ｐ＝ｅｘｐ｛−［（ｔ_０＋ｔ）^β −ｔ_０ ^β］／τ_ｐ｝Ａ（ｔ_０）ｇ_ｐ …（６）
   ただし、τ _ｐ は緩和時間を表す。
5. 請求項３又は４に記載の方法をコンピュータに実行させるコンピュータプログラム。