DKT Element

DKT要素、Discrete Kirhihoff Triangle要素

Last Update:2009年05月15日

概要

DKT要素(Discrete Kirchhoff Triangle)は薄板３角形要素の中で広く使われる要素である。薄板要素とはKirchhoffの仮定が成り立つような要素である。

薄板の解析において一般的な変位法を用いた有限要素法ではKirchhoffの仮定は一種の制約条件として働く。このときような問題を変位法で解析する場合、補間関数が要素間で C_1 連続を満たさなければならない。これにはHermaite型のような高次の補間関数を必要とするので手続きが非常に煩雑であるという欠点がある。DKT要素ではKirchhoffの仮定が要素の中の離散的な点においてのみ成り立つとして計算を行う。このため、撓みや回転角の C_1 連続性は考慮しなくてもよく、一般的な補間関数を用いることができるので容易であるという利点がある。

DKT要素を用いた薄板の解析は１９６０年代にDhatt[1]らによって提案された。それから１０年の後、Batoz[2]によって種々の薄板要素が数値実験によって比較され、DKT要素の有効性が示された。

離散Kirchhoff条件を四角形に応用したDKQ(Disctete Kirchhoff Quadrilateral)やさらに多角形に応用した要素も提案されている。

薄板理論

薄板には次の３つの仮定が成り立つ。

中立面は変形で伸縮されない。
変形前に中立面に垂直だった物質中の直線は、変形後も中立面に対して垂直で直線である。
面外剪断による歪エネルギーは曲げエネルギーに対して十分小さいので無視できる。

これらの仮定はKirchhoffの仮定(Kirchhoff's hypothesis)と呼ばれ、これを元に微小変形薄板要素が作られる。

薄板がｘｙ平面上にある場合に板がｚ方向へ撓んでいるとしよう。

ｘｙｚの変位は次のようになる。

$u=z\theta_y$

$v=-z\theta_x$

z=w

ここでは撓み(Deflection)であり、 $\theta_x$ 、 $\theta_y$ はそれぞれ中立面の垂線のx軸まわり、y軸周りの回転量である。

$\theta_x= \frac{\partial w}{\partial y}$

$\theta_y=-\frac{\partial w}{\partial x}$

歪と応力

この時工学歪は次のようになる。

$\epsilon_x=\frac{\partial u}{\partial x}=-z\frac{\partial^2 w}{\partial x^2}$

$\epsilon_y=\frac{\partial v}{\partial y}=-z\frac{\partial^2 w}{\partial y^2}$

$\gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}=-2z\frac{\partial^2 w}{\partial x\partial y}$

平面応力条件では応力-歪関係は次のようになる。

$\{\begin{array}\sigma_x\\ \sigma_y\\ \tau_{xy}\end{array}\}=\frac{E}{1-\nu^2} \[\begin{array} 1&&\nu&&\\ \nu&&1&&\\ && && && \frac{1-\nu}{2}\\ \end{array}\] \{\begin{array}\epsilon_x\\ \epsilon_y\\ \gamma_{xy}\end{array}\}$

歪を撓みによって置き換えると、

$\{\begin{array}\sigma_x\\ \sigma_y\\ \tau_{xy}\end{array}\}=\frac{E}{1-\nu^2} \[\begin{array} 1&&\nu&&\\ \nu&&1&&\\ && && && \frac{1-\nu}{2}\\ \end{array}\] \{\begin{array}-z\frac{\partial^2 w}{\partial x^2}\\-z\frac{\partial^2 w}{\partial x^2}\\ -2z\frac{\partial^2 w}{\partial x\partial y}\end{array}\}$

撓みと曲率

板の中立面の曲率をたわみや撓み角 $\theta$ と対応づけてみよう。

板の中立面上の点 $\b{p}$ は次のように表されるのだった

$\b{p}(x,y)=$\begin{array}x\\y\\w(x,y)\end{array}$$

これを x,y で微分することにより接ベクトルを $\b{r}_x,\b{r}_y$ を求めると

$\b{r}_x=\frac{\partial \b{p}}{\partial x}=$\begin{array}1\\0\\\frac{\partial w}{\partial x}\end{array}$$

$\b{r}_y=\frac{\partial \b{p}}{\partial y}=$\begin{array}0\\1\\\frac{\partial w}{\partial y}\end{array}$$

さて接ベクトル $\b{r}_x,\b{r}_y$ に直交するベクトル $\b{r}_3$ は外積を用いて次のように書ける

$\b{r}_3=\frac{\b{r}_x\times\b{r}_y}{||\b{r}_x\times\b{r}_y||}=\frac{1}{\sqrt{1+(\frac{\partial w}{\partial y})^2+(\frac{\partial w}{\partial y})^2}}$\begin{array}-\frac{\partial w}{\partial x}\\-\frac{\partial w}{\partial y}\\1\end{array}$\simeq$\begin{array}-\frac{\partial w}{\partial x}\\-\frac{\partial w}{\partial y}\\1\end{array}$$

さて曲面の曲率テンソル $\kappa$ は次のように定義される

$\kappa_{ij}=\frac{\partial \b{r}_i}{\partial \b{x}_j}\cdot \b{r}_3$

これから

$\kappa_{xx}=$\begin{array}0\\0\\\frac{\partial^2 w}{\partial x^2}\end{array}$\cdot$\begin{array}-\frac{\partial w}{\partial x}\\-\frac{\partial w}{\partial y}\\1\end{array}$=\frac{\partial^2 w}{\partial x^2}$

