超弾性モデルの比較

次の例では,自分の実験データに最も適したモデルの選ぶために,さまざまな超弾性モデルを短い説明を交えて比較する.

はじめに

特定の材料用に正しい超弾性モデルを選ぶ際は,材料の力学的挙動,利用可能な実験データ,意図する応用方法を含むいくつかの要因を考える必要がある.

実験データは適切な超弾性モデルを選ぶのに必須である.さまざまな荷重条件の下における材料の応力ひずみ反応を測定するために,単軸引張試験,圧縮試験,せん断試験等の力学的試験が行われる.

例えば軟組織の場合,比較的簡単で制御された実験設定が可能なため,実験データを得るために単軸引張試験が広く使われる.標本は一方向での制御された張力を受け,正確な測定と分析か可能である.単軸荷重条件は,応力とひずみが主に1つの軸に沿って確実に適用されるようにし,力学的挙動の解析を簡約化する.変位データは試験マシンのグリップ間の力の関数として記録される.

実験データについての詳細は超弾性モノグラフでご覧いただきたい.

このチュートリアルで使う追加の関数のために有限要素コンテキストをロードする.

有限要素パッケージをロードする:

単軸試験データの生成

関心のある超弾性モデルのをフィットするために,前述の材料挙動を記述する実験データが必要である.次のセクションでは,単軸引張試験が応力ひずみ関係を表すデータを集めると考えられる.実験室に行く必要がないように,実験データはバーチャルで生成される.したがって,特定の構造挙動によって起こる応力と伸張の間の関係を解析的に表す.

超弾性モノグラフで述べた通り,単軸応力はひずみエネルギー と伸張 と以下のように関係している:

バーチャル実験データを生成するためには,構造モデルを選ぶ必要がある.例えば,多項式関数の簡潔性を持ちながら明確な非線形挙動を示すヨーモデルを使う.明らかに構造モデルや材料パラメータによって実験データはもとより,次の比較結果もすべて変化する.

まず,ヨーひずみエネルギーを単軸応力式に挿入して以下を得る:

ここで は最初のひずみ不変量 である.

非圧縮制約条件 のため,第2,第3の主伸張はメインの主伸張と という関係を持つ.最初のひずみ不変量はメインの伸長 について以下で表すことができる:

以下では,実験データは実際のデータの不完全性を模倣して,ノイズを加えた応力・伸張関係として生成される.加えられたノイズは,実行された材料モデルのフィットがどのようにデータの質品質に影響されるかも示す.

この例では2つの異なるデータ範囲を考えるということも大切である.実験からのデータは実際のモデルフィットに使われるデータよりも,通常大きい範囲に広がる.これはなるべく正確なフィットを作成したいからであり,関心のある範囲のデータだけをフィットしそれ以上は使わないようにすることができる.

単軸試験データを生成する補助関数を作成する:
ヨー単軸試験データを生成する材料パラメータを定義する:
実験手順パラメータを定義する:
ヨーモデルの単軸引張応力データを生成する:
ノイズを生成する:
ランダムなノイズを可視化する:
生成したノイズを応力データに加える:

バーチャルの実験データは最大伸張の に達したが,モデルのフィットは制限された範囲で行われる.これは後でモデルの正確性がどの程度選んだデータ範囲に依存してるかを示すものである.実際,後で見るように,ある範囲で同様に動作するモデルでも,その範囲を超えると大きな差を示すことがある.

フィットのための範囲を制限する:
2つの異なる範囲の間の差を示す.赤い点のデータはモデルのフィットに使われ,灰色のデータは利用可能なデータである:

これをもう一度強調するが,ここでモデルフィットに使われた最大伸張の選択は任意である.実際のデータを使ってフィットを行う場合は,実際の応用場面で想定される伸張を考慮に入れ,それに応じたデータ範囲を選ばなければならない.

モデルのフィット

次のステップでは,単軸引張試験データに基づいて超弾性材料モデルフィットを生成する.通常データは一連の実験から取得するか,ここで行っているようにバーチャル材料データを作成するかする.このセクションではいくつかの非圧縮性超弾性材料モデルを同じデータにフィットする.

