特許 7382478 ポリマー材料の破壊靭性シミュレーション装置

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

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

(11)【特許番号】

(24)【登録日】2023-11-08

(45)【発行日】2023-11-16

(54)【発明の名称】ポリマー材料の破壊靭性シミュレーション装置

(51)【国際特許分類】

   G06F 30/20 20200101AFI20231109BHJP

   G01N 3/00 20060101ALI20231109BHJP

   G06F 113/26 20200101ALN20231109BHJP

【ＦＩ】

G06F30/20

G01N3/00 T

G06F113:26

【請求項の数】 11

(21)【出願番号】P 2022196272

(22)【出願日】2022-12-08

【審査請求日】2023-09-08

【早期審査対象出願】

(73)【特許権者】

【識別番号】000219602

【氏名又は名称】住友理工株式会社

(74)【代理人】

【識別番号】110000648

【氏名又は名称】弁理士法人あいち国際特許事務所

(72)【発明者】

【氏名】渡邊 直樹

(72)【発明者】

【氏名】矢島 高志

(72)【発明者】

【氏名】笠井 誠司

(72)【発明者】

【氏名】小山 敏幸

【審査官】堀井 啓明

(56)【参考文献】

【文献】特開２０２０－９５４９７（ＪＰ，Ａ）

【文献】特開２０１９－１０６１０８（ＪＰ，Ａ）

【文献】特開２０１６－８１１８３（ＪＰ，Ａ）

【文献】堀井 裕太 Yuta Horii Yuta Horii，フェーズフィールド微視的弾性論（反復法）を用いた破壊靱性の評価 Calculation of fracture toughness based on the iteration method of phase-field microelasticity theory，日本鉄鋼協会講演論文集「材料とプロセス」 第１６９回 ，第28巻，日本，脇本 眞也 一般社団法人日本鉄鋼協会，2015年03月01日，p.311

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

Ｇ０６Ｆ ３０／００－３０／２８

Ｇ０１Ｎ ３／００

Ｇ０６Ｆ １１３／２６

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

【請求項1】

マトリックスポリマーを含むポリマー材料の破壊靭性シミュレーション装置であって、
モデル定義領域における位置要素ごとに、前記マトリックスポリマーに対応するマトリックスポリマーモデルまたはクラックモデルにより構成されたシミュレーションモデルであって、前記マトリックスポリマーモデルの弾性率と前記クラックモデルの弾性率とが異なる弾性率となるように設定され、かつ、少なくとも前記マトリックスポリマーモデルの弾性率として応力および歪に応じた弾性率が設定された前記シミュレーションモデルを記憶する記憶部と、
前記シミュレーションモデルにおいて、複数の前記クラックモデルを連続させたクラックモデル群の位置および大きさを設定するクラックモデル群位置設定部と、
前記クラックモデル群位置設定部により前記クラックモデル群が設定された前記シミュレーションモデルにおいて、初期条件として前記シミュレーションモデル全体に対して所定応力または所定歪を付与し、前記シミュレーションモデルが力学的平衡状態となるように演算処理を行う演算処理部と、
前記シミュレーションモデルが前記力学的平衡状態となるときにおいて、弾性歪エネルギーに、前記初期条件の付与に伴って減少したポテンシャルエネルギーを加算することにより得られるギブスエネルギーを算出するギブスエネルギー算出部と、
前記クラックモデル群位置設定部により設定された前記クラックモデル群の大きさを表すクラックサイズと、前記ギブスエネルギー算出部により算出された前記ギブスエネルギーとの関係に基づいて、前記クラックサイズが増大した場合に前記ギブスエネルギーが減少する程度を表すエネルギー解放率を算出するエネルギー解放率算出部と、
を備える、ポリマー材料の破壊靭性シミュレーション装置。

【請求項2】

前記演算処理部は、
現在の前記シミュレーションモデルの前記位置要素のそれぞれにおける弾性率を現在の弾性率として、前記初期条件が付与された前記シミュレーションモデルが前記力学的平衡状態となるように演算処理を行うことにより、前記力学的平衡状態となるときの前記シミュレーションモデルの位置要素ごとの応力および歪を算出し、
前記マトリックスポリマーモデルの前記位置要素ごとの応力および歪に基づいて、当該マトリックスポリマーモデルの位置要素における変更後の弾性率を算出し、
前記マトリックスポリマーモデルの前記位置要素における弾性率を前記変更後の弾性率に変更して、再び前記シミュレーションモデルが前記力学的平衡状態となるように演算処理を行うことにより、前記力学的平衡状態となるときの前記シミュレーションモデルの前記位置要素ごとの応力および歪を算出し、
予め設定した安定状態に近づくまで、前記変更後の弾性率の算出、および、前記シミュレーションモデルの前記位置要素ごとの応力および歪の算出を繰り返し行うことにより、前記シミュレーションモデルが前記安定状態としての前記力学的平衡状態となるときの前記シミュレーションモデルの前記位置要素ごとの応力、歪および弾性率を算出する、請求項１に記載のポリマー材料の破壊靭性シミュレーション装置。

【請求項3】

前記安定状態は、前記弾性歪エネルギーのマイナス２倍が前記ポテンシャルエネルギーに一致する状態から所定値の範囲内の状態である、請求項２に記載のポリマー材料の破壊靭性シミュレーション装置。

【請求項4】

前記安定状態は、前記シミュレーションモデルの前記位置要素のすべての歪の変化が所定値以下となる状態である、請求項２に記載のポリマー材料の破壊靭性シミュレーション装置。

【請求項5】

前記マトリックスポリマーモデルは、応力－歪の曲線上の点と起点とを結ぶ直線の傾きを表すセカント弾性率に基づいて設定されている、請求項１または２に記載のポリマー材料の破壊靭性シミュレーション装置。

【請求項6】

前記クラックモデル群位置設定部は、操作者の入力により前記クラックモデル群の位置および大きさを設定するように構成されている、請求項１または２に記載のポリマー材料の破壊靭性シミュレーション装置。

【請求項7】

前記クラックモデル群位置設定部は、
演算処理により前記シミュレーションモデルの前記位置要素における応力に基づいて拡大する前記クラックモデル群の位置および大きさを決定し、
拡大後の前記クラックモデル群の位置および大きさを、次の前記シミュレーションモデルを構成する前記クラックモデル群の位置および大きさとして設定するように構成されている、請求項１または２に記載のポリマー材料の破壊靭性シミュレーション装置。

【請求項8】

前記ポリマー材料は、前記マトリックスポリマーにフィラーが分散されたポリマー複合材料であり、
前記シミュレーションモデルは、前記モデル定義領域の前記位置要素ごとに、前記マトリックスポリマーモデル、前記フィラーに対応するフィラーモデル、または前記クラックモデルにより構成されたシミュレーションモデルである、請求項１または２に記載のポリマー材料の破壊靭性シミュレーション装置。

【請求項9】

前記マトリックスポリマーモデルは、ポリマー母相モデル、前記ポリマー母相モデルと前記フィラーモデルとの間に位置するポリマー界面相モデル、により構成され、
前記ポリマー界面相モデルは、前記位置要素ごとに、前記ポリマー界面相モデルの弾性率を、所定の最小弾性率と所定の最大弾性率との間の割合を設定するためのフェーズフィールド変数を用いて定義されている、請求項８に記載のポリマー材料の破壊靭性シミュレーション装置。

【請求項10】

前記ポリマー界面相モデルの弾性率は、前記フェーズフィールド変数が前記ポリマー母相モデル側から前記フィラーモデル側にかけて連続的に変化する値となるように定義されている、請求項９に記載のポリマー材料の破壊靭性シミュレーション装置。

【請求項11】

前記弾性歪エネルギーは、前記シミュレーションモデルの前記位置要素ごとの応力および歪に基づいて算出され、
前記ポテンシャルエネルギーは、前記シミュレーションモデル全体における歪に基づいて算出される、請求項１または２に記載のポリマー材料の破壊靭性シミュレーション装置。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、ポリマー材料の破壊靭性シミュレーション装置に関する。

