超弾性

目次

はじめに

ゴムや発泡体のような材料は大きい変形にさらされてもまだ完全に弾性のままであり得る.つまり荷重が取り除かれると,変形は塑性変形のない完全に可逆となるのである.これらは超弾性材料またはグリーン弾性体と言われる.ゴムや発泡体以外にも,ゴムのような形態の生物組織やポリマーが超弾性材料のカテゴリに入る.亜弾性材料とは対照的に,超弾性材料は大きい変形や回転を受ける.

小さい変形の場合,変形の前後におけるオブジェクトの形状の差は,オブジェクトのサイズと比較して小さい.大きい変形の場合はもはやそうはいかず,変形を説明する必要がある.例えば,大きい変形は力が作用している面の形状に影響を与える場合がある.小さい変形の場合は,形状の変化はごくわずかなので無視される.

大きい変形に潜む数式は変形勾配で紹介してある.先に進む前に,このセクションに関係のある用語とその公式をまとめる:

ここまでで見てきた応力測定はコーシー応力 だけである.コーシー応力は空間と変形構成で定義された応力測定であり,最終的に計算したいものである.しかし超弾性材料モデルでは他の応力測定が使われる.これは,一般的に変形構成のモデルを定義することがずっと難しいからである.それよりも非変形材料構成のモデルを定義した方がよい.非変形構成で応力を計算したら,コーシー応力を得ることができる.これは後で示す.ここでは平衡方程式が力密度 と第一ピオラ・キルヒホッフ応力 を関連付ける.基準構成の運動方程式は以下で与えられる:

ここで はともに非変形構成の体積力密度と質量密度である.つまり,密度は材料の初期体積で割った質量のことである. は未知の変位ベクトルである. は[]で測定された第一ピオラ・キルヒホッフ応力であり非変形構成で定義されているため非変形面積あたりの力 である.第一ピオラ・キルヒホッフ応力はPK1と省略されることが多い.

初期構成で設定を表したいため,ここまで平衡方程式にピオラ・キルヒホッフ応力を使い,最終の変形構成における応力を求めたいならばコーシー応力を使う.しかしそれ以上のことが必要である.コーシー応力は現在の面積あたりの現在の力 である.

超弾性材料モデルを確立するための一般的な方法は,スカラーのひずみエネルギー密度関数 から始めることである.ひずみエネルギー密度関数 は実験から構築され,指定された材料を定義する.これは変形によって材料に保存されたエネルギーの測定である.単位非変形体積あたりのこのエネルギー密度関数はひずみ測定の関数として公式化される.それからエネルギー密度関数はそのひずみ測定および応力測定の結果について微分される.

構成挙動が現在の変形状態のみの関数として表せる場合,その材料は弾性と考えられる.この条件が真である場合,位置 における粒子の応力測定は現在の変形勾配 の関数である.変形勾配 は第一ピオラ・キルヒホッフ応力 と共役であり,以下のように書ける:

関連する変形率テンソル を持つ応力テンソル の二重縮約 は単位基準体積あたりの内部の機械的仕事率を記述する.したがってひずみエネルギー密度関数は以下で与えられる[17, c. 5.2, p. 142]:

構成方程式は,互換の応力測定とひずみ測定の間の関係を記述しなければならない[5, p. 498].第一ピオラ・キルヒホッフ応力 と変形勾配 はそのような互換の応力・ひずみのペアである.ひずみエネルギー密度関数は,異なるひずみ測定についての同等の形式で表すことができる.

立体の構成挙動は材料座標で与えられることが多く,固体力学ではラグランジュ記述が好まれる[16, p. 61].固体力学の設定では,現在の構成のある部分は通常簡単に分からないため,空間構成で表される応力テンソルを使うのは便利ではない.ほとんどの場合,基準構成または中間構成で表される応力テンソルを使った方が簡単である[16, p. 109, p. 146].エネルギー密度関数をPK1と変形勾配 で表すだけでなく,右コーシー・グリーン変形テンソル でも表すために,新しい疑似応力測定 を導入する. は[]で測定される第二ピオラ・キルヒホッフ応力である.第二ピオラ・キルヒホッフ応力PK2と略されることが多い.疑似応力とは,この応力が非変形構成でも変形構成でも定義されていないことを意味する.

応力場 がひずみ場 と共役に動作するように,応力場 はひずみ場 と共役に動作する[16, p. 155ff].言い換えると,グリーン・ラグランジュひずみ は第二ピオラ・キルヒホッフ応力 と互換なのである[17, c. 5.2, p. 143]:

他の応力や共役ひずみも存在するが,ここでは触れない.

上で説明したように,ひずみエネルギー密度関数は異なるひずみ測定の同等の形式で表すことができる.例えばエネルギー密度関数 はグリーン・ラグランジュひずみ で表すことができる.エネルギー密度関数のその他の表し方として,例えばコーシーテンソル ,微小ひずみ の不変式を使うというものがある.エネルギー密度関数がどのように表されるか(例えば )によって,エネルギー関数はその独立変数(例えば )について微分される.

ひずみエネルギー密度関数が変形勾配 で表され,そのひずみエネルギー密度関数の導関数が変形勾配 について取られると,PK1応力が得られる:

変形勾配 が第一ピオラ・キルヒホッフ応力 と互換だからである.

グリーン・ラグランジュひずみ ,コーシー変形テンソル ,またはひずみ測定についての微小ひずみ の関数である,ひずみエネルギー密度関数 は結果的に第二ピオラ・キルヒホッフ応力 となる.

ひずみ測定 は第二ピオラ・キルヒホッフ応力 と互換だからである.

エネルギー密度関数が第二ピオラ・キルヒホッフ応力 と互換の形式で表される場合,それが平衡方程式で使えるように第一ピオラ・キルヒホッフ応力に変換されなければならない:

結果の第二ピオラ・キルヒホッフ応力 は次の関係で第一ピオラ・キルヒホッフ応力 に変換することができる:

そして平衡方程式は以下になる:

ナンソンの公式は,応力がどのように互いに変換されるかを記述するものであり,以下のリストのようにまとめられる[5, p. 482]:

文献読むことについてさらに複雑になるのは,著者によって応力の定義が異なり,転置されたものを使う場合があるということである[11, ch2.2.3 p. 44].ここでは,例えば を使った場合,2つ目の下付き文字 は面法線ベクトルに相当する.

材料と空間構成の間の数量の変換は押し出し,引き戻しと呼ばれる[16, p. 82].

要約すると以下のようになる:

参考までに,微小ひずみ測定 の場合,応力測定 はおおよそ同じであり となる.この意味で,微小ひずみ もPK2 と互換である.

有限要素パッケージをロードし,$HistoryLengthを0に設定する:
変数を設定する:
直方体のメッシュを生成する:

サンブナン・キルヒホッフモデル

サンブナン・キルヒホッフモデルは亜弾性材料の簡単なモデルである.サンブナン・キルヒホッフモデルは現象モデルである.つまり,これは観察された挙動を記述するのである.このモデルは大きい変形の形態において欠損があるため,その目的は主に教育的なものである.欠損がある場合,サンブナン・キルヒホッフモデルは圧縮の下で柔らかくなるが,これは物理的なものではない.この場合は,カスタムの構成方程式を定式化する例として,サンブナン・キルヒホッフモデルを使う.

最も簡単な構造モデルは線形弾性である.ひずみで定義すると,エネルギー密度は以下で与えられる[14, p 13]:

ここで は微小ひずみテンソル, は第一ラメ定数, は剛性率 としても知られる第二定数である.ラメ定数は,別の材料パラメータのセクションで説明してあるようにヤング率およびポアソン比に関連付けることができる.は,線形弾性で対応するエネルギー密度関数である.

サンブナン・キルヒホッフモデルは,上のエネルギー密度関数を使うが,微小ひずみテンソル の代りにグリーン・ラグランジュひずみ を使う.

グリーン・ラグランジュひずみ の場合,エネルギー密度関数 は以下で与えられる[14, p 15]:

グリーン・ラグランジュひずみ についてのエネルギー密度関数の導関数は以下になる:

ここで は第二ピオラ・キルヒホッフ応力, は単位行列である.

超弾性材料の法則は非常に非線形であるため,解くのが難しい.最終的な変形の量によっては,線形化された方程式の集合を反復的に解くという標準的なアプローチで解を求めることは現実的ではない.これはWolfram言語の問題ではなく,すべてのPDE解法が直面するチャレンジである.このチャレンジに挑むために2つのワークフローが利用できる.どちらのアプローチも,物体に作用する力が徐々に増加する.この2つのアプローチとは以下である.

最初のアプローチは,変形の量に影響を及ぼすパラメータ を持つパラメトリック関数を作成するというものである.最初このパラメータ は0に設定し,方程式系の解を得る.次にこの解を, をわずかに増加させたパラメトリック系の次の評価の初期シードとして使う.今度はその解が,パラメータ が,システムの完全変形に同等の1になるまで同様に使われる.この方法は亜弾性モデルについてのセクションですでに使っている.

2つ目のアプローチは同じことを達成するために,物体に作用する力の量を徐々に増加させて時間積分を使う.時間積分を使う主な利点は,NDSolveには次の時間刻みを自動的に選択する適用型時間刻みメカニズムが装備されているという点である.パラメトリック力増加アプローチでは,次の刻み幅を手動で選ばなければならないが,疑似時間積分ではその必要がない.ここで0から任意に選んで終了時間の1まで方程式系を時間積分する.時間変数 はパラメータ の代りになり,物体に作用する力を徐々に増加させるために使われる.このアプローチはネオフック材料モデルについてのセクションで示す.以下の例では,パラメトリック力増加アプローチを示す.

例として単位立方体を使う.立方体を で固定する一方で,正の 方向に張力を掛ける.

材料パラメータでサンブナン・キルヒホッフモデルを指定する:
力の加わった変形 についてのパラメトリック関数を作成する:
最大変形 ,ステップ数 nsteps,初期変位ベクトルを指定する:

次では,前の解を使って現在の解のための非線形解法を開始する.

変位が増加するときの解の進行状況を監視する:

非線形解法が収束しない場合は,まず取るステップ数 nsteps を増やしてみる.ソルバが収束しないのにはさまざまな原因がある.要素数を増やすだけでも,要素数が少なければ見付けられた解が見付けられなくなる可能性がある.その場合の一つの解決策として,粗いメッシュを使って解を求め,それを初期シードとして細かいグリッド計算を行うというものがある.このアプローチは有限要素の最善の使い方である「メモリを多く必要とする偏微分方程式」で紹介してある.そこでは例に反復ソルバが使われているが,反復ソルバを指定するオプションを取り除くことによってデフォルトの線形ソルバを使うことができる.

超弾性材料モデルでは,変形が大きく可視なので,変形をスケールしたい場合が多い.この場合,変形またはもとの立体の"ScalingFactor"を1に設定することができる.

計算された変位に基づいて変形したメッシュを生成する.

以下の可視化は縮尺通りである.

変形した物体を赤で,変形していない物体をエッジフレームで可視化する:
変形した物体の体積を計算する:

体積はもとの体積より小さくなっていることが分かる.このトピックについての詳細は圧縮性のセクションで説明してある.

物体のひずみを計算する:

ひずみからコーシー応力を計算することは,線形弾性の場合とは多少異なる.線形弾性の場合,応力を復元するためにひずみを使うことができた.ここでは,構成方程式が第二ピオラ・キルヒホッフ応力とともに動作するため,コーシー応力に変換できるように変形勾配 も計算する必要がある.

物体のコーシー応力を計算する:

関心のある応力が,線形弾性材料の場合のように構成方程式で使われる応力と同じならば,そのときに限り変位は与える必要はない.

物体の第二ピオラ・キルヒホッフ応力を計算する:
物体の第一ピオラ・キルヒホッフ応力を計算する:

第一ピオラ・キルヒホッフ応力は対称ではない.したがって,使われるテンソルが対称でなければならないミーゼス応力の可視化等に使うことはできない.

コーシー応力のミーゼス応力を計算する:
ミーゼス関数を抽出する:
ミーゼス応力の最小・最大値を求める:

変形した物体上のミーゼス応力を可視化したい.

変形したメッシュ上のミーゼス応力の補間関数を作成する:
変形した物体上のミーゼス応力の等高線を表示する:

次に,サンブナン・キルヒホッフモデルで線形弾性が再生できるかどうかを調べる.微小ひずみ では,エネルギー密度関数 は以下で与えられる[14, p 13]:

は第一ラメ定数, は弾性係数 としても知られている第二定数である.ラメ定数はヤング率とポワソン比と関連付けることができる.

微小ひずみ についてのエネルギー密度関数 の導関数は以下である:

PDEモデルを設定するためにヘルパー関数を作成する:
材料パラメータを設定する:
微小ひずみ で等方性サンブナン・キルヒホッフモデルを使う:
変位を計算する:

次に線形弾性機能を使って同じモデルを計算する.応力が比較しやすいように,線形弾性材料で一般的な工学ひずみを使わない.パラメータ pars"EngineeringStrain"Falseと設定すると工学ひずみを使わないようにすることができる.これは線形材料法則と,真ひずみが使われる非線形材料法則を比較する場合におもしろいかもしれない.詳細はSolidMechanicsStrainのシンボルページに記述してある.

工学ひずみに使用を無効にする:
線形弾性モデルの解を計算する:
両方のモデルの最小・最大変位の差分を比較する:

サンブナン・キルヒホッフモデルには限界がある.このモデルは強力な圧縮に対する耐性が弱い.サンブナン・キルヒホッフの弾性物体が圧縮されるとき,最初復元力で反応する.圧縮の臨界閾値(単一軸の変形しない次元のほぼ58%)を超えると,復元力は減少する[15].つまり,サンブナン・キルヒホッフモデルは圧縮されると柔らかくなるのである.

ユーザ定義の材料法則を書く方法に関する情報を与えるために,実際の実装を見てみよう.ユーザ定義の材料法則(UMATと言われることもある)を加える例は亜弾性材料のセクションで示してある.サンブナン・キルヒホッフ材料モデルの実際の実装は以下で示す通りである.このモデルの構築方法は,他の構成モデルの構築の例として使うことができる.

サンブナン・キルヒホッフ材料モデルの実装:
MaterialModel[_, "Isotropic", "VenantKirchhoff", vars_, pars_, data__] :=           

Module[
    {X, dim, lambda, mu, EE, ee, rules, stressMatrix, strainMeasure},

    values = GetSolidMechanicsMaterialParameterList[vars, pars,
        {"LameParameter", "ShearModulus"}, <||>];
    If[ FailureQ[values], Return[ $Failed, Module]; ];
    {lambda, mu} = values;

    X = vars[[-1]];
    dim = Length[X];
    EE = Array[ee, {dim, dim}];
    stressMatrix = lambda * Tr[EE] * IdentityMatrix[dim] + 2 * mu * EE;

    (* get the strain measure *)
    strainMeasure = "StrainMeasure" /. data;
    rules = Thread[Rule[Flatten[EE], Flatten[strainMeasure]]];
    stressMatrix = stressMatrix /. rules;

    stressMatrix
]

新しい材料モデルを加える

ここまで最初の超弾性材料モデルを見てきたので,新しい材料モデルを加える方法を説明する.ここでは前に紹介したサンブナン・キルヒホッフモデルを使って,これをカスタムモデルとして取り入れる方法を示す.

材料モデルを書く:

材料モデルを作ったら,固体力学関数にその利用方法を指示する必要がある.これはパラメータ連想で適切なパラメータを指定することによって行う.

超弾性材料顧客モデルのパラメータを指定する:

"SolidMechanicsMaterialModelFunction"は新しい材料モデルを指定する.材料パラメータは,前述のサンブナン・キルヒホッフモデルの例と同じものである."StrainMeasure"には"GreenLagrange"ひずみを使う.利用できるひずみ測定についての詳細はひずみの理論セクションをご覧いただきたい.次にSolidMechanicsPDEComponentに,材料法則で指定した応力測定を告げる必要がある.材料法則に使われる応力測定は"ConstitutiveStressMeasure"である.非線形平衡方程式の応力測定は"EquilibriumStressMeasure"である.例えば,"FirstPiolaKirchhoff"応力で材料法則を指定したら,それに応じて"ConstitutiveStressMeasure"を設定するのである.デフォルトでは"OutputStressMeasure""Cauchy"応力であるが,これは変更することができる.応力と応力の間の変換は自動的に行われる.

ネオ・フックモデル

名前が示す通り,ネオ・フック材料モデルは,線形弾性のフックモデルを大きい変形に拡張したものである.圧縮性ネオ・フックモデル のエネルギー密度関数は以下で与えられる[5, p. 499]:

同じエネルギー密度関数に対する他の定式もいろいろ存在するが,通常以下のように不変量に基づいている.

ここで は第一不変量であり,であるがここでは追求しない.第二ピオラ・キルヒホッフ応力は以下で与えられる[5, p. 500]:

ユーザ定義の材料法則を書く方法に関する情報を与えるために,実際の実装を見てみよう.ユーザ定義の材料法則(UMATと言われることもある)を加える例は亜弾性材料のセクションで示してある.ネオ・フック材料法則の実際の実装はネオフックモデルのセクションで示す通りである.このモデルの構築方法は,他の構成モデルの構築の例として使うことができる.

ネオ・フックモデルは力学モデルである.つまり,これは材料の機械的原理から派生する.係数 は剛性率, は体積弾性率)が利用されることもある.

例として単位立方体を使う.立方体を で固定する一方で,正の 方向に張力を掛ける.物体を約40%引き伸ばしたい.

PDEモデルを設定するためにヘルパー関数を作成する:

ここで問題にしているモデルは圧縮性モデルである.つまり,変形した物体の体積は負荷の下で変化する可能性があるということである.圧縮性モデルは,後で説明する非圧縮性または微圧縮性のモデルと対照的である.実装された超弾性材料モデルのデフォルトは微圧縮性である.ここで説明するモデルを利用するためには,パラメータ"Compressibility"->"Compressible"を指定する必要がある.

材料パラメータで等方性ネオ・フックモデルを指定する:

まずパラメトリック力増加アプローチを試す.

ネオ・フックモデルのパラメトリック関数を設定する:
最終的な変位の1/100の変位を計算する:

境界変位での小さい最初の増加ですら収束できないことが分かる.収束解を得るために最初の変位をさらに減少させることは可能であるが,40%の伸張に向けてその解を反復すると,刻み幅が小さいままだと長時間かかる.別の方法として,PDEを時間依存PDEに変換して,初期条件と初期条件の導関数0を使ってPDEモデルを架空時間の1単位まで時間積分する.

架空時間の0から1まで時間依存PDEを解き,境界変位を徐々に増加する:

通常NDSolveやその関連関数で時間依存微分方程式の時間範囲を指定する場合,時間変数 ,開始時間,終了時間のリストとして指定する.上のNDSolveValueの呼出しでは,時間領域はで指定されている.ここで に設定されていることに注目してほしい.これは初期条件が与える時間から最終時間 まで時間積分が実行されるということであり,同時にNDSolveValueに結果の補間関数に最後の時間刻み だけを保存するように指示するものである.

後続のデータ処理をより効率的にするために,変位の時間依存補間関数を定常補間関数に変換する.この場合は,最後の時間刻みにだけしか関心がないので簡単にできる.InterpolatingFunctionオブジェクトはNDSolveが取ったさまざまな時間刻みにおけるすべてのノードの解をリストに保存する.したがって補間関数から値の最後の集合を抽出するだけである.これらの値は新しい補間関数にパッケージすることができる.

時間依存解を定常解に変換する:
変形したメッシュを生成する:
変形した物体の体積を計算する:

体積はもとの体積より大きいことが分かる.このトピックについての詳細は圧縮性のセクションに記載されている.それでも体積の増加は,全体的な体積と比べると小さい.

変形した物体を可視化する:
変形したメッシュ上で変位の補間関数を作成する:
変形したメッシュ上で総変位を可視化する:

上のグラフィックスではグラフィックスの境界ボックスが表示されているのであり,もとの変形していない立方体ではない.

物体のひずみを計算する:
物体のコーシー応力を計算する:
ミーゼス応力を計算する:
ミーゼス関数を抽出する:
ミーゼス応力の最小・最大値を求める:
変形したメッシュ上でミーゼス応力の補間関数を作成する:
変形した物体上のミーゼス応力の等高線を背後から見て表示する:

等高線は,使用された粗いメッシュに起因する数値的影響である.

ネオフックモデルは,超弾性組織の二軸引張試験アプリケーションモデルでも使われている.

圧縮性ネオフック材料モデルの実装:
MaterialModel[_, "Isotropic", "NeoHookeanCompressible", vars_, pars_, data__]:=
Module[
    {U, X, dim, lambda, mu, fp, f, idm, c, cinv, j, stressMatrix},

    values = GetSolidMechanicsMaterialParameterList[vars, pars,
        {"LameParameter", "ShearModulus"}, <||>];
    If[ FailureQ[values], Return[ $Failed, Module]; ];
    {lambda, mu} = values;

    dim = pars["EmbeddingDimension"];
    idm = IdentityMatrix[3];
    f = idm;
    f[[1 ;; dim, 1 ;; dim]] = DeformationGradient[vars, pars];

    (* in the plane strain case: f[[3,3]] = 1; *)
    If[ pars["SolidMechanicsModelForm"] == "PlaneStress",
        f[[3,3]] = Exp[-lambda/(lambda + 2 mu) Log[Det[f[[1;;dim, 1;;dim]]]]];
    ];

    j = Det[f];
    c = Transpose[f] . f;
    cinv = Inverse[c];
    stressMatrix = mu * (idm - cinv) + lambda * Log[j] * cinv;

    (* Second Piola-Kirchhoff stress *)
    Simplify[stressMatrix[[1;;dim, 1;;dim]]]
]

ひずみの不変量

ひずみエネルギー密度関数 は,変形勾配 ではなくコーシー・グリーンのひずみで表すことが多い.原理上は変形勾配 を使っても全く問題ない.しかし を使うことの欠点は,材料モデルを設計するときに直感的に理解しにくいということである.例えば,材料モデルの圧縮性は,エネルギー密度関数 を表すために変形勾配 ではなく異なる式を使ったときにうまく表せる.したがって はコーシー・グリーンのひずみ のようなひずみ測定,または変形勾配 から派生する不変量を使って表すことが一般的である[15, p. 11ff].

材料モデルは回転不変である.つまり,物体と負荷が回転してもモデルは同じ結果になるということである.超弾性材料モデルは,ひずみエネルギー密度関数が以下を満足するとき,そのときに限り回転不変である:

ここで は回転行列である.

ここで扱う2つ目の特性は,考慮する材料が等方性であるということである.材料特性がすべての方向において同じ場合,その材料は等方性と言う.超弾性材料モデルは,ひずみエネルギー密度関数が以下を満足するとき,そのときに限り等方性である:

ここで は回転行列である.その結果,回転不変であり等方性でもある材料は以下を満足する:

(任意の回転 および のとき)

特異値分解 を使うと,回転不変であり等方性でもある材料は以下を満足する:

ここで の3つの特異値 である. は主ストレッチと言われる.ひずみエネルギー密度関数 の特異値で定義することもできるが,特異値分解は計算的に高価なので,等方性材料には行われない.その代りに,コーシー・グリーンひずみテンソル の3つの等方性不変量 を使うというものがある.これらは特異値と同等であるが計算がしやすいのである.等方性不変量は以下のように定義できる:

ここで はそれぞれ第一,第二,第三の等方性不変量であり, は右コーシー・グリーンの変形テンソルである.文献では の定義方法が少なくとも2つある:

不変量は, の特異値である主ストレッチ を使って表すこともできる:

特別な場合も記しておきたい.ひずみがない場合,すべてのストレッチは となり,不変量はそれぞれ定数の となる.

要するに,ひずみエネルギー密度関数 はコーシー・グリーンひずみとひずみ不変量 のどちらで表すこともできる:

不変量は,回転に依存せず計算しやすいひずみのユニークな測定を与える.

実際,右コーシー・グリーンひずみ に関する不変量の導関数を知っておくのも便利である[16, p. 216]:

ここで は単位行列を表す.

このことから,不変量 に基づく第二ピオラ・キルヒホッフ応力 を構築することができる.この場合,等方性不変量 を使ってひずみエネルギー密度関数 を表す.第二ピオラ・キルヒホッフ応力 は以下の連鎖法則を使って計算できる:

