有限歪弾性ブリック要素¶
このブリック要素は,大きな変形弾性のためのいくつかの古典的な超弾性構成則を実装します.
有限歪み弾性に関するいくつかの復習¶
\(\Omega\) を基準配置とし, \(\Omega_t\) を弾性体の変形配置とします.そして, \(X \in \Omega\) に対して次のように変形を表します \(\Phi(x) = u(X) + X\) .ベクトルフィールド \(u\) は,初期位置に対する変位です.
Cauchy-Greenテンソルは次のように定義されます.
変形勾配テンソル(Green-Lagrange)は次の通りです.
(線形弾性の場合, \({\nabla u^T}{\nabla u}\) は無視されます).
ここで次の通りであり
2つのテンソル \(E\) と \(C\) は,有限ひずみ弾性構成則を記述するために使用されます.
主要な不変量および導関数¶
有限ひずみ弾性構成則の記述は,しばしば,変形テンソルの主要な不変量を必要とします.
\(i_1,i_2,i_3\) は次数 \(1,2\) と \(3\) の不変量です.
\(H\) 方向のテンソル \(E\) の不変量の導関数 \(H\) は次のとおりです.
ゆえに
また再喝すると
不変量の2次導関数は,次の式によって定義される4次テンソルとなります.
表記法 \(A:B\) は Frobenius 積を表します. この積 \(A:B = \displaystyle\sum_{ij}A_{ij}B_{ij}\) には以下の特性があります.
また,
この性質により,Cauchy-Greenテンソル不変量の関数として,特に汎用化されたBlatz-Koひずみエネルギーの場合は,構成的規則を表現することが可能になります.
弾性ポテンシャルエネルギーとその導関数¶
基準配置における応力は,第2種Piola-Kirchhoff応力テンソル \({\hat{\hat{\sigma}}} = \nabla\Phi^{-1}\sigma\nabla\Phi^{-t}~\det \nabla\Phi\) によって記述することができます. ここで,\(\sigma\) は,変形された構成 \(\Omega_t\) のCauchy応力テンソルです.超弾性の構成則は次のように与えられます.
ここで, \({W}\) は材料のひずみエネルギーの密度です.全ひずみエネルギーは次式となります.
また, \(v\) 方向のエネルギーの導関数は次のように書くことができます.
これは次式によります.
Aが対称であるとき, \({\hat{\hat{\sigma}}}\) の場合には, \(A:B = A:(B+B^T)/2\) となります.
別の方法として,静的平衡を考慮するものがあり,これは基準配置において以下のように書くことができます.
部分積分により,次のようになります.
接線行列¶
変位 \(u\) は固定されています.接線行列を得るためには, \(u\) を \(u+h\) に置き換えます.
線形部分に関して \(h\) を考慮します.ここで
それは \(v\) と \(h\) に関して対称であり.これは次のように書き換えることができます.
ここで \(\mathcal{A}\) は対称であり \(3\times3\times3\times3\) テンソルは \(\mathcal{A}_{ijkl} = ((\frac{\partial^2 W}{\partial E^2})_{ijkl} + (\frac{\partial^2 W}{\partial E^2})_{ijlk})/2\) により与えられます.
いくつかの古典的な構成則¶
線形化された:Saint-Venant Kirchhoff則(微小変形)
¶
3つのパラメータ Mooney-Rivlin法
¶
圧縮可能な材料です.
ここで \(c_1\), \(c_2\) と \(d_1\) は係数が与えられ,次の式になります.
その後
非圧縮性の材料です.
非圧縮性制約 \(i_3( C) = 1\) はLagrange乗数 \(p\) (圧力)で扱われます.
拘束: \(\sigma = -pI \Rightarrow {\hat{\hat{\sigma}}} = -p\nabla\Phi\nabla\Phi^{-T}\det\nabla\Phi\)
Ciarlet-Geymonat 法
¶
Lame 係数 \(\lambda, \mu\) は \(\max(0,\frac{\mu}{2}-\frac{\lambda}{4})<a<\frac{\mu}{2}\) のようになります.([ciarlet1988] を参照してください).
一般化されたBlatz-Ko法
¶
\(\frac{\partial}{\partial C} {W}(C) = \displaystyle\sum_{j}\frac{\partial W}{\partial i_j(C)} \frac{\partial i_j(C)}{\partial C}\) と \(\frac{\partial^2}{\partial C^2} {W}(C) = \displaystyle\sum_{j} \displaystyle\sum_{k} \frac{\partial^2 W}{\partial i_j(C) \partial i_k(C)} \frac{\partial i_k(C)}{\partial C} \otimes \frac{\partial i_j(C)}{\partial C} + \displaystyle\sum_{j} \frac{\partial W}{\partial i_j(C)} \frac{\partial^2 i_j(C)}{\partial C^2}\) であるため,Cauchy-Greenテンソル不変量に対して歪エネルギー関数の導関数を計算しなければなりません( \(\frac{\partial i_j}{\partial E}(C;H) = 2 \frac{\partial i_j}{\partial C}(C;H)\) であるため,以下の点 \(E\) に関して不変量導関数を計算する必要はありません.)
平面ひずみ超弾性
¶
以上のすべてのモデルは,ボリューム領域で有効です.対応する平面歪み2次元モデルは,応力テンソルと4次テンソル \(\mathcal{A}\) をその平面成分に制限することで得られます.
モデルに非線形弾性ブリック要素を追加する¶
この要素は有限ひずみ弾性問題を表しています.ファイル getfem/getfem_nonlinear_elasticity.h
と getfem/getfem_nonlinear_elasticity.cc
で定義されています.このブリック要素をモデルに追加する関数は以下の通りです.
ind = getfem::add_nonlinear_elasticity_brick
(md, mim, varname, AHL, dataname, region = -1);
AHL
は,考えられる超弾性法則を表す getfem::abstract_hyperelastic_law
型のオブジェクトです.次の中から選択する必要があります.
getfem::SaintVenant_Kirchhoff_hyperelastic_law AHL;
getfem::Ciarlet_Geymonat_hyperelastic_law AHL;
getfem::Mooney_Rivlin_hyperelastic_law AHL(compressible, neohookean);
getfem::plane_strain_hyperelastic_law AHL(pAHL);
getfem::generalized_Blatz_Ko_hyperelastic_law AHL;
Saint-Venant Kirchhoff則は,2つのラメ係数で定義された線形化された法則であり,Ciarlet Geymonat則は2つのラメ係数と追加の係数(\(\lambda, \mu, a\))で定義されます.
Mooney-Rivlin則は,2つのオプションのフラグを受け取ります.最初のものは,材料が圧縮可能かどうかを判定し(\(d_1 \neq 0\)),2番目のフラグは材料がNeo-Hookean(\(c_2 = 0\))かどうかを判定します.これらのフラグに応じて,1〜3個の係数が必要な場合があります.デフォルトでは,それは非圧縮性かつNeo-Hookeanでないと定義されているため,2つの物質係数(\(c_1\), \(c_2\))が必要です.この場合,大きなひずみの非圧縮条件で使用する必要があります.
平面ひずみ超弾性則は,パラメータとして超弾性則の指針をとり,2D平面ひずみ近似を実行します.
md
はモデル変数, mim
は積分法, varname
は文字列に変数が追加される変数の名前, dataname
は法則の係数を表すモデル内のデータの名前です(定数または有限要素法を記述することができる)また region
はその項が考慮される領域です(デフォルトではすべてのメッシュ).
tests
ディレクトリ内のプログラム nonlinear_elastostatic.cc
と interface/tests/matlab
ディレクトリ内のプログラム demo_nonlinear_elasticity.m
はこの要素の非圧縮性条件での使用の例です.
新しい超弾性構成則の追加は,歪エネルギー,応力テンソルおよび応力テンソルの導関数の表現を与えることにあることに留意してください.詳細は,ファイル getfem/getfem_nonlinear_elasticity.cc
を参照してください.特に,不変量およびそれらの導関数の表現が利用可能です.
Von MisesまたはTrescaの応力を計算する関数も利用できます.
VM = compute_Von_Mises_or_Tresca
(md, varname, AHL, dataname, mf_vm, VM, tresca)
これは,有限要素法 mf_vm 上のVon MisesまたはTresca応力の自由度のベクトルを返します. tresca
はブール値で,Tresca応力に対しては true
,Von Mises応力には false
となります.
モデルに有限ひずみの非圧縮性ブリック要素を追加する¶
このブリック要素は,大きなひずみ問題で非圧縮性条件を追加します
圧力を表すLagrange乗数が混合式で導入されます.このブリック要素をモデルに追加する関数は次の通りです.
ind = add_nonlinear_incompressibility_brick
(md, mim, varname, multname, region = -1)
ここで, md
はモデル, mim
は積分法, varname
は非圧縮性条件が追加されたモデルの変数, multanme
は圧力に対応する乗数変数(変数の有限要素法と乗数の間で少なくとも線形のLadyzhenskaja-Babuska-Brezzi極限条件が満たされていることに注意してください.) region
は,項が考慮されるメッシュ領域に対応するオプションのパラメータです(デフォルトではすべてのメッシュ).
高水準汎用構築バージョン¶
汎用弱形式言語 (GWFL) は, GetFEM で実装された超弾性ポテンシャルおよび構成規則にアクセスできるようにします.これにより,言語で直接使用することができます.たとえば,モデル内の汎用的な構築ブリック要素を使用するか,または特定の数量の補間(応力など)が可能です.
非線形弾性に役立つ非線形演算子の言語のリストです.
Det(M) % determinant of the matrix M
Trace(M) % trace of the matrix M
Matrix_i2(M) % second invariant of M (in 3D): (sqr(Trace(m)) - Trace(m*m))/2
Matrix_j1(M) % modified first invariant of M: Trace(m)pow(Det(m),-1/3).
Matrix_j2(M) % modified second invariant of M: Matrix_I2(m)*pow(Det(m),-2/3).
Right_Cauchy_Green(F) % F' * F
Left_Cauchy_Green(F) % F * F'
Green_Lagrangian(F) % (F'F - Id(meshdim))/2
Cauchy_stress_from_PK2(sigma, Grad_u) % (Id+Grad_u)*sigma*(I+Grad_u')/det(I+Grad_u)
ポテンシャルのリストです.
Saint_Venant_Kirchhoff_potential(Grad_u, [lambda; mu])
Plane_Strain_Saint_Venant_Kirchhoff_potential(Grad_u, [lambda; mu])
Generalized_Blatz_Ko_potential(Grad_u, [a;b;c;d;n])
Plane_Strain_Generalized_Blatz_Ko_potential(Grad_u, [a;b;c;d;n])
Ciarlet_Geymonat_potential(Grad_u, [lambda;mu;a])
Plane_Strain_Ciarlet_Geymonat_potential(Grad_u, [lambda;mu;a])
Incompressible_Mooney_Rivlin_potential(Grad_u, [c1;c2])
Plane_Strain_Incompressible_Mooney_Rivlin_potential(Grad_u, [c1;c2])
Compressible_Mooney_Rivlin_potential(Grad_u, [c1;c2;d1])
Plane_Strain_Compressible_Mooney_Rivlin_potential(Grad_u, [c1;c2;d1])
Incompressible_Neo_Hookean_potential(Grad_u, [c1])
Plane_Strain_Incompressible_Neo_Hookean_potential(Grad_u, [c1])
Compressible_Neo_Hookean_potential(Grad_u, [c1;d1])
Plane_Strain_Compressible_Neo_Hookean_potential(Grad_u, [c1;d1])
Compressible_Neo_Hookean_Bonet_potential(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Bonet_potential(Grad_u, [lambda;mu])
Compressible_Neo_Hookean_Ciarlet_potential(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Ciarlet_potential(Grad_u, [lambda;mu])
第2種Piola-Kirchhoff応力テンソルのリストです.
Saint_Venant_Kirchhoff_PK2(Grad_u, [lambda; mu])
Plane_Strain_Saint_Venant_Kirchhoff_PK2(Grad_u, [lambda; mu])
Generalized_Blatz_Ko_PK2(Grad_u, [a;b;c;d;n])
Plane_Strain_Generalized_Blatz_Ko_PK2(Grad_u, [a;b;c;d;n])
Ciarlet_Geymonat_PK2(Grad_u, [lambda;mu;a])
Plane_Strain_Ciarlet_Geymonat_PK2(Grad_u, [lambda;mu;a])
Incompressible_Mooney_Rivlin_PK2(Grad_u, [c1;c2])
Plane_Strain_Incompressible_Mooney_Rivlin_PK2(Grad_u, [c1;c2])
Compressible_Mooney_Rivlin_PK2(Grad_u, [c1;c2;d1])
Plane_Strain_Compressible_Mooney_Rivlin_PK2(Grad_u, [c1;c2;d1])
Incompressible_Neo_Hookean_PK2(Grad_u, [c1])
Plane_Strain_Incompressible_Neo_Hookean_PK2(Grad_u, [c1])
Compressible_Neo_Hookean_PK2(Grad_u, [c1;d1])
Plane_Strain_Compressible_Neo_Hookean_PK2(Grad_u, [c1;d1])
Compressible_Neo_Hookean_Bonet_PK2(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Bonet_PK2(Grad_u, [lambda;mu])
Compressible_Neo_Hookean_Ciarlet_PK2(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Ciarlet_PK2(Grad_u, [lambda;mu])
材料パラメータに関する導関数は,Saint Venant Kirchhoff超弾性法のためには個別に実装されていないことに注意してください.したがって,パラメータをモデルの他の変数に依存させることはできません(導関数は実装するのはそれほど複雑ではありませんが,現時点では古い実装のラッパーのみが記述されています).
モデルの結合は,弱形式レベルで行われることに注意してください.汎用的な方法では,ポテンシャルで問題を定義しないことをお勧めします.第1に,ポテンシャルの2次導関数(接線系を得るために必要)が非常に複雑で最適化不可能であり,第2に,ポテンシャルレベルで主結合を得ることができないという2つの理由があります.従って,ポテンシャルの使用は,ポテンシャルの実際の計算に限定されるべきです.
モデルまたはga_workspaceの変数 u
にSaint Venant-Kirchhoff超弾性項を追加するには,次の構築文字列を使用します.
"((Id(meshdim)+Grad_u)*(Saint_Venant_Kirchhoff_PK2(Grad_u,[lambda;mu]))):Grad_Test_u"
その場合,データ lambda
と mu
を model/ga_workspace で宣言する必要があることに注意してください.もちろん,いくつかのデータに応じて陽な定数や式で置き換えることも可能です.
非圧縮性のMooney-Rivlin則に関しては,非圧縮性の項を完成されなければなりません.たとえば,次の非圧縮ブリック要素について,
ind = add_finite_strain_incompressibility_brick(md, mim, varname, multname, region = -1);
このブリック要素は, p
が乗数で u
変位を表す変数であれば項 p*(1-Det(Id(meshdim)+Grad_u))
を与える単純なものです.
モデルに超弾性項を加えることは,次の関数によっても行うことができます.
ind = add_finite_strain_elasticity_brick(md, mim, lawname, varname, params,
region = size_type(-1));
ここで, md
はモデル, mim
は積分法, varname
は大きなひずみ変位を表すモデルの変数, lawname
は構成則名で Saint_Venant_Kirchhoff
, Generalized_Blatz_Ko
, Ciarlet_Geymonat
, Incompressible_Mooney_Rivlin
, Compressible_Mooney_Rivlin
, Incompressible_Neo_Hookean
, Compressible_Neo_Hookean
, Compressible_Neo_Hookean_Bonet
または Compressible_Neo_Hookean_Ciarlet
のうちのどれかです.params
は小ベクトルまたはベクトルフィールドとして定義された構成則のパラメータを表す文字列です.
Von Mises応力は次の関数で補間することができます
void compute_finite_strain_elasticity_Von_Mises(md, varname, lawname, params, mf_vm, VM,
rg=mesh_region::all_convexes());
ここで, md
はモデル, varname
は有限ひずみ変位を表すモデルの変数, lawname
は構成則名(前述のブリック要素を参照), params
は文字列補間が行われる法則のパラメータ ,mf_vm
は(不連続なのが好ましい)Lagrange有限要素法を表現し, VM
は補間された値が格納される model_real_plain_vector
型のベクトルを表します.