特定のモデルのひずみエネルギー を考慮することで,フィット過程の間の所望の構成モデルに対する材料パラメータが特定できる解析式が導ける.

各モデルについて,材料パラメータの特定過程は同じである.まず方程式(1)で示してあるように,単軸応力 を計算するために第1,第2の不変量についてひずみエネルギーを導出する.その後試験データを表す最適な材料パラメータを探すために数値フィットを行う.

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

ネオフックモデル

非圧縮性ネオフックモデルのひずみエネルギーは以下である:

ここで材料パラメータは である.

単軸試験では,方程式(2)からの応力は以下になる:

ひずみエネルギーを定義する:
材料パラメータを特定する:
試験データとフィットしたモデルを比較する:

パラメータ2つのムーニー・リブリンモデル

パラメータを2つ持つ非圧縮性ムーニー・リブリンモデルのひずみエネルギーは以下である:

ここで は2つの材料パラメータである.

一貫性のあるフィットには,ひずみエネルギーの正値性を満足するために追加の制約条件 も必要である.この制約条件によって数値的な材料安定性が確実になる.

単軸試験では,方程式(3)からの応力は以下になる:

ひずみエネルギーを定義する:
材料パラメータを特定する:
試験データとフィットしたモデルを比較する:

パラメータ5つのムーニー・リブリンモデル

パラメータを5つ持つ非圧縮性ムーニー・リブリンモデルのひずみエネルギーは以下である:

単軸試験では,方程式(4)からの応力は以下になる:

ひずみエネルギーを定義する:
材料パラメータを特定する:
試験データとフィットしたモデルを比較する:

パラメータ9つのムーニー・リブリンモデル

パラメータを9つ持つ非圧縮性ムーニー・リブリンモデルのひずみエネルギーは以下である:

単軸試験では,方程式(5)からの応力は以下になる:

ひずみエネルギーを定義する:
材料パラメータを特定する:
試験データとフィットしたモデルを比較する:

Arruda-Boyce

非圧縮性Arruda-Boyceモデルのひずみエネルギーは以下である:

ここで はモデルの2つの材料パラメータであり, はテイラー級数からの数値係数である:

単軸試験では,方程式(6)からの応力は以下になる:

ひずみエネルギーを定義する:
材料パラメータを特定する:
試験データとフィットしたモデルを比較する:

Gent

非圧縮性Gentモデルのひずみエネルギーは以下である:

ここで は2つの材料係数である.

単軸試験では,方程式(7)からの応力は以下になる:

ひずみエネルギーを定義する:
材料パラメータを特定する:
試験データとフィットしたモデルを比較する:

ヨー

非圧縮性ヨーモデルのひずみエネルギーは以下である:

単軸試験では,方程式(8)からの応力は以下になる:

ひずみエネルギーを定義する:
材料パラメータを特定する:
試験データとフィットしたモデルを比較する:

解析比較

次のセクションではフィットしたさまざまなモデルを比較する.フィットの範囲に焦点を当てると,ヨー,パラメータ5つのムーニー・リブリン,パラメータ9つのムーニー・リブリンがその多項式特性のため利用可能なデータに最もよくフィットしているように見える.通常のムーニー・リブリン,Arruda-Boyce,Gentはネオフックモデルのように,線形のような挙動に還元されているように見える.

フィット範囲における実験データとフィットしたモデルをプロットする:

モデルの挙動は厳密にフィット範囲に限られる.伸張 を超えてモデルを使うと,大部分のモデルは残りの実験データから逸れていくが,パラメータ5つのムーニー・リブリンとヨーモデルは想定される応力範囲に近いままである.モデルはデータの小さい部分 にだけ特にフィットされており,その伸張を超えた挙動を示し,フィット範囲で与えられる妥当性を超えてモデルを使わないよう注意を喚起している.全体のデータ範囲へのフィットは異なる結果となる.

フィットしたモデルをデータ範囲外まで拡張する:

有限要素法の単軸試験

次のセクションではいくつかの有限要素法解析を実装して,モデル間の数値的差を示す.以下の解析は超弾性挙動を示す軟組織の標本で実行された単軸引張試験である.材料を記述するために,以前フィットしたモデルを使う.

3Dモデル