不変量の導関数を使うと式を以下のように簡約できる:

圧縮性

材料の記述における重要な要素の一つとして,その材料が変形したときの体積挙動がある.通常は,圧縮性,非圧縮性,微圧縮性のいずれかの材料モデルに分類する.圧縮性材料では,変形した物体の体積は変化してもよい.非圧縮性材料では,変形の下での体積変化はゼロである.この非圧縮性は非圧縮性についての数学定義で表され,で与えられる.

ある意味,非圧縮性材料モデルは,であることを要求することによって非圧縮性が強制される,ポワソン比が の線形弾性材料と同等である.ポワソン比 は,荷重方向と垂直の方向における材料挙動を表す.ポワソン比は のようなストレッチによる体積の相対変化を表すことができる.体積変化がゼロである非圧縮性材料では,ポワソン比は でなければならない.次にはこれが無限体積弾性率 を導き,モデルを解くときに非常に硬くなる.硬いモデルを解く数値的困難を避けるために,微圧縮性モデルが導入された.

非圧縮性材料をモデル化するために,ひずみエネルギー密度関数 を定積(一定の体積)部分 と体積部分 に分ける:

ここで である.定積部分 は体積の変化のない変形を表す.体積部分 は非圧縮性をモデル化する.実際,体積部分は非圧縮性を強制するために使われる.は体積比 にだけ依存するように構築されている.を強制するモデルは非圧縮性である.

として与えられる[16, p. 230].しかし,非圧縮性のために を設定したので, となる. はひずみエネルギー密度関数 が定積変形にのみ依存することを示すために使われる.

完全非圧縮性材料の場合,体積部分 は定義されないが,定数を持つひずみエネルギーで仮定される.非圧縮性を強制するために,静水圧 を導入する:

以下の連鎖法則

を使って を微分すると が与えられる.次に を使うと,以下で表すことのできる第二ピオラ・キルヒホッフ応力 になる:

上の微分では定積部分は について微分され,2項目のみが に依存する.

を強制するとモデルの硬さが大きくなるため,数値的に問題となり得る.さらに,追加の自由度が必要になる.したがって厳密に を満足するのではなく,その条件を少し緩めて を目標にする.これが微圧縮性モデルと呼ばれるもので,まずまずの妥協点と言える.これで体積部分は以下で示される:

非圧縮性条件 を強制するために,静水圧 を以下で表すことができる:

ここで は材料の体積弾性率である. の実験値が利用できない場合,によって剛性率 に関連付けることができる. はスケール因子であり,通常の範囲にある.

微圧縮性材料モデルの体積弾性率の実装:
k = EstimateBulkModulus[vars, pars, <|s -> 10^4, "ShearModulus" -> mu|>];
p = - k * (j - 1);

微圧縮性材料モデルについての第二ピオラ・キルヒホッフ応力は以下である:

超弾性モデルのコレクション

ここで示すモデルは,現象論的モデル,機構論的モデル,ハイブリッドモデルに分類できる.材料の現象論的記述とはモデルが実験データになるべく一致するように作成されているが物理的根拠は持たないということである.これに対して機構論的モデルは,物理的過程に基づいており,実験データに近づけようとするものである.ハイブリッドモデルはこの両方の面からなる.

現象論的モデル

機構論的モデル

ハイブリッドモデル

ここで提示する超弾性モデルにはすべてモデルパラメータがある.剛性率のように,中には"ShearModulus"のような文字列のパラメータとして指定されるものもある.名前を持たないとか方程式の一部とかの他のパラメータは 等の形式的な記号パラメータとして与えられる.表記的アルファベット文字についてのガイドページを見ると形式文字の作成方法が分かる.ワークフローに関して言うと形式文字は以下の例からコピーすることができる.

超弾性材料モデルの比較は「超弾性モデルの比較」の応用例で見ることができる.

ムーニー・リブリンモデル

一般化されたリブリン超弾性モデルは顕著な非線形変形を見せるゴムのような材料の機械的動作を記述するために広く使われている.このモデルは,大きなひずみを受けたり,等方性,均質性,微圧縮性のような特性を持つ材料に適している.これは生体力学,軟組織力学,エラストマー工学を含む分野でよく使われ,現象論的モデルと見なされている.

一般化されたリブリンモデル[18]のひずみエネルギー密度関数 は以下で与えられる:

ここで は適用された応力の等容反応に関連した材料係数を表す.係数 となるように常に0に設定される.

方程式(1)の最後の項は圧縮性ムーニー・リブリンモデルに関連している.ここで は応力適用による体積応答に関連する.Wolfram言語の現行バージョンは圧縮性のバージョンを提供しないので に設定される.

モデルを線形弾性と一致させるために,体積弾性率を に設定し,剛性率を に設定することも必要である.

材料係数 の指数 の特別な組合せはよく知られるモデルに導かれる.例えば,のいずれか2つのパラメータを持つムーニー・リブリンモデルはネオフックモデルになる.最もよく使われるムーニー・リブリンモデルは2個,5個,9個のパラメータを持つタイプである.

非圧縮性ムーニー・リブリンモデル

方程式(2)のムーニー・リブリンモデルのひずみエネルギーは最初と2番目のひずみ不変量の線形結合である.これはメルビンムーニー[19]によって提唱されてから,不変量についてロナルドリブリン[18]によって形式化された.最もよく使われるのはパラメータを2個,5個,9個持つムーニー・リブリンモデルであり,使われない係数はゼロに等しい.パラメータを9個持つひずみエネルギー関数は以下になる:

5個のパラメータを持つモデルでは係数 のみが使われる.ひずみエネルギーは以下になる:

2個のパラメータの場合は係数 のみが使われる.ひずみエネルギーは以下になる:

材料パラメータは通常 および 等の名前が付けられる.

等容部分は,等容ひずみ不変量 のみに依存する.非圧縮性制約条件 のため,最初の等容不変量は になる.

等容部分は以下のように応力に寄与している:

ここで である.

第二ピオラ・キルヒホッフ応力 は以下のように表される:

完全非圧縮性ムーニー・リブリンモデルは,現在Wolfram言語に実装されていない.

微圧縮性ムーニー・リブリンモデル

9パラメータの微圧縮性ムーニー・リブリンモデルの第二ピオラ・キルヒホッフ応力 は以下のようになる:

2パラメータと5パラメータの場合も,ピオラ・キルヒホッフ応力は同様の構造になる.使われない係数はゼロである.

ムーニー・リブリンモデルのパラメータ によって剛性率 と関連付けられる.

実験的値が利用できない場合に体積弾性率 を設定するために, によって剛性率 と関連付けることができる. はスケール因子であり通常の範囲である.つまり,剛性率を使って を調整するのである.

モデルパラメータ s はオプショナルでありデフォルトのになる.パラメータ"BulkModulus"が与えられると はその値に設定される.

次に微圧縮性ムーニー・リブリンモデルの使い方の例を挙げる.

モデルパラメータを持つ,等方性の2パラメータのムーニー・リブリンモデルを指定する:
モデルパラメータを持つ,等方性の5パラメータのムーニー・リブリンモデルを指定する:
モデルパラメータを持つ,等方性の9パラメータのムーニー・リブリンモデルを指定する:
強制荷重 についてのパラメトリック関数を作成する:
最大荷重 ,ステップ数 nsteps,初期変位ベクトルを指定する:

前の解を使って現在の解のための非線形解法を開始する.

変位が増加するときの解の進行状況を監視する:
計算した変位に基づいて変形メッシュを生成する:

以下の可視化はスケールされている.

変形した物体を赤で,変形していない物体をエッジフレームで可視化する:
変形した物体の体積を計算する:

微圧縮性のムーニー・リブリン材料モデルの実際の実装は以下の通りである.このモデルの構築方法は,他の構造モデルの構築の例として使うことができる.

ムーニー・リブリン材料モデルの実装:
MaterialModel[_, "Isotropic", "MooneyRivlin", vars_, parsIn_, data__] :=
Module[
    {dim, f, idm, c, cinv, j, pars, stressMatrix, p, i1, i1iso, C10, C01, C11,         C20, C02, C30, C03, C21, C12, k, ciso, mu, s, i2iso, PK2Iso, idg, temp,         idgInverse, idgInverseDot, idgJacobian},

    pars = parsIn;
    If[ KeyExistsQ[pars, Subscript[c, 1]],
        pars[Subscript[c, 1, 0]] = pars[Subscript[c, 1]];
    ];
    If[ KeyExistsQ[pars, Subscript[c, 2]],
        pars[Subscript[c, 0, 1]] = pars[Subscript[c, 2]];
    ];

    values = GetSolidMechanicsMaterialParameterList[vars, pars, {
        Subscript[c, 1, 0], Subscript[c, 0, 1], Subscript[c, 1, 1],
        Subscript[c, 2, 0], Subscript[c, 0, 2],
        Subscript[c, 3, 0], Subscript[c, 0, 3],
        Subscript[c, 2, 1], Subscript[c, 1, 2]
        }, <|
        Subscript[c, 1, 1] -> 0, Subscript[c, 2, 0] -> 0, Subscript[c, 0, 2] -> 0,         Subscript[c, 3, 0] -> 0, Subscript[c, 0, 3] -> 0, Subscript[c, 2, 1] -> 0,         Subscript[c, 1, 2] -> 0
        |>
    ];

    If[ FailureQ[values], Return[ $Failed, Module]; ];
    {C10, C01, C11, C20, C02, C30, C03, C21, C12} = values;

    dim = pars["EmbeddingDimension"];
    idm = IdentityMatrix[3];
    f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];
    idg = pars["InelasticDeformationGradient"];
    temp = Through[idg[vars, pars]];
    temp = ExtendedDeformationGradient /@ temp;

    idgInverse = Inverse /@ temp;
    idgInverseDot = Dot @@ idgInverse;
    idgJacobian = (Times @@ Det /@ temp);
    f = f . idgInverseDot;
    j = Det[f];
    c = Transpose[f] . f;
    cinv = Inverse[c];
    i1 = Tr[c];
    i1iso = i1 * j^(-2/3);
    ciso = c * j^(-2/3);
    i2iso = (Tr[ciso]^2 - Tr[ciso . ciso])/2;

    PK2Iso = 2 ((C10 + 2 C20 (i1iso - 3) + C11 (i2iso - 3) +
        2 C30 (i1iso - 3)^2 + 2 C21 (i1iso - 3) (i2iso - 3) +
        C12 (i2iso - 3)^2 + i1iso (C01 + C11 (i1iso - 3) +
        2 C02 (i2iso - 3) C21 (i1iso - 3)^2 + 2 C12 (i1iso - 3) (i2iso - 3) +
        2 C03 (i2iso - 3)^2)) idm - (C01 + C11 (i1iso - 3) +
        2 C02 (i2iso - 3)) ciso);

    mu = 2 (C10 + C01);
    k = EstimateBulkModulus[vars, pars, <|s -> 10^3, "ShearModulus" -> mu|>];
    p = - k (j - 1);

    If[pars["SolidMechanicsModelForm"] == "PlaneStress",
        j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));
        p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);
    ];

    stressMatrix = PK2Iso - p * j * cinv;

    (* pull-back *)
    stressMatrix = idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

    (* returns second Piola-Kirchhoff stress *)
    stressMatrix[[1 ;; dim, 1 ;; dim]]
]

ヨーモデル

ヨーモデルはゴムのような材料の弾性特性の現象論的記述を表す.つまり,このモデルは観察された挙動を記述するのである.このモデルは微圧縮性材料および非圧縮性材料に使われる.Wolfram言語に実装されているヨーモデルは,微圧縮性モデルである.

リブリンのゴム弾性理論[16, p. 238]から,ひずみエネルギー関数をひずみ不変量の無限ベキ級数として表すことが可能である.ヨーは であると仮定し,第一ひずみ不変量 の立方体であると提唱した.この仮定は がいかにゼロに近いかということを示すいくつかの実験に基づいている[1, p. 245].

ヨーモデルのひずみエネルギー密度関数 は以下で与えられる:

ここで として与えられる第一ひずみ不変量である.

ヨーモデルは という仮定から導かれているので,このモデルには圧縮性のバージョンはない.つまり であり,非圧縮性をモデル化する.

