Total Lagrange

Total-Lagrange法による幾何学的非線形弾性問題の有限要素法解析

Last Update:2009年05月15日

支配方程式の弱形式化

公称応力を用いた平衡方程式

運動量保存の法則から任意の領域について

$\left.\frac{\partial}{\partial{t}}\right|_{\b{X}}\int_v \rho {\b{v}} dv = \int_v \rho {\b{g}} dv + \int_s {\b{t}} ds$

が成り立つ。 dv=JdV と、変形後の表面力 $\b{t}$ に表面の変化率をかけた表面力 $\tilde{\b{t}}=\frac{ds}{dS}\b{t}$ を用いて、次のように領域を変形前の領域に戻すことができる。

$\left.\frac{\partial}{\partial{t}}\right|_{\b{X}}\int_V \rho {\b{v}} JdV = \int_V \rho {\b{g}} Jdv + \int_S \tilde{\b{t}} dS$

また、 $\rho J=\rho_0$ と $\Pi^T\cdot\b{N}=\tilde{t}$ で定義される、第一Piola-Kirhihoffテンソル（公称応力） $\Pi$ を用いて、

$\left.\frac{\partial}{\partial{t}}\right|_{\b{X}}\int_V \rho_0 {\b{v}} dV = \int_V \rho_0 {\b{g}} dv + \int_S N\cdot\Pi dS$

第１項では積分領域が時間によって変化しなく、また $\rho_0$ も時間不変であるから、物質微分を積分の中にいれることが可能である。また、第３項にガウスの発散定理を適応して、

$\int_V\rho_0\left.\frac{\partial\b{v}}{\partial{t}}\right|_{\b{X}} dV = \int_V \rho_0 {\b{g}} dv + \int_V \nabla_X \cdot\Pi dV$

となる。この式が任意の領域において成り立つので、積分を取り除くことができる。

公称応力を用いた平衡方程式
$\rho_0\left.\frac{\partial\b{v}}{\partial{t}}\right\|_{\b{X}} = \rho_0 {\b{g}} + \nabla_X \cdot\Pi$

弱形式化した平衡方程式

上の平衡方程式を弱形式化する。

上式の両辺に固定境界条件が設定されている面 S_1 上で０となる任意の関数 $\delta \b{u}$ をかけて領域内で積分すると、

$\int_V \delta \b{u} \cdot \rho_0\left.\frac{\partial\b{v}}{\partial{t}}\right|_{\b{X}} dV = \int_V \delta\b{u}\cdot\rho_0 {\b{g}}dV + \int_V\delta \b{u}\cdot(\nabla_X \cdot\Pi)dV$

第三項について変形を行う。

$\delta \b{u}\cdot(\nabla_X \cdot\Pi)=\delta u_j\frac{\partial\Pi_{ij}}{\partial X_i}=\frac{\partial}{\partial X_i}(\delta u_j\Pi_{ij})-\frac{\partial\delta u_j}{\partial X_i}\Pi_{ij}=\nabla_X\cdot(\Pi\cdot\delta \b{u})-\Pi^T:\delta F$

これを上式に代入して、

$\int_V \delta \b{u} \cdot \rho_0\left.\frac{\partial\b{v}}{\partial{t}}\right|_{\b{X}} dV = \int_V \delta\b{u}\cdot\rho_0 {\b{g}}dV + \int_V \nabla_X\cdot(\Pi\cdot\delta \b{u})dV-\int_V\Pi^T:\delta F dV$

右辺第二項にガウスの発散定理を適応して、

$\int_V \delta \b{u} \cdot \rho_0\left.\frac{\partial\b{v}}{\partial{t}}\right|_{\b{X}} dV = \int_V \delta\b{u}\cdot\rho_0 {\b{g}}dV + \int_S N\cdot(\Pi\cdot\delta \b{u})dS-\int_V\Pi^T:\delta F dV$

右辺第二項は次の関係を用いる

$N\cdot\Pi=\tilde{\b{t}}$

次のような弱形式化した平衡方程式が得られる。

弱形式化した平衡方程式２（第一Piola-Kirhihoff応力）
$\int_V \delta \b{u} \cdot \rho_0\left.\frac{\partial\b{v}}{\partial{t}}\right\|_{\b{X}} dV+\int_V\Pi^T:\delta F dV = \int_V \delta\b{u}\cdot\rho_0 {\b{g}}dV + \int_S \tilde{\b{t}}\cdot\delta \b{u}dS$

右辺第三項について次のような変換を行う。