単軸引張試験は材料を力学的に試験するための標準化された手順である.ここでは矩形の標本を一端で引っ張る.まずこのセクションでは完全な3Dモデルを示す.計算の複雑性を減少させるために,後続のセクションでは2Dモデルを使う.

変数定義:
大きさの定義:
境界点定義:
指定されたクランプ領域を持つ境界メッシュを作成して可視化する:
平面メッシュを作成して可視化する:
3次元メッシュを作成して可視化する:
メッシュ内のノード数を求める:

一つの次元(厚さ)が他の次元よりもずっと小さいこの位相構造の形状に対するメッシュを生成するためには,多数の要素が必要である.これは主に小さい(厚さ)層に十分な数の要素が必要で,要素の大きさが構造の大きさに比べて非常に小さくならなければならないためである.

厚さの方向に少なくとも3層の要素を使う理由は大きく分けて2つある.一つは厚さ方向の応力勾配を正確に捉えるためであり,もう一つは張力を受けた構造の幾何学的挙動を捉えるためである.平面応力式は単一の要素からできた層の方が適している場合がある.さらに全体的な構造は想定よりも硬い結果になる可能性がある.

z座標に沿った3つの異なる理想層をハイライトするメッシュを可視化する:
ネオフック材料モデルの材料特性を分析する:
パラメトリックソルバを作成する:
境界の牽引力を定義する:
解の変数を定義する:
荷重 をゆっくり増加させることで変位について解く:
変形したメッシュを計算する:
ひずみを計算する:
変形したメッシュ上にひずみ場をマップする:
可視化のために変形したメッシュにもう一度メッシュを施す:
変形したメッシュ上の結果のひずみ場を表示する:

この標本は120%まで変形している.しかし,右側は固定境界として表されていなかった.ここでは単に境界に力を加えただけで端が大きく変形したことを示している.

変形した端の結果のひずみ場を表示する:

実験手順,クランプの使い方をよりよく表すために,固定領域のモデル化を変える.標本の固定領域は,クランプ内で固定されたままであり,それが引っ張られるのであるから,ほぼ変形しない状態のままになるはずである.この挙動を達成するために,クランプのモデルとしてクランプが適用される領域で,標本の材料特性をずっと硬くする.こうすることで標本の中央に変形を集中させ,より正確な結果となる.

また,要素の3つの線を厚さで使用して正確な結果を得る必要がある.これでモデルの大きさが大きくなる.モデルの最初の表現は2次元で得られる.

2Dでの比較

問題を2Dに還元することによって,いくつかの材料モデルを比較してその差を表示することが素早くできる.ここで上の3D一軸引張試験を2Dで再実行する.一つ異なるのは,領域を拡張して標本を引っ張るクランプの特別な領域を作ったことである.クランプはより硬く,その挙動は少なくとも2倍の大きさの硬度係数を持つネオフックモデルで再生される.この手順で端の変形を減少させ標本全体の解析を拡張することができる.

まず,すべての材料モデルのアプリケーションで使えるメッシュと境界変数を定義する.またいくつかの補助関数を定義して,解析されるそれぞれのモデルについて同じ境界問題,対応するパラメトリックソルバ,解法を定義する.

変数定義:
マーカー定義:
マーカーを含む2Dメッシュで再びメッシュを作成する:

以下では,どのシミュレーションにも同じ境界条件が使われる.

左側の固定条件を定義する:
最大牽引力を定義する:
牽引力に対する境界条件を定義する:

以下のシミュレーションでは,考慮されるそれぞれのモデルに対するいくつかの情報を抽出することができる.各解析について,パラメトリック解法の同じステップにおける応力・ひずみ関係を抽出する.

すべてのシミュレーションにおいて,クランプ領域は高剛性のネオフック材料としてモデル化する.これによって,先端の変形が減少し,標本のひずみに集中することができる.

パラメトリックソルバのステップ数を定義する:
クランプに異なる材料を実装するための補助関数を定義する:
パラメトリックソルバを実装するための補助関数を定義する:
形状の面積を計算する:
PDEを解き,決まったステップにおける変位,応力,ひずみを抽出するための補助関数を定義する:

ネオフックモデル