$\kappa_{yy}=$\begin{array}0\\0\\\frac{\partial^2 w}{\partial y^2}\end{array}$\cdot$\begin{array}-\frac{\partial w}{\partial x}\\-\frac{\partial w}{\partial y}\\1\end{array}$=\frac{\partial^2 w}{\partial y^2}$

$\kappa_{xy}=\kappa_{yx}=$\begin{array}0\\0\\\frac{\partial^2 w}{\partial x\partial y}\end{array}$\cdot$\begin{array}-\frac{\partial w}{\partial x}\\-\frac{\partial w}{\partial y}\\1\end{array}$=\frac{\partial^2 w}{\partial x\partial y}$

ここで曲率テンソル $\kappa$ を次のようにベクトル表記する。

$\kappa=\{\begin{array}\kappa_x\\ \kappa_y\\ 2\kappa_{xy}\end{array}\} =\{\begin{array}\frac{\partial^2 w}{\partial x^2}\\\frac{\partial^2 w}{\partial x^2}\\ 2\frac{\partial^2 w}{\partial x\partial y}\end{array}\} =\{\begin{array}\frac{\partial \theta_y}{\partial x}\\ -\frac{\partial \theta_x}{\partial y}\\ \frac{\partial \theta_y}{\partial y}-\frac{\partial \theta_x}{\partial x}\end{array}\}$

歪エネルギー

さてここで曲げによる歪エネルギー U_b を求めてみよう。

$U=\frac{1}{2}\int_{\Omega}\int^{\frac{t}{2}}_{-\frac{t}{2}}\{\begin{array}\epsilon_x\\ \epsilon_y\\ \gamma_{xy}\end{array}\}^T\{\begin{array}\sigma_x\\ \sigma_y\\ \tau_{xy}\end{array}\}dzd\Omega\\\qquad\qquad= \frac{1}{2}\int_{\Omega}\int^{\frac{t}{2}}_{-\frac{t}{2}}\{\begin{array}-z\frac{\partial^2 w}{\partial x^2}\\-z\frac{\partial^2 w}{\partial x^2}\\ -2z\frac{\partial^2 w}{\partial x\partial y}\end{array}\}^T\frac{E}{1-\nu^2} \[\begin{array} 1&&\nu&&\\ \nu&&1&&\\ && && && \frac{1-\nu}{2}\\ \end{array}\] \{\begin{array}-z\frac{\partial^2 w}{\partial x^2}\\-z\frac{\partial^2 w}{\partial x^2}\\ -2z\frac{\partial^2 w}{\partial x\partial y}\end{array}\}dzd\Omega$

ここで $D=\int^{t/2}_{-t/2}\frac{z^2E}{1-\nu^2}dz=\frac{E}{1-\nu^2}\[\frac{1}{3}z^3\]^{t/2}_{-t/2}=\frac{Et^3}{12(1-\nu^2)}$ とおき

$D_b=D\[\begin{array} 1&&\nu&&\\ \nu&&1&&\\ && && && \frac{1-\nu}{2}\\ \end{array}\]$ とおくと次のように書ける。

$U=\frac{1}{2}\int_{\Omega}\{\begin{array}\frac{\partial^2 w}{\partial x^2}\\\frac{\partial^2 w}{\partial x^2}\\ 2\frac{\partial^2 w}{\partial x\partial y}\end{array}\}^TD_b\{\begin{array}\frac{\partial^2 w}{\partial x^2}\\\frac{\partial^2 w}{\partial x^2}\\ 2\frac{\partial^2 w}{\partial x\partial y}\end{array}\}d\Omega$

さらに変形後の中立面の曲率 $\kappa$ ベクトルを使うと歪エネルギーは次のように書ける

薄板の歪エネルギー
$U=\frac{1}{2}\int_{\Omega}\b{\kappa}D_b\b{\kappa}d\Omega$

これはKirchhoff条件を仮定しない一般的な問題の曲げ歪エネルギー U_b にあたる。つまり薄板の曲げによる歪エネルギーは、剪断歪エネルギーが０で曲げ歪エネルギーだけである。

曲げモーメント、捩じりモーメント

行列表記をやめると次のように書ける。

$\sigma_x=\frac{-zE}{1-\nu^2}\[\frac{\partial^2 w}{\partial x^2}+\nu\frac{\partial^2 w}{\partial y^2}\]$

$\sigma_y=\frac{-zE}{1-\nu^2}\[\nu\frac{\partial^2 w}{\partial x^2}+\frac{\partial^2 w}{\partial y^2}\]$

$\tau_{xy}=\frac{-2zE}{1-\nu^2}$\frac{1-\nu}{2}$$\frac{\partial^2 w}{\partial x\partial y}$$

曲げモーメント、捩りモーメントは次のとおり

$M_x=\int^{\frac{t}{2}}_{-\frac{t}{2}}\sigma_x\cdot zdz=\frac{-Et^3}{12(1-\nu^2)}\[\frac{\partial^2 w}{\partial x^2}+\nu\frac{\partial^2 w}{\partial y^2}\]=D\[-\frac{\partial^2 w}{\partial x^2}-\nu\frac{\partial^2 w}{\partial y^2}\]$