【背景技術】

【0002】

特許文献１には、マトリックスゴムの破壊現象を再現するためのシミュレーションについて記載されている。このシミュレーションに用いられるゴム材料モデルは、マトリックスゴムをモデル化した複数個のゴム要素と、ゴム要素間を剥離可能に連結する第１粘着要素とを用いて定義されている。また、ゴム材料モデルは、ゴム要素および第１粘着要素に加えて、フィラーをモデル化したフィラー要素と、ゴム要素とフィラー要素との間を剥離可能に連結する第２粘着要素とを用いて定義されている。

【0003】

特許文献２には、物性値が異なる複数の高分子成分を含むゴム材料、すなわちポリマーブレンド材料の破壊現象をシミュレーションすることについて記載されている。ポリマーブレンド材料モデルは、第１高分子成分を第１粗視化粒子で表現した第１高分子鎖モデル、および、第２高分子成分を第２粗視化粒子で表現した第２高分子鎖モデルを用いて定義され、さらに、第１高分子鎖モデルと第２高分子鎖モデルとの間に界面が形成されるように定義されている。そして、シミュレーションにより、ポリマーブレンド材料の変形を計算して、第１粗視化粒子間の相互作用を定義するための第１相互作用ポテンシャル、および、第２粗視化粒子間の相互作用を定義するための第２相互作用ポテンシャルに基づいて、界面での破壊現象を計算している。このようにして、第１高分子鎖モデルおよび第２高分子鎖モデルが存在しない空孔（ボイド）が生じる現象が計算される。

【先行技術文献】

【特許文献】

【0004】

【文献】特開２０１８－２５４８９号公報

【文献】特開２０２０－８６８３５号公報

【発明の概要】

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

【0005】

材料の破壊特性として、破壊靭性特性が知られている。破壊靭性特性の指標の１つとして、エネルギー解放率がある。エネルギー解放率は、クラックサイズの進展に際して失われるエネルギーを、クラックサイズで除したものと定義される。つまり、エネルギー解放率は、クラックサイズが増大した場合にギブスエネルギーが減少する程度を表す。

【0006】

シミュレーションモデルにおいてクラックサイズを設定した状態で、ギブスエネルギーを取得することで、エネルギー解放率を得ることができる。しかし、破壊靭性特性を評価する対象を、マトリックスポリマーを含むポリマー材料とする場合には、シミュレーションが容易ではない。つまり、ポリマー材料のエネルギー解放率を高精度に得ることができない。

【0007】

エネルギー解放率が小さいほど、破壊靭性が高い。破壊靭性が高いほど、クラックが進展しにくい。ポリマー材料のシミュレーションにおいて、エネルギー解放率、すなわちクラックの進展度合いを把握することができれば、現存するポリマー材料の破壊靭性を把握することができるようになる。さらには、将来において、クラックが進展しにくいポリマー材料の指針を提供することもできるようになる。

【0008】

本発明は、かかる背景に鑑みてなされたものであり、ポリマー材料のエネルギー解放率を高精度に算出することにより、ポリマー材料の破壊靭性特性を評価することができるポリマー材料の破壊靭性シミュレーション装置を提供しようとするものである。

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

【0009】

本発明の一態様は、マトリックスポリマーを含むポリマー材料の破壊靭性シミュレーション装置であって、
モデル定義領域における位置要素ごとに、前記マトリックスポリマーに対応するマトリックスポリマーモデルまたはクラックモデルにより構成されたシミュレーションモデルであって、前記マトリックスポリマーモデルの弾性率と前記クラックモデルの弾性率とが異なる弾性率となるように設定され、かつ、少なくとも前記マトリックスポリマーモデルの弾性率として応力および歪に応じた弾性率が設定された前記シミュレーションモデルを記憶する記憶部と、
前記シミュレーションモデルにおいて、複数の前記クラックモデルを連続させたクラックモデル群の位置および大きさを設定するクラックモデル群位置設定部と、
前記クラックモデル群位置設定部により前記クラックモデル群が設定された前記シミュレーションモデルにおいて、初期条件として前記シミュレーションモデル全体に対して所定応力または所定歪を付与し、前記シミュレーションモデルが力学的平衡状態となるように演算処理を行う演算処理部と、
前記シミュレーションモデルが前記力学的平衡状態となるときにおいて、弾性歪エネルギーに、前記初期条件の付与に伴って減少したポテンシャルエネルギーを加算することにより得られるギブスエネルギーを算出するギブスエネルギー算出部と、
前記クラックモデル群位置設定部により設定された前記クラックモデル群の大きさを表すクラックサイズと、前記ギブスエネルギー算出部により算出された前記ギブスエネルギーとの関係に基づいて、前記クラックサイズが増大した場合に前記ギブスエネルギーが減少する程度を表すエネルギー解放率を算出するエネルギー解放率算出部と、
を備える、ポリマー材料の破壊靭性シミュレーション装置にある。

【発明の効果】

【0010】

上記態様によれば、シミュレーションモデルは、モデル定義領域における位置要素ごとに、マトリックスポリマーに対応するマトリックスポリマーモデルまたはクラックモデルにより構成されている。そして、マトリックスポリマーモデルおよびクラックモデルは、それぞれ、応力および歪に応じた弾性率が設定されている。つまり、シミュレーションモデルは、位置要素ごとに、応力および歪に応じた弾性率が設定されている。

【0011】

ここで、マトリックスポリマーは、応力に対する歪の関係を表す弾性率が線形の関係を有していない。そこで、上記のように、シミュレーションモデルに応力または歪を付与した場合に、マトリックスポリマーモデルの弾性率を、シミュレーションモデルにおける位置要素ごとの応力や歪に応じて異なる値に設定できる。つまり、マトリックスポリマーは、応力および歪に応じて弾性率が変化することから、シミュレーションモデルは、このようなマトリックスポリマーの弾性率の変化現象を適切に捉えることができるモデルとなる。

【0012】

このように位置要素ごとに応力および歪に応じた弾性率が設定されたシミュレーションモデルを用いて、演算処理部が演算処理を行っている。具体的には、演算処理部は、初期条件としてシミュレーションモデル全体に対して所定応力または所定歪を付与して、シミュレーションモデルが力学的平衡状態となるように演算処理を行う。

【0013】

そして、ギブスエネルギー算出部が、シミュレーションモデルが力学的平衡状態となるときにおいて、弾性歪エネルギーに、初期条件の付与に伴って減少したポテンシャルエネルギーを加算することにより得られるギブスエネルギーを算出する。従って、算出されるギブスエネルギーは、高精度な値となる。

【0014】

エネルギー解放率算出部が、クラックサイズとギブスエネルギーとの関係に基づいて、エネルギー解放率を算出している。このようにして算出されたエネルギー解放率は、高精度となる。

【0015】

以上より、ポリマー材料のエネルギー解放率を高精度に算出することにより、ポリマー材料の破壊靭性特性を評価することができるポリマー材料の破壊靭性シミュレーション装置を提供することができる。

【図面の簡単な説明】

【0016】

【図1】実施形態１における破壊靭性シミュレーション装置の構成を示す図。

【図2】シミュレーション対象Ｐとしての、フィラーおよびマトリックスポリマーを含むポリマー材料の断面画像を示す図。

【図3】図１に示す破壊靭性シミュレーション装置を構成するベースモデル作成処理部により作成されたベースモデルＭ１を示す図であって、フィラーモデルＭａ、ポリマー母相モデルＭｂ、ポリマー界面相モデルＭｃを定義したモデルを示す図。

【図4】ベースモデルＭ１の拡大模式図（左図）、および、当該拡大図における位置要素ごとのモデル種を示す図（右図）。

【図5】ポリマー母相モデルＭｂおよびポリマー界面相モデルＭｃを含むマトリックスポリマーモデルのセカント弾性率を示す図。

【図6】ベースモデルＭ１にクラックモデルＭｄを追加したシミュレーションモデルＭ２を示す図。

【図7】シミュレーションモデルＭ２の拡大模式図（左図）、および、当該拡大模式図における位置要素ごとのモデル種を示す図。