非圧縮性ヨーモデル

ヨーモデルの定積部分は以下である:

定積部分は,定積ひずみ不変量 だけに依存する.非圧縮条件 により,第一定積不変量は である.

定積部分は,以下のように応力に寄与する:

ここで である.

第二ピオラ・キルヒホッフ応力 は以下で表される:

完全非圧縮ヨーモデルは現在Wolfram言語に実装されていない.

微圧縮性ヨーモデル

微圧縮性ヨーモデルの第二ピオラ・キルヒホッフ応力 は以下となる:

第一ひずみ不変量の定積部分は で表される.

ヨーモデルを利用するためには,のモデルパラメータを指定する必要がある.あとは の値を指定することである.

ヨーモデルパラメータ によって剛性率 に関連付けられる.体積弾性率 を設定するのに,利用できる実験値がない場合, によって剛性率 に関連付けることができる. はスケール因子であり,通常の範囲である.つまり,剛性率を使って を調整するのである.

モデルパラメータsはオプショナルであり,デフォルトのになる.パラメータ"ShearModulus"が設定されており,同時にモデルパラメータc1が与えられていない場合,に設定される.パラメータ"BulkModulus"が与えられると, はその値に設定される.

次に微圧縮性ヨーモデルの使い方の例を示す.

モデルパラメータを使って等方性ヨーモデルを指定する:
強制荷重 に対するパラメトリック関数を作成する:
最大荷重 ,ステップ数 nsteps,初期変位ベクトルを指定する:

前の解を使って現在の解のための非線形解法を開始する.

変位が増加する間の解の進捗状況を監視する:
計算された変位に基づいて変形されたメッシュを生成する:

次の可視化は原寸に比例している.

変形した物体を赤で,変形していない物体を辺のフレームとして可視化する:
変形した物体の体積を計算する:

微圧縮性ヨー材料モデルの実際の実装は以下の通りである.このモデルの構築方法は,他の構成モデルを構築する際の例として使うことができる.

ヨーの構成モデルを使った層状血管」の応用例ではヨー材料モデルを利用している.

ヨー材料モデルの実装:
MaterialModel[_, "Isotropic", "Yeoh", vars_, parsIn_, data__] :=
Module[
    {pars, dim, f, idm, c, cinv, j, stressMatrix, p, i1, i1iso, c1, c2, c3, mu, k,
    PK2Iso, idg, temp, idgInverse, idgInverseDot, idgJacobian},

    pars = parsIn;
    If[ KeyExistsQ[pars, "ShearModulus"],
        pars[Subscript[c, 1]] = pars["ShearModulus"]/2;
    ];

    values = GetSolidMechanicsMaterialParameterList[vars, pars, {
        Subscript[c, 1], Subscript[c, 2], Subscript[c, 3]}, <||>];

    If[ FailureQ[values], Return[ $Failed, Module]; ];
    {c1, c2, c3} = values;

    dim = pars["EmbeddingDimension"];
    idm = IdentityMatrix[3];
    f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];
    idg = pars["InelasticDeformationGradient"];
    temp = Through[idg[vars, pars]];
    temp = ExtendedDeformationGradient /@ temp;

    idgInverse = Inverse /@ temp;
    idgInverseDot = Dot @@ idgInverse;
    idgJacobian = (Times @@ Det /@ temp);
    f = f . idgInverseDot;
    j = Det[f];
    c = Transpose[f] . f;
    cinv = Inverse[c];
    i1 = Tr[c];
    i1iso = i1 * j^(-2/3);

    PK2Iso = 2 * (c1 + 2*c2*(i1iso - 3) + 3*c3*(i1iso - 3)^2)*IdentityMatrix[3];
    mu = 2 * c1;

    k = EstimateBulkModulus[vars, pars, <|s -> 10^4, "ShearModulus" -> mu|>];
    p = - k * (j - 1);

    If[pars["SolidMechanicsModelForm"] == "PlaneStress",
        j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));
        p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);
    ];

    stressMatrix = PK2Iso - p * j * cinv;

    (* pull-back *)
    stressMatrix =idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

    (* returns second Piola-Kirchhoff stress *)
    Simplify[stressMatrix[[1 ;; dim, 1 ;; dim]]]
]

ネオフックモデル

ネオフックモデルは線形弾性の拡張として前のセクションで紹介した.次では非圧縮性と微圧縮性のバージョンを説明する.またネオフックモデルはムーニー・リブリンモデルの特殊形であるが,機構論的モデルと考えられる.

非圧縮性ネオフックモデル

非圧縮性ネオフックモデルのひずみエネルギー関数にはネオフックモデルの等容項しか含まれない.

等容部分は,等容ひずみ不変量 のみに依存する.非圧縮性制約条件 があるため,最初の等容不変量は である.

等容部分は以下のように応力に寄与している:

ここで である.

第二ピオラ・ キルヒホッフ応力 は以下となる:

完全非圧縮性ネオフックモデルは現在Wolfram言語に実装されていない.

微圧縮性ネオフックモデル

微圧縮性ネオフックモデルの第二ピオラ・キルヒホッフ応力 は以下となる:

ネオフックモデルは圧縮性形式および微圧縮性形式で利用できるが,デフォルトは微圧縮性形式である.圧縮性ネオフックモデルを使うためには,モデルパラメータ μ とパラメータ "Compressibility"->"NearlyIncompressibile"を指定することが必要である.

体積弾性率 を指定するのに,利用できる実験値がない場合, によって剛性率 に関連付けることができる. はスケール因子であり,通常の範囲である.つまり,剛性率を使って を調整するのである.パラメータ"BulkModulus"が与えられると はその値に設定される.

次は微圧縮性ネオフックモデルの使い方の例である.

材料パラメータを持つ等方性微圧縮性ネオフックモデルを指定する:
ネオフックモデルのパラメータ関数を設定する:
最大荷重 ,ステップ数 nsteps,初期変位ベクトルを指定する:
変位が増加するときの解の進行状況を監視する:
計算した変位に基づいて変形メッシュを生成する:
変形した物体を赤で,変形していない物体をエッジフレームで可視化する:
変形した物体の体積を計算する:
微圧縮性ネオフック材料モデルの実装:
MaterialModel[_, "Isotropic", "NeoHookeanNearlyIncompressible", vars_, pars_, data__] :=
Module[
    {U, X, dim, lambda, mu, fp, f, idm, c, cinv, j, PK2Iso, stressMatrix, idg,
    temp, idgInverse, idgInverseDot},

    If[ KeyExistsQ[pars, Subscript[c, 1]],
        mu = 2 * pars[Subscript[c, 1]];
        pars["ShearModulus"] = mu;
    ,(* else *)
        mu = GetSolidMechanicsMaterialParameter[vars, pars, "ShearModulus"];
    ];
    If[ FailureQ[mu], Return[ $Failed, Module]; ];

    dim = pars["EmbeddingDimension"];
    idm = IdentityMatrix[3];
    f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];
    idg = pars["InelasticDeformationGradient"];
    temp = Through[idg[vars, pars]];
    temp = ExtendedDeformationGradient /@ temp;

    idgInverse = Inverse /@ temp;
    idgInverseDot = Dot @@ idgInverse;
    idgJacobian = (Times @@ Det /@ temp);

    f = f . idgInverseDot;
    j = Det[f];
    c = Transpose[f] . f;
    cinv = Inverse[c];
    PK2Iso = mu * idm;

    k = EstimateBulkModulus[vars, pars, <|s -> 10^2|>];
    p = - k (j - 1);

    If[pars["SolidMechanicsModelForm"] == "PlaneStress",
        j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));
        p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);
    ];

    stressMatrix = PK2Iso - p * j * cinv;

    (* pull-back *)
    stressMatrix = idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

    (* Second Piola-Kirchhoff stress *)
    Simplify[stressMatrix[[1;;dim, 1;;dim]]]
]

Arruda-Boyceモデル

Arruda-Boyceモデルは高分子物質を記述するために使われる機構論的モデルである.このモデルは 個のポリマー鎖を含む代表的な基本セルの統計力学に基づいている.例えば直方体の基本セルは,中心からそれぞれの角に1個ずつの合計8個のポリマー鎖を持つ.これらのセルは有限要素メッシュには関係していないが,モデルの理論的導出のための基本をなす.

361.gif

直方体基本セルにおけるArruda-Boyceモデルの8鎖ネットワーク.

Arruda-Boyce [20]モデルのひずみエネルギー密度関数 は以下で与えられる:

ここで はそれぞれ単鎖のセグメント数とネットワークのポリマー鎖数を表し, はボルツマン定数, は絶対温度(単位は[K])を表す.鎖の伸長 は次のように最初のひずみ不変量に関連している:

ランジュバン関数 は以下で与えられる:

-1 は逆ランジュバン関数を表し,以下のように鎖の伸長 の関数である:

ランジュバン統計を使うことで,鎖のエントロピーに基づく鎖伸び切り限界が説明できる.

逆ランジュバン関数を計算することは,計算的に高価である.ひずみエネルギー関数を計算的により効率的にしたバージョンは,逆ランジュバン関数を打ち切られたテイラー級数[20]の5項で近似することで作成することができる.これでひずみエネルギー密度関数は次で与えられる:

このバージョンのひずみエネルギー密度関数も,最初のひずみ不変量の関数である方程式(3)で を置換している.

このモデルには,2つの材料係数が必要である.一つは材料係数 ,もう一つは以下のように鎖ネットワークのロッキング伸長に関連した鎖係数 β である:

はテイラー展開からの数値係数である:

非圧縮性Arruda-Boyceモデル

ひずみエネルギー関数は次のようになる:

等容部分は,等容ひずみ不変量 のみに依存する.非圧縮性制約条件 があるため,最初の等容不変量は である.

等容部分は以下のように応力に寄与している:

ここで である.

不変量のセクションで説明した通り,線形弾性との一貫性のために,初期剛性率は以下を導く条件 を満足する必要がある:

添字表記 は,初期状態においてひずみがゼロの場合張力 を持つことを表している.このことから となることが分かる.

第二ピオラ・キルヒホッフ応力 は次のように表せる:

完全非圧縮Arruda-Boyceモデルは現在Wolfram言語に実装されていない.

微圧縮性Arruda-Boyceモデル

微圧縮性Arruda-Boyceモデルの第二ピオラ・キルヒホッフ応力 は以下である:

最初のひずみ不変量の等容部分は で表される.

Arruda-Boyceモデルを利用するためには,2つの材料係数が必要である.一つは材料係数 ,もう一つは のように鎖ネットワークのロッキング伸長に関連した鎖係数 である.

次に微圧縮性Arruda-Boyceモデルの使い方の例を示す.

モデルパラメータを使って等方性Arruda-Boyceモデルを指定する:
強制荷重 に対するパラメトリック関数を作成する:
最大荷重 ,ステップ数 nsteps,初期変位ベクトルを指定する:

次では,前の解を使って現在の解のための非線形解法を開始する.

変位が増加するときの解の進行状況を監視する:
計算した変位に基づいて変形メッシュを生成する:

以下の可視化はスケールされている.

変形した物体を赤で,変形していない物体をエッジフレームで可視化する:
変形した物体の体積を計算する:

微圧縮性Arruda-Boyce材料モデルの実際の実装は以下の通りである.このモデルの構築方法は,他の構成モデルを構築する際の例として使うことができる.

Arruda-Boyce材料モデルの実装:
MaterialModel[_, "Isotropic", "ArrudaBoyce", vars_, pars_, data__] :=
Module[
    {dim, f, idm, c, cinv, j, stressMatrix, p, i1, i1iso, c1, lambdam, k, mu, s,     alpha, beta, i1isopower, pp, PK2Iso, idg, temp, idgInverse, idgInverseDot,
    idgJacobian},

    values = GetSolidMechanicsMaterialParameterList[vars, pars,
        {Subscript[c, 1],Subscript[λ, m]}, <||>];
    If[ FailureQ[values], Return[ $Failed, Module]; ];
    {c1, lambdam} = values;

    dim = pars["EmbeddingDimension"];
    idm = IdentityMatrix[3];
    f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];
    idg = pars["InelasticDeformationGradient"];
    temp = Through[idg[vars, pars]];
    temp = ExtendedDeformationGradient /@ temp;

    idgInverse = Inverse /@ temp;
    idgInverseDot = Dot @@ idgInverse;
    idgJacobian = (Times @@ Det /@ temp);

    f = f . idgInverseDot;
    j = Det[f];
    c = Transpose[f] . f;
    cinv = Inverse[c];
    i1 = Tr[c];
    i1iso = i1 *j^(-2/3);

    alpha = {1/2, 1/20, 11/1050, 19/7000, 519/673750};
    pp = {0, 1, 2, 3, 4};
    beta = Power[1/lambdam^2, pp];
    i1isopower = Power[i1iso, pp];

    PK2Iso = 2 * c1 * ((alpha * beta) . i1isopower) * idm;