$\b{\Pi}^T : \delta \b{F} = (\b{S}\cdot\b{F}^T)^T : \delta \b{F} = \tr(\b{S}\cdot\b{F}^T\cdot\delta \b{F} )\\ \qquad\qquad\qquad\qquad= \tr$\b{S}\cdot\frac{1}{2}(\b{F}^T\cdot\delta \b{F}+\delta\b{F}^T\cdot\b{F}) $ = \tr$\b{S}\cdot\delta(\frac{1}{2}\b{F}^T\cdot\b{F}) $ = \b{S}:\delta(\frac{1}{2}\b{F}^T\cdot\b{F}) = \b{S}:\delta\b{E}$

弱形式化した平衡方程式は次のように書くことができる。

弱形式化した平衡方程式２（第二Piola-Kirhihoff応力）
$\int_V \delta \b{u} \cdot \rho_0\left.\frac{\partial\b{v}}{\partial{t}}\right\|_{\b{X}} dV+\int_V\b{S}:\delta\b{E}dV = \int_V \delta\b{u}\cdot\rho_0 {\b{g}}dV + \int_S \tilde{\b{t}}\cdot\delta \b{u}dS$

の弱形式化した平衡方程式が２つ出てきたが、 $S,\delta E$ は対称であるので後の形式で弱形式化した平衡方程式のほうが簡単であることが多い。Total-Lagrange法では通常後の弱形式化した平衡方程式を使う。

支配方程式の有限要素法離散化

$\int_V \b{S}:\delta\b{E}dV=Q^a_i\delta u^a_i$

のように表されるとき Q^a_i を内力ベクトルと呼ぶ

以下内力ベクトルを求める。

$E_{ij}=\frac{1}{2}(\frac{\partial u_i}{\partial X_j}+\frac{\partial u_j}{\partial X_i}+\frac{\partial u_h}{\partial X_i}\frac{\partial u_h}{\partial X_j})$

$\delta E_{ij}=\frac{1}{2}(\frac{\partial \delta u_i}{\partial X_j}+\frac{\partial \delta u_j}{\partial X_i}+\frac{\partial \delta u_h}{\partial X_i}\frac{\partial u_h}{\partial X_j}+\frac{\partial u_h}{\partial X_i}\frac{\partial \delta u_h}{\partial X_j})$

非対称化した仮想歪テンソルの導入

の対称性を生かして計算を簡単にすることができる。

ここで $\delta \bar{E}$ を次のように定義する。

$\delta \bar{E}=\b{F}^T\cdot\delta\b{F}$

明らかに、 $\delta \bar{\b{E}}$ の対称成分は $\delta \b{E}$ である。これをテンソルとして書くと、

$\delta \b{E}=\frac{1}{2}(\delta\bar{\b{E}}+\delta\bar{\b{E}}^T)$

つまり、ここで第二Piola-Kirhihoff応力 $\b{S}$ の対称性とテンソルの内積の性質を利用して、

$\b{S}:\delta\b{E}=\b{S}:\frac{1}{2}(\delta\bar{\b{E}}+\delta\bar{\b{E}}^T)=\frac{1}{2}(\b{S}:\delta\bar{\b{E}}+\b{S}:\delta\bar{\b{E}}^T) = \frac{1}{2}(\b{S}:\delta\bar{\b{E}}+\b{S}^T:\delta\bar{\b{E}}^T) = \frac{1}{2}(\b{S}:\delta\bar{\b{E}}+\b{S}:\delta\bar{\b{E}})=\b{S}:\delta\bar{\b{E}}$

$\b{S}:\delta\bar{\b{E}}$ の中の項の数は $\b{S}:\delta\b{E}$ の数の半分である。 $\delta\bar{\b{E}}$ を用いることで計算が簡単になったことがわかる。

内挿関数の導入

$\delta\bar{E}$ を成分で書くと次のとおり

$\delta \bar{E}_{ij}=F_{hi}\delta F_{hj}=(\delta_{hi}+Z_{hi})F_{hj}=F_{ij}+F_{hj}Z_{hi}=\frac{\partial \delta u_i}{\partial X_j} + \frac{\partial \delta u_h}{\partial X_j}\frac{\partial u_h}{\partial X_i}$

ここで $\delta u$ の添え字がになるように添え字を入れえると次のようになる。

$\delta \bar{E}_{gh}=\frac{\partial \delta u_i}{\partial X_h}\delta_{ig} + \frac{\partial \delta u_i}{\partial X_h}\frac{\partial u_i}{\partial X_g}$