$M_y=\int^{\frac{t}{2}}_{-\frac{t}{2}}\sigma_y\cdot zdz=\frac{-Et^3}{12(1-\nu^2)}\[\nu\frac{\partial^2 w}{\partial x^2}+\frac{\partial^2 w}{\partial y^2}\]=D\[-\nu\frac{\partial^2 w}{\partial x^2}-\frac{\partial^2 w}{\partial y^2}\]$

$M_{xy}=\int^{\frac{t}{2}}_{-\frac{t}{2}}\tau_{xy}\cdot zdz=\frac{-2Et^3}{12(1-\nu^2)}$\frac{1-\nu}{2}$$\frac{\partial^2 w}{\partial x\partial y}$=D$\frac{1-\nu}{2}$$-\frac{\partial^2 w}{\partial x\partial y}$$

１次元梁要素の変位法による解析

１次元梁の問題の問題を通して、梁問題を変位法で解析する場合、撓みとその勾配が要素間で連続でなければならないことを示しておく。１次元梁問題ではこのような補間は比較的容易であるが、シェルの問題になると複雑な補間が必要であり、実用的ではない。

後に示すようにDKTは辺上で撓みを３次補間しているのと等価であり、ある種のアナロジーが成り立つ。

両端固定梁

次のような両端が変位と回転が固定されているような梁の問題を考えよう

w(0)=w(l)=w'(0)=w'(l)=0

を曲げ剛性、つまり k=EI （但しはヤング率は断面２次モーメント）とするとのように書ける。

$M=k\frac{\partial^2 w}{\partial x^2}$

$\frac{\partial^2 M}{\partial x^2}=f$

よって、強形式の支配方程式は