【図8】図１に示す破壊靭性シミュレーション装置を構成するベースモデル作成処理部による処理を示すフローチャート。

【図9】図１に示す破壊靭性シミュレーション装置を構成するクラックモデル群位置設定部による処理を示すフローチャート。

【図10】図１に示す破壊靭性シミュレーション装置を構成する演算処理部による処理を示すフローチャート。

【図11】演算処理部による処理において、シミュレーションモデルＭ２に対して所定応力または所定歪を付与した状態を示す図。

【図12】図１に示す破壊靭性シミュレーション装置を構成するエネルギー解放率算出部において算出されるエネルギー解放率を説明する図。

【図13】実施形態２における破壊靭性シミュレーション装置の構成を示す図。

【図14】図１３に示すクラックモデル群位置設定部により、新たなクラックモデル群を自動的に設定する過程を示す図。

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

【0017】

（実施形態１）
１－１．破壊靭性シミュレーション装置１の概要
破壊靭性シミュレーション装置１は、ポリマー材料をシミュレーション対象とする。シミュレーション対象であるポリマー材料は、少なくとも、ゴムなどのマトリックスポリマーを含む。ポリマー材料は、例えば、タイヤ、ダンパ、靴、スポーツ用品など、種々の分野で使用される。ポリマー材料に含まれるマトリックスポリマーは、１種としても良いし、複数種類としても良い。

【0018】

シミュレーション対象のポリマー材料は、マトリックスポリマーに、目的に応じたフィラーが分散されたポリマー複合材料としても良い。また、ポリマー複合材料は、凝集したフィラーが不均一に分布している。フィラーには、補強を目的としたもの、導電性を目的としたものなどがある。なお、シミュレーション対象のポリマー材料は、フィラーを含まないものとしても良い。

【0019】

例えば、ポリマー複合材料の例としてのタイヤは、マトリックスポリマーに、カーボンブラックなどの充填剤としてのナノフィラーが分散されることにより、補強された加硫ゴムが用いられる。

【0020】

そして、シミュレーション対象のポリマー材料がクラックを含む場合において、ポリマー材料に荷重がかけられることに伴って、クラックが進展することがある。ここで、クラックは、ポリマー材料の使用によって発生する空乏領域のみならず、ポリマー材料の使用初期（製品完成時）において予め存在する空乏領域を含む。

【0021】

そこで、破壊靭性シミュレーション装置１は、ポリマー材料の破壊靭性を評価するために、クラックサイズの進展に際して失われるエネルギー量を表すエネルギー解放率を得ることを目的とする。以下に詳細に説明するが、破壊靭性シミュレーション装置１が算出するエネルギー解放率は、クラックサイズが増大した場合にギブスエネルギーが減少する程度を表す。

【0022】

１－２．破壊靭性シミュレーション装置１の構成の概要
破壊靭性シミュレーション装置１の構成の概要について、図１を参照して説明する。本実施形態においては、シミュレーション対象Ｐのポリマー材料は、マトリックスポリマーにフィラーが分散されたポリマー複合材料である場合を例にあげる。

【0023】

破壊靭性シミュレーション装置１は、ベースモデル作成処理部２、クラックモデル群位置設定部３、記憶部４、演算処理部５、ギブスエネルギー算出部６、エネルギー解放率算出部７を備える。各部の概要について説明する。なお、各部の詳細は後述する。

【0024】

ベースモデル作成処理部２は、シミュレーションモデルＭ２のベースとなるベースモデルＭ１を作成する。ベースモデル作成処理部２は、モデル定義領域における位置要素ごとに、マトリックスポリマーに対応するマトリックスポリマーモデルＭｂ，Ｍｃ、または、フィラーに対応するフィラーモデルＭａにより構成されたベースモデルＭ１を作成する。ベースモデルＭ１は、マトリックスポリマーモデルＭｂ，Ｍｃの弾性率として、位置要素ごとの応力および歪に応じた弾性率が設定されている。また、ベースモデルＭ１は、フィラーモデルＭａの弾性率として、予め設定された弾性率が設定されている。

【0025】

クラックモデル群位置設定部３は、ベースモデルＭ１に、クラックに対応するクラックモデルＭｄを付加することにより、シミュレーションモデルＭ２を作成する。特に、クラックモデル群位置設定部３は、シミュレーションモデルＭ２において、複数のクラックモデルＭｄを連続させたクラックモデル群の位置および大きさを設定する。クラックモデルＭｄの弾性率は、予め設定された弾性率が設定されている。

【0026】

記憶部４は、ベースモデル作成処理部２およびクラックモデル群位置設定部３により作成されたシミュレーションモデルＭ２を記憶する。つまり、記憶部４は、モデル定義領域における位置要素ごとに、フィラーモデルＭａ、マトリックスポリマーモデルＭｂ，ＭｃまたはクラックモデルＭｄにより構成されたシミュレーションモデルＭ２を記憶する。シミュレーションモデルＭ２は、フィラーモデルＭａの弾性率、マトリックスポリマーモデルＭｂ，Ｍｃの弾性率、および、クラックモデルＭｄの弾性率が異なる弾性率となるように設定されている。さらに、マトリックスポリマーモデルＭｂ，Ｍｃの弾性率は、位置要素ごとに応力および歪に応じた弾性率が設定されている。

【0027】

演算処理部５は、シミュレーションモデルＭ２において、初期条件としてシミュレーションモデルＭ２全体に対して所定応力または所定歪を付与し、シミュレーションモデルＭ２が力学的平衡状態となるように演算処理を行う。本実施形態においては、演算処理部５は、シミュレーションモデルＭ２が安定状態に近づくまで、シミュレーションモデルＭ２の弾性率の変更と、力学的平衡状態となるような演算処理とを繰り返し行う。

【0028】

ギブスエネルギー算出部６は、シミュレーションモデルＭ２が力学的平衡状態となるときにおいて、弾性歪エネルギーに、初期条件の付与に伴って減少したポテンシャルエネルギーを加算することにより得られるギブスエネルギーを算出する。

【0029】

エネルギー解放率算出部７は、クラックモデル群の大きさを表すクラックサイズと、ギブスエネルギーとの関係に基づいて、クラックサイズが増大した場合にギブスエネルギーが減少する程度を表すエネルギー解放率を算出する。上述したように、エネルギー解放率が、ポリマー材料の破壊靭性を表す指標となる。エネルギー解放率が小さいほど、破壊靭性が高いことを表す。つまり、破壊靭性が高いほど、クラックが進展しにくいことになる。

【0030】

１－３．ベースモデルＭ１
ベースモデル作成処理部２により作成されるベースモデルＭ１について、図２～図５を参照して説明する。ベースモデル作成処理部２は、図２に示すようなシミュレーション対象Ｐのポリマー材料の断面画像に対応するベースモデルＭ１、すなわち図３に示すようなベースモデルＭ１を作成する。

【0031】

図２に示すように、シミュレーション対象Ｐのポリマー材料として、マトリックスポリマーＰｂに、フィラーＰａが不均一に分散されたものを例に挙げる。図２は、実際のポリマー材料の断面画像である。このように、実際のポリマー材料は、フィラーＰａが均一に分散されることはなく、不均一に分散されている。シミュレーション対象Ｐのポリマー材料の断面画像は、フィラーＰａおよびマトリックスポリマーＰｂにより表されている。

【0032】

フィラーＰａは、例えば、カーボンブラックなどである。フィラーＰａは、カーボンブラックの他に、種々の機能を発揮するための充填剤を適用することができる。また、フィラーＰａは、ナノフィラーが凝集することにより構成されている。従って、凝集されたフィラーＰａの塊は、形も大きさも異なる。マトリックスポリマーＰｂは、ゴムやエラストマーであって、例えば、スチレンブタジエンゴム（ＳＢＲ）などを適用する。

【0033】