ここで $\delta u$ ,が次のように形状関数Nを用いて要素 V_e 内で離散化されているとする。

$\delta u_i= N^p \delta u^p_i$ 、 u_i=N^r u^r_i

$\delta \bar{E}_{gh} = \frac{\partial N^p \delta u^p_i}{\partial X_h}\delta_{ig} + \frac{\partial N^p \delta u^p_i}{\partial X_h}\frac{\partial N^r u^r_i}{\partial X_g}=\[B_{gh}\]^p_i \delta u^p_i$

ここで $\[B_{gh}\]^p_i$ は変位と歪を関係づける行列であり、

$\[B_{gh}\]^p_i=\frac{\partial N^p}{\partial X_h}\delta_{ig} + \frac{\partial N^p}{\partial X_h}\frac{\partial N^r}{\partial X_g} u^r_i=\frac{\partial N^p}{\partial X_h}(\delta_{ig} + \frac{\partial u_i}{\partial X_g})=\frac{\partial N^p}{\partial X_h}F_{ig}$

Total Lagarange法における変位と歪を対応づける行列
$\[B_{gh}\]^p_i=\frac{\partial N^p}{\partial X_h}\delta_{ig} + \frac{\partial N^p}{\partial X_h}\frac{\partial N^r}{\partial X_g} u^r_i=\frac{\partial N^p}{\partial X_h}F_{ig}$

とおいた。これを用いると要素内内力は次のようにかける。

$\{Q_e\}^p_i=\int_V_eS_{gh}\[B_{gh}\]^p_idV$

内力はこれを要素ごとに足し合わせたものであり、次のようになる。

内力ベクトル
${Q}^a_i=\sum_e\{Q_e\}^p_i$ (但し要素内節点番号pの節点番号はa) $\{Q_e\}^p_i=\int_V_eS_{gh}\[B_{gh}\]^p_idV$

$\delta E_{gh}=\frac{\partial E_{gh}}{\partial u^p_i}\delta u^p_i=\frac{1}{2}$\[B_{gh}\]^p_i+\[B_{hg}\]^p_i$\delta u^p_i$

ここで、 $\delta u^p_i$ は任意だから消去して次の関係を得る。

$\frac{\partial E_{gh}}{\partial u_i^p}=\frac{1}{2}$\[B_{gh}\]^p_i+\[B_{hg}\]^p_i$$

接線剛性行列

F=Q を満たすような解を求めなければならないが、非線形有限要素法の場合、 ${Q}^a_i$ は変位に対して非線形である。また ${F}^a_i$ も非線形である場合が多い

1次以上の連立方程式は計算機で解くことは不可能であるので、この非線形連立方程式を線形化して反復的に解くことで近似的にこの方程式を解く。この線形化のやり方には数々の解き方があるがもっとも良く使われるのがNewton-Raphson法である。

Newton-Raphson法

上の関係は ${Q}^a_i={F}^a_i$ となる。

ここで次のように残差ベクトル ${R}^p_i$ を定義する。

${F}^a_i-{Q}^a_i={R}^a_i$

$\b{Q}$ 、 $\b{F}$ は節点変位 u^a_j の関数であるから、残差についても同じであるので、残差は ${R}^a_i({u}^b_j)$ のように書ける。

ここで、残差が０つまり ${R}^a_i({u}^b_j)=0$ の時、 ${u}^b_j$ は解の一つである。

さて、ある解となる変位 ${u_n}^b_j$ の周りで外力を1次の項までTayler展開したとすると、

${R}^a_i({u}^b_j)\sime{R}^a_i({u_n}^b_j)+\frac{\partial {R}^a_i({u_n}^b_j)}{\partial u^b_j}\Delta u^b_j$

と近似できる。

この近似から求められる、残差が０となるような変位 ${u_{n+1}}^a_j$ は上式の右辺を０として

${u_{n+1}}^a_j={u_n}^a_j+\Delta {u_n}^a_j$

ここで増分量 $\Delta {u_n}^a_j$ は上式の右辺を０とした式から求められ、連立一次方程式

${R}^a_i({u_n}^b_j)+\frac{\partial {R}^a_i({u_n}^b_j)}{\partial u^b_j}\Delta {u_n}^b_j=0$ つまり、 $\[-\frac{\partial {R}^a_i({u_n}^b_j)}{\partial u^b_j}\]\Delta {u_n}^b_j={R}^a_i({u_n}^b_j)$