$\{\begin{array}k\frac{\partial^4 w}{\partial x^4}=f\end{array}\\w(0)=w(l)=w'(0)=w'(l)=0$

となる。強形式ではは４回微分可能でなければならない。

弱形式化

この問題は次のように弱形式化される

$find\; v\in W\; st\; \int_0^l k w'' v'' dx=\int_0^l f v dx, \;\;\forall v\in W$

ここでは２階微分が可能で、２階微分がリーマン可積分、つまり２階微分が区分連続であること。それに加えて、境界条件を満たす集合である、つまり

$W=\{\;v\;|\;v\in C[0,l],\; v'\in C[0,l],\; v''\in PC[0,l],\; v(0)=v(l)=v'(0)=v'(l)=0 \;\}$

強形式の式が成り立てば弱形式の式が成り立つことを示す。

強形式の式に任意の $v\in W$ を乗じて、区間 [0,l] で積分すると、

$\int_0^l k \frac{\partial^4 w}{\partial x^4} vdx=\int_0^l f v dx \\\;\;\leftrightarrow \int_0^l \frac{\partial }{\partial x}(k \frac{\partial^3 w}{\partial x^3} v)dx-\int_0^l \frac{\partial}{\partial x}(k\frac{\partial^2 w}{\partial x^2}\frac{\partial v}{\partial x}) dx+\int_0^l k\frac{\partial^2 w}{\partial x^2}\frac{\partial^2 v}{\partial x^2}\ dx=\int_0^l f v dx$

左辺第１，２項は v(0)=v(l)=v'(0)=v'(l)=0 から０になる。よって

$\int_0^l k\frac{\partial^2 w}{\partial x^2}\frac{\partial^2 v}{\partial x^2}\ dx=\int_0^l f v dx$

が成り立つ。

補間関数

区間 I=[0,l] を分割して I_i=[a,b] となっていることにしよう。つまり、要素の中で、は C^2 級関数であるとすると、の２階微分は区間連続である。さて、が連続であるためには、節点を区間要素の両端に置けば、隣合う要素同士で値を共有するので良かった。同様にv'が連続であるためには勾配を表す節点を区間要素の両端に置けばよい。

要素 I_i の中で補間関数の２階微分が定数であるとしよう。つまり、要素間で補間関数は３次関数となる。要素 I_i の４つの節点(両端における変位と変位勾配)に対応する補間関数 N^1,N^2,N^3,N^4 は次のようになる。

$N_i^1(a)=1,\;N_i^1(b)=0,\;N_i^1'(a)=0,\;N_i^1'(b)=0$

$N_i^2(a)=0,\;N_i^2(b)=0,\;N_i^2'(a)=1,\;N_i^2'(b)=0$

$N_i^3(a)=0,\;N_i^3(b)=1,\;N_i^3'(a)=0,\;N_i^3'(b)=0$

$N_i^4(a)=0,\;N_i^4(b)=0,\;N_i^4'(a)=0,\;N_i^4'(b)=1$

DKT板曲げ要素

薄板の歪エネルギーは面の曲率のみを用いて書くことができることがわかった。面の曲率は回転の微分から求められる。面内の回転の分布がわかればエネルギーが求まるということがわかる。

x軸周りの回転、y軸周りの回転が次のように３角形の頂点と辺の中点における回転によって２次補間されるとする

$\b{\theta}=\sum^6_{i=1} N^i \b{\theta}^b$

さて頂点における回転量が変数として与えられているので、辺の中点におけるの回転量がわかれば回転量の３角形内の分布を求めることができる。分布が求まれば微分をとって曲率がもとまる。その曲率をもとにエネルギー、剛性行列がもとまる。

DKT要素作成の流れ
頂点における撓み、ｘｙ軸周りの回転 ↓ 頂点および辺の中点におけるｘｙ軸周りの回転 ↓ ｘｙ軸周りの回転の補間 ↓ 回転の微分をとって曲率を求める ↓ エネルギーを計算し、剛性行列を求める。

以下に示すように撓みは要素の頂点および辺の中点で C_1 級連続である。しかし、中点以外の辺上では連続性は考慮されない。よって厳密にいえばDKT要素は非適合要素の部類に入る。DKT法はエネルギーを得るために直接変位を使う変位法と違い、回転角を得るためだけに撓みを使い、回転角から導かれる曲率をもとにエネルギーを計算しているので、このような変位の不連続性があってもあまり気にする必要はない。

辺の中点における回転量

下図のように、三角形abcがあった場合に辺bcの中点dがあるとする。

辺bcの中点dにおける回転量 $\b{\theta}^d$ を辺の頂点の撓み w^b,w^c と回転 $\b{\theta}^b,\b{\theta}^c$ によって求める方法を考えよう。これには以下の３つの条件を使う。

$\b{\theta}^d\cdot \b{s}^d=\frac{1}{2}(\b{\theta}^b\cdot \b{s}^d+\b{\theta}^c\cdot \b{s}^d)$
辺の中点における辺に添った方向の撓みの微分は辺上で変位が３次関数と仮定したときの値に等しい。
３角形要素の頂点及び辺の中点でKirhihoff条件を満たす(Disctete Kirhihoff条件)

１番目の条件

1番目の条件は辺に添ったベクトル周りの回転において中点の値がが頂点の値から線形補間されることを意味している。

中点において、辺に添ったベクトル周りの回転についてKirchhoff条件から回転量を定めようとすると撓みの辺と垂直な方向への微分が必要になる。これは要素内部の撓みの情報が必要であることを意味している。すると辺を介して向かい合う要素同士で中点における、辺と垂直な方向の微分の連続性を満たすことができなくなる。微分の連続性を保つためには、このように辺を介して向かい合う要素が両方とも共有している辺両端の情報から回転を決めることが必要不可欠となってくる。これは次に求める、辺に添ったベクトル周りの回転だけでなく、辺に垂直なベクトル周りの回転についてもいえる。

さて、１番目の条件が仮定された時、中点における回転を求める際に,辺の法線周りの回転の中点での値 $\b{\theta}^d\cdot\b{n}^d$ がわかれば、座標変換を使って $\b{\theta}=(\b{\theta}^d\cdot\b{s}^d)\b{s}^d+(\b{\theta}^d\cdot\b{n}^d)\b{n}^d$ というように中点での回転を求めることができる。 $\b{\theta}^d\cdot\b{n}^d$ を頂点の回転と撓みから求めるために２番目と３番目の条件を使う。

２番目の条件

2番目の条件について考えてみよう。回転が面内で２次補間されていることから、辺上の回転は２次関数である。座標変換すると辺の法線周りの回転も２次関数である。辺に添った方向での撓みの微分に関してKirchhoff条件が成り立つ場合、ｗの微分が辺に垂直なベクトル周りの回転になるということであるから、撓みは辺上で３次関数になることを意味している。辺上の撓みをとりあえず求めてみよう。辺上の撓みの分布がわかれば３番目の条件から $\b{\theta}^d\cdot{s}^d$ を求めることができる。

辺上の撓みの分布が3次関数であることが仮定されている。３次関数を定めるためには何か４つ条件があればよい（例えば４つの点における値など）。今３角形の頂点における値から辺の中点における値を求めたかったから、４つの条件を辺の両端における値とその微分としよう。

辺の長さを l_d とするとき、頂点で０頂点で l_d をとるような媒介変数を用いての辺における値 $\bar{w}^d(s)$ を表す。これらの条件は次のように書ける。

$\bar{w}^d(0)=w^b$ $\bar{w}^d(l_d)=w^c$

$\frac{\partial \bar{w}^d}{\partial s}(0)=\left.\frac{\partial w}{\partial \b{s}^d}\right|_b=(\nabla w)^b\cdot\b{s}^d$ $\frac{\partial \bar{w}^d}{\partial s}(l_d)=\left.\frac{\partial w}{\partial\b{s}^d\right|_c=(\nabla w)^c\cdot\b{s}^d$

と l_d-s に関する対称性に注意すると、簡単な計算から w^d(s) は次のように書ける

$\bar{w}^d(s)=\frac{(l_d-s)^2\{3s+(l_d-s)\}}{l_d^3}w^b+\frac{s^2\{3(l_d-s)+s\}}{l_d^3}w^c+\frac{s(l_d-s)^2}{l_d^2}\left.\frac{\partial w}{\partial \b{s}^d}\right|_b+\frac{s^2(l_d-s)}{l_d^2}\left.\frac{\partial w}{\partial \b{s}^d}\right|_c$

上の微分を計算すると

$\frac{\partial \bar{w}^d}{\partial s}(s)=\frac{-6s(l_d-s)}{l_d^3}w^b+\frac{6(l_d-s)s}{l_d^3}w^c-\frac{(l_d-s)\{2s-(l_d-s)\}}{l_d^2}\left.\frac{\partial w}{\partial \b{s}^d}\right|_b-\frac{s\{2(l_d-s)-s\}}{l_d^2}\left.\frac{\partial w}{\partial \b{s}^d}\right|_c$

上式に $s=\frac{l_d}{2}$ を代入して辺の中点での辺に添った方向の微分を求める。

$\left.\frac{\partial w}{\partial \b{s}^d}\right|_d=\frac{\partial \bar{w}^d}{\partial s}$\frac{l_d}{2}$=-\frac{3}{2l_d}w^b+\frac{3}{2l_d}w^c-\frac{1}{4}\left.\frac{\partial w}{\partial \b{s}^d}\right|_b-\frac{1}{4}\left.\frac{\partial w}{\partial \b{s}^d}\right|_c$

書き換えると

$(\nabla w)^d\cdot \b{s}^d=-\frac{3}{2l_d}w^b+\frac{3}{2l_d}w^c-\frac{1}{4}(\nabla w)^b\cdot\b{s}^d-\frac{1}{4}(\nabla w)^c\cdot\b{s}^d$

さて、辺の中点上での辺に添った方向 $\b{s}^d$ のの微分が上のように表されることがわかった。これと３番目の条件を用いて $\theta^d_s$ を求めてみよう。

３番目の条件

２番目の条件はKirhihoff条件を離散的に満たすからDiscrete Kirhihoff条件（離散Kirhihoff条件）と呼ばれる。

以下の式が頂点と辺の中点で成り立つ。

$\theta_x=\frac{\partial w}{\partial y}\qquad\qquad\theta_y=-\frac{\partial w}{\partial x}$

これを用いて撓みの勾配を計算すると

$\nabla w=\frac{\partial w}{\partial x}\b{e}_x+\frac{\partial w}{\partial y}\b{e}_y=-\theta_y\b{e}_x+\theta_x\b{e}_y$

$(\nabla w)\cdot \b{s}=-\theta_y(\b{e}_x\cdot\b{s})+\theta_x(\b{e}_y\cdot\b{s})=\theta_y (-s_x)+\theta_x s_y=\theta_y n_y+\theta_x n_x=\b{\theta}\cdot\b{n}$

これを上の式に代入すると

$\b{\theta}^d\cdot \b{n}^d=-\frac{3}{2l_d}w^b+\frac{3}{2l_d}w^c-\frac{1}{4}\b{\theta}^b\cdot\b{n}^d-\frac{1}{4}\b{\theta}^c\cdot\b{n}^d$

となり、法線方向の回転が得られる。

中点での回転の値

$\b{\theta}^d\cdot \b{n}^d=-\frac{3}{2l_d}w^b+\frac{3}{2l_d}w^c-\frac{1}{4}\b{\theta}^b\cdot\b{n}^d-\frac{1}{4}\b{\theta}^c\cdot\b{n}^d$

$\b{\theta}^d\cdot \b{s}^d=\frac{1}{2}(\b{\theta}^b\cdot \b{s}^d+\b{\theta}^c\cdot \b{s}^d)$

であった。ここから $\b{\theta}^d$ を求めよう。 $\b{s}^d$ と $\b{n}^d$ は直交しているから

$\b{\theta}^d=(\b{\theta}^d\cdot \b{n}^d)\b{n}^d+(\b{\theta}^d\cdot \b{s}^d)\b{s}^d\\\qquad\qquad=(-\frac{3}{2l_d}w^b+\frac{3}{2l_d}w^c-\frac{1}{4}\b{\theta}^b\cdot\b{n}^d-\frac{1}{4}\b{\theta}^c\cdot\b{n}^d)\b{n}^d+(\frac{1}{2}\b{\theta}^b\cdot \b{s}^d+\frac{1}{2}\b{\theta}^c\cdot \b{s}^d)\b{s}^d$

これで中点における回転量が求まったわけであるが、成分を書き出すことで、もう少し簡単な形に整理してみよう。

$\theta^d_x=(-\frac{3n^d_x}{2l_d})w^b+(\frac{3n^d_x}{2l_d})w^c+(-\frac{{n^d_x}^2}{4}+\frac{{s^d_x}^2}{2})(\theta^b_x+\theta^c_x)+(-\frac{n^d_yn^d_x}{4}+\frac{s^d_ys^d_x}{2})(\theta^b_y+\theta^c_y)$

$\theta^d_y=(-\frac{3n^d_y}{2l_d})w^b+(\frac{3n^d_y}{2l_d})w^c+(-\frac{n^d_xn^d_y}{4}+\frac{s^d_xs^d_y}{2})(\theta^b_x+\theta^c_x)+(-\frac{{n^d_y}^2}{4}+\frac{{s^d_y}^2}{2})(\theta^b_y+\theta^c_y)$

更に辺の法線ベクトル接線ベクトルを頂点座標と対応づけよう。

$\b{s}^d=(\b{p}^c-\b{p}^b)/l_d=$\begin{array}x_{cb}\\y_{cb}\end{array}$/l_d\quad\rightarrow\quad s^d_x=x_{cb}/l_d,\quad s^d_y=y_{cb}/l_d$

$\b{s}^d\bot\b{n}^d\qquad\rightarrow\qquad\b{n}^d=$\begin{array}y_{cb}\\-x_{cb}\end{array}$/l_d\quad\rightarrow\quad n^d_x=y_{cb}/l_d,\quad n^d_y=-x_{cb}/l_d$

を代入すると

$\theta^d_x=(-\frac{3y_{cb}}{2l_d^2})w^b+(\frac{3y_{cb}}{2l_d^2})w^c+(-\frac{{y_{cb}}^2}{4l^2_d}+\frac{{x_{cb}}^2}{2l^2_d})(\theta^b_x+\theta^c_x)+(\frac{3x_{cb}y_{cb}}{4l^2_d})(\theta^b_y+\theta^c_y)$

$\theta^d_y=(\frac{3x_{cb}}{2l_d^2})w^b+(-\frac{3x_{cb}}{2l_d^2})w^c+(\frac{3x_{cb}y_{cb}}{4l^2_d})(\theta^b_x+\theta^c_x)+(-\frac{{x_{cb}}^2}{4l^2_d}+\frac{{y_{cb}}^2}{2l^2_d})(\theta^b_y+\theta^c_y)$

さらに次のように変数を定める

$a_d=-x_{cb}/l_d^2$

$b_d=\frac{3}{4}x_{cb}y_{cb}/l^2_d$

$c_d=(\frac{{x_{cb}}^2}{4}-\frac{{y_{cb}}^2}{2})/l^2_d=\frac{1}{4}-\frac{3}{4}\frac{y^2_{cb}}{l^2_d}$

$d_d=-y_{cb}/l_d^2$

$e_d=(\frac{{y_{cb}}^2}{4}-\frac{{x_{cb}}^2}{2})/l^2_d=\frac{3}{4}\frac{y^2_{cb}}{l_d^2}-\frac{1}{2}$

これらを用いると最終的に辺の中点における回転量は次のとおり。

$\theta^d_x=+1.5d_dw^b-1.5d_dw^c-e_d(\theta^b_x+\theta^c_x)+b_d(\theta^b_y+\theta^c_y)$

$\theta^d_y=-1.5a_dw^b+1.5a_dw^c+b_d(\theta^b_x+\theta^c_x)-c_d(\theta^b_y+\theta^c_y)$

節点の撓みとｘｙ軸回転から曲率を求める行列

撓み角の２次補間

$\theta_x=N^1\theta^1_x + N^2\theta^2_x + N^3\theta^3_x + N^4\theta^4_x + N^5\theta^5_x + N^6\theta^6_x\\\qquad\qquad=N^1\theta^1_x + N^2\theta^2_x + N^3\theta^3_x\\\qquad\qquad\qquad+N^4\{+1.5d_4w^2-1.5d_4w^3-e_4(\theta^2_x+\theta^3_x)+b_4(\theta^2_y+\theta^3_y)\}\\\qquad\qquad\qquad+N^5\{+1.5d_5w^3-1.5d_5w^1-e_5(\theta^3_x+\theta^1_x)+b_5(\theta^3_y+\theta^1_y)\}\\\qquad\qquad\qquad+N^6\{+1.5d_6w^1-1.5d_6w^2-e_6(\theta^1_x+\theta^2_x)+b_6(\theta^1_y+\theta^2_y)\}\\\qquad\qquad=-1.5(N^5d_5-N^6d_6)w^1+(N^1-N^5e_5-N^6e_6)\theta^1_x+(N^5b_5+N^6b_6)\theta^1_y\\\qquad\qquad\qquad-1.5(N^6d_6-N^4d_4)w^2+(N^2-N^6e_6-N^4e_4)\theta^2_x+(N^6b_6+N^4b_4)\theta^2_y\\\qquad\qquad\qquad-1.5(N^4d_4-N^5d_5)w^3+(N^3-N^4e_4-N^5e_5)\theta^3_x+(N^4b_4+N^5b_5)\theta^3_y$

$\theta_y$ についても同様に次のように補間関数を用いて表される。

$\theta_y=N^1\theta^1_y + N^2\theta^2_y + N^3\theta^3_y + N^4\theta^4_y + N^5\theta^5_y + N^6\theta^6_y\\\qquad\qquad=N^1\theta^1_y + N^2\theta^2_y + N^3\theta^3_y\\\qquad\qquad\qquad+N^4\{-1.5a_4w^2+1.5a_4w^3+b_4(\theta^2_x+\theta^3_x)-c_4(\theta^2_y+\theta^3_y)\}\\\qquad\qquad\qquad+N^5\{-1.5a_5w^3+1.5a_5w^1+b_5(\theta^3_x+\theta^1_x)-c_5(\theta^3_y+\theta^1_y)\}\\\qquad\qquad\qquad+N^6\{-1.5a_6w^1+1.5a_6w^2+b_6(\theta^1_x+\theta^2_x)-c_6(\theta^1_y+\theta^2_y)\}\\\qquad\qquad=+1.5(N^5a_5-N^6a_6)w^1+(N^5b_5+N^6b_6)\theta^1_x+(N^1-N^5c_5-N^6c_6)\theta^1_y\\\qquad\qquad\qquad+1.5(N^6a_6-N^4a_4)w^2+(N^6b_6+N^4b_4)\theta^2_x+(N^2-N^6c_6-N^4c_4)\theta^2_y\\\qquad\qquad\qquad+1.5(N^4a_4-N^5a_5)w^3+(N^4b_4+N^5b_5)\theta^3_x+(N^3-N^4c_4-N^5c_5)\theta^3_y$

行列、ベクトル表記

次のように未知数を並べたベクトル $\b{r}$ を作る

$\b{r}^T=\{w^1,w^2,w^3,\theta^1_x,\theta^1_y,\theta^2_x,\theta^2_y,\theta^3_x,\theta^3_y\}$

$\theta_x,\theta_y$ は行列 $\b{H}_x,\b{H}_y$ を用いて次のように書ける

$\theta_x=\b{H}_x\b{r}\qquad\qquad \theta_y=\b{H}_y\b{r}$

但し行列 $\b{H}_x,\b{H}_y$ は次のとおり

$\b{H}^T_x=\{\begin{array}-1.5(N^5d_5-N^6d_6)\\ -1.5(N^6d_6-N^4d_4)\\-1.5(N^4d_4-N^5d_5)\\ (N^1-N^5e_5-N^6e_6) \\(N^5b_5+N^6b_6)\\ (N^2-N^6e_6-N^4e_4)\\ (N^6b_6+N^4b_4)\\ (N^3-N^4e_4-N^5e_5)\\ (N^4b_4+N^5b_5)\end{array}\} \qquad\qquad \b{H}^T_y=\{\begin{array}+1.5(N^5a_5-N^6a_6)\\ +1.5(N^6a_6-N^4a_4)\\ +1.5(N^4a_4-N^5a_5)\\ (N^5b_5+N^6b_6)\\ (N^1-N^5c_5-N^6c_6)\\ (N^6b_6+N^4b_4)\\ (N^2-N^6c_6-N^4c_4)\\ (N^4b_4+N^5b_5)\\ (N^3-N^4c_4-N^5c_5) \end{array}\}$

線形形状関数の導入

補間関数 $N^i\quad(i=1,\ldots,6)$ はサブパラメトリック２次補間であるから１次補間関数 $L_i\quad(i=1,\ldots,3)$ を使って次のように表すことができる。

$N^1=L_1(2L_1-1) \qquad\qquad N^2=L_2(2L_2-1) \qquad\qquad N^3=L_3(2L_3-1)$

$N^4=4L_2L_3 \qquad\qquad N^5=4L_3L_1 \qquad\qquad N^6=4L_1L_2 \qquad\qquad$

これを代入すると、 $\b{H}_x,\b{H}_y$ はそれぞれ次のように書くことができる。

$\b{H}^T_x=\{\begin{array}-6(L_1L_3d_5-L_1L_2d_6)\\ -6(L_1L_2d_6-L_2L_3d_4)\\ -6(L_2L_3d_4-L_1L_3d_5)\\ L_1(2L_1-1)-4L_1L_3e_5-4L_1L_2e_6\\ 4(L_1L_3b_5+L_1L_2b_6)\\ L_2(2L_2-1)-4L_1L_2e_6-4L_2L_3e_4\\ 4(L_1L_2b_6+L_2L_3b_4)\\ L_3(2L_3-1)-4L_2L_3e_4-4L_1L_3e_5\\ 4(L_2L_3b_4+L_1L_3b_5)\end{array}\}\qquad\qquad\b{H}^T_y=\{\begin{array}+6(L_1L_3a_5-L_1L_2a_6)\\ +6(L_1L_2a_6-L_2L_3a_4)\\ +6(L_2L_3a_4-L_1L_3a_5)\\ 4(L_1L_3b_5+L_1L_2b_6)\\ L_1(2L_1-1)-4L_1L_3c_5-4L_1L_2c_6\\ 4(L_1L_2b_6+L_2L_3b_4)\\ L_2(2L_2-1)-4L_1L_2c_6-4L_2L_3c_4\\ 4(L_2L_3b_4+L_1L_3b_5)\\ L_3(2L_3-1)-4L_2L_3c_4-4L_1L_3c_5 \end{array}\}$

撓み角の微分

さて $\b{H}_x,\b{H}_y$ の微分を計算して、微分された角度を求める行列を作ろう。 $\b{H}_x,\b{H}_y$ を直接ｘやｙで微分するのではなく、形状関数で微分しておくてと簡単である。ここは一旦形状関数での微分を計算する。

$\frac{\partial L_1}{\partial L_2}=-1\qquad\qquad\frac{\partial L_1}{\partial L_3}=-1$

に注意して、計算を行うと以下のとおり

$\frac{\partial H_x^T}{\partial L_2}=\{\begin{array}-(L_2-L_1)t_6+L_3t_5\\+(L_2-L_1)t_6+L_3t_4\\-L_3(t_4+t_5)\\-1+(L_2-L_1)r_6+L_3r_5\\-(L_2-L_1)q_6-L_3q_5\\1+(L_2-L_1)r_6-L_3r_4\\-(L_2-L_1)q_6+L_3q_4\\-L_3(r_4-r_5)\\L_3(q_4-q_5)\\\end{array}\}\qquad\qquad\frac{\partial H_x^T}{\partial L_3}=\{\begin{array}+(L_3-L_1)t_5-L_2t_6\\+L_2(t_4+t_6)\\-L_2t_4-(L_3-L_1)t_5\\-1+L_2r_6+(L_3-L_1)r_5\\-L_2q_6-(L_3-L_1)q_5\\-L_2(r_4-r_6)\\L_2(q_4-q_6)\\1-L_2r_4+(L_3-L_1)r_5\\L_2q_4-(L_3-L_1)q_5\end{array}\}$

$\frac{\partial H_y^T}{\partial L_2}=\{\begin{array}+(L_2-L_1)p_6-L_3p_5\\-(L_2-L_1)p_6-L_3p_4\\+L_3(p_4+p_5)\\-(L_2-L_1)q_6-L_3q_5\\2-6L_1-(L_2-L_1)r_6-L_3r_5\\-(L_2-L_1)q_6+L_3q_4\\-2+6L_2-(L_2-L_1)r_6+L_3r_4\\L_3(q_4-q_5)\\L_3(r_4-r_5)\end{array}\}\qquad\qquad\frac{\partial H_y^T}{\partial L_3}=\{\begin{array}-(L_3-L_1)p_5+L_2p_6\\-L_2(p_4+p_6)\\+L_2p_4+(L_3-L_1)p_5\\-L_2q_6-(L_3-L_1)q_5\\2-6L_1-L_2r_6-(L_3-L_1)r_5\\L_2(q_4-q_6)\\L_2(r_4-r_6)\\L_2q_4-(L_3-L_1)q_5\\-2+6L_3+L_2r_4-(L_3-L_1)r_5\end{array}\}$

但し、

$p_d=-6x_{cb}/l^2_d$

$q_d=3x_{cb}y_{cb}/l^2_d$

$r_d=3y^2_{cb}/l^2_d$

$t_d=-6y_{cb}/l^2_d$

中立面の曲率

ここで中立面の曲率 $\kappa$ テンソルをベクトル表記したものは次のように表されるのであった。

$\{\begin{array}\kappa_x\\ \kappa_y\\ \kappa_{xy}\end{array}\}=\{\begin{array}\frac{\partial \theta_y}{\partial x}\\-\frac{\partial \theta_x}{\partial y}\\\frac{\partial \theta_y}{\partial y}-\frac{\partial \theta_x}{\partial x}\end{array}\}$

$\hat{H}^T_x,\hat{H}^T_y$ を利用して、上は $\theta_x,\theta_y$ の微分を計算すると曲率ベクトルは次のとおりになる。

$\{\begin{array}\kappa_x\\ \kappa_y\\ \kappa_{xy}\end{array}\}=\{\begin{array}\frac{\partial L_2}{\partial x}\frac{\partial H_y}{\partial L_2}+\frac{\partial L_3}{\partial x}\frac{\partial H_y}{\partial L_3}\\ -\frac{\partial L_2}{\partial y}\frac{\partial H_x}{\partial L_2}-\frac{\partial L_3}{\partial y}\frac{\partial H_x}{\partial L_3}\\ \frac{\partial L_2}{\partial y}\frac{\partial H_y}{\partial L_2}+\frac{\partial L_3}{\partial y}\frac{\partial H_y}{\partial L_3}-\frac{\partial L_2}{\partial x}\frac{\partial H_x}{\partial L_2}-\frac{\partial L_3}{\partial x}\frac{\partial H_x}{\partial L_3}\end{array}\}\cdot\b{r}$

これが頂点の撓みとxy回転から、面内の点における曲率を求める行列 $\b{B}$ である。

剛性行列

さて、DKT要素は離散点においてKirchhoff条件が成り立つと仮定して作った要素であった。離散点以外ではKirchhoff条件が成り立たつかどうかは分からない。Kirchhoff条件が成り立たない板の歪エネルギーを計算するためにはMindlin-Reissnerの厚板理論を使って剪断エネルギーも考慮しなければならない。しかし、ここでもう一つ仮定をおく。

DKT要素内の歪エネルギーは剪断エネルギーを無視できる

こではエネルギーを計算するときのみ、要素内全体でKirchhoff仮定が成り立つとしたということである。

歪エネルギーは曲げエネルギ U_b と等しく、を用いて次のように書ける

$U=U_b=\frac{1}{2}\int_{\Omega}\{\begin{array}\kappa_x\\ \kappa_y\\ \kappa_{xy}\end{array}\}^T\b{D}_b\{\begin{array}\kappa_x\\ \kappa_y\\ \kappa_{xy}\end{array}\}d\Omega=\frac{1}{2}\sum_e\int_{\Omega}_e\{\begin{array}\kappa_x\\ \kappa_y\\ \kappa_{xy}\end{array}\}^T\b{D}_b\{\begin{array}\kappa_x\\ \kappa_y\\ \kappa_{xy}\end{array}\}d\Omega\\\qquad\qquad=\frac{1}{2}\sum_e\int_{\Omega_e}\b{r}^T\b{B}^T\b{D}_b\b{B}\b{r}\Omega=\frac{1}{2}\b{r}^T\[\sum_e\int_{\Omega_e}\b{B}^T\b{D}_b\b{B}\Omega\]\b{r}$

明らかに

$\b{K}_e=\int_{\Omega_e}\b{B}^T\b{D}_b\b{B}\Omega$

が要素剛性行列で

$\b{K}=\sum_e \b{K}_e$

が全体剛性行列である。

エネルギーの停留条件 $\delta U=0$ から釣り合いの位置では次を満たす

$\delta\b{r}^T\b{K}\b{r}=\delta\b{r}^T\cdot\b{F}$

ここで $\b{F}$ は $\b{r}$ に対応する右辺ベクトルである。

被積分関数は二次関数なので３点による数値積分で十分である。

[2]によると、積分点１，２，３は辺の中点にあり重みは1/3づつにすればよい。

積分点番号	積分点の位置	重み
1	$(0,\frac{1}{2},\frac{1}{2})$	1/3
2	$(\frac{1}{2},0,\frac{1}{2})$	1/3
3	$(\frac{1}{2},\frac{1}{2},0)$	1/3

非線形有限要素法のためのテンソル解析の基礎	久田俊明　著
Numerical Methods in Structural Mechanics	Zdenek Bittnar, Jiri Sejnoha著

Made by Nobuyuki UMETANI 梅谷　信行

n.umetani@gmail.com