ベースモデルＭ１は、図２に示すようなシミュレーション対象Ｐのポリマー材料の断面画像に対応するモデルである。本実施形態においては、ベースモデルＭ１は、例えば、メゾスケールのモデル定義領域、すなわち、１０ｎｍ～１００μｍを一辺とするモデル定義領域を有する。ベースモデルＭ１をメゾスケールの定義領域を有するモデルとすることにより、上述した目的である破壊靭性特性を評価することができる。

【0034】

ベースモデルＭ１は、モデル定義領域における位置要素ごとに、フィラーＰａに対応するフィラーモデルＭａ、マトリックスポリマーＰｂに対応するマトリックスポリマーモデルＭｂ，Ｍｃにより構成されている。本実施形態においては、マトリックスポリマーＰｂに対応するマトリックスポリマーモデルＭｂ，Ｍｃを、ポリマー母相モデルＭｂとポリマー界面相モデルＭｃとに分けて設定する。従って、ベースモデルＭ１は、モデル定義領域における位置要素ごとに、フィラーモデルＭａ、ポリマー母相モデルＭｂ、および、ポリマー界面相モデルＭｃのいずれか１つのモデルにより構成されている。

【0035】

図４の左側に、フィラーモデルＭａ、ポリマー母相モデルＭｂ、ポリマー界面相モデルＭｃの拡大模式図を示し、図４の右側に、拡大模式図において各位置要素のモデル種を示す。図４の右側の図において、記号Ａ，Ｂ，Ｃは、順に、フィラーモデルＭａ、ポリマー母相モデルＭｂ、ポリマー界面相モデルＭｃに相当する。

【0036】

フィラーモデルＭａは、図３に示すように、ナノフィラーが凝集することにより構成されている。従って、個々のフィラーモデルＭａは、形も大きさも異なる。

【0037】

ポリマー母相モデルＭｂは、マトリックスポリマー単層として存在する領域に位置する。従って、ポリマー母相モデルＭｂは、マトリックスポリマー単層の性質を有する。ポリマー母相モデルＭｂは、フィラーモデルＭａから僅かに離れた位置に位置する。また、ポリマー母相モデルＭｂは、Ｏｇｄｅｎモデルを適用する。

【0038】

ポリマー界面相モデルＭｃは、ポリマー母相モデルＭｂとフィラーモデルＭａとの間に位置する。特に、ポリマー界面相モデルＭｃは、フィラーモデルＭａの表面に位置する。そこで、ポリマー界面相モデルＭｃは、フィラーモデルＭａの表面からの厚みにより表される。つまり、ポリマー界面相モデルＭｃの位置は、フィラーモデルＭａの表面から、予め設定されたポリマー界面相モデルＭｃの厚み情報に基づいて設定される。ここで、ポリマー界面相モデルＭｃの厚みは、予め設定することができ、マトリックスポリマーＰｂの材質やフィラーＰａの材質などによって適宜設定することができる。

【0039】

さらに、各モデルＭａ，Ｍｂ，Ｍｃは、それぞれに応じた初期弾性率が設定されている。フィラーモデルＭａの弾性率と、マトリックスポリマーモデルＭｂ，Ｍｃの弾性率とは、異なる弾性率となるように設定されている。

【0040】

本実施形態においては、フィラーモデルＭａの弾性率は、一定値となるように設定されている。例えば、フィラーモデルＭａの弾性率は、２０ＭＰａなどとする。ただし、フィラーモデルＭａの弾性率は、この値に限られるものではない。

【0041】

ポリマー母相モデルＭｂの弾性率Ｃｍ（ｒ）は、位置ｒに応じて異なる値に設定される。具体的には、式（１）に示すように、ポリマー母相モデルＭｂの弾性率Ｃｍ（ｒ）は、位置要素ごとに、応力および歪に応じた弾性率が設定されている。式（１）において、ｒは、モデル定義領域における位置を表す。

【0042】

【数1】

【0043】

ポリマー母相モデルＭｂの弾性率Ｃｍ（ｒ）は、図５に示すように、応力－歪の曲線上の点と起点とを結ぶ直線の傾きを表すセカント弾性率Ｃ^０（ｒ）となるように設定される。ベースモデルＭ１においては、ポリマー母相モデルＭｂの初期弾性率Ｃ^０（ｒ）は、図５に示す応力－歪の曲線における起点からの立ち上がりの傾きに一致する。

【0044】

ポリマー界面相モデルＭｃの弾性率Ｃinter（ｒ）は、ポリマー母相モデルＭｂと同様に位置ｒに応じて異なる値に設定しても良いし、一定値となるように設定しても良い。

【0045】

ポリマー界面相モデルＭｃの弾性率Ｃinter（ｒ）は、位置ｒに応じて異なる値に設定する場合には、式（２）に示すように、位置ｒの関数であるフェーズフィールド変数φ（ｒ）を用いて定義される。

【0046】

【数2】

【0047】

フェーズフィールド変数φ（ｒ）は、ベースモデルＭ１の位置要素ごとの弾性率Ｃ（ｒ）を、所定の最小弾性率Ｃ^０（ｒ）と所定の最大弾性率Ｃ^１との間の割合を設定するための変数である。所定の最小弾性率Ｃ^０（ｒ）は、ポリマー母相モデルＭｂの弾性率Ｃｍ（ｒ）に相当する弾性率Ｃ^０（ｒ）である。つまり、所定の最小弾性率Ｃ^０（ｒ）は、図５に示すセカント弾性率から得られる。所定の最大弾性率Ｃ^１は、マトリックスポリマーＰｂとして取り得る最大弾性率としても良いし、フィラーモデルＭａの弾性率としても良い。

【0048】

フェーズフィールド変数φ（ｒ）は、０～１の値をとる。つまり、フェーズフィールド変数φ（ｒ）が０の場合には、ポリマー界面相モデルＭｂの弾性率Ｃｍ（ｒ）は、所定の最小弾性率Ｃ^０（ｒ）となる。フェーズフィールド変数φ（ｒ）が１の場合には、ポリマー界面相モデルＭｂの弾性率Ｃｍ（ｒ）は、所定の最大弾性率Ｃ^１となり、

【0049】

フェーズフィールド変数φ（ｒ）は、位置ｒに応じて異なる値に設定することができる。例えば、フェーズフィールド変数φ（ｒ）は、フィラーモデルＭａからの距離に応じて異なる値に設定することができる。具体的には、フェーズフィールド変数φ（ｒ）は、フィラーモデルＭａに近いほど１に近い値とし、フィラーモデルＭａから遠いほど０に近い値とする。

【0050】

特に、フェーズフィールド変数φ（ｒ）を、ポリマー母相モデルＭｂ側からフィラーモデルＭａ側にかけて連続的に変化する値となるように定義すると良い。この場合、ポリマー界面相モデルＭｃの弾性率Ｃinter（ｒ）は、ポリマー母相モデルＭｂからフィラーモデルＭａ側にかけて連続的に変化する。

【0051】

また、フェーズフィールド変数φ（ｒ）は、位置ｒによらず、一定値とすることもできる。例えば、フェーズフィールド変数φ（ｒ）を、０．５などに設定する。

【0052】

また、ポリマー界面相モデルＭｃの弾性率Ｃinter（ｒ）は、フェーズフィールド変数φ（ｒ）を用いた関数とするのではなく、一定値とすることもできる。例えば、１０ＭＰａなどとする。ポリマー界面相モデルＭｃの弾性率は、例えば、フィラーモデルＭａの弾性率の中間程度の値などに設定する。

【0053】

１－４．シミュレーションモデルＭ２
シミュレーションモデルＭ２について、図６および図７を参照して説明する。シミュレーションモデルＭ２は、上述したように、クラックモデル群位置設定部３により作成され、記憶部４に記憶される。

【0054】

シミュレーションモデルＭ２は、ベースモデルＭ１に、クラックに対応するクラックモデルＭｄを付加することにより作成される。図６および図７に示すように、シミュレーションモデルＭ２においては、複数のクラックモデルＭｄを連続させたクラックモデル群が設定されている。クラックモデル群は、複数のクラックモデルＭｄを連続させて構成することにより、位置および大きさが設定される。

【0055】