を解くことで得られる。

さて、変位 ${u_n}^b_j$ が解に十分近ければTaylor展開した２次以降の項の影響は小さくなり、 ${u_{n+1}}^b_j$ は解のよい近似になっているので、 ${u_{n+1}}^a_j$ は ${u_n}^a_j$ よりも解に近づいていると考えられる。

ここで再び ${u_{n+1}}^b_j$ の周りで同じように1次のTaylor展開を行い解をさらに近づけることができる。反復的にこのような操作を行うことで必要とされる精度まで解を計算することをNewton-Raphson法と呼ぶ

Newton-Raphson法のアルゴリズム

接線剛性行列

$[K_e]^{p,q}_{i,j} =\frac{\partial \{Q_e\}^p_i}{\partial u^q_j}$

$\frac{\partial \{Q_e\}^p_i}{\partial u^q_j} = \int_V_e\frac{\partial S_{gh}}{\partial u^q_j}\[B_{gh}\]^p_i + S_{gh}\frac{\partial \[B_{gh}\]^p_i}{\partial u^q_j}dV=[{K_L}_e]^{p,q}_{i,j}+[{K_{NL}}_e]^{p,q}_{i,j}$

3式の第一項は歪を固定した時の応力の変位による増分に対する剛性であり、初期変位行列(initial displacement matrix)、
３式の第二項は応力を固定した時の歪の変位にようる増分に対する剛性であり、初期応力行列(initial stress matrix)、あるいは幾何剛性行列(geometric stiffness matrix)と呼ばれる。

初期応力行列

$[{K_{NL}}_e]^{p,q}_{i,j} = \int_V_e S_{gh}\frac{\partial \[B_{gh}\]^p_i}{\partial u^q_j}dV$

を求めてみよう.

歪-変位関係行列 $\[B_{gh}\]^p_i$ は

$\[B_{gh}\]^p_i=\frac{\partial N^p}{\partial X_h}\delta_{ig} + \frac{\partial N^p}{\partial X_h}\frac{\partial N^r}{\partial X_g} u^r_i$

であった、これを節点変位 u^q_j について微分をとってみると、第一項は消えて次のようになる。

$\frac{\partial \[B_{gh}\]^p_i}{\partial u^q_j} = \frac{\partial N^p}{\partial X_h}\frac{\partial N^r}{\partial X_g}\delta_{ij}\delta_{rq} = \frac{\partial N^p}{\partial X_h}\frac{\partial N^q}{\partial X_g}\delta_{ij}$

よって初期応力行列は次のとおり

$[{K_{NL}}_e]^{p,q}_{i,j}=\int_V_eS_{gh}\frac{\partial N^p}{\partial X_h}\frac{\partial N^q}{\partial X_g}\delta_{ij}dV$

初期変位行列

$[{K_L}_e]^{p,q}_{i,j} = \int_V_e\frac{\partial S_{gh}}{\partial u^q_j}\[B_{gh}\]^p_idV$

を求めてみよう．

まず、 $\frac{\partial S_{gh}}{\partial u^q_j}$ を計算する。

$d\b{S}=\b{C}:d\b{E}$

の関係があるとする。これを成分で書くと、

$dS_{gh}=C_{ghef}dE_{ef}$

ここで応力は歪の関数つまり変位の関数であるから,に対する全微分は変位についての偏微分により、次のように書ける。

$\frac{\partial S_{gh}}{\partial u^q_j}du^q_j=C_{ghef}\frac{\partial E_{ef}}{\partial u^q_j}du^q_j$

これが任意の du^q_j について成り立つから、両辺から消去して、

$\frac{\partial S_{gh}}{\partial u^q_j}=C_{ghef}\frac{\partial E_{ef}}{\partial u^q_j}$

である。これを更に変更して変位と関係づけよう

$\frac{\partial S_{gh}}{\partial u^q_j}=C_{ghef}\frac{\partial E_{ef}}{\partial u^q_j}=C_{ghef}\frac{1}{2}$\[B_{ef}\]^q_j+\[B_{fe}\]^q_j$=\frac{1}{2}C_{ghef}\[B_{ef}\]^q_j + \frac{1}{2}C_{ghfe}\[B_{ef}\]^q_j\\ \qquad\qquad\qquad =\frac{1}{2}(C_{ghef}+C_{ghfe})\[B_{ef}\]^q_j = \bar{C}_{ghef}\[B_{ef}\]^q_j$

第３式から第４式への変形は第二項について添え字 e,f を入れ替えた。また、

