(57)【要約】【課題】分離型有限要素法による非圧縮Mooney-Rivlin超弾性体の変形解析を、より高速化高精度化して実用化する。【解決手段】構成方程式に低減不変量を導入して、変形前の不定圧力の初期値を零にする。非圧縮Mooney-Rivlin超弾性体の基礎方程式を、分離型有限要素法であるGSMAC有限要素法により離散化する。Newmark-β法の考えにより、速度について時間を２次の精度で近似する。各時刻において不定圧力ｐが∇・ｖ＝０を満足するように決める。各時刻の状態値に基づいて予測子の初期値を決め、予測子が収束するまで、同時緩和法により予測子の反復計算を行う。収束した予測子の値を次の時刻の状態値とする。この計算を繰り返すことにより、超弾性体の変形を解析する。【選択図】図１
2005338974.tif 000002
【特許請求の範囲】
【請求項１】
固体の基礎方程式系を分離型有限要素法により離散化して、固体の変形を解析する固体解析方法であって、前記基礎方程式系のうちの構成方程式を、圧力項とその他の項に分解して、前記圧力にGSMAC有限要素法を適用して解くことを特徴とする固体解析方法。
【請求項２】
速度を時間の２次の精度で近似するとともに、各時刻において速度の発散が零であることを不定圧力が満足する半離散化式に基づき、運動方程式を複数の式に分離することで、非圧縮超弾性体の基礎方程式系を分離型有限要素法により離散化した離散化式に従って非圧縮超弾性体の変形を解析する超弾性体解析方法であって、非圧縮超弾性体の状態の初期値を入力し、各時刻の非圧縮超弾性体の状態値に基づいて前記離散化式の予測子の初期値を決め、同時緩和法により前記予測子が収束するまで反復計算を行い、収束した予測子の値を次の時刻の非圧縮超弾性体の状態値とする計算を繰り返すことにより時間進行することを特徴とする超弾性体解析方法。
【請求項３】
前記分離型有限要素法としてGSMAC有限要素法を用いることを特徴とする請求項２記載の超弾性体解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分離型有限要素法による弾性体解析方法に関し、特に、非圧縮Mooney-Rivlin超弾性体を分離型有限要素法により解析する弾性体解析方法に関する。
【背景技術】
【０００２】
超弾性体とは、構成方程式において高次ひずみの項が存在するために、応力とひずみの間に非線形な関係がある固体である。複雑な鎖状分子構造を有する高分子材料は、超弾性体に属する。超弾性体は、非圧縮性物質であるとみなすことができる。
【０００３】
非圧縮性物質の解析では、非圧縮条件から求めた拘束条件と運動方程式を連立して解く必要がある。超弾性体解析の場合には、非圧縮条件を密度の式に代入することで、変形こう配の第３不変量が１であるという条件を導出し、これを変形の非圧縮拘束条件とする。超弾性体の有限要素法解析では、混合型有限要素法（mixed finite element method）が用いられる。混合型解法を混合補間法ということもある。混合型有限要素法とは、非圧縮拘束条件と運動方程式の連立方程式を直接的に解いて近似解を求める解法である。
【０００４】
混合型有限要素法として、Lagrange未定乗数法（Lagrange multiplier method）やペナルティ法（Penalty method）がよく用いられる。Lagrange未定乗数法では、圧力を非圧縮拘束条件に対する未定乗数として扱うのに対して、ペナルティ法では、微小な圧縮性を表すペナルティ係数を導入することで、運動方程式における圧力項を消去する。
【０００５】
しかし、混合型有限要素法では、離散化式の行列計算において直接法以外の方法を用いると安定性が悪くなるという問題がある。その理由は、非圧縮拘束条件に圧力項が現れないために、難散化式における全体係数行列の対角項に零が含まれるからである。また、未知量である変位と圧力に同じ次数の補間関数を用いると、非圧縮拘束条件を満足する解が得られないという問題もある。このことは、線形問題に限りinf-sup条件（Babuska-Brezzi-Kikuchi条件、LBB条件）を用いて証明される。変位の補間関数の次数より圧力の補間関数の次数を下げて、未知量の配置を工夫することで、その問題点を解決することができる。しかし、変位と圧力の補間関数の組合せが悪いと、解が得られない場合もある。混合型有限要素法でも動解析は可能であるが、低容量・高速計算に適していない。
【０００６】
ところで、非圧縮性流体解析の場合には、非圧縮条件を連続の式に代入することで、速度の発散が零という条件を導出し、これを流れの非圧縮拘束条件とする。非圧縮性流体の有限要素法解析では、混合型以外に、分離型有限要素法（decoupled finite element method）が用いられる。分離型有限要素法とは、運動方程式を時間についてのみ離散化した半離散化式から、非圧縮拘束条件を利用して分離式を求め、それを解いて近似解を求める間接的な解法である。このとき、分離式中に、圧力または修正速度ポテンシャルのPoisson方程式が導入される。このPoisson方程式を解くことによって、非圧縮拘束条件を満足する解を得ることができる。
【０００７】
流体の分離型有限要素法は、非圧縮性流体の差分法で発展してきたMAC法（Marker And Cell method）、SMAC法（Simplified MAC method）、HSMAC法（Highly Simplified MAC method）に代表されるMAC系解法を、有限要素法に拡張したものである。差分法の分野では、MAC系解法によって、微小な圧縮性を有する流体も解析されている。MAC法、SMAC法を拡張した分離型有限要素法では、離散化したPoisson方程式の行列計算方法として、共役こう配法などの反復法がよく用いられる。実現象を対象とする大規模な計算になると、この行列計算には膨大な計算機記憶容量と計算時間が必要になるという問題点が生じる。
【０００８】
分離型有限要素法では、解くべき離散化Poisson方程式の係数行列の大きさは、混合型有限要素法で解くべき方程式の係数行列の大きさに比べて小さくなり、係数行列の対角項は正の値となる。また、inf-sup条件を考慮する必要がないので、速度と圧力に同じ次数の補間関数を用いることも可能である。さらに、各時刻において非圧縮拘束条件がある程度の精度で満足されると、MAC系解法の特徴であるサイクル誤差自己調整法の原理により、時間進行に伴い誤差が累積することはない。これらは、混合型有限要素法に対する利点である。
【０００９】
HSMAC法を拡張した分離型有限要素法として、GSMAC有限要素法（Generalized-Simplified MAC-FEM）がある。この手法の大きな特徴は、HSMAC法に導入されている優対角近似による離散化Poisson方程式の簡易計算を利用した速度と圧力の同時緩和法を、非構造なメッシュで利用できることである。速度と圧力の同時緩和法を用いると、離散化Poisson方程式を解くための全体係数行列を作成する手順を省くことができる。これにより、複雑な形状を有する問題の大規模計算を目指した離散化Poisson方程式の低容量・高速計算が可能となる。
【００１０】
固体の有限要素法解析において、混合型有限要素法が用いられる解析例は多く見られるが、分離型有限要素法が用いられた解析例は見られない。固体解析に分離型有限要素法を適用するための本質的な困難はない。そこで、超弾性体を混合型有限要素法で解析する場合の問題点を解消するために、本発明者は、非特許文献１において、非圧縮Mooney-Rivlin超弾性体をGSMAC有限要素法で解析することを提案した。分離型有限要素法としてGSMAC有限要素法を用いる。流体解析で用いられる分離型有限要素法を固体解析に適用するためには、Euler表示の基礎方程式にLagrange法を導入し、振動問題を高精度に解析できるNewmark-β法を導入する。固体である非圧縮超弾性体を流体の計算手法である分離型有限要素法で解析できる。
【００１１】
この有限要素法の計算では、変位を得るために、初期時刻における剛性行列が必要となる。従来のアルゴリズムでは、初期時刻における剛性行列をGaussの数値積分法によって計算している。低容量・高速計算を考慮すると、GSMAC有限要素法で導入されている解析的表示を用いて剛性行列を計算するアルゴリズムに変更する必要がある。
【非特許文献１】橋本学、棚橋隆彦：「非圧縮Mooney-Rivlin超弾性体のGSMAC有限要素法に関する研究（流体と固体の連成有限要素法解析に向けて）」、計算工学会講演論文集、Vol.8(2003年5月)．
【発明の開示】
【発明が解決しようとする課題】
【００１２】
しかし、非特許文献１において提案した解析方法を、より高精度・低容量で高速計算ができるようにして実用化するためには、以下のような課題を解決する必要がある。有限要素法の形状関数として、双１次関数を使用しているが、より高精度な計算を行うために、高次補間の形状関数も使用可能にする必要がある。例えば、双１次形状関数を用いた変位型有限要素法では、曲げ変形を行うと得られる解が固くなること（Shear locking）が知られている。高次補間の形状関数を用いると、このようなことは起こらない。本発明は、分離型有限要素法による非圧縮Mooney-Rivlin超弾性体の解析を、より高速化高精度化して実用化することを目的とする。
【課題を解決するための手段】
【００１３】
上記の課題を解決するために、本発明では、固体の基礎方程式系を分離型有限要素法により離散化して、固体の変形を解析する固体解析方法を、基礎方程式系のうちの構成方程式を、圧力項とその他の項に分解して、圧力にGSMAC有限要素法を適用して解く方法とした。
【００１４】
また、超弾性体解析方法を、超弾性体の基礎方程式を分離型有限要素法により離散化し、Newmark-β法により速度について時間を２次の精度で近似し、各時刻において不定圧力ｐが∇・ｖ＝０を満足するように決め、各時刻の状態値に基づいて予測子の初期値を決め、予測子が収束するまで同時緩和法により予測子の反復計算を行い、収束した予測子の値を次の時刻の状態値とする計算を繰り返すことにより、超弾性体の変形を解析する方法とした。この解析計算においては、変形前の不定圧力の初期値を有しないようにするために、低減不変量を導入したMooney-Rivlin体の構成方程式を用いる。
【発明の効果】
【００１５】
本発明では、上記のように構成したことにより、非圧縮超弾性体の変形を、より高速に高精度で解析でき、十分実用的に使用できる。
【発明を実施するための最良の形態】
【００１６】
以下、本発明を実施するための最良の形態について、図１〜図22を参照しながら詳細に説明する。
【実施例】
【００１７】
本発明の実施例は、超弾性体の基礎方程式を分離型有限要素法により離散化し、Newmark-β法により速度について時間を２次の精度で近似し、各時刻において不定圧力ｐが∇・ｖ＝０を満足するように決め、各時刻の状態値に基づいて予測子の初期値を決め、予測子が収束するまで同時緩和法により予測子の反復計算を行い、収束した予測子の値を次の時刻の状態値とする計算を繰り返すことにより、超弾性体の変形を解析する超弾性体解析方法である。
【００１８】
図１は、本発明の実施例における超弾性体解析方法の概略流れ図である。図２は、超弾性体解析方法において前提とする点の動きと非圧縮条件の説明図である。図３は、超弾性体解析方法で用いるGSMAC有限要素法の流れ図である。図４は、超弾性体解析方法で用いる同時緩和法の流れ図である。図５は、平面ひずみ状態（引張変形）の解析モデルの説明図である。図６は、引張変形解析モデルの境界条件を示す図である。図７は、引張変形解析計算に使用した共通解析条件のデータを示す表である。図８は、引張変形解析計算に使用したケース別解析条件のデータを示す表である。図９は、引張変形前の解析メッシュを示す図である。図10は、引張変形ひずみの経時変化を示す図である。図11は、引張変形解析領域の全面積の経時変化を示す図である。図12は、引張変形後の解析メッシュを示す図である。図13は、引張変形ひずみの静的状態における相対誤差と厳密値を示す表である。
【００１９】
図14は、平面ひずみ状態（せん断変形）の解析モデルの説明図である。図15は、せん断変形解析モデルの境界条件を示す図である。図16は、せん断変形解析計算に使用した共通解析条件のデータを示す表である。図17は、せん断変形解析計算に使用したケース別解析条件のデータを示す表である。図18は、せん断変形前の解析メッシュを示す図である。図19は、せん断変形ひずみの経時変化を示す図である。図20は、せん断変形解析領域の全面積の経時変化を示す図である。図21は、せん断変形後の解析メッシュを示す図である。図22は、相対誤差と厳密値を示す表である。
【００２０】
本発明の実施例における超弾性体解析方法の原理と手順を説明する。例として、非圧縮Mooney-Rivlin超弾性体を解析対象とするが、一般の固体に適用可能である。本方法は一般化すれば、固体の基礎方程式系を分離型有限要素法により離散化して、固体の変形を解析する固体解析方法である。基礎方程式系のうちの構成方程式を、圧力項とその他の項に分解して、圧力にGSMAC有限要素法を適用して解く。
【００２１】
最初に、図１を参照しながら、超弾性体解析方法の概要を説明する。高分子材料や生体軟組織などの非圧縮Mooney-Rivlin超弾性体の基礎方程式を立てる。その際に、低減不変量（reduced invariant）を導入して、無変形状態で不定圧力が初期値を持たないようにする。この基礎方程式を、分離型有限要素法であるGSMAC有限要素法により離散化する。その際、弾性振動問題を精度よく解析できるようにするために、Newmark-β法の考えを導入して、速度に関しては時間を２次の精度にして近似する。また、時間の離散化である半離散化式を立てる際に、各時刻において、不定圧力ｐが∇・ｖ＝０を満足するように決める。運動方程式を分離して、固体要素の離散化式を立てる。各時刻の状態値に基づいて予測子の初期値を決め、予測子が収束するまで、同時緩和法により予測子の反復計算を行う。収束した予測子の値を次の時刻の状態値とする。この計算を繰り返すことにより、超弾性体の変形を解析する。
【００２２】
次に、非圧縮Mooney-Rivlin超弾性体の基礎方程式について説明する。基礎方程式は、以下に示すように、密度の式、Cauchyの運動方程式、非圧縮超弾性体の構成方程式、弾性ポテンシャル関数定義式、低減不変量定義式、変形−変形こう配関係式、変形こう配−変位こう配関係式である。
detＦ＝１                                                    （１ａ）
ρ(d2ｕ/dt2)＝∇・Ｔ＋ρｂ                                   （１ｂ）
Ｔ＝−ｐＩ＋２{(∂Ｗ/∂IIIB)(IIIB)I
＋((∂Ｗ/∂IB)＋(∂Ｗ/∂IIB)IB)Ｂ
−(∂Ｗ/∂IIB)Ｂ・Ｂ}                                      （１ｃ）
Ｗ＝ｃ1({IB}−３)＋ｃ2({IIB}−３)                            （１ｄ）
{IB}＝IB/(IIIB)1/3                                           （１ｅ）
{IIB}＝IIB/(IIIB)2/3                                         （１ｆ）
Ｂ＝Ｆ・ＦT                                                  （１ｇ）
Ｆ＝Ｉ＋ｕ∇←x                                              （１ｈ）
ただし、Ｆは変形こう配、ρは密度、ｕは変位、d/dtは実質時間微分、∇は空間座標のナブラ、ＴはCauchy応力、ｂは単位質量当たりの体積力、ｐは不定圧力、Ｉは恒等テンソル、Ｂは左Cauchy-Green変形、Ｗは弾性ポテンシャル関数である。∇←xは、物質座標に関する後形（右形）ナブラである。
【００２３】
IB，IIB，IIIBはそれぞれ、
IB＝trＢ                                                     （２ａ）
IIB＝（1/2）{(trＢ)2−tr(Ｂ・Ｂ)}                            （２ｂ）
IIIB＝detＢ                                                  （２ｂ）
で定義される不変量(invariant)、{IB},{IIB}は低減不変量(reduced invariant)、ｃ1，ｃ2はMooney-Rivlin定数、∇Xは物質座標のナブラである。空間座標における位置ｘと物質座標における位置Ｘの間には
ｘ＝Ｘ+ｕ                                                    （３）
の関係がある。また、運動の表現方法としてLagrange法が用いられる。
【００２４】
基礎方程式を
∇・ｖ＝０                                                   （４ａ）
ρ(dｖ/dt)＝−∇ｐr＋∇・τ＋ρｂ                            （４ｂ）
dｕ/dｔ＝ｖ                                                  （４ｃ）
ｐr＝ｐ＋(2/3)[ｃ1(trＢ)＋ｃ2{(trＢ)2‐tr(Ｂ・Ｂ)}]          （４ｄ）
τ＝２{ｃ1＋ｃ2(trＢ)}Ｂ‐２ｃ2Ｂ・Ｂ                        （４ｅ）
Ｂ＝Ｆ・ＦT                                                  （４ｆ）
Ｆ＝Ｉ＋ｕ∇←x                                              （４ｇ）
のように変形する。ただし、ｖは速度である。超弾性体の混合型有限要素法解析では、非圧縮拘束条件としてdetＦ＝１またはIIIB＝１の条件が一般に用いられる。しかし、分離型有限要素法によって非圧縮超弾性体の基礎方程式を離散化するので、非圧縮流体解析で用いられる式(４ａ)を用いる。
【００２５】
次に、基礎方程式を分離型有限要素法によって離散化する方法を説明する。時間進行において重要となることは２つある。１つ目は、弾性振動問題を解析できるように速度に対して時間２次精度であるNewmark-β法の考えを導入することである。２つ目は、時刻ｔ＝ｔn+1におけるｖn+1＝ｖ(ｘ,ｔn+1)の分布と物質点の位置ｘn+1＝ｘ(Ｘ,ｔn+1)を用いて定義される
∇n+1・ｖn+1≡∇・{ｖ(ｘ,ｔn+1)}|x=xn+1
が零であることを満足するように不定圧力ｐn+1を決めることである(図２参照)。
【００２６】
Newmark-β法の考えを導入したMAC系解法の時間進行方法を説明する。ｂn+1を既知として、基礎方程式をNewmark-β法に基づいて時間的に離散化すると、
∇n+1・ｖn+1＝０                                              （５ａ）
ρ[ｖn+1]＝−∇n+1ｐrn+1＋∇n+1・τn+1＋ρｂn+1               （５ｂ）
ｖn+1＝ｖn＋(Δｔ){γ[ｖn+1]＋(１−γ)[ｖn]}                  （５ｃ）
ｕn+1＝ｕn＋(Δｔ)ｖn＋(Δｔ)2{β[ｖn+1]＋((1/2)−β)[ｖn]}   （５ｄ）
ｐn+1＝ｐrn+1−(2/3)[ｃ1(trＢn+1)＋ｃ2{(trＢn+1)2−tr(Ｂn+1・Ｂn+1)}]
（５ｅ）
τn+1＝２{ｃ1＋ｃ2(trＢn+1)}Ｂn+1−２ｃ2Ｂn+1・Ｂn+1          （５ｆ）
Ｂn+1＝Ｆn+1・(Ｆn+1)T                                        （５ｇ）
Ｆn+1＝Ｉ＋ｕn+1∇←x                                         （５ｈ）
となる。ただし、[ｖ]はdｖ/dtを意味し、γ＝1/2，β＝1/6の場合には線形加速度法に対応する。
【００２７】
式(５)を近似的に満足するための陽的な多段階反復計算方法を考える。
ｖ(0)＝ｖn＋(Δｔ)[ｖn]                                      （６ａ）
ｕ(0)＝ｕn＋(Δｔ)ｖn＋(Δｔ)2[ｖn]/2                        （６ｂ）
を計算し、ｐr(0)＝ｐnrとおく。ｍ＋１回目の反復計算では、速度ｖn+1の予測子ｖ(m+1)が条件∇(m+1)・ｖ(m+1)＝０を満足すると、ｐrn+1の予測子pr(m+1)が求まるように、運動方程式の半離散化式を定める。
【００２８】
すると、
Ｆ(m)＝Ｉ＋ｕ(m)∇←x                                          （７ａ）
Ｂ(m)＝Ｆ(m)・(Ｆ(m))T                                         （７ｂ）
τ(m)＝２{ｃ1＋ｃ2(trＢ(m))}Ｂ(m)−２ｃ2Ｂ(m)・Ｂ(m)           （７ｃ）
∇(m)・ｖ(m+1)＝０                                             （７ｄ）
ρ[ｖ(m+1)]＝−∇(m)ｐr(m+1)＋∇(m)・τ(m)＋ρｂn+1            （７ｅ）
ｖ(m+1)＝ｖ(0)＋γ(Δｔ)([ｖ(m+1)]−[ｖn])                     （７ｆ）
ｕ(m+1)＝ｕ(0)＋β(Δｔ)2([ｖ(m+1)]−[ｖn])                    （７ｇ）
となる。
【００２９】
さらに、MAC系解法の時間進行より、式(７ｄ)，式(７ｅ)，式(７ｆ)を
（予測子ステップの式）
ρ{[ｖ]}＝−∇(m)ｐr(m)＋∇(m)・τ(m)＋ρｂn+1                （８ａ）
{ｖ}＝ｖ(0)＋γ(Δｔ)({[ｖ]}−[ｖn])                          （８ｂ）
（Poisson方程式）
(∇(m))2φp＝ρ∇(m)・{ｖ}                                    （８ｃ）
ρ[ｖ(m+1)]＝ρ{[ｖ]}−∇(m)φp/(γ(Δｔ))                    （８ｄ）
ρｖ(m+1)＝ρ{ｖ}−∇(m)φp                                   （８ｅ）
ｐr(m+1)＝ｐr(m)＋φp/(γ(Δｔ))                              （８ｆ）
のように分離する。ただし、修正速度ポテンシャルφpは
φp＝γδｐr(Δｔ)                                            （９）
で表される。
【００３０】
そして、収束判定条件
|ｖ(m+1)−ｖ(m)|≒０                                          （10）
を満足するときに
[ｖn+1]≒[ｖ(m+1)]
ｖn+1≒ｖ(m+1)
ｐrn+1≒ｐr(m+1)
ｕn+1≒ｕ(m+1)                                                （11）
とする。このように時間進行した後に、式(５ｅ)よりｐn+1を求める。
【００３１】
時間的な離散化式をGalerkin有限要素法で離散化する。
ｖα(0)＝ｖαn＋(Δｔ)[ｖαn]                                  （12ａ）
ｕα(0)＝ｕαn＋(Δｔ)ｖαn＋(Δｔ)2[ｖαn]/2                  （12ｂ）
を計算し、ｐre(0)＝ｐrenとおく。ｍ＋１回目の計算における有限要素Ωe(m)または基準配置の有限要素Ωe0における離散化式を示す。
【００３２】
（予測ステップの式）
ρ[Ｍ]αβ(m){[ｖ]}β＝Ｃα(m)ｐre(m)−２(ｃ1＋２ｃ2)Ｓα
−２[(ｃ1＋ｃ2)(trＳαβL)Ｉ＋ｃ2{２ＳαβL−(ＳαβL)T}]・ｕβ(m)
−４ｃ2ＳγβαN1・(ｕβ(m)ｕγ(m))T
−２ｃ2(ＳαβγN1−ＳγβαN1)(ｕβ(m)・ｕγ(m))
＋２ｃ2(ＳαβγN1＋ＳγβαN1)・(ｕβ(m)ｕγ(m))
−２ｃ2ＳαβγδN2{ｕβ(m)(ｕγ(m)・ｕδ(m))−(ｕβ(m)・ｕγ(m))ｕδ(m)}
＋ρ[Ｍ]αβ(m)ｂβn+1＋∫Γe0Ｎa{ｔ}(m)dΓ0                    （13ａ）
{ｖ}α＝ｖα(0)＋γ(Δｔ)({[ｖ]}α−[ｖ]αn)                      （13ｂ）
【００３３】
（Poisson方程式）
∫Γe(m)(∂φp/∂ｎ(m))dΓ(m)＝ρＣβ(m)・{ｖ}β                  （13ｃ）
（修正子ステップの式）
ρ[Ｍ]αβ(m)[ｖ]β(m+1)＝ρ[Ｍ]αβ(m){[ｖ]}β＋Ｃα(m)(φpe/γ(Δｔ))
−∫Γe(m)Ｎα(φp/γ(Δｔ))ｎ(m)dΓ(m)                         （13ｄ）
ρ[Ｍ]αβ(m)ｖβ(m+1)＝ρ[Ｍ]αβ(m){ｖ}β＋Ｃα(m)φpe
−∫Γe(m)Ｎα(φpｎ(m))dΓ(m)                                  （13ｅ）
ｐre(m+1)＝ｐre(m)＋(φpe/γ(Δｔ))                               （13ｆ）
ｕα(m+1)＝ｕα(0)＋β(Δｔ)2([ｖ]α(m+1)−[ｖ]αn)               （13ｇ）
【００３４】
ただし、Ｎαを双１次形状関数とすると、係数行列は
（質量行列）
Ｍαβ(m)＝∫Ωe(m)ＮαＮβdΩ(m)                                 （14ａ）
（勾配行列）
Ｃα(m)＝∫Ωe(m)∇(m)ＮαdΩ(m)                                  （14ｂ）
（剛性行列）
Ｓα＝∫Ωe0∇xＮαdΩ0                                           （14ｃ）
ＳαβL＝∫Ωe0∇xＮα∇xＮβdΩ0                                 （14ｄ）
ＳαβγN1＝∫Ωe0∇xＮα(∇xＮβ・∇xＮγ)dΩ0                    （14ｅ）
ＳαβγδN2＝∫Ωe0(∇xＮα・∇xＮβ)(∇xＮγ・∇xＮδ)dΩ0       （14ｆ）
で表され、[Ｍ]αβ(m)は集中質量行列である。また、Γe(m)またはΓe0は要素境界であり、{ｔ}(m)は公称応力である。
【００３５】
さらに、式(13ｃ)には速度と圧力の同時緩和法を利用した簡易計算を行う。そして、収束条件|ｖα(m+1)−ｖα(m)|＜ε2を満足するとき、
[ｖαn+1]≒[ｖα(m+1)]
ｖαn+1≒ｖα(m+1)
ｐren+1≒ｐre(m+1)
ｕαn+1≒ｕα(m+1)                                                （16）
とする。
【００３６】
このように時間進行した後に、
ｐen+1＝ｐren+1−(２ｃ1＋４ｃ2)
−(2/(3Ωe0))[２(ｃ1＋４ｃ2)Ｓα・ｕαn+1
＋(ｃ1＋４ｃ2)(trＳαβL)(ｕαn+1・ｕβn+1)
＋４ｃ2ＳαβL:(ｕαn+1ｕβn+1−ｕβn+1ｕαn+1)
＋４ｃ2ＳαβγN1・{ｕαn+1(ｕβn+1・ｕγn+1)−ｕγn+1(ｕβn+1・ｕαn+1)}
＋ｃ2ＳαβγδN2・{(ｕαn+1・ｕβn+1)・(ｕγn+1・ｕδn+1)
−(ｕαn+1・ｕγn+1)・(ｕβn+1・ｕδn+1)}]                       （17）
よりpen+1を求める。以上の数値計算法を含む全体の流れ図を図３に示し、速度圧力の同時緩和法による反復計算の流れ図を図４に示す。ただし、図４における<∇・ｖ>eは速度の発散の要素平均値を意味し、ε1は非圧縮拘束条件の収束判定基準である。
【００３７】
次に、解析例の平面ひずみ問題として、引張変形と単純せん断変形を説明する。応力に粘性応力項を加えて粘弾性体として扱い、ひずみ遅延を起こさせることによって静的状態を実現する。そして、静的な状態における厳密解との比較を行う。静的な状態の判定には|ｖn+1−ｖn|＜ε3を用いる。粘弾性体のモデルとして、式(４ｅ)を
τ＝２{ｃ1＋ｃ2(trＢ)}Ｂ‐２ｃ2Ｂ・Ｂ＋λt{４(ｃ1＋ｃ2)Ｄ}      （18）
に変更したVoigtモデルの構成方程式を用いる。ここで、λtは１次元の場合における遅延時間に相当するパラメータであり、Ｄは変形速度であり、
Ｄ＝(1/2)(∇ｖ＋ｖ∇←)                                         （19）
で表される。
【００３８】
図５のような平面ひずみ状態の解析モデルを考える。ただし、図５においてｘ＝０を対称面と仮定する。モデルの実線部分は基準配置、点線部分は現配置を表している。境界条件を図６に示すように与える。図６に示されるＷは厚さであり、平面ひずみ状態を仮定するために変形において一定である。ｌとｈの厳密値は
(ｌ/Ｌ)4−(Ｐ0/(２(ｃ1＋ｃ2)ＨＷ))(ｌ/Ｌ)3−１＝０            (20ａ)
ｈ/Ｈ＝Ｌ/ｌ                                                  (20ｂ)
より求めることができる。図５の引張変形におけるｌとｈについて検証する。図７と図８に、計算に使用した解析条件のデータを示す。ただし、図７において、代表長さはＬr＝Ｈとし、代表速度はＵr＝√(ｃ1＋ｃ2)/ρとする。減衰効果とΔｔの制限から考えて、λt＝0.01sを用いる。また、図９に計算に使用したメッシュを示す。計算の初期条件として、[ｖ]＝０，ｖ＝０，ｕ＝０，ｐ＝０(ｐr＝２ｃ1＋４ｃ2)を用いる。
【００３９】
図10に、計算から求めたGreen-LagrangeひずみＥの履歴を示す。ただし、Green-Lagrangeひずみは
Ｅ＝(1/2)(ＦT・Ｆ−Ｉ)                                          (21)
で表される。図10を見ると、ある時刻からひずみが一定値になることがわかる。これは、静的な状態を実現できたことを意味する。図11に、解析領域の全面積の時刻歴を示す。図11を見ると、全面積の変化は1.0×10-5％未満であることがわかる。全面積がほとんど変化しないことから、本解析では非圧縮拘束条件∇・ｖ＝０を十分に満足している。図12に、静的な状態とみなされる時刻におけるメッシュを示す。図12を見ると、メッシュが崩壊することなく、整然と配置されていることがわかる。図13に静的状態における相対誤差を厳密値と共に示す。図13を見ると、近似値は厳密値と良好に一致することがわかる。
【００４０】
次に、単純せん断変形の解析例について説明する。図14のような平面ひずみ状態の解析モデルを考える。モデルの実線部分は基準配置、点線部分は現配置を表している。境界条件を図15に示すように与える。角度が２γ0の応力状態が厳密値である。図14に示されるせん断変形におけるｕ0について検証する。代表長さはＬr＝Ｌ、代表速度はＵr＝√((ｃ1＋ｃ2)/ρ)とする。図16と図17に、計算に使用した解析条件のデータを示す。また、図18に、計算に使用したメッシュを示す。計算の初期条件として[ｖ]＝０，ｖ＝０，ｕ＝０，ｐ＝０(ｐr＝２ｃ1＋４ｃ2)を用いる。
【００４１】
図19に、計算から求めたGreen-Lagrangeひずみの履歴Ｅを示す。図19を見ると、静的な状態を実現できたことがわかる。図20に、解析領域の全面積の時刻歴を示す。図20を見ると、全面積の変化が大きいCase２でさえ、その変化は2.0×10-5％未満であることがわかる。よって、非圧縮拘束条件を十分に満足している。図21に、静的な状態とみなすことができる時刻におけるメッシュを示す。図21を見ると、メッシュが崩壊することなく、整然と配置されていることがわかる。図22に静的な状態における近似値と相対誤差を厳密値と共に示す。図22を見ると、近似値は厳密値と良好に一致することがわかる。
【００４２】
上記のように、本発明の実施例では、超弾性体解析方法を、超弾性体の基礎方程式を分離型有限要素法であるGSMAC有限要素法により離散化し、Newmark-β法により速度について時間を２次の精度で近似し、各時刻において不定圧力ｐが∇・ｖ＝０を満足するように決め、各時刻の状態値に基づいて予測子の初期値を決め、予測子が収束するまで同時緩和法により予測子の反復計算を行い、収束した予測子の値を次の時刻の状態値とする計算を繰り返すことにより、超弾性体の変形を解析する方法としたので、超弾性体の変形を十分に高速に、高精度で解析できる。
【産業上の利用可能性】
【００４３】
本発明の超弾性体の分離型有限要素法は、医療などの分野に応用可能である。高齢社会の急速な移行に伴って、福祉・医療分野が重要視されている。生体内の力学現象を解明することは、病的原因の解明やより効果的な手術方法の確立につながる。生体内の固体は軟組織を形成するものが多く、血管壁や心臓壁や眼球網膜など、超弾性体とみなすことができるものも多い。よって、超弾性体の有限要素法解析の重要性が増大する。
【００４４】
非定常な流体力が作用することによってゴム状構造物が振動する問題は工学上多く見られる。低容量・高速計算に適した動解析手法である本発明の分離型有限要素法は、このような問題を解析するのに適している。
【図面の簡単な説明】
【００４５】
【図１】本発明の実施例における超弾性体解析方法の概略流れ図、
【図２】超弾性体解析方法において前提とする点の動きと非圧縮条件の説明図、
【図３】超弾性体解析方法で用いるGSMAC有限要素法の流れ図、
【図４】超弾性体解析方法で用いる同時緩和法の流れ図、
【図５】平面ひずみ状態（引張変形）の解析モデルの説明図、
【図６】引張変形解析モデルの境界条件を示す図、
【図７】引張変形解析計算に使用した共通解析条件のデータを示す表、
【図８】引張変形解析計算に使用したケース別解析条件のデータを示す表、
【図９】引張変形前の解析メッシュを示す図、
【図１０】引張変形ひずみの経時変化を示す図、
【図１１】引張変形解析領域の全面積の経時変化を示す図、
【図１２】引張変形後の解析メッシュを示す図、
【図１３】引張変形ひずみの静的状態における相対誤差と厳密値を示す表、
【図１４】平面ひずみ状態（せん断変形）の解析モデルの説明図、
【図１５】せん断変形解析モデルの境界条件を示す図、
【図１６】せん断変形解析計算に使用した共通解析条件のデータを示す表、
【図１７】せん断変形解析計算に使用したケース別解析条件のデータを示す表、
【図１８】せん断変形前の解析メッシュを示す図、
【図１９】せん断変形ひずみの経時変化を示す図、
【図２０】せん断変形解析領域の全面積の経時変化を示す図、
【図２１】せん断変形後の解析メッシュを示す図、
【図２２】せん断変形解析計算の相対誤差と厳密値を示す表である。
【図１】
2005338974.tif 000003
【図２】
2005338974.tif 000004
【図３】
2005338974.tif 000005
【図４】
2005338974.tif 000006
【図５】
2005338974.tif 000007
【図６】
2005338974.tif 000008
【図７】
2005338974.tif 000009
【図８】
2005338974.tif 000010
【図９】
2005338974.tif 000011
【図１０】
2005338974.tif 000012
【図１１】
2005338974.tif 000013
【図１２】
2005338974.tif 000014
【図１３】
2005338974.tif 000015
【図１４】
2005338974.tif 000016
【図１５】
2005338974.tif 000017
【図１６】
2005338974.tif 000018
【図１７】
2005338974.tif 000019
【図１８】
2005338974.tif 000020
【図１９】
2005338974.tif 000021
【図２０】
2005338974.tif 000022
【図２１】
2005338974.tif 000023
【図２２】
2005338974.tif 000024