クラックモデル群は、例えば、図６に示すように、縦横比の大きな細長い楕円形状、一方向に長さを有する線形状（例えば、直線形状）などのように設定されている。ここで、クラックモデル群は、長手方向がシミュレーションモデルＭ２に付与される所定応力または所定歪の方向に対して直交する方向となるように設定されている。図６および図７においては、上下方向に所定応力または所定歪が付与されるため、クラックモデル群の長手方向は、左右方向に延在するように設定されている。そして、クラックモデルＭｄは、弾性率が０に設定されている。

【0056】

従って、シミュレーションモデルＭ２は、モデル定義領域における位置要素ごとに、フィラーモデルＭａ、ポリマー母相モデルＭｂ、ポリマー界面相モデルＭｃ、および、クラックモデルＭｄのいずれか１つのモデルにより構成されている。さらに、シミュレーションモデルＭ２は、フィラーモデルＭａの弾性率、マトリックスポリマーモデルＭｂ，Ｍｃの弾性率、クラックモデルＭｄの弾性率がそれぞれ異なる弾性率となるように設定されている。特に、少なくともポリマー母相モデルＭｂの弾性率は、応力および歪に応じた弾性率が設定されている。上述したように、ポリマー界面相モデルＭｃの弾性率は、一定値としても良いし、位置に応じて異なる弾性率となるように設定しても良い。

【0057】

１－５．ベースモデル作成処理部２の処理
ベースモデル作成処理部２によるベースモデルＭ１の作成処理について、図８を参照して説明する。ベースモデル作成処理部２は、ベースモデルＭ１の作成種を選択する（Ｓ１１）。ベースモデルＭ１の作成種としては、画像利用の方法、画像参照の方法、自由作成の方法から選択可能とされている。

【0058】

Ｓ１１において画像利用の方法が選択された場合には（Ｓ１２：画像利用）、まず、図２に示すようなシミュレーション対象Ｐのポリマー材料の断面画像を取得する（Ｓ２１）。断面画像は、フィラーＰａおよびマトリックスポリマーＰｂにより表されている。シミュレーション対象Ｐのポリマー材料の断面画像は、モデル定義領域に対応するサイズを有する。

【0059】

続いて、断面画像に基づいてベースモデルＭ１を作成する。まず、図３に示すように、断面画像に基づいて、フィラーＰａに対応するフィラーモデルＭａの位置を設定する（Ｓ２２）。例えば、図２に示す断面画像を２値化することにより、フィラーモデルＭａの位置を特定することができる。

【0060】

続いて、マトリックスポリマーＰｂに対応するマトリックスポリマーモデルを、ポリマー母相モデルＭｂとポリマー界面相モデルＭｃとに分けて設定する。まずは、ポリマー界面相モデルＭｃの位置を設定する（Ｓ２３）。ポリマー界面相モデルＭｃの位置は、予め設定されたポリマー界面相モデルＭｃの厚み情報に基づいて、フィラーモデルＭａの表面から厚み情報に相当する厚みの領域に設定される。

【0061】

続いて、マトリックスポリマーモデルの残りの領域を、ポリマー母相モデルＭｂの位置として設定する（Ｓ２４）。このようにして、ベースモデルＭ１が作成される。つまり、ベースモデルＭ１は、図３に示すように、モデル定義領域の位置要素ごとに、フィラーモデルＭａ、ポリマー母相モデルＭｂ、および、ポリマー界面相モデルＭｃのいずれか１つのモデルにより構成される。

【0062】

また、Ｓ１１において画像参照の方法が選択された場合には（Ｓ１２：画像参照）、まず、図２に示すようなシミュレーションの対象Ｐのポリマー材料の断面画像を取得する（Ｓ３１）。断面画像は、フィラーＰａおよびマトリックスポリマーＰｂにより表されている。

【0063】

続いて、断面画像に基づいてベースモデルＭ１を作成する。まず、図３に示すように、断面画像に基づいて、フィラーＰａに対応するフィラーモデルＭａの位置を設定する（Ｓ３２）。例えば、図２に示す断面画像を２値化することにより、フィラーモデルＭａの位置を特定することができる。続いて、操作者の入力により、フィラーモデルＭａの位置を修正する（Ｓ３３）。例えば、フィラーモデルＭａが存在しない位置要素をフィラーモデルＭａとして追加したり、フィラーモデルＭａとしての位置要素をフィラーモデルＭａでない状態にしたりする。

【0064】

続いて、マトリックスポリマーＰｂに対応するマトリックスポリマーモデルを、ポリマー母相モデルＭｂとポリマー界面相モデルＭｃとに分けて設定する。まず、ポリマー界面相モデルＭｃの位置を設定する（Ｓ３４）。ポリマー界面相モデルＭｃの位置は、予め設定されたポリマー界面相モデルＭｃの厚み情報に基づいて、フィラーモデルＭａの表面から厚み情報に相当する厚みの領域に設定される。

【0065】

続いて、マトリックスポリマーモデルの残りの領域を、ポリマー母相モデルＭｂの位置として設定する（Ｓ３５）。このようにして、ベースモデルＭ１が作成される。つまり、ベースモデルＭ１は、図３に示すように、モデル定義領域の位置要素ごとに、フィラーモデルＭａ、ポリマー母相モデルＭｂ、および、ポリマー界面相モデルＭｃのいずれか１つのモデルにより構成される。

【0066】

また、Ｓ１１において自由作成の方法が選択された場合には（Ｓ１２：自由作成）、まず、操作者の入力により、モデル定義領域において、フィラーモデルＭａの位置を設定する（Ｓ４１）。続いて、マトリックスポリマーＰｂに対応するマトリックスポリマーモデルを、ポリマー母相モデルＭｂとポリマー界面相モデルＭｃとに分けて設定する。まず、ポリマー界面相モデルＭｃの位置を設定する（Ｓ４２）。ポリマー界面相モデルＭｃの位置は、予め設定されたポリマー界面相モデルＭｃの厚み情報に基づいて、フィラーモデルＭａの表面から厚み情報に相当する厚みの領域に設定される。続いて、マトリックスポリマーモデルの残りの領域を、ポリマー母相モデルＭｂの位置として設定する（Ｓ４３）。

【0067】

このようにして、ベースモデルＭ１が作成される。つまり、ベースモデルＭ１は、図４に示すように、モデル定義領域の位置要素ごとに、フィラーモデルＭａ、ポリマー母相モデルＭｂ、および、ポリマー界面相モデルＭｃのいずれか１つのモデルにより構成される。

【0068】

１－６．クラックモデル群位置設定部３の処理
クラックモデル群位置設定部３による処理について図９を参照して説明する。クラックモデル群位置設定部３は、ベースモデル作成処理部２により作成されたベースモデルＭ１を取得する（Ｓ５１）。

【0069】

続いて、操作者の入力により、ベースモデルＭ１に対して、クラックモデルＭｄを付加する。このとき、複数のクラックモデルＭｄを連続させたクラックモデル群を設定する。このようにして、クラックモデル群の位置および大きさが設定される（Ｓ５２）。つまり、クラックモデル群位置設定部３は、操作者の入力によりクラックモデル群の位置および大きさを設定するように構成されている。

【0070】

このようにして、ベースモデルＭ１に対してクラックモデル群が付加設定されることにより、図６に示すようなシミュレーションモデルＭ２が作成される。作成されたシミュレーションモデルＭ２は、記憶部４に記憶される。

【0071】

１－７．演算処理部５の処理
演算処理部５による処理について図１０および図１１を参照して説明する。演算処理部５は、シミュレーションモデルＭ２を用いて、所定応力または所定歪を付与した場合に、力学的平衡状態となるように演算処理を行う。そして、演算処理部５は、力学的平衡状態になるときにおけるシミュレーションモデルＭ２の位置要素ごとの応力および歪を算出する。演算処理部５による処理は、具体的には、以下のように行われる。

【0072】