$\bar{C}_{ghef}=\frac{1}{2}(C_{ghef}+C_{ghfe})$

とおいた。と $\bar{C}$ は基本的に別物で両者を取り違えてよくバグの原因になるので注意が必要である。

これらを用いると初期変位行列は次のようになる。

$[{K_L}_e]^{p,q}_{i,j}=\int_V_e\frac{\partial S_{gh}}{\partial u^q_j}\[B_{gh}\]^p_idV=\int_V_e\bar{C}_{ghef}\[B_{ef}\]^q_j\[B_{gh}\]^p_idV$

接線剛性行列

以上をまとめると次のようになる。

接線剛性行列は以下のように要素接線剛性行列の和として計算される。

接線剛性行列
$[K]^{a,b}_{i,j}=\sum_e[K_e]^{p,q}_{i,j}$ (但し、はの要素内節点番号) $[K_e]^{p,q}_{i,j}=\int_V_e\bar{C}_{ghef}\[B_{ef}\]^q_j\[B_{gh}\]^p_idV + \int_V_eS_{gh}\frac{\partial N^p}{\partial X_h}\frac{\partial N^q}{\partial X_g}\delta_{ij}dV$

実装に関して

要素剛性行列、要素内残差ベクトルを作るプログラムの流れ

1次補間の要素を使っても被積分関数が複雑なのでよっぽど高速化に興味がない限り解析積分を使わずに数値積分を用いる。

変位勾配テンソルを求める
変位勾配テンソルからGreen-Lagrange歪を求める
Green-Lagrange歪から物質の構成式を用いて第二Piola-Kirchhoff応力を求める
変位勾配テンソルから歪-変位関係行列を求める
第二Piola-Kirhihoff応力と歪-変位関係行列から要素節点等価内力ベクトル $\{Q_e\}$ を求める
物質の構成式から求まる $\bar{C}$ と第二Piola-Kirchhoff応力と歪-変位関係行列から要素剛性行列を得る。
要素節点等価外力ベクトル $\{F_e\}$ から要素節点等価外力ベクトル $\{Q_e\}$ を引くことで要素節点等価残差ベクトル $\{R_e\}$ を得る。

変位勾配テンソルを求めるルーティン

$Z=\b{u}\otimes\nabla=\frac{\partial u_i}{\partial X_j}\b{e}_i\otimes\b{e}_j$

$Z_{ij}=\frac{\partial u_i}{\partial X_j}=\frac{\partial N^a}{\partial X_j}u^a_i$

ここで、次のようにプログラムの中で変数が対応づけられているとすると、

ndim:空間の次元の数、nnoel:要素あたりの節点数
dudx[i][j]= $Z_{ij}$ 、edisp[a][i] = 、dndx[a][i] = $\frac{\partial N^a}{\partial X_i}$

変位勾配テンソルを求めるルーティンは次のとおり

double dudx[ndim][ndim];
{  // 変位勾配テンソルを求めるルーティン
   for(unsigned int i=0;i<ndim*ndim;i++){ *(&dudx[0][0]+i)=0.0; }
   for(unsigned int inoel=0;inoel<nnoel;inoel++){
      for(unsigned int idim=0;idim<ndim;idim++){
      for(unsigned int jdim=0;jdim<ndim;jdim++){
         dudx[idim][jdim] += edisp[inoel][idim]*dndx[inoel][jdim];
      }
      }
   }
}

変位勾配テンソルからGreen-Lagrange歪を求めるルーティンの実装例

$E=\frac{1}{2}(Z+Z^T+Z^T\cdot Z)$ である。

つまり、 $E_{ij}=\frac{1}{2}(Z_{ij}+Z_{ji}+Z_{ki}Z_{kj})$

ここでプログラムの中で次のように変数が対応付けれているとすると、

ndim:空間の次元の数、nnoel:要素あたりの節点数
dudx[i][j]= $Z_{ij}$ 、strain[i][j]= $E_{ij}$

変位勾配テンソルからGreen-Lagrange歪をもとめるルーティンの実装例は以下のとおり

double strain[ndim][ndim];
for(unsigned int idim=0;idim<ndim;idim++){
for(unsigned int jdim=0;jdim<ndim;jdim++){
   strain[idim][jdim] = 0.5*( dudx[idim][jdim] + dudx[jdim][idim] );
   for(unsigned int kdim=0;kdim<ndim;kdim++){
      strain[idim][jdim] += 0.5*dudx[kdim][idim]*dudx[kdim][jdim];
   }
}
}