ネオフック材料モデルの材料特性を設定する:
パラメトリックソルバを作成する:
荷重 をゆっくりと増加させて変位について解く:
変形したメッシュを計算する:

方向への荷重と標本の形状のため,最もおもしろいひずみ効果はx方向の変形を示す のひずみ成分で観察することができる.同じ可視効果はEquivalentStrain関数を使って得ることができる.

ひずみ場を計算する:
牽引方向のひずみをプロットする:

ヨーモデル

ヨー材料モデルの材料特性を設定する:
パラメトリックソルバを作成する:
荷重 をゆっくり増加させることで変位について解く:
変形したメッシュを計算する:
ひずみ場を計算する:
牽引方向のひずみをプロットする:

パラメータ2つのムーニー・リブリンモデル

パラメータを2つ持つムーニー・リブリン材料モデルの材料特性を設定する:
パラメトリックソルバを作成する:
荷重 をゆっくり増加させることで変位について解く:
変形したメッシュを計算する:
ひずみ場を計算する:
牽引方向のひずみをプロットする:

パラメータ5つのムーニー・リブリンモデル

パラメータを5つ持つムーニー・リブリン材料モデルの材料特性を設定する:
パラメトリックソルバを作成する:
荷重 をゆっくり増加させることで変位について解く:
変形したメッシュを計算する:
ひずみ場を計算する:
牽引方向のひずみをプロットする:

Gentモデル

Gent材料モデルの材料特性を設定する:
パラメトリックソルバを作成する:
荷重 をゆっくり増加させることで変位について解く:
変形したメッシュを計算する:
ひずみ場を計算する:
牽引方向のひずみをプロットする:

Arruda-Boyce

Arruda-Boyce材料モデルの材料特性を設定する:
パラメトリックソルバを作成する:
荷重 をゆっくり増加させることで変位について解く:
変形したメッシュを計算する:
ひずみ場を計算する:
牽引方向のひずみをプロットする:

アニメーション

すべてのモデル名のリストを生成する:
すべての変位履歴のリストを生成する:
等高線プロットの上方境界ひずみと下方境界ひずみを求める:
境界を拡大して大きいプロットにする:
色付きマップのために最大ひずみをスケールする規則を定義する:
応力・ひずみプロットの境界を定義する:
各フレームで変形した標本を更新するための補助関数を定義する:
各モデルのグラフィックスを生成する補助関数を定義する:
各フレームについての行を生成する補助関数を作成する:
各フレームを生成する補助関数を定義する:
アニメーションでのフレーム数を定義する:
フレームのリストを生成する:
アニメーションを表示する:

まとめ

数値結果ではモデル間のいくつかの違いが示される.最も顕著な違いはフィットの範囲内と範囲外でモデルを比較すると現れる.例えばGentモデルとムーニーリブリンモデルはフィット範囲内では似ているが,5パラメータを持つムーニーリブリン,9パラメータを持つムーニーリブリン,ヨーの各モデルのように,フィット範囲外で発散し始める.

超弾性モデルはすべて非線形性を示すが,その強度は材料パラメータと同様にモデルによって異なる.挙動がほぼ直線として描かれるネオフックモデルやムーニーリブリンモデル等のモデルは,平均的な意味で実験データにはっきりフィットするが,逐一では解析的比較で示すように差は大きい.このような場合,所望のモデルは実験データ範囲だけではなく,アプリケーションで想定される応力ひずみ範囲にもに重点を置いてフィットする必要がある.

他の挙動の差はヨーのような非線形モデルよりもネオフックモデルのようなより線形モデルの標本の初期硬度に見られる.材料モデルには線形に近い領域と非常に非線形となる領域がある.例えば,ヨーモデルの場合,初期の変形がネオフックモデルの変形よりも大きいアニメーションでよく見ることができる.一旦モデルが非常に非線形な領域に近付くと,モデルの硬度が増すためヨーモデルの変形は小さくなる.ネオフックモデルは同じ応力に対して大きい変形を見せる.

実験手順で通常行われるように,これも力を変位条件としてではなく境界条件として適用する理由である.さらにこのような境界荷重では異なるモデル間の標本の異なる変形を得ることができる.応力の差についてのグラフィカルな外見によって,変形の差はもっと簡単に理解できる.

用語集

参考文献

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