図１０に示すように、演算処理部５は、記憶部４に記憶されているシミュレーションモデルＭ２を取得する（Ｓ６１）。シミュレーションモデルＭ２は、位置要素ごとに、フィラーモデルＭａ、ポリマー母相モデルＭｂ、ポリマー界面相モデルＭｃ、および、クラックモデルＭｄの少なくとも１つのモデルより構成されている。

【0073】

続いて、シミュレーションモデルＭ２全体に対して初期条件を付与する（Ｓ６２）。初期条件は、図１１に示すように、シミュレーションモデルＭ２の上下端を上下に引っ張る方向に、所定応力を付与する。ただし、初期条件として、所定応力に替えて、上下方向の所定歪を付与するようにしても良い。

【0074】

続いて、現在のシミュレーションモデルＭ２の位置要素のそれぞれにおける弾性率を現在の弾性率として、初期条件が付与されたシミュレーションモデルＭ２が力学的平衡状態となるように演算処理を行う（Ｓ６３）。最初の段階において、現在のシミュレーションモデルＭ２の位置要素のそれぞれにおける弾性率は、フィラーモデルＭａ、ポリマー母相モデルＭｂ、ポリマー界面相モデルＭｃ、クラックモデルＭｄのそれぞれにおいて初期設定された弾性率となる。そして、力学的平衡状態となるように演算処理を行うことにより、力学的平衡状態となるときのシミュレーションモデルＭ２の位置要素ごとの応力および歪を算出する。

【0075】

力学的平衡状態となるように演算処理を行うことについて、具体的には、数値解析により、式（３）に示す力学的平衡方程式（力の釣り合い方程式）を解くことにより行われる。ここで、物体が静止している状態においては、物体の内部の応力が釣り合っていることから、式（３）を満たすことになる。

【0076】

【数3】

【0077】

続いて、力学的平衡状態となるときのシミュレーションモデルＭ２が、予め設定した安定状態であるか否かを判定する（Ｓ６４）。安定状態とは、ビリアル定理に基づいて、シミュレーションモデルＭ２全体の弾性歪エネルギーＥのマイナス２倍がポテンシャルエネルギーＰに一致する状態から所定値の範囲内の状態である。最も安定な状態とは、式（４）を満たす場合である。

【0078】

【数4】

【0079】

ここで、弾性歪エネルギーＥは、シミュレーションモデルＭ２の位置要素ごとの応力および歪に基づいて算出される。また、ポテンシャルエネルギーＰは、シミュレーションモデルＭ２全体における応力および歪に基づいて算出される。

【0080】

また、安定状態は、ビリアル定理に基づく判定の他に、次のようにすることもできる。安定状態は、シミュレーションモデルＭ２の位置要素のすべての歪の変化が所定値以下となる状態とすることもできる。

【0081】

Ｓ６４の判定において、予め設定した安定状態に到達したと判定されれば（Ｓ６５：Ｙｅｓ）、演算処理部５による処理を終了する。

【0082】

一方、まだ安定状態に到達していないと判定された場合には（Ｓ６５：Ｎｏ）、ポリマー母相モデルＭｂの位置要素における変更後の弾性率を算出する（Ｓ６６）。このとき、ポリマー母相モデルＭｂの位置要素ごとの弾性率は、当該位置要素における応力および歪に対応するセカント弾性率とする。セカント弾性率は、図５に示す応力－歪の曲線上の点と起点とを結ぶ直線の傾きである。また、ポリマー界面相モデルＭｃは、一定値でない場合には、変更後のポリマー母相モデルＭｂの弾性率に応じて変更する。

【0083】

再び、Ｓ６３に戻り、マトリックスポリマーモデルＭｂ，Ｍｃの位置要素における弾性率を変更後の弾性率に変更して、再びシミュレーションモデルＭ２が力学的平衡状態となるように演算処理を行う（Ｓ６３）。そして、力学的平衡状態となるときのシミュレーションモデルＭ２の位置要素ごとの応力および歪を算出する。

【0084】

そして、再び、安定状態であるか否かを判定し（Ｓ６５）、安定状態となるまで、弾性率を変更して、力学的平衡状態となるように演算処理を行うことを繰り返す。

【0085】

このように、シミュレーションモデルＭ２が予め設定した安定状態に近づくまで、収束計算を行う。つまり、シミュレーションモデルＭ２が安定状態に近づくまで、変更後の弾性率の算出、および、シミュレーションモデルＭ２の位置要素ごとの応力および歪の算出を繰り返し行う。繰り返し処理により、最終的には、シミュレーションモデルＭ２が安定状態としての力学的平衡状態となるときのシミュレーションモデルＭ２の位置要素ごとの応力、歪および弾性率を算出する。

【0086】

１－８．ギブスエネルギー算出部６の処理
ギブスエネルギー算出部６の処理について説明する。ギブスエネルギーＧは、シミュレーションモデルＭ２が安定状態に近づいたときの力学的平衡状態になるときにおいて、弾性歪エネルギーＥに、初期条件の付与に伴って減少したポテンシャルエネルギーＰを加算することにより得られる。ギブスエネルギーＧは、式（５）により得られる。

【0087】

【数5】

【0088】

式（５）は、式（６）のように表すことができる。式（６）において、右辺第１項が、弾性歪エネルギーＥであり、右辺第２項が、ポテンシャルエネルギーＰである。

【0089】

【数6】

【0090】

１－９．エネルギー解放率算出部７の処理
エネルギー解放率算出部７の処理について図１２を参照して説明する。図１２のグラフは、横軸をクラックサイズａとし、縦軸をギブスエネルギーＧとする。クラックサイズａは、図６に示すように、クラックモデル群の所定方向長さである。本実施形態においては、図１１に示すようにシミュレーションモデルＭ２に所定応力または所定歪を付与する場合において、クラックサイズａは、所定応力または所定歪の付与方向に対して直交する方向の長さとする。

【0091】

上述したように、クラックモデル群位置設定部３によりクラックモデル群の位置および大きさが設定されることにより、シミュレーションモデルＭ２が作成される。このシミュレーションモデルＭ２におけるギブスエネルギーＧが、演算処理部５の処理およびギブスエネルギー算出部６の処理により取得できる。従って、１種のクラックモデル群を設定することにより、１個のギブスエネルギーＧが取得できる。図１２において、クラックサイズａが最も小さな値の場合のギブスエネルギーＧをプロットする。

【0092】

続いて、クラックモデル群位置設定部３により、クラックサイズａを大きくしたシミュレーションモデルＭ２を作成して、当該シミュレーションモデルＭ２におけるギブスエネルギーＧを取得する。このようにして、複数のクラックサイズａのそれぞれについてのギブスエネルギーＧを取得する。そうすると、図１２に示すように、複数のプロット点を得ることができる。

【0093】

図１２から分かるように、クラックサイズａとギブスエネルギーＧとの関係は、通常、クラックサイズａが増大するほど、ギブスエネルギーＧが減少する関係を有する。シミュレーションが適切に行われている場合は、この関係を満たす。図１２に示すように、クラックサイズａに対するギブスエネルギーＧの関係は、線形の関係を有する。

【0094】

エネルギー解放率ｇは、式（７）の関係を有する。つまり、エネルギー解放率ｇは、クラックサイズａが増大した場合に、ギブスエネルギーＧが減少する程度を表す。図１２において、プロットした点を結ぶ直線の傾きを符号反転した値が、エネルギー解放率ｇとなる。

【0095】

【数7】

【0096】

そこで、複数のクラックサイズａのそれぞれについてのギブスエネルギーＧを取得すると、複数のプロット点に基づいて一次近似直線を算出する。一次近似直線は、図１２の破線に相当する。そして、一次近似直線の傾きを符号反転させた値を、エネルギー解放率ｇとして算出する。このようにして、対象のベースモデルＭ１にクラックモデルが付加された場合におけるエネルギー解放率ｇを得ることができる。

【0097】