歪-変位関係行列を作るルーティンの実装例

$F_{ij}=Z_{ij}+\delta_{ij}$

ここでプログラム中で次のように変数が関係付けられているとすると、

ndim:空間の次元の数、nnoel:要素あたりの節点数
dudx[i][j]= $Z_{ij}$ 、f_mat[i][j]= $F_{ij}$ 、disp2stain[p][i][g][h]= $\[B_{gh}\]^p_i$ 、dndx[a][i] = $\frac{\partial N^a}{\partial X_i}$

歪-変位関係行列を作るルーティンの実装例は次のようになる。

double disp2strain[nnoel][ndim][ndim][ndim];
{  // 歪-変位関係行列を作るルーティン
   double f_mat[ndim][ndim];   // 変形勾配テンソル
   // 変形勾配テンソルFを求める
   for(unsigned int idim=0;idim<ndim;idim++){
      for(unsigned int jdim=0;jdim<ndim;jdim++){
         f_mat[idim][jdim] = dudx[idim][jdim];
      }
      f_mat[idim][idim] += 1.0;
   }
   // 変形勾配テンソルから歪-変位関係行列を作る
   for(unsigned int pnoel=0;pnoel<nnoel;pnoel++){
   for(unsigned int idim=0;idim<ndim;idim++){
      for(unsigned int gdim=0;gdim<ndim;gdim++){
      for(unsigned int hdim=0;hdim<ndim;hdim++){
         disp2strain[pnoel][idim][gdim][hdim] = dndx[pnoel][hdim]*f_mat[idim][gdim];
      }
      }
   }
   }
}

第二Piola-Kirhihoff応力と歪-変位関係行列から要素節点等価内力ベクトル $\{Q_e\}$ を求めるルーティンの実装例

$\{Q_e\}^p_i=\int_V_eS_{gh}\[B_{gh}\]^p_idV$

これは数値積分を用いて次のように書くことができる。

$\{Q_e\}^p_i\sime\sum_{\alpha}\omega_{\alpha}S_{gh}\[B_{gh}\]^p_iJ_{\alpha}$

ここで、次のようにプログラムの中で変数が対応づけられているとすると、

ndim:空間の次元の数、nnoel:要素あたりの節点数
detwei= $\omega_{\alpha}J_{\alpha}$
eforce_in[p][i]= $\{Q_e\}^p_i$ 、stress[g][h]= $S_{gh}$ 、

$\[B_{gh}\]^p_i$

要素節点等価内力ベクトルを求めるルーティンの実装例は以下のとおり

// 内力を求める
for(unsigned int pnoel=0;pnoel<nnoel;pnoel++){
for(unsigned int idim=0;idim<ndim;idim++){
   for(unsigned int gdim=0;gdim<ndim;gdim++){
   for(unsigned int hdim=0;hdim<ndim;hdim++){
      eforce_in[pnoel][idim] += detwei*stress[gdim][hdim]*disp2strain[pnoel][idim][gdim][hdim];
   }
   }
}
}

埋め込み座標を用いた定式化

埋め込み座標系における歪み

埋め込み座標を用いて、Total Lagrange法の定式化をしてみよう．

点の曲座標 r_i に沿ったベクトルを G_i とする．また、点が変形後にの位置にあるとして、曲座標 r_i に沿った接ベクトルを g_i とする．つまり、

$G_i=\frac{\partial X}{\partial r_i}$ 、 $\b{g}_i=\frac{\partial x}{\partial r_i}$

となる．これらは共変基底ベクトル(covariant base vector)と呼ばれる． G_i が変形後 g_i となったので、変形勾配テンソルを用いると、

$\b{g}_i=F\cdot G_i$

が成り立つ.よって、

埋め込み座標系での変形勾配テンソル
$F=\b{g}_i\otimes G^i$

となることがわかる．ここで G^i は反変基底ベクトル(contravariant base vector)と呼ばれる、

$G^i\cdot G_j = \delta^i_j$

を満たすベクトルである．これらを用いると、

$C=F^T\cdot F=(G^j\otimes g_j)\cdot(g_i\otimes G^i)=$\b{g}_i\cdot\b{g}_j$G^i\otimes G^j=g_{ij}G^i\otimes G^j$

$I=(G_i\cdot I\cdot G_j)G^i\otimes G^j=G_{ij}G^i\otimes G^j$

$E=\frac{1}{2}(C-I)=\frac{1}{2}(g_{ij}-G_{ij})G^i\otimes G^j$