mu = 2 c1 ((pp + 1) * Power[3, pp] * alpha) . beta;
    k = EstimateBulkModulus[vars, pars, <|s -> 10^2, "ShearModulus" -> mu|>];
    p = - k (j - 1);

    If[pars["SolidMechanicsModelForm"] == "PlaneStress",
        j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));
        p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);
    ];

    stressMatrix = PK2Iso - p * j * cinv;

    (* pull-back *)
    stressMatrix =idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

    stressMatrix[[1 ;; dim, 1 ;; dim]]
]

Gentモデル

Gent超弾性モデルは,材料が限界状態に達する点 の最大値が存在するという仮定に基づくハイブリッドモデルである.この状態は で表され,分子鎖の最大伸長率に関連すると仮定される.

Gent [21]モデルのひずみエネルギー密度関数 は以下で与えられる:

ここで は最初のひずみ不変量である.Gentモデルでは,剛性率 と極限値 の2つの材料パラメータが必要である.極限値は で与えられる.ここで はロッキング伸長と言われる.に近づくにつれて特異性に近づくため はロッキング伸長と呼ばれるのである.

非圧縮性Gentモデル

ひずみエネルギー関数は以下である:

等容部分は,等容ひずみ不変量 のみに依存する.非圧縮性制約条件 があるため,最初の等容不変量は である.

等容部分は以下のように応力に寄与している:

ここで である.

第二ピオラ・キルヒホッフ応力 は以下で表される:

完全非圧縮Gentモデルは現在Wolfram言語に実装されていない.

微圧縮性Gentモデル

微圧縮性Gentモデルの第二ピオラ・キルヒホッフ応力 は以下となる:

最初のひずみ不変量の等容部分は で表される.

Gentモデルを利用するためには,モデルパラメータ を指定する必要がある.

次に微圧縮性Gentデルの使い方の例を示す.

モデルパラメータを使って等方性Gentモデルを指定する:
強制荷重 に対するパラメトリック関数を作成する:
最大荷重 ,ステップ数 nsteps,初期変位ベクトルを指定する:

次では,前の解を使って現在の解のための非線形解法を開始する.

変位が増加するときの解の進行状況を監視する:
計算した変位に基づいて変形メッシュを生成する:

以下の可視化はスケールされている.

変形した物体を赤で,変形していない物体をエッジフレームで可視化する:
変形した物体の体積を計算する:

微圧縮性Gent材料モデルの実際の実装は以下の通りである.このモデルの構築方法は,他の構成モデルを構築する際の例として使うことができる.

Gent材料モデルの実装:
MaterialModel[_, "Isotropic", "Gent", vars_, pars_, data__] :=
Module[
    {dim, f, idm, c, cinv, j, stressMatrix, p, i1, i1iso, Jm, mu, k, PK2Iso,
    idg, temp, idgInverse, idgInverseDot, idgJacobian},

    values = GetSolidMechanicsMaterialParameterList[vars, pars,
        {"ShearModulus", Subscript[J, m]}, <||>];
    If[ FailureQ[values], Return[ $Failed, Module]; ];
    {mu, Jm} = values;

    dim = pars["EmbeddingDimension"];
    idm = IdentityMatrix[3];
    f = ExtendedDeformationGradient[DeformationGradient[vars, pars] ];
    idg = pars["InelasticDeformationGradient"];
    temp = Through[idg[vars, pars]];
    temp = ExtendedDeformationGradient /@ temp;

    idgInverse = Inverse /@ temp;
    idgInverseDot = Dot @@ idgInverse;
    idgJacobian = (Times @@ Det /@ temp);

    f = f . idgInverseDot;
    j = Det[f];
    c = Transpose[f] . f;
    cinv = Inverse[c];
    i1 = Tr[c];
    i1iso = i1 *j^(-2/3);
    PK2Iso = (mu * Jm/(Jm - i1iso + 3)) idm;

    k = EstimateBulkModulus[vars, pars, <|s -> 10^2, "ShearModulus" -> mu|>];
    p = - k (j - 1);

If[pars["SolidMechanicsModelForm"] == "PlaneStress",
        j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));
        p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);
    ];

    stressMatrix = PK2Iso - p*j *cinv;

    (* pull-back *)
    stressMatrix =idgJacobian*idgInverseDot.stressMatrix.Transpose[idgInverseDot];

    stressMatrix[[1 ;; dim, 1 ;; dim]]
]

平面ひずみモデルと平面応力モデル

線形弾性の場合,平面ひずみと平面応力の形式がある.平面ひずみモデルと平面応力モデルの形式は超弾性モデルに変更することもできる.線形弾性平面モデルには超弾性の場合にも該当する注意事項がある.詳細は「線形弾性」セクションをご覧いただきたい.

平面ひずみ

平面ひずみモデルの最初の仮定は, 方向には変位がないということである:

この仮定は 方向を含むすべてのひずみ成分はゼロであることを示す.変形勾配 の定義を見る:

次が得られる:

還元された二次元変形勾配 は次のように表すことができる:

変形勾配 の完全なものと還元されたものを記号的に生成する:

代数的にはこれらは同じである.例えば以下になる:

変形勾配 の完全なものと還元されたものの行列式を調べる:
変形勾配 の完全なものと還元されたものの逆行列を調べる:

この類似性は,左・右コーシーテンソル等の他のコンストラクトにも該当する.つまり,平面ひずみでは以下の形式を持つように変形勾配を設定することだけが必要なのである:

3Dの場合と同様に計算が実行され,最後には2x2の部分成分が抽出される.

平面応力

平面応力の場合,物事はもっと複雑になる.平面応力状態は3つの主応力の一つがゼロという特徴がある.任意であるがゼロ応力の方向に垂直な2つの座標軸を選ぶことによって,応力状態は以下の構造を持つ:

応力は還元された2x2行列 で表すことができる.ゼロ応力が 方向ならば応力テンソルは以下になる:

複雑なのは,平面応力条件を満足することである.これを行う一般的な方法は,厚さの変化として,平面外の変形を介して平面応力条件を満足することである.このために変形勾配に 成分を導入して の形式を以下にする:

運動学から 成分を決定することは不可能である.これは問題の構造モデル[28]から決定しなければならない.つまり,すべての超弾性モデルについてこれを一般的に行える方法はないのである.次では圧縮性モデルと微圧縮性モデルの導出を見る.

圧縮性ネオフックモデルの平面応力

[5, p.502]に従って,三次元的関係から始める:

以下のように,必要条件 を課す:

平面応力状態では厚さは小さいので, 項は1に近い.これによって 関連項をテイラー級数で表すことができる:

簡約化を使う:

最終結果:

微圧縮性モデルの平面応力

微圧縮性モデルの場合は異なるアプローチを使う.幸いにもこのアプローチはすべての微圧縮性モデルに使うことができる.

等容部分と体積部分に分解した第二ピオラ・キルヒホッフ応力から始める:

等容項は以下で与えられる:

微圧縮性材料の体積部分は以下で与えられる:

ここでの目的は,平面応力状態 が満足するような および の式を求めることである.

ここで平面応力条件 を課すために,成分となるように0である必要がある.これにより以下が導ける:

ゆえに は以下のように表すことができる:

ここで非圧縮制約条件が暗示的に考慮される.

の形式に基づくと,コーシーグリーンテンソル は以下の形式を持つ:

の対称性のため,の4つの成分だけを求める必要がある.さらに非圧縮制約条件によって,成分を以下のように残りの成分の関数として表すことができる:

以下が成り立つ:

ここで:

非圧縮制約条件

および

のため,以下が得られる:

したがって次のようになる:

ここで の成分で表す:

等方性材料ではでなければならない.は利用可能である. は次のように表せる:

ヤコビアン は変形勾配の還元形式で以下のように修正される:

修正された を使うと,第二ピオラ・キルヒホッフ応力は平面応力で表すことができる:

微圧縮性モデルのコードでは,これは以下で表せる:

微圧縮性超弾性モデルについての平面応力の実装:
If[pars["SolidMechanicsModelForm"] == "PlaneStress",
    j = Det[f] (1/(f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]]));
    p = PK2Iso[[1, 1]] /((f[[1, 1]]*f[[2, 2]] - f[[1, 2]]*f[[2, 1]])^2);
];

平面ひずみと平面応力の比較

平面応力または平面ひずみのどちらのモデルを選ぶかは,パラメータ pars の"SolidMechanicsModelForm"連想キーで設定することができる.デフォルトの2D固体力学モデル形式は"PlaneStress"である.

次の例では,左側が固定され右側が引っ張られる正方形領域についての平面応力と平面ひずみの比較を示す.

形状の2D変数とメッシュを設定する:
材料特性とモデルパラメータを指定する:

線形弾性ひずみ」のセクションで述べたように,平面ひずみモデルの高さは物体の縦・横と比較して必ずしも大きい必要はない.唯一の必要条件は 方向の終点が 平面の制約条件ではないということである.

平面応力モデルの形式を指定する:
平面ひずみモデルの形式を指定する:
最大荷重 100000の2Dモデルを生成する関数を作成する:
平面応力の例を解く:
平面ひずみの例を解く:
計算した変位に基づいて変形したメッシュを作成する:
平面応力について変形した物体を赤で,平面ひずみについては青で,変形していない物体を輪郭線として可視化する:
軸上で における水平変位を比較する:

変形したメッシュは2つの異なるモデルについて異なる変位を示している.平面ひずみのメッシュは, 方向の変形に対する暗示的な制約条件により,伸長方向に沿って減少した変位を示す平面応力よりも硬く見える.

超弾性モデルのキャリブレーション

超弾性材料モデルを使うためには,モデルパラメータを文献から取るか実験データにフィットするかする必要がある.適切な超弾性モデルを選ぶためには,実験データが重要である.一軸引張,圧縮,せん断等の力学的試験を行ってさまざまな荷重条件下での材料の応力ひずみ応答を測定する.これらの試験から得られたデータは,特定の超弾性材料モデルに対する適切なモデルパラメータを求めるために使うことができる.最後のステップで,ある超弾性材料モデルが実際の実験データをどの程度うまく表しているかを評価することができる.これらのステップはすべてWolfram言語で行うことができ,それについてはこのセクションで説明する.

超弾性組織の二軸引張試験」の応用モデルも利用できる.

一軸試験のデータフィット

一軸引張試験は,比較的簡単で制御された実験設定のため,軟組織やゴム状の材料の試験のために広く使われている.この試験では標本は単一方法への制御された引張の対象となり,正確な測定と解析が可能となる.一軸荷重条件は応力とひずみが主に一つの軸に沿って適用されることを確実にし,力学的挙動の解析を簡単にする.この試験は,荷重方向に沿って軟組織が構造的に一様で等方であるという仮定に基づいている.この仮定は現実では完璧に真だとは言えないかもしれないが,多くの軟組織のタイプに対して妥当な近似を提供する.

変位データは,試験機械のグリップ間の力の関数として記録される.一般にいくつかの曲線上のモデルをフィットするためには,複数の測定サイクルが繰り返される.

507.gif

上の図は,2つのクランプが標本を引っ張る一軸引張試験の典型的な実験設定を示している.

力変位データは,応力ひずみ曲線を計算するために使われる.応力は力 および標本の横断面積 に関連する:

一方,ひずみは標本の長さ の変動に関連する:

軟質材料では,一般に変位が大きいため,ひずみではなく伸長が使われる:

応力ひずみ曲線は によって応力ひずみ曲線に関連している.実験データから,応力ひずみ,より正確に言うと応力伸長曲線が得られる.

次に一軸引張試験の実験データから構成モデルの調整までの進め方を説明する一般的な手順を示す.