１－１０．作用効果
破壊靭性シミュレーション装置１は、
モデル定義領域における位置要素ごとに、マトリックスポリマーＰｂに対応するマトリックスポリマーモデルＭｂ，ＭｃまたはクラックモデルＭｄにより構成されたシミュレーションモデルＭ２であって、マトリックスポリマーモデルＭｂ，Ｍｃの弾性率とクラックモデルＭｄの弾性率とが異なる弾性率となるように設定され、かつ、少なくともマトリックスポリマーモデルＭｂ，Ｍｃの弾性率として応力および歪に応じた弾性率が設定されたシミュレーションモデルＭ２を記憶する記憶部４と、
シミュレーションモデルＭ２において、複数のクラックモデルＭｄを連続させたクラックモデル群の位置および大きさを設定するクラックモデル群位置設定部３と、
クラックモデル群位置設定部３によりクラックモデル群が設定されたシミュレーションモデルＭ２において、初期条件としてシミュレーションモデルＭ２全体に対して所定応力または所定歪を付与し、シミュレーションモデルＭ２が力学的平衡状態となるように演算処理を行う演算処理部５と、
シミュレーションモデルＭ２が力学的平衡状態となるときにおいて、弾性歪エネルギーＥに、初期条件の付与に伴って減少したポテンシャルエネルギーＰを加算することにより得られるギブスエネルギーＧを算出するギブスエネルギー算出部６と、
クラックモデル群位置設定部３により設定されたクラックモデル群の大きさを表すクラックサイズａと、ギブスエネルギー算出部６により算出されたギブスエネルギーＧとの関係に基づいて、クラックサイズａが増大した場合にギブスエネルギーＧが減少する程度を表すエネルギー解放率ｇを算出するエネルギー解放率算出部７と、を備える。

【0098】

上記の破壊靭性シミュレーション装置１によれば、シミュレーションモデルＭ２は、モデル定義領域における位置要素ごとに、マトリックスポリマーＰｂに対応するマトリックスポリマーモデルＭｂ，ＭｃまたはクラックモデルＭｄにより構成されている。そして、マトリックスポリマーモデルＭｂ，ＭｃおよびクラックモデルＭｄは、それぞれ、応力および歪に応じた弾性率が設定されている。つまり、シミュレーションモデルＭ２は、位置要素ごとに、応力および歪に応じた弾性率が設定されている。

【0099】

ここで、マトリックスポリマーＰｂは、応力に対する歪の関係を表す弾性率が線形の関係を有していない。そこで、上記のように、シミュレーションモデルＭ２に応力または歪を付与した場合に、マトリックスポリマーモデルＭｂ，Ｍｃの弾性率を、シミュレーションモデルＭ２における位置要素ごとの応力や歪に応じて異なる値に設定できる。つまり、マトリックスポリマーＰｂは、応力および歪に応じて弾性率が変化することから、シミュレーションモデルＭ２は、このようなマトリックスポリマーＰｂの弾性率の変化現象を適切に捉えることができるモデルとなる。

【0100】

このように位置要素ごとに応力および歪に応じた弾性率が設定されたシミュレーションモデルＭ２を用いて、演算処理部５が演算処理を行っている。具体的には、演算処理部５は、初期条件としてシミュレーションモデルＭ２全体に対して所定応力または所定歪を付与して、シミュレーションモデルＭ２が力学的平衡状態となるように演算処理を行う。

【0101】

そして、ギブスエネルギー算出部６が、シミュレーションモデルＭ２が力学的平衡状態となるときにおいて、弾性歪エネルギーに、初期条件の付与に伴って減少したポテンシャルエネルギーＰを加算することにより得られるギブスエネルギーＧを算出する。従って、算出されるギブスエネルギーＧは、高精度な値となる。

【0102】

エネルギー解放率算出部７が、クラックサイズａとギブスエネルギーＧとの関係に基づいて、エネルギー解放率ｇを算出している。このようにして算出されたエネルギー解放率ｇは、高精度となる。

【0103】

以上より、ポリマー材料のエネルギー解放率ｇを高精度に算出することができる。結果として、ポリマー材料の破壊靭性特性を適切に評価することができる。

【0104】

また、演算処理部５は、
現在のシミュレーションモデルＭ２の位置要素のそれぞれにおける弾性率を現在の弾性率として、初期条件が付与されたシミュレーションモデルＭ２が力学的平衡状態となるように演算処理を行うことにより、力学的平衡状態となるときのシミュレーションモデルＭ２の位置要素ごとの応力および歪を算出し、
マトリックスポリマーモデルＭｂ，Ｍｃの位置要素ごとの応力および歪に基づいて、当該マトリックスポリマーモデルＭｂ，Ｍｃの位置要素における変更後の弾性率を算出し、
マトリックスポリマーモデルＭｂ，Ｍｃの位置要素における弾性率を変更後の弾性率に変更して、再びシミュレーションモデルＭ２が力学的平衡状態となるように演算処理を行うことにより、力学的平衡状態となるときのシミュレーションモデルＭ２の位置要素ごとの応力および歪を算出し、
予め設定した安定状態に近づくまで、変更後の弾性率の算出、および、シミュレーションモデルＭ２の位置要素ごとの応力および歪の算出を繰り返し行うことにより、シミュレーションモデルＭ２が安定状態としての力学的平衡状態となるときのシミュレーションモデルＭ２の位置要素ごとの応力、歪および弾性率を算出するようにしている。

【0105】

つまり、演算処理部は、シミュレーションモデルＭ２が安定状態に近づくまで、繰り返し、弾性率を変更して、変更された弾性率を用いて力学的平衡状態となるように演算処理を行っている。ここで、マトリックスポリマーは、応力や歪が変化すると、弾性率が変化する材料である。

【0106】

そのため、シミュレーションモデルＭ２において初期条件としての弾性率を用いて力学的平衡状態になったとしても、マトリックスポリマーモデルＭｂ，Ｍｃの各位置要素においては、応力や歪が変化することに伴って、弾性率も変化することになってしまう。

【0107】

そこで、初期条件としての弾性率を用いて力学的平衡状態になったとしても、変更後の応力や歪に対応する弾性率に変化させて、再び力学的平衡状態となるように演算処理を行い、この処理を繰り返し行っている。つまり、安定状態になるまで収束計算を行っている。この繰り返し計算を行うことにより、マトリックスポリマーの各位置要素における弾性率、応力、歪の関係が、安定した状態とすることができる。そして、演算処理部５により、安定状態におけるシミュレーションモデルＭ２の位置要素ごとの応力および歪を用いて、エネルギー解放率ｇを算出するため、高精度なエネルギー解放率ｇを得ることができる。

【0108】

また、安定状態の１つの例は、弾性歪エネルギーのマイナス２倍がポテンシャルエネルギーに一致する状態から所定値の範囲内の状態である。いわゆる、ビリアル定理に基づく状態である。このように、ビリアル定理を用いて、安定状態を定義することにより、演算処理部５が算出するシミュレーションモデルＭ２の位置要素ごとの応力および歪は、高精度なものとなる。

【0109】

また、安定状態の他の例は、シミュレーションモデルＭ２の位置要素のすべての歪の変化が所定値以下となる状態である。このように、歪の変化を用いたとしても、演算処理部５が算出するシミュレーションモデルＭ２の位置要素ごとの応力および歪は、高精度とすることができる。

【0110】

また、マトリックスポリマーモデルＭｂ，Ｍｃは、図５に示すような、応力－歪の曲線上の点と起点とを結ぶ直線の傾きを表すセカント弾性率に基づいて設定されている。セカント弾性率も用いることにより、演算処理部５における繰り返し演算を高速にすることができる。

【0111】

また、クラックモデル群位置設定部３は、操作者の入力によりクラックモデル群の位置および大きさを設定するように構成されている。操作者が自由にクラックモデル群を設定することができるため、自由度が高い。

【0112】

また、ポリマー材料は、マトリックスポリマーＰｂにフィラーＰａが分散されたポリマー複合材料である。シミュレーションモデルＭ２は、モデル定義領域の位置要素ごとに、マトリックスポリマーモデルＭｂ，Ｍｃ、フィラーＰａに対応するフィラーモデルＭａ、またはクラックモデルＭｄにより構成されたシミュレーションモデルＭ２である。

【0113】