とかける

埋め込み座標系でのGreen-Lagrange歪
$E=\frac{1}{2}(g_{ij}-G_{ij})G^i\otimes G^j$

仮想歪みエネルギー式

埋め込み座標を使って、仮想歪みエネルギー式を具体的に計算してみよう．

第二Piola-Kirchhoff応力と、Green-Lagrange歪みの変分を埋め込み座標系で成分表記すると次のようになる．

$S=S^{ij}G_i\otimes G_j$

$\delta E=\delta E_{ij}G^i\otimes G^j$

成分表示された応力と仮想歪みを使った仮想歪みエネルギーは以下のとおり．

$S:\delta E=(S^{ij}G_i\otimes G_j) : (\delta E_{kl}G^k\otimes G^l) = S^{ij}\delta E_{kl} (G_i\cdot G^k)(G_j\cdot G^l)=S^{ij}\delta E_{ij}$

$\delta E = \frac{1}{2}\delta g_{ij}G^i\otimes G^j = \frac{1}{2}$\delta \b{g}_i\cdot \b{g}_j+\b{g}_i\cdot\delta \b{g}_j$G^i\otimes G^j = \frac{1}{2}$\frac{\partial\delta\b{u}}{\partial r_i}\cdot \b{g}_j+\b{g}_i\cdot\frac{\partial\delta\b{u}}{\partial r_j}$G^i\otimes G^j$

ここで、次のような非対称化した仮想歪み $\delta\bar{E}$ を考える．

$\delta\bar{E} = $\b{g}_i\cdot\frac{\partial\delta\b{u}}{\partial r_j}$G^i\otimes G^j$

これを用いると、

$S:\delta E=S:\delta\bar{E}$

となる．ここで、変位が u_i=N^p u^p_i のように離散化されているとすると、

$\delta\bar{E}_{gh} = \b{g}_g\cdot\frac{\partial\delta\b{u}}{\partial r_h} = \b{g}_g_i\frac{\partial\delta \b{u}_i}{\partial r_h} = \b{g}_g_i\frac{\partial N^p \delta u^p_i}{\partial r_h} = $\b{g}_g_i\frac{\partial N^p}{\partial r_h}$ \delta u^p_i = \[B_{gh}\]^p_i \delta u^p_i$

のように書ける．ここで、 $\[B_{gh}\]^p_i$ は変位と歪みを関係づけるマトリックスで

変位と歪みを関係づける行列
$\[B_{gh}\]^p_i = \b{g}_g_i\frac{\partial N^p}{\partial r_h}$

である．これを用いると要素あたりの仮想歪みエネルギーの式は

$\int_{V_e} S:\delta E dV = \int_{V_e}S:\delta\bar{E} dV = \int_{V_e}S^{gh}\[B_{gh}\]^p_i dV u^p_i$

となるので、節点等価内力 $\{Q_e\}^p_i$ は

節点等価内力
$\{Q_e\}^p_i = \int_{V_e}S^{gh}\[B_{gh}\]^p_i dV$

となる．

接線剛性

要素接線剛性行列は要素節点等価内力を変位で微分したものである．

$[K_e]^{p,q}_{i,j} =\frac{\partial \{Q_e\}^p_i}{\partial u^q_j}$

以下、この２項について計算する．

初期応力行列

$[{K_{NL}}_e]^{p,q}_{i,j} = \int_V_e S_{gh}\frac{\partial \[B_{gh}\]^p_i}{\partial u^q_j}dV$

を求めてみよう.

$\frac{\partial \[B_{gh}\]^p_i}{\partial u^q_j} = \frac{\partial}{\partial u^q_j}$\b{g}_g_i\;\frac{\partial N^p}{\partial r_h}$ = $\frac{\partial}{\partial u^q_j}\frac{\partial (x_i + u_i)}{\partial r_g}$\frac{\partial N^p}{\partial r_h} = $\frac{\partial}{\partial u^q_j}\frac{\partial u_i}{\partial r_g}$\frac{\partial N^p}{\partial r_h}$

u_j=N^q u^q_j

と補間されているとすると、

$$\frac{\partial}{\partial u^q_j}\frac{\partial u_i}{\partial r_g}$\frac{\partial N^p}{\partial r_h} = $\frac{\partial}{\partial u^q_j}\frac{\partial N^q u^q_i}{\partial r_g}$\frac{\partial N^p}{\partial r_h} = \frac{\partial N^p}{\partial r_h}\frac{\partial N^q}{\partial r_g}\delta_{ij}$