超弾性モデルの場合,材料挙動は第二ピオラ・キルヒホッフ応力 はコーシー・グリーンひずみテンソル)を持つひずみエネルギー で表される.コーシー応力は抵抗関係 は変形勾配)を介して第二ピオラ・キルヒホッフ応力に関連している.

主伸長間の関係を簡約するために,を意味する完全非圧縮性材料の仮説から始める.実際のFEM解析では,材料モデルの微圧縮性公式を使う.しかし,パラメータ推定では,非圧縮性のでも微圧縮性でも未知のパラメータは同じである.

単純な一軸引張下の完全非圧縮性モデルを仮定すると,非圧縮性制約条件 は主伸長間の関係を導く.引張方向に沿って となり,他の主伸長は以下のようになる:

次に,非圧縮性材料の第二ピオラ・キルヒホッフ応力テンソルは以下である:

ここで上で述べたように である.

一軸引張では,均一応力状態はおよびに減少する.適用された伸長 に関する方程式(4)の伸長関係を思い出すと, となる.これは自明ではなく,このチュートリアルの範疇を超える.関心のある方は[16]をご覧いただきたい. 応力成分は,以下のように,適用された伸長 に関するひずみエネルギー導関数に関連させることができる:

は一軸引張コーシー応力である.超弾性材料モデルの材料パラメータを見付けることは,特定の材料モデルのひずみエネルギー関数 への実験データの最適フィットを求める最適化となる.

一軸応力を計算するための補助関数を定義する:

次では一般のフィット過程を示す.[22]からのデータが超弾性構成モデルでフィットされる.フィットの結果は,構成モデルパラメータ定義について超弾性フレームワーク内で直接使うことができるため,例は以下のように定式化される.

実験データをロードする:
実験データをプロットする:

方程式(5)を使って一軸応力を実験データにフィットする.

例としてヨーモデルを使うが,これは他の材料モデルと同じように動作する.モデルパラメータ を持つヨーモデルのひずみエネルギー密度関数は以下で与えられる:

ひずみエネルギー密度関数を定義する:

非圧縮性ヨーモデルに対する方程式(6)の一軸応力成分は以下である:

一軸応力成分を計算する:
ヨーモデルをフィットする:
テストデータとフィットしたモデルの比較を示す:

上のプロットから分かるように,ヨー材料モデルはデータ範囲全体でうまくフィットすることはできない.これは珍しいことではない.これに対処する方法の一つに,フィットモデルの範囲を狭めるというものがある.

制限された範囲でヨーモデルをフィットする:
テストデータとフィットしたモデルの比較を示す:

フィットしたモデルのパラメータは超弾性材料シミュレーションの設定で直接使えるようになった.

変数を定義する:
メッシュを定義する:
フィット結果に基づいてヨーモデルを指定する:
強制荷重 に対するパラメトリック関数を作成する:
最大荷重 ,ステップ数 nsteps,初期変位ベクトルを指定する:
変位が増加するときの解の進行状況を監視する:
計算した変位に基づいて変形メッシュを生成する:
変形した物体を赤で,変形していない物体をエッジフレームで可視化する:

乗法分解

乗法分解は,熱膨張のような物理現象を超弾性材料の変形に取り込むために使われるアプローチである.

まず変形勾配を弾性部分 と非弾性部分 に分割する:

弾性部分は材料の説明と構造挙動に関連しており,非弾性部分は熱膨張のような追加の物理現象を捉える.「非弾性」という名前は構造挙動に対応する弾性に対比したものである.構造挙動によらない挙動はすべて非弾性によるものと考えられる.ここでいう非弾性とは,それがなくなってもずっと残る効果があるという意味ではない.

他の物理現象を考慮しなくてよい場合, であり以下が得られる:

どちらの場合でも はこれまでに見てきた変形勾配の役割を果たしている.混乱しないように,非弾性特性 がある場合, を全変形勾配と呼ぶ.

全変形勾配 は基準配置()を現配置()に投影する.分解のため,中間配置()はに投影する非弾性変形勾配によって表され,弾性変形勾配はに投影する.

560.gif

全変形勾配 は基準配置 を現配置 に投影する.分解により,中間配置 が存在する.非弾性変形勾配 に投影し,弾性変形勾配 に投影する.

材料モデルには弾性変形勾配が必要であり,以下のように定義できる:

弾性第二ピオラ・キルヒホッフ応力 は,ひずみエネルギー導関数に関連する現配置でまだ表される:

は超弾性材料挙動を表すひずみエネルギー密度関数であり, は全変形を表すコーシーグリーンテンソルである:

第一ピオラ・キルヒホッフ応力は,現配置から見られるように,基準配置で常に定義される.この場合,中間配置が現配置に対する基準配置であるため,PK1は中間配置で定義されるということである.これを弾性第一ピオラ・キルヒホッフ応力という:

平衡方程式を設定するためには,第一ピオラ・キルヒホッフ応力 を得るために,弾性第一ピオラ・キルヒホッフ応力 は基準配置に引き戻される必要がある.基準配置で表される は非弾性特性[30]の原因ともなる. は全第二ピオラ・基準配置でから得ることができ,次のように表される:

これで標準の が使われる.

非弾性変形勾配は"InelasticDeformationGradient"パラメータで指定することができる.

"InelasticDeformationGradient"を指定する:

マルチフィジックス結合

モデルについて複数の非弾性の寄与を考えなければならない場合がある.熱効果,弾塑性,ポリマー膨潤,軟組織成長[29,30,31]等の物理現象を加えることができる.このような場合,乗法分解が数回実行される.

変形勾配を弾性成分と非弾性成分に分割することを繰り返すことができる.その結果いくつかの非弾性成分を得る:

この場合,正しくPDEを書くために[30],いくつかの中間配置と引戻し操作が必要となる.非弾性変形勾配が適用される順序も関係がある.

いくつかの"InelasticDeformationGradient"を指定する:
超弾性材料モデルの乗法変形勾配の実装:
idg = pars["InelasticDeformationGradient"];
temp = Through[idg[vars, pars]];
temp = ExtendedDeformationGradient /@ temp;

idgInverse = Inverse /@ temp;
idgInverseDot = Dot @@ idgInverse;
idgJacobian = (Times @@ Det /@ temp);
f = f . idgInverseDot;

熱および物質移動等のいくつかの物理分野での応用例は「吸湿膨張」の応用モデルで利用できる.

熱弾性

熱機械結合問題を記述するために,非弾性変形勾配を使って温度効果 を取り入れる[29].熱の関与は純粋に体積測定によるものであると想定されることが多い:

ここで は熱ひずみであり,これは固体力学モノグラフの「熱弾性」セクションで説明したものと全く同じである.

ThermalInelasticDeformationGradientの実装:

ThermalInelasticDeformationGradientを使うためには,線形弾性の場合と同じように,熱ひずみを計算するのに必要なパラメータすべてを指定しなければならない.

結合アプローチ

一般的に,熱膨張等の他の物理現象を伴う固体力学のマルチフィジックス結合のアプローチには2通りの方法がある.一つはまず熱問題を解いてから,その解を固体力学アプリケーションの入力として使う逐次解法である.このアプローチの利点は,2つの方程式系の大きさが結合した系よりも小さいままなので計算的なコストが少なくて済むという点である.これは方程式が解かれる速さと使われる要素の数に影響を及ぼす.

2つ目のアプローチは温度と変形を同時に解く完全結合アプローチである.このアプローチは,変形が測定可能な方法で温度分布にも影響を及ぼす場合,より正確な解となることがある.

逐次結合アプローチ

逐次的に解かれる熱機械問題では,熱問題は非弾性変形勾配に対する熱膨張関係によって使われる温度場を得ることによって別に解かれる.

588.gif

逐次解法では,一つの現象に関連した一つの問題が先に解かれ,その解が2つ目の物理分野の入力として使われる.上の図は固体力学分野の熱荷重として使われる熱問題の処理を示している.

上で述べたように,逐次アプローチは計算コストを抑えるが,対処すべき複雑性が導入されることがある.超弾性で現れる可能性のある大きい変位のため,一つの重要な面は,逐次解法の繰り返しのたびに温度等の関連する物理量の解を変形した領域にマップし,マップし返すという点である.

次で説明するように,一般に完全結合アプローチの方がより正確な結果になり問題の式を簡約する.

完全結合アプローチ

完全結合アプローチでは固体力学PDEと熱PDEが同時に解かれる.単独の方程式系が設定されて解かれる.完全結合の設定では,各部分系が互いに影響を及ぼす.この特定の場合,変位は温度分布に影響を及ぼし,温度分布は変位に影響を及ぼす.つまり結合は両方向に働く.

熱方程式の最も簡単な形式は以下で与えられる:

は熱伝導率, は熱源を表す.熱方程式についての詳細は「熱移動」モノグラフを一読されたい.

熱方程式は中間配置で表される.熱方程式を固体力学PDEと一緒に解くためには,熱方程式は固体力学方程式と同じ参照枠に収まるように基準配置に戻されなければならない.これには熱方程式の発散演算子と勾配演算子の両方を引戻す必要があり[29, pp 521],参照配置で表される以下の定常状態熱方程式が導かれる.

ここで は熱非弾性変形勾配,である. 項は発散の前に置かれるため,この方程式は有限要素ソルバには向かないのである.詳細は「PDEの係数形式」のセクションで説明してある.発散内部でヤコビアン を引き出す必要がある.これは が発散内で持つ効果を補う別の項を加えることで行える.

これは適切な係数形式であり,実装する方程式は以下になる:

完全結合アプローチの例

次のセクションでは完全結合の2D熱機械分析を示す.

簡単なゴムの正方形領域が左側で固定され,右側で引っ張られている.この機械問題は上面で高温を,底面で対流熱流束を表す境界条件で記述される熱問題とつながっている.

600.gif

機械問題の境界条件.左側は固定され右側は引っ張られる.熱問題では上面が ,対流境界の熱移動係数は である.

関連する異なる現象のため,PDEは空間位置に依存する変位と温度の関数で表される.

固体力学と熱移動PDEに対する変数を宣言する:
矩形領域のメッシュを生成する:
シリコンゴムに対する材料パラメータを設定する:
熱パラメータを設定する:

熱非弾性変形勾配は両方の現象にかかわっており,「熱弾性」のセクションで説明するようにThermalInelasticDeformationGradient関数によって記述される.

PDEの2つの集合が別々に設定されてから完全結合パラメトリックソルバで結合される.これにより荷重乗数 を介して境界条件を徐々に増加させることが可能になる.

機械の記述についての材料パラメータを設定する:
最大荷重を に指定する:
牽引境界条件を設定する:
固定境界条件を設定する:
機械PDE集合を設定する:

次に熱部分を設定する.

熱非弾性変形勾配を設定する:
熱非弾性ヤコビアンと変形勾配の転置を設定する:
熱の記述に対する材料パラメータを設定する:
補償項を設定する:
定常状態熱方程式を設定する:
温度境界条件を設定する:
対流境界条件を設定する:
温度PDE集合を設定する:

最終的な結合PDEには機械と熱の両方の成分が含まれるため,解のベクトルには変位 および温度 場が含まれる.

大域的変数を定義する:

ここでパラメトリックソルバを設定する.熱膨張と外部荷重値が乗数パラメータ に接続される.

パラメータ を持つ強制境界条件に対するパラメトリック関数を作成する:

このパラメトリック関数によって,乗数 を0から1まで徐々に大きくすることができるため,境界における辺の牽引力と熱交換の両方が増加させられる.乗数は二次関数的に増加する.これは線形増加の法則に関して初期ステップで境界条件に課せられる合計ステップ数を減らすのに役立つ.二次関数的増加で収束の問題が発生した場合は,線形増加に戻ることをお勧めする.

乗数の最大値を設定する:
初期の解ベクトルのステップ数 nsteps を指定する:
荷重因子 が二次関数的に増加するときの解法の進行状況を監視する:
変位の結果を抽出する:
温度の結果を抽出する:
計算した変位に基づいて変形メッシュを生成する:
変形した物体を赤で,変形しなかった物体を黒で可視化する:
変形したメッシュ上に温度場をマップする:
温度の結果を可視化する:

