Pressure Stiffness

圧力剛性行列

Last Update:2009年05月06日

概要

圧力剛性とは，非線形有限要素法において，圧力境界条件の変位による厳密な微分によって生じる接線剛性のことである．この接線剛性行列を考慮することで，Newton-Raphson法による増分計算を少ない反復数で収束させる効果がある．

圧力境界条件

外力による仮想仕事の変分 $\delta\b{W}_{ext}$ とすると，

$\delta\b{W}_{ext} = \int_s \delta\b{u} \cdot \b{f} ds = \int_s \delta\b{u} \cdot (p\b{n}) ds$

要素節点等価外力

要素の中で，要素内の節点の番号づけ $\{p,q\}$ が，グローバルな節点の番号 $\{a,b\}$ に対応しているとする．

また，形状関数によって，位置 $\b{x}$ と仮想変位 $\delta\b{u}$ が離散化されているとする.つまり，

$\b{x}=N^p\b{x}^p,\;\;\; \delta\b{u}=N^q\delta\b{u}^q$

であるとする．このとき，仮想変位によるエネルギーの変分を次のように表すことができる．

$\delta\b{W}_{ext} = \int_s \delta\b{u} \cdot (p\b{n}) ds\\\;\;\; = \sum_e^n \int_{s_e} \delta\b{u} \cdot (p\b{n}) ds = \sum_e^n \int_{s_e} N^p\delta\b{u}^p \cdot (p\b{n}) ds = \delta u^a_i \sum_e^n \int_s p N^p n_i ds\\\;\;\; = \delta u^a_i\sum_e^nF^p_i\\\;\;\; = \delta\b{u}^a_i \b{F}^a_i = \{\delta\b{u}\}\cdot\{\b{F}\}$

ここで， ${F_e}^p_i$ は要素内節点等価外力ベクトルであり，

${F_e}^p_i=\int_{s_e} p N^p n_i ds$

である．この ${F_e}^p_i$ を具体的に各要素について計算してみよう．

アイソパラメトリック要素の場合

アイソパラメトリック要素の場合は次のように単位法線ベクトルが表される．

$\b{n}=\frac{\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2}}{\left|\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2}\right|}$

また微小面積は次のようになる．

$ds=\left|\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2}\right|dr_1dr_2$

これを用いると要素内等価節点外力は次のとおり．

${F_e}^p_i=\int_{s_e}p N^p \frac{ (\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2})_i}{\left|\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2}\right|} \left|\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2}\right|dr_1dr_2\\ \;\;\; = \int_{s_e}p N^p(\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2})_idr_1dr_2$

具体的にアイソパラメトリック要素の形状が四角形の場合と三角形の場合でこれを計算してみよう．

四角形の場合

$r_1=[-1,1],\;\;r_2=[-1,1]$

${F_e}^p_i = \int_{-1}^1\int_{-1}^1 p N^p(\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2})_idr_1dr_2$

三角形の場合

$r_1=[0,1-r_2],\;\;r_2=[0,1]$

${F_e}^p_i = \int_{0}^1\int_{0}^{1-r_2} p N^p(\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2})_idr_1dr_2$

三角形一次要素の場合

三角形一次要素の場合は，法線が面内で一定である．

$\b{a}=(\b{x}^1-\b{x}^0)\times (\b{x}^2-\b{x}^0)$

とおくと，面積 $\Delta$ は

$\Delta = \frac{1}{2}|\b{a}|$

と書ける．また，面の単位法線ベクトル $\b{n}$ は

$\b{n}=\frac{\b{a}}{|\b{a}|}$

となる．これらを用いると

${F_e}^p_i = \int_{s_e} p N^p n_i ds = p n_i\int_{s_e} N^p ds = \frac{1}{3}pn_i\Delta = \frac{1}{6}pa_i\\\;\;\; = \frac{1}{6}p\b{a}\cdot\b{e}_i$

となる．

圧力剛性行列

圧力剛性行列(pressuer stiffness matrix)は圧力による等価節点外力の節点変位による微分で計算できる．

$K^{ab}_{ij}=\frac{\partial F^a_i}{\partial u^b_j} = \frac{\partial}{\partial u^b_j}\{\sum_n^e {F_e}^p_i\} = \sum_e^n \frac{\partial {F_e}^p_i}{\partial u^q_j} = \sum_n^e {K_e}^{pq}_{ij}$

となる． ${K_e}^{pq}_{ij}$ は要素剛性行列で，

${K_e}^{pq}_{ij} = \frac{\partial {F_e}^p_i}{\partial u^q_j}$

である．これを具体的に計算すると以下のようになる．

アイソパラメトリック要素の場合

圧力はUpdate Lagrange的に取り扱う．

$\frac{\partial\b{x}}{\partial r_m}=\frac{\partial N^s\b{x}^s}{\partial r_m}=\frac{\partial N^s}{\partial r_m}\b{x}^s$

