|TOP Page|


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

要素節点等価外力

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

また,形状関数Nによって,位置\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}

となる.こちらも座標の添え字に対して反対称(Skew Symmetric)

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

であることがわかる.

参考

Links

東京大学新領域創成科学研究科、久田・杉浦研究室、渡邊先生の講義資料
http://www.sml.k.u-tokyo.ac.jp/members/nabe/

Books

非線形有限要素法のためのテンソル解析の基礎 久田俊明 著


Made by Nobuyuki UMETANI  梅谷 信行
n.umetani@gmail.com