材料は等方性であり,牽引力だけが右側すべてについての水平変位定数を持つ水平方向の拡張を導く.しかし垂直方向の温度勾配により,下側がわずかな圧縮効果を示し上側の拡張が増大していることで結果的に曲がったような効果をもたらしている.

全ひずみ

グリーン・ラグランジュひずみテンソル は全ひずみを表す.乗法分解のため,これは全コーシー・グリーンひずみテンソル 内部の非弾性変形と弾性変形の両方に寄与する.

全ひずみを計算する:
最小ひずみと最大ひずみを計算する:
全ひずみのいくつかの成分を可視化する:

非弾性ひずみ

全グリーン・ラグランジュひずみテンソルには弾性および非弾性の寄与が含まれている:

この例における非弾性の寄与は,熱変形勾配 を介して以下のように完全に指定される:

これで全ひずみテンソルとは別に熱ひずみを計算し可視化することができる.

非弾性コーシー・グリーンひずみテンソルを計算する:
非弾性熱グリーン・ラグランジュひずみテンソルを計算する:

次に熱ひずみ が計算できる.非弾性熱変形勾配は温度 と荷重乗数 の関数として定義される.これらは実際の値で置き換える必要がある.材料が等方性ですべての対角成分は同じなので,熱ひずみの対角成分の一つを考えるだけで充分である.

熱ひずみを計算する:

変形していないメッシュや変形したメッシュ上の温度結果を置換することは可能であるが,プロット領域はそれに応じて選ばなければならない.deformedTemperatureを置換したため,deformedMesh上にプロットする必要がある.

熱ひずみを可視化する:

弾性ひずみ

弾性の寄与を抽出するために,全コーシー・グリーンひずみから始めて方程式(7)の定義を反転する:

これでグリーン・ラグランジュひずみテンソルの定義を使って以下のように弾性の寄与を計算することができる:

非弾性コーシー・グリーンひずみテンソルのパラメータを置換する:

全コーシー・グリーンひずみテンソル を計算するために,について を解き以下を得る:

全コーシー・グリーンひずみテンソル を計算する:
弾性コーシー・グリーンひずみテンソル を計算する:
弾性グリーン・ラグランジュひずみテンソル を計算する:
弾性ひずみ のいくつかの成分を可視化する:

等容ひずみと体積ひずみの分割

ひずみを等容に分割したものは偏差ひずみとしても知られているが,体積成分は材料が外部または内部の荷重にどのように応答するかということの洞察を提供する変形を強調する一般的な方法である.等容ひずみは体積の変化なしで形状の変化を伴う変形の一部を指す.体積ひずみは材料の圧縮や膨張による体積の変化を表す.

この例では2つの異なる現象を伴っている.これらは体積熱効果と荷重牽引で,結果として体積寄与と等容寄与の両方になる.

等容ひずみへの体積ひずみの関係内部で弾性寄与と非弾性寄与を分割するのに前述の解釈を使うことができる.全コーシー・グリーンひずみテンソルも次のように分解できる:

ここで および はそれぞれ体積および等容のコーシー・グリーンひずみテンソルを表す.

コーシー・グリーンひずみテンソルの等容寄与は次のように表される:

ここで である.次に体積寄与は以下になる:

または,それをテンソルの体積部分として直接表すこともできる:

ここで は単位行列である.ヤコビアン指数の分母はテンソルの次元(通常空間関係では二次テンソル)に関連していることも重要である.

弾性/非弾性ヤコビアンを使って,弾性または非弾性のコーシー・グリーンひずみテンソルの等容寄与と体積寄与を判別するために同じ関係を使うことができる.

この例では非弾性寄与はThermalInelasticDeformationGradientで表され,その体積寄与と等容寄与を調べることができる.

非弾性ヤコビアンの最終結果を置き換える:
領域次元を定義する:
等容非弾性コーシー・グリーンひずみテンソルを計算する:
等容非弾性グリーン・ラグランジュひずみテンソルを計算する:
体積非弾性コーシー・グリーンひずみテンソルを計算する:
体積非弾性グリーン・ラグランジュひずみテンソルを計算する:
等容非弾性グリーン・ラグランジュひずみテンソルの平均を計算する:

熱ひずみは純粋な体積膨張なので,予想通り,等容寄与は数値的にゼロである.実際,非弾性体積寄与は等方性であり完全に前述の熱ひずみに相当する.

体積非弾性グリーン・ラグランジュひずみテンソルを可視化する:

応力

機械的応力の評価は,材料が組み合された効果にどのように反応するかを理解する重要な一面である.変形勾配の乗法分解により,結果のコーシー応力は機械的荷重および熱荷重の両方を含む平衡方程式を満足する.

平衡方程式は変形設定または基準設定で表されるので,熱場に関連する中間領域にマップすることは物理的な意味を持たない.特に応力は基準設定で表される.

コーシー応力を計算する:
ミーゼス応力を計算する:

乗数 と温度 を代入すると,式は簡単な補間関数ではなく,変形していないメッシュ上で評価すると単独の補間関数を構築する複合式になる.

ミーゼス応力の結果を代入し,補間関数を構築する:
ミーゼス応力を変形したメッシュ上にマップする:
変形したメッシュ状のミーゼス応力を単位[kPa]で可視化する:

左上側の応力はその他の部分と比べて高い.これは上部領域の温度による大きい体積膨張と対称的な固定制約条件による.この集中にも関わらず,応力はほぼ一定であり境界荷重を満足する範囲にある.

複合材料構成モデル

複合材料とは複数の材料からなる.さまざまな構成モデルで記述される複合材料領域の設定を説明するために,基本の正方形を3つの領域に分け,それぞれを異なる構成超弾性モデルで説明する.

変数を定義する:
次元を定義する:
材料の接触面を区別した内面を含む境界を定義する:

各領域に対して異なる構成モデルを導入するために,ElementMarkerを使う.

以下の公式は,各層の間が完全に接着しているという暗黙の仮定に依存する.このモデルは完全接着として一つの材料領域から別の材料領域へと荷重を移動させることを観察するのが重量である.応力は力学的平衡により接触面全体で連続的な実体となる一方,応力は異なる構成モデルに依存する可能性がある.

マーカーの定義:
材料の接触面を区別した内面を含む境界を定義する:
メッシュを生成する:
メッシュと領域を可視化する:
左側の固定された条件を定義する:
牽引力についての境界条件を定義する:

異なる構成モデルを適用するために,ヘルパー関数を使う.この関数の目的は,ElementMarkerによって各成分が実際の構成モデルに送られる構成モデルを返すことである.

この関数はparsで与えられるすべてのモデルを抽出することで動作する.与えられたモデルのそれぞれに対して構成モデルが作成される.それから区分条件が設定され,材料モデルの数を照合する.モデルは1から始まる連続番号であるElementMarkerを送る.したがって与えられる最初のモデルはElementMarker==1を持つことになる.

複合材料ヘルパー関数を定義する:

次の例では,3つの異なる材料を3つの異なるサブドメインArruda-Boyce,ヨー,ネオフックに適用する.このためにまず大域的な pars を生成してからサブモデルパラメータを生成する.最後のステップでこれらを大域的 pars に保存する.

大域的材料パラメータを定義する:
ネオフックモデルの材料パラメータを定義する:
ヨーモデルの材料パラメータを定義する:
Arruda Boyceモデルの材料パラメータを定義する:
サブモデルを設定する:
パラメトリック関数を設定する:
最大牽引力を定義する:
ステップ数を定義する:
変位を初期化する:
解く:
変形メッシュを生成する:
変形オブジェクトを表示する:

物体の変形から分かる通り,異なる構成モデルが使ってあるため,材料の接触面全体にわたり変形の非連続性が見られる.

横等方性の超弾性材料

横等方性材料については,線形弾性の場合として固体力学モノグラフで紹介した.このセクションではその概念を超弾性材料に拡張し,それをどのように繊維強化材料等に使うことができるかを示す.

超弾性材料は材料が変形によって得たエネルギーをどのように保存するかを記述する,ひずみエネルギー密度関数 によって特徴づけられる.完全等方材料の場合,このひずみエネルギー密度関数 は右コーシー・グリーンひずみテンソル の3つの主不変量 の関数である.

横等方性材料の場合,材料を強化する繊維等の材料対称性を含める必要がある.基準設定の繊維の方向,つまり一般的に対称軸に接する単位ベクトルを で表す.現在の設定では類似した単位ベクトルを定義し で表す.これら2つのベクトルの関係は以下で与えられる:

ここで は変形勾配, は方向 への繊維の伸長である.これを強調するために繊維の伸長を考慮し,材料行列は考慮しない.

横等方性超弾性材料のひずみエネルギー密度関数を記述するために,スペンサー[24]はそれぞれ繊維方向と右コーシー・グリーンひずみテンソル に依存する2つの擬似不変量 を導入した:

最初の疑似不変量 は方向 の伸長の平方に等しい:

これらの擬似不変量の右コーシー・グリーンひずみテンソル についての導関数を知ると便利である:

ここではテンソル積であり,関数TensorProductを使って計算することができる.

ここで連鎖規則を使って,横等方性超弾性材料の第二ピオラ・キルヒホッフ応力テンソル についての方程式を導出する:

不変量と擬似不変量の導関数を使う:

これは,等方性超弾性固体の第二ピオラ・キルヒホッフ応力テンソルの定式を最後の2項で拡張する.この定式の性質を考えると,ひずみエネルギー密度関数 を材料の等方性応答を記述する関数と横等方性応答を記述する関数の2つに分けることが可能である.これにより以下の関係が得られる:

ここで は主不変量 の関数,は擬似不変量 の関数である.項 は強化モデルと言われることもある.

標準的な強化材料モデル

このセクションでは繊維強化超弾性材料のモデル化の方法を説明する.標準的な強化モデルは,強化モデル の最も簡単な例であろう.これは簡単なためかなり人気が高く,繊維強化材料や生物組織のモデル化に使われる.例えば,脳幹の弾性応答のモデル化[23]に使うことができる.ひずみエネルギー密度関数 は以下の形式である:

ここで μ は剛性率,は第一擬似不変量,は基準設定の繊維方向の単位ベクトル, は左コーシー・グリーンひずみテンソルである.

無次元パラメータ は繊維剛性を示す. の値が大きいほど,繊維剛性は大きい.生体力学ではこの値は通常10代の低い値になる.例えば,生後4週間の子豚の脳幹の剛性は10 [23]くらいである.ならば,は基本のひずみエネルギー関数を回復する.

次に標準的な強化モデルの操作方法の例を示す.この例は圧縮性ネオフックモデルのひずみエネルギー密度関数に基づいており,結果として次のひずみエネルギー密度関数になる:

ここで λ はラメパラメータ, は変形勾配である.

材料パラメータを線形弾性の場合で考えると,最適なのは実験を行うことである.他にも,各方向の材料パラメータを推定しようとする均質化シミュレーションを使う方法もある.もっと簡単なものとして,行列の材料特性を横方向に使って,ひずみエネルギー密度関数 の等方性部分と繊維のパラメータを軸方向,つまり の横等方性部分に使うというものもある.

横等方性材料の仕組みの概要は[16]にある.繊維強化材料に焦点を当てた詳しい概要は[24]をご覧いただきたい.

繊維強化材料

例として繊維がどのように変形に影響するかを調べる.で維持されておりまっすぐな繊維を持つ2D矩形領域を使うが,で規定の圧力を加える.繊維がどのように変形に影響するかを調べるために,さまざまな繊維剛性およびさまざまな方向を調べる.

変数を定義する:
領域を定義する:
メッシュを生成し可視化する:
基本パラメータを定義する:
標準的な強化材料モデルの繊維方向および繊維剛性を定義する:

および の値はパラメータとして保持される.これは,それらがどのように変形に影響するかを調べるためである.

モデルを設定する:
パラメータ を持つパラメトリックプロットを設定する:
繊維の方向と剛性の集合を定義する:
規定の繊維剛性 および角度 においてパラメトリック関数を評価し,かかる時間を測定する:

シミュレーションの結果は保存されない.パラメトリック関数が計算結果をキャッシュに残すから保存する必要はないのである.同じ入力値を使うと,キャッシュされた結果が取り出され返される.

繊維の剛性 および角度 でのメッシュの変形を可視化する:
メッシュがどのように変形したかを計算する関数を設定する:
変形したメッシュおよびもとの領域(うすい灰色)を可視化する:

変形が繊維強化材料の繊維にどのような影響を与えるかを調べるためには,現在の設定 での繊維の方向が必要である.これは基本設定 における繊維方向に変形勾配 を適用すること,で計算できる.変形勾配 であり,関数DeformationGradientで計算できる.

を計算するために関数を設定する:
を変形したメッシュ上に投影する関数を設定する:
の値および繊維変位 によって,変形した繊維設定,初期の設定,最終的な設定をプロットする関数を作成する:
の値と繊維の変位 についての変位プロットを生成する:
高忠実度の可視化を得るために,ラスタ化処理をコメントアウトする:
フレームをアニメーション化する:

このプロットは背景にもとの繊維変位(灰色)を示し,変形した物体を新しい繊維方向で示している.

標準的な強化材料モデルを使って簡単な例を解くことができたが,ここで繊維特性がどのように結果に影響を及ぼすかを調べることができる.この材料は確かに繊維強化材料のように振舞っていることが分かる.のとき,変形が繊維方向に依存しないのは予想通りである.では標準的な等方性ネオフックモデルだからである.また,適用された圧力に垂直な繊維の場合,材料の伸長は繊維が変形と平行な場合よりも大きい.次に,繊維が適用された圧力に垂直であっても変形は繊維剛性に依存することも分かった.圧力を加え始めると,繊維はもはやまっすぐではなくなり,適用された圧力に垂直になるためである.

についての変形を示す:

繊維が補強の役割を果たしていることが分かる.適用された圧力に対して繊維が垂直の場合,材料の伸長はずっと大きい.

のときの変形を示す:

先の場合よりも差は小さいが,適用された圧力に対して繊維が平行の場合でも繊維が変形に影響を及ぼしていることが分かる.

曲がった繊維

上の例で,繊維強化材料の使い方を示した.簡単なものにするために,2次元に限定し繊維方向 が定数であると仮定した.次のステップとして,このような制約なしで同じような問題を考える.唯一の仮定は,材料が均質,つまり材料パラメータが一定ということである.一辺で保持され,別の辺に規定の圧力が適用される単位立方体を考える.

変数を設定する:
領域を作成する:
領域のメッシュを生成する:
繊維剛性の特定の値 と非定数の繊維方向角度 に対するパラメータを更新する:
立方体の中の繊維を可視化する:
モデルを設定する:
パラメトリック圧力 に対するパラメトリックソルバを設定する:
選ばれた圧力に対する変位を計算し,かかった時間を監視する:
変形した物体を可視化する:

結果の変形は一様ではないことが分かる.これは材料が非等方性,特に横等方性という性質を持つ結果である.次に現在の設定において変形したメッシュと繊維を示す.

変形したメッシュを生成し可視化する:

次に現在の設定における繊維を可視化する.

現在の設定における繊維方向 を計算する.これは (Fは変形勾配 )に等しい.

もとのメッシュの を計算する:
変形したメッシュに を投影する:
現在の設定における繊維を可視化する:

上の可視化では,メッシュの要素があまりに変形しているので,ToElementMesh[ToBoundaryMesh[deformedMesh]]を使って変形したメッシュのメッシュを再び生成する必要がある.これは可視化のためだけに必要な手順である.

結果の変形は,繊維が曲がっているため,前の例のものほど単純ではない.このプロットでは,結果の変形は繊維方向に少しねじれていることが分かる.

曲がった繊維を含む例も解くことができた.

非伸縮繊維

材料内の繊維は,変形の際に長さが変化しなければ非伸縮と考えられる.つまり伸長 ということである.という制約条件は,材料の微圧縮性を強制するのに使われる条件と類似した方法で強制される.ひずみエネルギーのラグランジュの未定乗数 をひずみエネルギー密度関数に加える:

ここで に依存しない.これは伸長が である,つまり第4不変量 に等しいということである.したがって,ひずみエネルギー密度関数 に依存しないのである.

完全非伸縮性は現在Wolfram言語には実装されていない.しかし の代りに だけを達成しようとすることで微圧縮性アプローチを模倣することはできる.これは以下のように設定して行う:

ここで は繊維の材料剛性を決定する.ひずみエネルギー密度関数は以下の形式を持つ:

ひずみエネルギー関数の最後の項は,繊維の方向における変形にペナルティを課す.ひずみエネルギー関数 に依存しないならば,標準強化材料モデルを得ることになる.つまり標準強化モデルは材料を繊維に沿って非伸縮にする方法として見ることができる.

次では,十分剛性のある標準強化材料モデルを使ってほとんど非伸縮の材料を得る方法を示す.

領域を設定しメッシュを生成する:
変数と材料パラメータを設定する:
繊維の方向と繊維の剛性を設定する:
モデルを設定する:
方程式を解いて変位を計算する:
変形した物体を表示する:
ひずみを計算する:
応力とそのノルムを計算する:
変形した物体の応力を表示する:
PlotRangeをAllに設定して変形した物体の応力を表示する:
変形勾配 を計算する:
繊維の伸長を計算する:
繊維の伸長をプロットする.一つは概要でありもう一つはPlotRangeをAllに設定したものである:

繊維の伸長 はほぼあらゆる場所でほぼ1に等しいことが分かる.唯一の例外は対角要素付近の領域であるが,そこでも伸長は小さい.伸長が大きい領域は応力が大きい領域に対応している.

対角の伸長をプロットし,2つ目はPlotRangeをAllに設定したものをプロットする:

伸長が対角に沿ってすぐに消失しているのが分かる.また線上に一様にジャンプが現れているのも分かる.これはメッシュに対応しており,細かいメッシュを取るとジャンプが多くなる.

2つの繊維のグループからなる材料

ここまでは1つの繊維のグループで強化された材料だけを考えてきた.より複雑な材料を考える場合,この仮定は制約的でありすぎる可能性がある.例えば,血管壁の層は主に等方性マトリックス材料であるエラスチン,2つの繊維のグループ,コラーゲンが対称的ならせんに配列されて構成されている.

このような種類の材料は,その対称性の中心となる軸がないため,もう横等方性ではない.しかし,一般的なアプローチは1つの繊維のグループの材料の場合と似ている.

まず,単位ベクトル で標準設定の繊維方向を記述し,単位ベクトル で現在の設定における繊維方向を表す.ここで2つの擬似不変量と類似のものが加えられるが,今回は2つ目の繊維のグループに対して行う:

右コーシー・グリーンひずみテンソル に関するこれらの擬似不変量の導関数は,の導関数に類似している:

しかし,2つの繊維のグループ間の相互作用を含む必要がある可能性があるため,これらの新しい擬似不変量を加えるだけでは不十分かもしれない.したがって以下も加える:

擬似不変量 は標準設定における2つの繊維のグループ間の角度で決まり,変形の影響は受けない.このような理由で,これは省略することができる.

右コーシー・グリーンひずみテンソルに関する の導関数は以下に等しい:

ここで連鎖規則を使って第二ピオラ・キルヒホッフ応力テンソルの公式を書くことができる:

ならば,2つの繊維のグループは直交であり,材料は繊維に垂直な平面および繊維が存在する平面について直交異方性である.

ここまで見てきた通り,横等方性材料をモデル化するために使ったアプローチを拡張して,2つの繊維のグループで強化される材料をモデル化することができる.これは血管等の生物材料または織り繊維によって強化されるファイバーグラスのような材料のモデル化にも使うことができる.

参考文献

1.  Backstrom, G; Simple Deformation and Vibration; GB Publishing; 2006

2.  The Efficient Engineer, https://www.youtube.com/watch?v=aQf6Q8t1FQE; retrieved 2021/4/27

3.  Cook, R; et. al.; Concepts and Applications of Finite Element Analysis; Wiley 2002; 4th Edition

4.  Bhatti, M.; Fundamental Finite Element Analysis and Application; Wiley 2005

5.  Bhatti, M.; Advanced Topics in Finite Element Analysis of Structures; Wiley 2006

6.  Brand, M.; Grundlagen FEM mit Solidworks; Vieweg+Teuber 2011; ISBN: 978-3-8348-1306-0

7.  https://en.wikipedia.org/wiki/Modal_analysis_using_FEM, retrieved 01/06/2021

8.  Nasdala, L.; FEM-Formelsammlung Statik und Dynamik; Vieweg+Teuber 2010

9.  Alawadhi, E. M.; Finite Element Simulations Using ANSYS; CRC Press 2010

10.  Constantinescu A., Korsunsky A.; Elasticity with Mathematica; Cambridge University Press 2007

11.  Bower, A.; Applied Mechanics of Solids; CRC Press 2010; ISBN-13: 978-1439802472; (http://solidmechanics.org/contents.php)

12.  The Efficient Engineer, https://www.youtube.com/watch?v=xkbQnBAOFEg; retrieved 2021/9/13

13.  Kelly, P.A.; Mechanics Lecture Notes: An introduction to Solid Mechanics. (http://homepages.engineering.auckland.ac.nz/~pkel015/SolidMechanicsBooks/index.html)

14.  Moosbrugger, C. (Editor); Atlas of Stress-strain Curves; ASM International 2002; ISBN 0-87170-739-X (https://books.google.com/books?id=up5KS9fd_pkC)

15.  Sifakis, E., Barbič, J.; Finite Element Method Simulation of 3D Deformable Solids; M&C Publishers 2016; ISBN 9781627054423

16.  Holzapfel, A.; Nonlinear Solid Mechanics; Wiley 2020; ISBN 978-0-471-82319-3

17.  Bonet, J., Wood, R.; Nonlinear Continuum Mechanics for Finite Element Analysis; Cambridge University Press 1997; ISBN 0-521-57272-X

18.  Mooney, M. Large elastic deformations of isotropic materials. Journal of Applied Physics 1940. 11: 582592. https://doi.org/10.1063/1.1712836

19.  Rivlin, R. S. Large elastic deformations of isotropic materials. Philosophical Transactions of the Royal Society of London. Mathematical and Physical Sciences 1948. 241: pg 379397. https://doi.org/10.1098/rsta.1948.0024

20.  Arruda, M. Boyce, M. A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. Journal of the Mechanics and Physics of Solids. 1993. Volume 41, Issue 2. Pages 389-412,. https://doi.org/10.1016/0022-5096(93)90013-6

21.  Gent, A. N. A New Constitutive Relation for Rubber. Rubber Chemistry and Technology. 1996. 69 (1): 5961. https://doi.org/10.5254/1.3538357

22.  Treolar, L.R.G.; Stress-strain data for vulcanised rubber under various types of deformation; Transactions of the Faraday Society, 1943

23.  Ning, X., Zhu, Q., Lanir, Y., Margulies, S. S.: A Transversely Isotropic Viscoelastic Constitutive Equation for Brainstem Undergoing Finite Deformation. Journal of Biomechanical Engineering, 2006. 128: pg 925-933. DOI: 10.1115/1.2354208

24.  Spencer, A. J. M.: Deformations of fibre-reinforced materials, 1972

25.  Nekliudova, E. A., Semenov, A. S., Melnikov, B. E., Semenov, S. G.: Experimental research and finite element analysis of elastic and strength properties of fiberglass composite material. Magazine of Civil Engineering 3 (47), 2014

26.  Otero F., Oller S., Martinez X., Salomón O.: Numerical homogenization for composite materials analysis. Comparison with other micro mechanical formulations. Composite Structures. 2015 Apr 1;122:405-16.

27.  Charalambakis, N.: Homogenization techniques and micromechanics: A survey and perspectives. 2010: 030803.

28.  Pascon, J. P. (2019). Large deformation analysis of plane-stress hyperelastic problems via triangular membrane finite elements. International Journal of Advanced Structural Engineering, vol 3, pp. 331-350, https://doi.org/10.1007/s40091-019-00234-w.

29.  Yosibash, Z. Weiss, D. Hartmann, S. (2014). High-order FEMs for thermo-hyperelasticity at finite strains. Computers and Mathematics with Applications, vol 67, pp. 477-496, https://doi.org/10.1016/j.camwa.2013.11.003

30.  Wriggers, P. Nonlinear finite element methods; Springer 2008

31.  Gasser, T. C. Vascular Biomechanics: Concepts, Models and applications; Springer 2021 https://doi.org/10.1007/978-3-030-70966-2