(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】6385113
(24)【登録日】2018年8月17日
(45)【発行日】2018年9月5日
(54)【発明の名称】高分子モデルの緩和弾性率を解析する装置、方法及びコンピュータプログラム。
(51)【国際特許分類】
G06F 17/11 20060101AFI20180827BHJP
G06F 19/00 20180101ALI20180827BHJP
G01N 3/00 20060101ALI20180827BHJP
【FI】
G06F17/11
G06F19/00 110
G01N3/00 K
【請求項の数】5
【全頁数】11
(21)【出願番号】特願2014-78650(P2014-78650)
(22)【出願日】2014年4月7日
(65)【公開番号】特開2015-201009(P2015-201009A)
(43)【公開日】2015年11月12日
【審査請求日】2016年12月27日
(73)【特許権者】
【識別番号】000003148
【氏名又は名称】東洋ゴム工業株式会社
(74)【代理人】
【識別番号】110000729
【氏名又は名称】特許業務法人 ユニアス国際特許事務所
(72)【発明者】
【氏名】高橋 宏幸
(72)【発明者】
【氏名】日野 理
【審査官】
清木 泰
(56)【参考文献】
【文献】
特開2013−250242(JP,A)
【文献】
特開2003−345780(JP,A)
【文献】
河端昌也,森山史朗,會田裕昌,ETFEフィルムの粘弾性挙動について,[online],一般社団法人日本膜構造協会,2005年,all 8 pages,[平成29年12月21日検索],インターネット,URL,http://www.makukouzou.or.jp/article/pdf/2005/2005-01.pdf
【文献】
Yohji KAWASAKI, Hiroshi WATANABE, Takashi UNEYAMA,A Note for Kohlrausch-Williams-Watts Relaxation Function,Journal of the Society of Rheology Japan,日本,The Society of Rheology Japan,2011年,Vol:39, No:3,Pages:127-131
【文献】
日野理,高橋宏幸,分子動力学による高分子粘弾性の解析法に関する研究,日本ゴム協会2013年年次大会研究発表講演会講演要旨,一般社団法人日本ゴム協会,2013年 5月23日,Page:87
(58)【調査した分野】(Int.Cl.,DB名)
G06F17/10−17/18
G06F19/00
G06F17/50
G01N 3/00− 3/62
G01N17/00−19/10
(57)【特許請求の範囲】
【請求項1】
高分子モデルをN個のグループに分割するグループ設定部と、
予め設定された解析条件を用いて平衡状態における前記高分子モデルの分子動力学計算を行い、前記グループ単位での応力の時間発展データを算出し、前記グループ単位の緩和弾性率Gji(t)[tは時間、i,jはグループに付けた番号を表す]を算出する緩和弾性率算出部と、
時刻tにおけるグループ毎の緩和弾性率Gji(t)を要素とするN×N個の行列A(t)を生成する行列生成部と、
終了時刻(t0+t)における行列A(t0+t)と固有ベクトルgpの積を、初期時刻t0における行列A(t0)と指数関数で表される固有値λpと固有ベクトルgpとの積で表した固有値方程式を解いて最大N組の固有値λp及び固有ベクトルgpを取得し、取得した固有値λp及び固有ベクトルgpを用いて行列A(t)を近似する固有値問題解決部と、を備え、
前記指数関数として、KWW(Kohlausch-Williams-Watts)関数を用いることを特徴とする高分子モデルの緩和弾性率を解析する装置。
【請求項2】
前記固有値方程式は、次の式(6)で表される請求項1に記載の装置。
A(t0+t)gp=exp{−[(t0+t)β−t0β]/τp}A(t0)gp …(6)
ただし、τpは緩和時間を表す。
【請求項3】
コンピュータが実行する方法であって、
高分子モデルをN個のグループに分割するステップと、
予め設定された解析条件を用いて平衡状態における前記高分子モデルの分子動力学計算を行い、前記グループ単位での応力の時間発展データを算出し、前記グループ単位の緩和弾性率Gji(t)[tは時間、i,jはグループに付けた番号を表す]を算出するステップと、
時刻tにおけるグループ毎の緩和弾性率Gji(t)を要素とするN×N個の行列A(t)を生成するステップと、
終了時刻(t0+t)における行列A(t0+t)と固有ベクトルgpの積を、初期時刻t0における行列A(t0)と指数関数で表される固有値λpと固有ベクトルgpとの積で表した固有値方程式を解いて最大N組の固有値λp及び固有ベクトルgpを取得し、取得した固有値λp及び固有ベクトルgpを用いて行列A(t)を近似するステップと、を含み、
前記指数関数として、KWW(Kohlausch-Williams-Watts)関数を用いることを特徴とする高分子モデルの緩和弾性率を解析する方法。
【請求項4】
前記固有値方程式は、次の式(6)で表される請求項3に記載の方法。
A(t0+t)gp=exp{−[(t0+t)β−t0β]/τp}A(t0)gp …(6)
ただし、τpは緩和時間を表す。
【請求項5】
請求項3又は4に記載の方法をコンピュータに実行させるコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、粘弾性(複素弾性率)を算出するために有用な、高分子モデルの緩和弾性率を解析する装置、方法及びコンピュータプログラムに関する。
【背景技術】
【0002】
高分子モデルなどの粘弾性材料の複素弾性率(貯蔵弾性率、損失弾性率、損失正接)は、緩和弾性率を表す式をフーリエ変換することで得られることが知られている。緩和弾性率は、グリーン・久保公式によりシステム全体の応力に基づいて算出される。緩和弾性率を表現するための級数近似として、例えば特許文献1のように、プロニー級数近似を行うことが知られている。特許文献1では、各級数の係数及び緩和時間の決定方法について開示しているようである。
【0003】
さらに、進んだ解析方法として、本発明の発明者らは、非特許文献1において、システムをいくつかのグループに分割しておき、緩和弾性率を級数近似したときの各係数に対する各グループの寄与を算出する可能に、固有値問題を解くことを述べている。
【先行技術文献】
【特許文献】
【0004】
【特許文献1】特開2003−345780号公報
【非特許文献】
【0005】
【非特許文献1】分子動力学による高分子粘弾性の解析法に関する研究,日野 理、高橋 宏幸,日本ゴム協会2013年年次大会講演要旨集,87頁
【発明の概要】
【発明が解決しようとする課題】
【0006】
高分子モデルの緩和弾性率を解析するにあたり、上記特許文献1の方法では、システム全体の緩和弾性率を表現する各級数の緩和時間τ
pと係数を最小二乗法などによりフィティングしているため、フーリエ変換により複素弾性率を適切に得ることができる。しかし、緩和弾性率を級数近似したときの各係数に対する、システムを構成する各グループの寄与を算出することはできない。
【0007】
一方、非特許文献1の方法では、固有値問題を解くことで得られた解を用いれば、緩和弾性率を級数近似したときの各係数に対する、システムを構成する各グループの寄与を算出する可能である。しかし、級数近似した結果をフーリエ変換して粘弾性(複素弾性率)を算出し、算出結果をプロットした場合に、
図6Aに示すように、本来では有るはずのない凹凸部がプロット結果に表れてしまう。このような凹凸部は、固有値問題に対応可能に近似する課程において人為的に生じたと考えられる。事実と異なり、このような間違ったデータに基づき設計を進めると不具合が生じるおそれがある。
【0008】
本発明は、このような課題に着目してなされたものであって、その目的は、緩和弾性率を級数近似したときの各係数に対する、システムを構成する各グループの寄与を算出可能であり、且つ、フーリエ変換による適切な複素弾性率を算出可能に、高分子モデルの緩和弾性率を解析する装置、方法及びコンピュータプログラムを提供することである。
【課題を解決するための手段】
【0009】
本発明は、上記目的を達成するために、次のような手段を講じている。
【0010】
すなわち、本発明の高分子モデルの緩和弾性率を解析する装置は、高分子モデルをN個のグループに分割するグループ設定部と、予め設定された解析条件を用いて平衡状態における前記高分子モデルの分子動力学計算を行い、前記グループ単位での応力の時間発展データを算出し、前記グループ単位の緩和弾性率G
ji(t)[tは時間、i,jはグループに付けた番号を表す]を算出する緩和弾性率算出部と、時刻tにおけるグループ毎の応力G
ji(t)を要素とするN×N個の行列A(t)を生成する行列生成部と、終了時刻(t
0+t)における行列A(t
0+t)を、初期時刻t
0における行列A(t
0)と指数関数と固有ベクトルg
pを含むプロニー級数形式で表した固有値方程式を解き、行列A(t)を近似する固有値問題解決部と、を備え、前記指数関数として、exp{−[(t
0+t)
β+t
0β]/τ
p}で表されるKWW(Kohlausch-Williams-Watts)関数を用いることを特徴とする。
ただし、τ
pは緩和時間を表す。
【0011】
本発明の高分子モデルの緩和弾性率を解析する方法は、コンピュータが実行する方法であって、高分子モデルをN個のグループに分割するステップと、予め設定された解析条件を用いて平衡状態における前記高分子モデルの分子動力学計算を行い、前記グループ単位での応力の時間発展データを算出し、前記グループ単位の緩和弾性率G
ji(t)[tは時間、i,jはグループに付けた番号を表す]を算出するステップと、時刻tにおけるグループ毎の応力G
ji(t)を要素とするN×N個の行列A(t)を生成するステップと、終了時刻(t
0+t)における行列A(t
0+t)を、初期時刻t
0における行列A(t
0)と指数関数と固有ベクトルg
pを含むプロニー級数形式で表した固有値方程式を解き、行列A(t)を近似するステップと、を含み、前記指数関数として、exp{−[(t
0+t)
β+t
0β]/τ
p}で表されるKWW(Kohlausch-Williams-Watts)関数を用いることを特徴とする。
ただし、τ
pは緩和時間を表す。
【0012】
このように、周波数空間において、従来の指数関数よりも広いKWW関数を用いるので、得られた近似式をフーリエ変換して複素弾性率を算出した場合に、従来の指数関数で生じていたデータの欠けが生じることなく、複素弾性率を適切に得ることが可能となる。したがって、複素弾性率を適切に得ることができる点と、級数近似したときの各成分がどういう緩和時間に寄与しているかを表現できる点とを両立できる。
【0013】
具体的な前記固有値方程式は、次の式(6)で表される。
A(t
0+t)g
p=exp{−[(t
0+t)
β+t
0β]/τ
p}A(t
0)g
p …(6)
【0014】
KWW関数のパラメータβを適切な値に設定するためには、前記KWW関数のパラメータβは、次の式(9)を満足する値に設定されることが好ましい。
G(t)=G(t
0)exp{−[(t+t
0)
β+t
0β]/2t
0} …(9)
【0015】
本発明は、上記方法を構成するステップをコンピュータに実行させるプログラムとして特定可能である。
【図面の簡単な説明】
【0016】
【
図1】本発明の高分子モデルの緩和弾性率を解析する装置を示すブロック図。
【
図2】装置が実行する解析処理ルーチンを示すフローチャート。
【
図3B】
図3Aに示すモデルを用いた分子動力学計算により得られた緩和弾性率を示す図。
【
図4】応力緩和を表現する指数緩和に関する説明図。
【
図5】装置が実行する固有値解析処理ルーチンを示すフローチャート。
【
図6A】従来法の級数近似をフーリエ変換して算出した複素弾性率を示す図。
【
図6B】本発明の級数近似をフーリエ変換して算出した複素弾性率を示す図。
【発明を実施するための形態】
【0017】
以下、本発明の一実施形態を、図面を参照して説明する。
【0018】
[高分子モデルの緩和弾性率を解析する装置]
本実施形態の装置は、粘弾性(複素弾性率)を算出するために有用な、高分子モデルの緩和弾性率を解析する装置である。具体的には、装置は、分子動力学計算を行い、時間と共に変化する応力を算出し、当該応力に基づいて緩和弾性率を算出する。得られる緩和弾性率は離散的なデータであるため、プロニー級数近似及び固有値問題を解き、時間的に連続的なデータを得る。得られた結果をフーリエ変換すれば、粘弾性、すなわち複素弾性率(貯蔵弾性率、損失弾性率、損失正接)を得ることが可能となる。
【0019】
図1に示すように、装置1は、グループ設定部10と、緩和弾性率算出部11と、行列生成部12と、固有値問題解決部13と、β決定部14と、を有する。これら各部10〜14は、CPU、メモリ、各種インターフェイス等を備えたパソコン等の情報処理装置において予め記憶されている
図2に示す解析処理ルーチンをCPUが実行することによりソフトウェア及びハードウェアが協働して実現される。
【0020】
図1に示すグループ設定部10は、メモリに記憶されている高分子モデルをN個のグループに分割する設定を行う(
図3A参照)。
図3Aの例は、500個の粒子で構成される複数の分子鎖を、N=10個のグループに分割した図である。なお、メモリには、キーボードやマウス等の既知の操作部を介してユーザからの操作を受け付け、高分子モデルに関する情報の設定、緩和応力を計算するための分子動力学計算に用いる各種解析条件が設定されている。
【0021】
図1に示す緩和弾性率算出部11は、予め設定された解析条件を用いて平衡状態における高分子モデルの分子動力学計算を行い、グループ単位での応力σ
i(t)の時間発展データを算出し、グループ単位の緩和弾性率G
ji(t)[tは時間、i,jはグループに付けた番号を表す]を算出する。
【0022】
緩和弾性率G(t)は、グリーン・久保公式に基づき次の式(1)で算出可能である。
【数1】
Vは体積、Tは温度、σはシステム(系)全体の応力、S
全時間は全解析時間における単位時間の数を表す。
ここで、システム全体の応力σ(t)は、次の式(2)のように、グループ毎の部分応力σ
i(t)として表せる。
【数2】
上記2式は、次の式(3),(4)のように相関行列として表現できる。
【数3】
【0023】
図3Bは、
図3Aに示すモデルを用いて分子動力学計算により算出した緩和弾性率を示す。横軸が時間[τ]を示し、縦軸が弾性率[ε/σ
3]を示す。
【0024】
図1に示す行列生成部12は、時刻tにおけるグループ毎の応力G
ji(t)を要素とするN×N個の行列A(t)を生成する。行列A(t)は、次の式(5)で表される。
【数4】
【0025】
図1に示す固有値問題解決部13は、終了時刻(t
0+t)における行列A(t
0+t)
と固有ベクトルgpの積を、初期時刻t
0における行列A(t
0)と指数関数
で表される固有値λpと固有ベクトルg
pとの積で表した固有値方程式を解き、行列A(t)を近似する。上記固有値方程式は、次の式(6)〜(8)で表される。
【数5】
ここで、λ
pは固有値、g
pは固有ベクトル、τ
pは緩和時間である。
【0026】
exp{−[(t
0+t)
β−t
0β]/τ
p}は、拡張された指数関数(Stretched exponential function)と呼ばれるKWW(Kohlausch-Williams-Watts)関数である。KWW関数は、通常の指数関数に比べて減衰が遅い関数であり、フーリエ変換した周波数領域において広いスペクトルを持つ特徴を有する。
【0027】
KWW関数のパラメータβは1よりも小さい値であるが、初期時刻t
0と終了時刻(t
0+t)に応じた適切な値を採用するためには、次の式(9)を満足する値に設定されるのが好ましい。
図1に示すβ決定部14は、次の式(9)を数値的に解き、同式(9)を満足するβを決定する。
【数6】
パラメータβは小さいほど減衰が遅くなり、周波数領域でスペクトルが広くなる。本実施形態の例では、0.5〜0.7の値が好ましかった。なお、本実施形態では、β決定部14を設けてβを動的に決定しているが、パラメータβを予め決定しておき、メモリに設定しておいてもよい。
【0028】
具体的に、固有値問題解決部13は、式(6)に示す固有値問題を解くにあたり、長時間部分から順次緩和解析を行い、全領域の応力緩和を、個々の緩和解析結果の和で表現する(プロニー級数表示)。すなわち、
図4に示すように、緩和時間(全体)は、緩和解析(1)〜(3)の総和で再現される。固有値問題の解法は、非特許文献1と同じであるので、詳細は述べないが、簡単に説明すれば、
図5のフローに基づく。すなわち、ステップST11において、長時間部分(終端緩和付近)に初期時刻t
0及び終了時刻(t
0+t)を設定し、固有値問題解決部13が、ステップST12において式(6)の固有値問題を解き、最大N組の固有値λ
pと固有ベクトルg
pを得る。次のステップST13において、式(10)でA(t)を近似した結果をA’(t)とする。
【数7】
ここで、緩和弾性率は、次の式(13)で近似される。
【数8】
なお、グループiのg
pへの寄与は、式(13)の第i行の要素の和で定義できる。
【0029】
次のステップST14において、全時間を近似できたか否かを判定する。近似できていない場合には、ステップST15において、A(t)−A’(t)を新たなA(t)に更新し、ステップST11の処理に戻る。
【0030】
[高分子モデルの緩和弾性率を解析する方法]
上記装置1を用いて緩和弾性率を解析する方法について、
図2を用いて説明する。
まず、ステップST1において、グループ設定部10が、高分子モデルをN個のグループに分割する。次のステップST2において、緩和弾性率算出部11が、予め設定された解析条件を用いて平衡状態における高分子モデルの分子動力学計算を行い、グループ単位での応力の時間発展データを算出し、グループ単位の緩和弾性率を算出する。次のステップST3において、行列生成部12が、時刻tにおけるグループ毎の応力G
ji(t)を要素とするN×N個の行列A(t)を生成する。次のステップST4において、固有値問題解決部13が、終了時刻(t
0+t)における行列A(t
0+t)を、初期時刻t
0における行列A(t
0)と指数関数と固有ベクトルg
pを含むプロニー級数形式で表した固有値方程式を解き、行列A(t)を近似する。
【0031】
本発明の効果を示すために、従来(非特許文献1)との比較を示す。従来では、次の式(14)を用いて固有値問題を解いていた。
【数9】
従来の方法で得られた級数近似の結果をフーリエ変換すると、
図6Aに示すように、損失正接(tanδ)において、本来では有るはずのない凹凸部(図中にて円で囲んで示す)がプロット結果に表れてしまう。これは、応力緩和を指数関数の足し算で表現するにあたり、関数の緩和時間軸における分布に粗密があると周波数空間において強度がない部分ができやすく、その結果、スペクトルが無い部分ができてしまい、従来方法では凹凸部が生じてしまうと考えられる。従来のように一般的な指数関数を用いた緩和解析法では、上記凹凸部のような欠けが生じやすい、緩和弾性率のような実時間の値では問題がないが、フーリエ変換で得られる複素弾性率のような周波数空間では問題となる。
【0032】
一方、本発明のようにKWW関数を用いれば、
図6Bに示すように、上記凹凸部は出現しない。すなわち、本発明は、特許文献1に記載の利点(複素弾性率を適切に得ることができる点)と、非特許文献1に記載の利点(級数近似したときの各成分がどういう緩和時間に寄与しているかを表現できる点)とを両立する解析方法であることが分かる。
【0033】
以上のように、本実施形態の高分子モデルの緩和弾性率を解析する装置は、高分子モデルをN個のグループに分割するグループ設定部10と、予め設定された解析条件を用いて平衡状態における高分子モデルの分子動力学計算を行い、グループ単位での応力の時間発展データを算出し、グループ単位の緩和弾性率G
ji(t)[tは時間、i,jはグループに付けた番号を表す]を算出する緩和弾性率算出部11と、時刻tにおけるグループ毎の応力G
ji(t)を要素とするN×N個の行列A(t)を生成する行列生成部12と、終了時刻(t
0+t)における行列A(t
0+t)を、初期時刻t
0における行列A(t
0)と指数関数と固有ベクトルg
pを含むプロニー級数形式で表した固有値方程式を解き、行列A(t)を近似する固有値問題解決部13と、を有する。前記指数関数として、exp{−[(t
0+t)
β+t
0β]/τ
p}で表されるKWW(Kohlausch-Williams-Watts)関数を用いる。
【0034】
本実施形態の高分子モデルの緩和弾性率を解析する方法は、コンピュータが実行する方法であって、高分子モデルをN個のグループに分割するステップ(ST1)と、予め設定された解析条件を用いて平衡状態における高分子モデルの分子動力学計算を行い、グループ単位での応力の時間発展データを算出し、グループ単位の緩和弾性率G
ji(t)[tは時間、i,jはグループに付けた番号を表す]を算出するステップ(ST2)と、時刻tにおけるグループ毎の応力G
ji(t)を要素とするN×N個の行列A(t)を生成するステップ(ST3)と、終了時刻(t
0+t)における行列A(t
0+t)を、初期時刻t
0における行列A(t
0)と指数関数と固有ベクトルg
pを含むプロニー級数形式で表した固有値方程式を解き、行列A(t)を近似するステップ(ST4)と、を含む。前記指数関数として、exp{−[(t
0+t)
β+t
0β]/τ
p}で表されるKWW(Kohlausch-Williams-Watts)関数を用いる。
【0035】
このように、周波数空間において、従来の指数関数よりも広いKWW関数を用いるので、得られた近似式をフーリエ変換して複素弾性率を算出した場合に、従来の指数関数で生じていたデータの欠けが生じることなく、複素弾性率を適切に得ることが可能となる。したがって、複素弾性率を適切に得ることができる点と、級数近似したときの各成分がどういう緩和時間に寄与しているかを表現できる点とを両立できる。
【0036】
本実施形態では、具体的に、前記固有値方程式は、次の式(6)で表される。
A(t
0+t)g
p=exp{−[(t
0+t)
β+t
0β]/τ
p}A(t
0)g
p …(6)
【0037】
本実施形態では、前記KWW関数のパラメータβは、次の式(9)を満足する値に設定している。
G(t)=G(t
0)exp{−[(t+t
0)
β+t
0β]/2t
0} …(9)
これによれば、KWW関数のパラメータβを適切な値に設定可能となる。
【0038】
本実施形態に係るコンピュータプログラムは、上記方法を構成する各ステップをコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。言い換えると、上記方法を使用しているとも言える。
【0039】
以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。
【0040】
上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。
【符号の説明】
【0041】
10…グループ設定部
11…緩和弾性率算出部
12…行列生成部
13…固有値問題解決部
14…β決定部