となる．よって、

$[{K_{NL}}_e]^{p,q}_{i,j} = \int_V_e S_{gh}\frac{\partial \[B_{gh}\]^p_i}{\partial u^q_j}dV = \int_V_e S_{gh}\frac{\partial N^q}{\partial r_g}\frac{\partial N^p}{\partial r_h}\delta_{ij}dV$

となる．

初期変位行列

$[{K_L}_e]^{p,q}_{i,j} = \int_V_e\frac{\partial S_{gh}}{\partial u^q_j}\[B_{gh}\]^p_idV$

を求めてみよう．

まず、 $\frac{\partial S^{gh}}{\partial u^q_j}$ を計算する．

ここで、微小な歪みの増分に対する応力の応答が次のように表されているとする．

dS=C:dE

ここで、は４階のテンソルで

$C=C^{ijkl}G_i\otimes G_j\otimes G_k\otimes G_l$

と成分表示できるとする．

$\frac{\partial S}{\partial E_{ef}}=\frac{\partial}{\partial E_{ef}}$C^{ghef}G_g\otimes G_h\otimes G_e\otimes G_f: E_{ef}G^e\otimes G^f$ = C^{ghef}G_g\otimes G_h$

$\frac{\partial S^{gh}}{\partial u^q_j} = \frac{\partial }{\partial u^q_j} $G^g\cdot S\cdot G^h$ = G^g\cdot$\frac{\partial S}{\partial u^q_j}$\cdot G^h = G^g\cdot$\frac{\partial S}{\partial E_{ef}}\frac{\partial E_{ef}}{\partial u^q_j}$\cdot G^h = G^g\cdot$\frac{\partial S}{\partial E_{ef}}$\cdot G^h \;\frac{1}{2}$[B_{ef}]^q_j+[B_{fe}]^q_j$\\\;\;\; = G^g\cdot$C^{ghef}G_g\otimes G_h$\cdot G^h \;\frac{1}{2}$[B_{ef}]^q_j+[B_{fe}]^q_j$ = C^{ghef}\;\frac{1}{2}$[B_{ef}]^q_j+[B_{fe}]^q_j$\\\;\;\; = \frac{1}{2}(C^{ghef}+C^{ghfe})[B_{ef}]^q_j = \bar{C}^{ghef}[B_{ef}]^q_j$

これらを代入すると、

$[{K_L}_e]^{p,q}_{i,j} = \int_V_e\frac{\partial S_{gh}}{\partial u^q_j}\[B_{gh}\]^p_idV=\int_V_e\bar{C}^{ghef}\[B_{ef}\]^q_j\[B_{gh}\]^p_idV$

のようになる．

以上をまとめると、埋め込み座標を用いた要素接線剛性行列は次のとおり

接線剛性行列
$[K_e]^{p,q}_{i,j} = \int_V_e\bar{C}^{ghef}\[B_{gh}\]^p_i\[B_{ef}\]^q_j + S^{gh}\frac{\partial N^q}{\partial r_g}\frac{\partial N^p}{\partial r_h}\delta_{ij}dV$

参考にしたもの

Books

非線形有限要素法のためのテンソル解析の基礎	久田俊明　著
Nonlinear Theory of Elasticity: Applications in Biomechanics	Taber著

Made by Nobuyuki UMETANI 梅谷　信行

n.umetani@gmail.com

Total Lagrange

目次

支配方程式の弱形式化

公称応力を用いた平衡方程式

問題設定

応力境界条件について

弱形式化した平衡方程式

支配方程式の有限要素法離散化

非対称化した仮想歪テンソルの導入

内挿関数の導入

接線剛性行列

Newton-Raphson法

Newton-Raphson法のアルゴリズム

接線剛性行列

初期応力行列

初期変位行列

接線剛性行列

実装に関して

要素剛性行列、要素内残差ベクトルを作るプログラムの流れ

変位勾配テンソルを求めるルーティン

変位勾配テンソルからGreen-Lagrange歪を求めるルーティンの実装例

歪-変位関係行列を作るルーティンの実装例

第二Piola-Kirhihoff応力と歪-変位関係行列から要素節点等価内力ベクトル $\{Q_e\}$ を求めるルーティンの実装例

埋め込み座標を用いた定式化

埋め込み座標系における歪み

仮想歪みエネルギー式

接線剛性

初期応力行列

初期変位行列

参考にしたもの

Links

Books