$\b{x}=\b{X}+\b{u}$

$\frac{\partial }{\partial u^q_j}$\frac{\partial\b{x}}{\partial r_m}$=\frac{\partial }{\partial u^q_j}$\frac{\partial N^s}{\partial r_m}\b{x}^s$ = \frac{\partial N^s}{\partial r_m}\frac{\partial x_k^s\b{e}_k}{\partial u^q_j} = \frac{\partial N^s}{\partial r_m}\delta_{sq}\delta_{kj}\b{e}_k = \frac{\partial N^q}{\partial r_m}\b{e}_j$

${K_e}^{pq}_{ij} = \frac{\partial }{\partial u^q_j}{F_e}^p_i = \frac{\partial }{\partial u^q_j}\int_{s_e}p N^p(\frac{\partial\b{x}}{\partial r_1}\times\frac{\partial\b{x}}{\partial r_2})_idr_1dr_2\\\;\;\;\; = \int_{s_e}p N^p \{ \frac{\partial }{\partial u^q_j}$\frac{\partial\b{x}}{\partial r_1}$\times\frac{\partial\b{x}}{\partial r_2} + \frac{\partial\b{x}}{\partial r_1}\times\frac{\partial }{\partial u^q_j}$\frac{\partial\b{x}}{\partial r_2}$\}_idr_1dr_2\\\;\;\;\; = \int_{s_e}p N^p \{ $ \frac{\partial N^q \b{e}_j}{\partial r_1}$\times\frac{\partial\b{x}}{\partial r_2} + \frac{\partial\b{x}}{\partial r_1}\times$\frac{\partial N^q\b{e}_j}{\partial r_2} $\}\cdot\b{e}_i dr_1dr_2\\\;\;\;\; = \int_{s_e}p N^p \{ \frac{\partial N^q}{\partial r_1}$\b{e}_j\times\frac{\partial\b{x}}{\partial r_2}$\cdot\b{e}_i + \frac{\partial N^q}{\partial r_2}$\frac{\partial\b{x}}{\partial r_1}\times\b{e}_j $\cdot\b{e}_i\} dr_1dr_2\\\;\;\;\; = \int_{s_e}p N^p $\b{e}_i\times\b{e}_j$\cdot$\frac{\partial N^q}{\partial r_1}\frac{\partial\b{x}}{\partial r_2}-\frac{\partial N^q}{\partial r_2}\frac{\partial\b{x}}{\partial r_1}$dr_1dr_2$

座標の添え字に対して反対称(Skew Symmetric)

${K_e}^{pq}_{ij}=-{K_e}^{pq}_{ji}$

であることがわかる．

三角形一次要素の場合

$\frac{\partial\b{x}^p}{\partial u^q_j} = \delta_{pq}\b{e}_j$

$\b{a}=(\b{x}^1-\b{x}^0)\times (\b{x}^2-\b{x}^0)$

$\{\begin{array}{l}\frac{\partial \b{a}}{\partial \b{u}^0_j} = \frac{\partial}{\partial \b{u}^0_j}\{(\b{x}^1-\b{x}^0)\times (\b{x}^2-\b{x}^0)\} = (-\b{e}_j)\times (\b{x}^2-\b{x}^0) + (\b{x}^1-\b{x}^0)\times(-\b{e}_j) = \b{e}_j\times(\b{x}^1-\b{x}^2)\\ \frac{\partial \b{a}}{\partial \b{u}^1_j} = \frac{\partial}{\partial \b{u}^1_j}\{(\b{x}^1-\b{x}^0)\times (\b{x}^2-\b{x}^0)\} = \b{e}_j\times(\b{x}^2-\b{x}^0)\\ \frac{\partial \b{a}}{\partial \b{u}^2_j} = \frac{\partial}{\partial \b{u}^2_j}\{(\b{x}^1-\b{x}^0)\times (\b{x}^2-\b{x}^0)\} = (\b{x}^1-\b{x}^0)\times \b{e}_j = \b{e}_j\times (\b{x}^0-\b{x}^1)\end{array}$

${K_e}^{pq}_{ij} = \frac{\partial {F_e}^p_i}{\partial u^q_j} = \frac{\partial}{\partial u^q_j}$\frac{1}{6}p\b{a}\cdot\b{e}_i$ = \frac{1}{6}p\frac{\partial\b{a}}{\partial u^q_j}\cdot\b{e}_i$

これを用いると，

$\{\begin{array}{l} {K_e}^{p0}_{ij} = \frac{1}{6}p(\b{e}_i\times\b{e}_j)\cdot(\b{x}^1-\b{x}^2)\\ {K_e}^{p1}_{ij} = \frac{1}{6}p(\b{e}_i\times\b{e}_j)\cdot(\b{x}^2-\b{x}^0)\\ {K_e}^{p2}_{ij} = \frac{1}{6}p(\b{e}_i\times\b{e}_j)\cdot(\b{x}^0-\b{x}^1)\end{array}$