一般に、フィラーＰａが分散されたポリマー複合材料は、破壊靭性特性を評価することが容易ではない。しかし、本実施形態を適用することにより、容易に破壊靭性特性を評価することができる。

【0114】

また、マトリックスポリマーモデルＭｂ，Ｍｃは、ポリマー母相モデルＭｂ、ポリマー母相モデルＭｂとフィラーモデルＭａとの間に位置するポリマー界面相モデルＭｃ、により構成されるようにしている。ポリマー界面相モデルＭｃは、位置要素ごとに、ポリマー界面相モデルＭｃの弾性率Ｃinter（ｒ）を、所定の最小弾性率Ｃ^０（ｒ）と所定の最大弾性率Ｃ^１（ｒ）との間の割合を設定するためのフェーズフィールド変数φ（ｒ）を用いて定義されている。これにより、より高精度に破壊靭性特性を評価することができる。ただし、ポリマー界面相モデルＭｃの弾性率Ｃinter（ｒ）は、フェーズフィールド変数φ（ｒ）を用いずに一定値とすることもできる。

【0115】

また、ポリマー界面相モデルＭｃの弾性率Ｃinter（ｒ）は、フェーズフィールド変数φ（ｒ）がポリマー母相モデルＭｂ側からフィラーモデルＭａ側にかけて連続的に変化する値となるように定義されると良い。これにより、より高精度に破壊靭性特性を評価することができる。

【0116】

また、弾性歪エネルギーＥは、シミュレーションモデルＭ２の位置要素ごとの応力および歪に基づいて算出され、ポテンシャルエネルギーＰは、シミュレーションモデルＭ２全体における歪に基づいて算出される。

【0117】

（実施形態２）
２－１．破壊靭性シミュレーション装置１００の構成
実施形態２の破壊靭性シミュレーション装置１００の構成について、図１３を参照して説明する。実施形態１の破壊靭性シミュレーション装置１と同一構成については、同一符号を付して説明を省略する。

【0118】

実施形態２の破壊靭性シミュレーション装置１００は、実施形態１に対して、クラックモデル群位置設定部１０３が相違する。図１に示す実施形態１の破壊靭性シミュレーション装置１を構成するクラックモデル群位置設定部３は、操作者の入力によってクラックモデル群の位置および大きさが設定されるように構成されていた。一方、実施形態２のクラックモデル群位置設定部１０３は、初回においては、操作者の入力によってクラックモデル群の位置および大きさが設定されるように構成されている。しかし、クラックモデル群位置設定部１０３は、２回目以降においては、操作者の入力によらずに、自動的に、クラックモデル群の位置および大きさが設定されるように構成されている。

【0119】

ここで、実際のポリマー材料において、既存クラックの周辺において応力が大きな部位が新たなクラックとして進展する。そこで、シミュレーションモデルＭ２において、新たなクラックモデル群の設定についても、同様の考え方を用いて行うこととする。

【0120】

クラックモデル群位置設定部１０３は、演算処理部５により算出されたシミュレーションモデルＭ２の位置要素ごとの応力の大きさに基づいて、クラックモデルが進展する部位を判定して、新たなクラックモデル群を設定することとする。従って、操作者が初期条件としてクラックモデル群を入力しさえすれば、繰り返しクラックモデル群を変更して、演算処理部５の処理、および、ギブスエネルギー算出部６の処理が、自動的に繰り返し行われる。その結果、自動的に、エネルギー解放率算出部７がエネルギー解放率ｇを算出するために必要な、クラックサイズａとギブスエネルギーＧとを表す複数のプロット点を取得することができる。つまり、エネルギー解放率ｇの算出が容易となる。

【0121】

２－２．クラックモデル群の設定
クラックモデル群位置設定部１０３による新たなクラックモデル群の設定について、図１４を参照して説明する。

【0122】

図１４の（ａ）（ｂ）（ｃ）（ｄ）は、ポリマー母相モデルＭｂの一部およびクラックモデルＭｄの一部を図示し、ポリマー母相モデルＭｂ中に、クラックモデルＭｄの進展過程を示す。図１４（ａ）に示すように、ポリマー母相モデルＭｂ１の内部に、連続するクラックモデルＭｄ１により構成されたクラックモデル群が存在する場合を初期状態とする。

【0123】

演算処理部５による処理により、シミュレーションモデルＭ２全体に所定応力または所定歪が付与された状態において、シミュレーションモデルＭ２の位置要素ごとの応力が算出される。図１４（ｂ）に示すように、例えば、黒印で示す部位σｍａｘが、応力が最大となる部位となる。図１４（ｂ）に示すように、一般に、クラックモデルＭｄ１が連続するクラックモデル群の端部付近が、大きな応力がかかる部位となる。

【0124】

続いて、図１４（ｃ）に示すように、クラックモデル群位置設定部１０３が、シミュレーションモデルＭ２の位置要素ごとの応力に基づいて、クラックモデルＭｄ１の進展する領域ΔＭｄを決定する。つまり、拡大するクラックモデル群の位置および大きさを決定する。例えば、位置要素ごとの応力が所定値以上であるか否かを判定し、応力が所定値以上となる位置要素を、クラックモデルＭｄ１が進展する要素と決定するようにしても良い。また、位置要素ごとの応力が応力最大値σｍａｘに対して所定割合（例えば、８０％）以上となる位置要素を、クラックモデルＭｄ１が進展する要素と決定するようにしても良い。

【0125】

そうすると、図１４（ｄ）に示すように、既存クラックモデルＭｄ１によるクラックモデル群に、進展領域ΔＭｄを追加することにより、新たなクラックモデルＭｄ２によるクラックモデル群が決定される。つまり、新たなポリマー母相モデルＭｂ２の内部に、新たな連続するクラックモデルＭｄ２により構成されたクラックモデル群が位置する。そして、新たなシミュレーションモデルＭ２が自動的に作成される。

【0126】

２－３．作用効果
上述したように、クラックモデル群位置設定部１０３は、演算処理によりシミュレーションモデルＭ２の位置要素における応力に基づいて拡大するクラックモデル群の位置および大きさを決定し、拡大後のクラックモデル群の位置および大きさを、次のシミュレーションモデルＭ２を構成するクラックモデル群の位置および大きさとして設定するように構成されている。

【0127】

このように、クラックモデル群位置設定部１０３が、自動的に次のクラックモデル群を決定することができるため、容易にエネルギー解放率ｇを得ることができる。

【符号の説明】

【0128】

１，１００ 破壊靭性シミュレーション装置
２ ベースモデル作成処理部
３，１０３ クラックモデル群位置設定部
４ 記憶部
５ 演算処理部
６ ギブスエネルギー算出部
７ エネルギー解放率算出部
Ｐａ フィラー
Ｐｂ マトリックスポリマー
Ｍ１ ベースモデル
Ｍ２ シミュレーションモデル
Ｍａ フィラーモデル
Ｍｂ ポリマー母相モデル（マトリックスポリマーモデル）
Ｍｃ ポリマー界面相モデル（マトリックスポリマーモデル）
Ｍｄ クラックモデル
Ｅ 弾性歪エネルギー
Ｐ ポテンシャルエネルギー
Ｇ ギブスエネルギー
ａ クラックサイズ
ｇ エネルギー解放率

【要約】

【課題】ポリマー材料のエネルギー解放率を高精度に算出することにより、ポリマー材料の破壊靭性特性を評価することができるポリマー材料の破壊靭性シミュレーション装置を提供すること。
【解決手段】破壊靭性シミュレーション装置１は、マトリックスポリマーモデルＭｂ，ＭｃまたはクラックモデルＭｄにより構成されたシミュレーションモデルＭ２全体に対して、初期条件として所定応力または所定歪を付与し、シミュレーションモデルＭ２が力学的平衡状態となるように演算処理を行う演算処理部５と、力学的平衡状態となるときのギブスエネルギーＧを算出するギブスエネルギー算出部６と、クラックサイズａが増大した場合にギブスエネルギーＧが減少する程度を表すエネルギー解放率ｇを算出するエネルギー解放率算出部７とを備える。
【選択図】図１

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】