FEM for 2D thin-walled beam

２次元薄肉梁要素の有限要素法解析

左から線形解析曲げ剛性(1)、非線形解析(曲げ剛性1)、非線形解析(曲げ剛性0.1)、非線形解析(曲げ剛性0.01)

Last Update:2009年10月25日

概要

２次元的な領域で変形する薄肉梁要素について解説する。梁要素とは柱や棒のように細長い構造物を、体積を持たない１次元的な分割によって空間離散化する方法である。つまり要素は断面積や断面に次モーメントなどの量をパラメータとしてもつ。また薄肉とは薄い構造物で成り立つ、Kirhhoff-Loveの仮定が成り立つということである。この仮定は分かり易く言えば、変形前に中立面に垂直な物体の中に埋め込まれたベクトルは変形後でも中立面に垂直ということである。薄肉構造物は面外せん断歪みを持たないので、歪みエネルギーは純粋に曲げによる歪みエネルギーだけとなる。曲げのエネルギーは撓みの2回微分で書くことが出来る。このことが薄肉梁の有限要素法において一般的にエルミート補間が必要であることを解説する。

最終的に求めたいのは非線形の薄肉梁要素であるが、幾つか段階を踏むと理解しやすい。2次元梁要素へのもっとも簡単なステップは、線形曲げ要素である。これは梁をｘ軸におき、変形を微小なｙ軸方向のたわみに限定した場合における問題である。次のステップは線形梁要素である。これは線形曲げ要素に梁の軸方向の剛性を加えて、梁の軸方向の変形も取り扱えるようにしたものである。非線形の定式化は、これらを大回転へ発展させて、梁のエネルギーを回転に依存しない形で定義してその微分を計算することで定式化を行う。

　静的問題では歪みポテンシャルエネルギーを導出し、その変分としてつりあいの式を導き出す。また動的問題では歪みポテンシャルエネルギーと力学的エネルギーを導出してハミルトンの原理を用いて運動方程式を導き出す。

静的線形曲げ問題

梁問題を変位法で解析する場合、撓みとその勾配が要素間で連続でなければならないことを示し、もっと単純な例であるエルミート補間について説明する。この補間において梁曲げ要素の定式化を解説する。

両端固定梁

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

w(0)=w(D)=w'(0)=w'(D)=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(D)=w'(0)=w'(D)=0$

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

弱形式化

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

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

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

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

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

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

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

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

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

が成り立つ。

補間関数

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

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

$N^1(0)=1,\;N^1(1)=0,\;\frac{\partial N^1(0)}{\partial\xi}=0,\;\frac{\partial N^1(1)}{\partial\xi}=0\;\; \Right\;\; N^1(\xi) = 1-3\xi^2+2\xi^3$

$N^2(0)=0,\;N^2(1)=0,\;\frac{\partial N^2(0)}{\partial\xi}=1,\;\frac{\partial N^2(1)}{\partial\xi}=0\;\; \Right\;\; N^2(\xi) = \xi-2\xi^2+\xi^3$

$N^3(0)=0,\;N^3(1)=1,\;\frac{\partial N^3(0)}{\partial\xi}=0,\;\frac{\partial N^3(1)}{\partial\xi}=0\;\; \Right\;\; N^3(\xi) = 3\xi^2-2\xi^3$

$N^4(0)=0,\;N^4(1)=0,\;\frac{\partial N^4(0)}{\partial\xi}=0,\;\frac{\partial N^4(1)}{\partial\xi}=1\;\; \Right\;\; N^4(\xi) = -\xi^2+\xi^3$

この時、要素 I_i=[a_i,b_i] の中の撓みの場は次のようにアイソパラメトリック補間される

$x=a_i(1-\xi)+b_i\xi\;\;\;\right\;\;\;dx = (b_i-a_i)d\xi$

$w = N^1 w^1 + N^2 \frac{\partial w^1}{\partial\xi} + N^3 w^2 + N^4 \frac{\partial w^2}{\partial\xi}\\\;\; = N^1 w^1 + N^2 (b_i-a_i)\frac{\partial w^1}{\partial x} + N^3 w^2 + N^4 (b_i-a_i)\frac{\partial w^2}{\partial x}$

曲げエネルギー

曲げによる歪みエネルギー W^b は次のように書くことができる。

$W^b= \frac{1}{2}\int_0^D EI$\frac{\partial^2 w}{\partial x^2}$^2dx = \sum_e U_e$

但し、 W^b_e は以下のような要素内の歪みエネルギーである。

$W^b_e = \frac{1}{2}\int_{a_i}^{b_i} EI$\frac{\partial^2 w}{\partial x^2}$^2dx$

ここで曲げ剛性を要素内一定であるとして積分の外にくくりだす。さらに b_i-a_i=L として表記を簡単にする。 $dx = ld\xi$ となる。これをもちいて積分と微分の変数変換を行って

$W^b_e = \frac{EI}{2}\int_0^1$\frac{1}{L^2}\frac{\partial^2 w}{\partial \xi^2}$^2Ld\xi = \frac{EI}{2L^3}\int_0^1$\frac{\partial^2 w}{\partial \xi^2}$^2d\xi$

ここで、各補間関数の２回微分は次のとおりであった。

$\{\begin{array}{llll}\frac{\partial^2 N^1(\xi)}{\partial \xi^2} &= +6(2\xi-1),\;\;& \;\; \frac{\partial^2 N^2(\xi)}{\partial \xi^2} &= +2(3\xi-2)\\ \frac{\partial^2 N^3(\xi)}{\partial \xi^2} &= -6(2\xi-1),\;\; & \;\; \frac{\partial^2 N^4(\xi)}{\partial \xi^2} &= +2(3\xi-1)\end{array}$

これを代入して

$w = N^1 w^1 + N^2 L\frac{\partial w^1}{\partial x} + N^3 w^2 + N^4 L\frac{\partial w^2}{\partial x}\\\Right\;\; \frac{\partial^2 w}{\partial \xi^2} = \frac{\partial^2 N^1}{\partial\xi^2} w^1 + \frac{\partial^2 N^2}{\partial\xi^2} L\frac{\partial w^1}{\partial x} + \frac{\partial^2 N^3}{\partial\xi^2}w^2 + \frac{\partial^2 N^4}{\partial\xi^2} L\frac{\partial w^2}{\partial x}\\\Right\;\; \frac{\partial^2 w}{\partial \xi^2} = +6(2\xi-1)w^1 +2(3\xi-2)L\frac{\partial w^1}{\partial x} -6(2\xi-1)w^2 +2(3\xi-1) L\frac{\partial w^2}{\partial x}\\\Right\;\; \frac{\partial^2 w}{\partial \xi^2} = \{\begin{array}+6(2\xi-1)\\ +2(3\xi-2)L\\ -6(2\xi-1)\\ +2(3\xi-1)L\end{array}\}^T \{\begin{array}w^1\\ \frac{\partial w^1}{\partial x}\\ w^2\\ \frac{\partial w^2}{\partial x}\end{array}\}\\\Right\;\;\frac{\partial^2 w}{\partial \xi^2} = \b{B}^T\b{U}_e$

但し, $\b{B}$ は頂点での変位とその１次微分と一般座標 $\xi$ における変位の２回微分を結びつける行列である。

$\b{B}^T = \{\begin{array}+6(2\xi-1),\;\; +2(3\xi-2)L,\;\; -6(2\xi-1),\;\; +2(3\xi-1)L\end{array}\}$

また、 $\b{U}_e$ は要素内未知変数ベクトルである。

$\b{U}^T_e = \{\begin{array}w^1,\;\; \frac{\partial w^1}{\partial x},\;\; w^2,\;\; \frac{\partial w^2}{\partial x}\end{array}\}$

これを用いると、曲げエネルギーは次のようになる。

$W^b_e = \frac{EI}{2L^3}\int_0^1$\frac{\partial^2 w}{\partial \xi^2}$^2d\xi = \frac{1}{2}\b{U}^T_e\left[ \frac{EI}{L^3}\int_0^1\b{B}\b{B}^T d\xi\right] \b{U}_e = \frac{1}{2}\b{U}_e^T K^b_e \b{U}_e$

但し、 K^b_e は次のとおり。

$K^b_e = \frac{EI}{L^3}\int_0^1\b{B}\b{B}^T d\xi = \frac{EI}{L^3}\[\begin{array}12 & 6L & -12 & 6L \\ 6L & 4L^2 & -6L & 2L^2\\ -12 & -6L & 12 & -6L \\ 6L & 2L^2 & -6L & 4L^2\end{array}\]$

これらを用いると全体の歪みエネルギー W^b は次のように書くことができる。

$W^b = \frac{1}{2}\b{U}^T\b{K}^b\b{U}$

ここで、 K^b は要素剛性行列 $\b{K}^b_e$ をマージした全体剛性行列であり、 $\b{U}$ は未知ベクトルである。

さて歪みエネルギーの変分を考えて、つりあいの式は次のようになる。

$\b{K}^b\b{U} = \b{F}$

動的線形曲げ要素

力学的エネルギー T^b は次のように書くことができる

$T^b = \frac{1}{2}\int_0^D \rho A$\frac{\partial w}{\partial t}$^2 dx = \sum_e\frac{1}{2}\int_{a_i}^{b_i} \rho A$\frac{\partial w}{\partial t}$^2 dx = \sum_eT^b_e$

$T^b_e = \frac{1}{2}\int_{a_i}^{b_i} \rho A$\frac{\partial w}{\partial t}$^2 dx = \frac{1}{2}\rho AL\int_0^1 $\frac{\partial w}{\partial t}$^2d\xi$

$w = N^1 w^1 + N^2 L\frac{\partial w^1}{\partial x} + N^3 w^2 + N^4 L\frac{\partial w^2}{\partial x}\\\;\;\Right \dot{w} = N^1 \dot{w}^1 + N^2 L\frac{\partial \dot{w}^1}{\partial x} + N^3 \dot{w}^2 + N^4 L\frac{\partial \dot{w}^2}{\partial x}\\\;\;\Right \dot{w} = \b{N}^T \dot{U}_e$

ここで $\b{N},\;\b{\dot{U}}_e$ は次のようにおいた

$\b{N}^T = \{1-3\xi^2+2\xi^3,\;\; \xi-2\xi^2+\xi^3,\;\; 3\xi^2-2\xi^3,\;\; -\xi^2+\xi^3 \}$

$\b{\dot{U}}_e^T = \{\dot{w}^1,\;\; \frac{\partial \dot{w}^1}{\partial x},\;\; \dot{w}^2,\;\;\frac{\partial \dot{w}^2}{\partial x}\}$

これを用いると運動エネルギーは次のように書くことができる。

$T^b_e = \frac{1}{2}\rho AL\int_0^1 $\frac{\partial w}{\partial t}$^2d\xi = \frac{1}{2}\dot{U}^T_e\[\rho AL\int_0^1 N N^T d\xi\]\dot{U}_e = \frac{1}{2}\dot{U}^T_\b{M}_e\dot{U}_e$

$M_e = \frac{\rho AL}{420}\left[\begin{array} 156 & 22L & 54 & -13L\\ 22L & 4L^2 & 13L & -3L^2\\ 54 & 13L & 156 & -22L\\ -13L & -3L^2 & -22L & 4L^2 \end{array}\right]$

これらを用いると全体の力学的エネルギーは次のように書くことができる。

$T^b = \frac{1}{2}\dot{\b{U}}^T\b{M}\dot{\b{U}}$

ここで、は要素剛性行列 $\b{M}_e$ をマージした全体剛性行列であり、 $\b{\dot{U}}$ は未知ベクトルの速度である。

さて、歪みポテンシャルエネルギー W^b と力学的エネルギーは次のように書けるということになる。

$W^b = \frac{1}{2}\b{U}^T\b{K}\b{U},\;\;\;\;T^b = \frac{1}{2}\b{\dot{U}}^T\b{M}\b{\dot{U}}$

ハミルトンの原理を使うと運動方程式は次のように書ける。

$\b{M}\b{\ddot{U}} + \b{K}\b{U} = \b{F}$

２次元線形梁要素

梁曲げの問題では梁の軸方向の歪みも加味して考える。簡単のため、最初にこの梁がｘ軸に平衡であるとして、ｘ方向の変位に対する剛性行列と質量行列を求める。

剛性行列

軸方向歪みによるエネルギー W^m は次のように書ける

$W^m = \frac{1}{2}\int_0^L AE \left(\frac{\partial u}{\partial x}\right)^2 dx = \sum_e\frac{1}{2}\int_{a_i}^{b_i} AE \left(\frac{\partial u}{\partial x}\right)^2 dx = \sum_e W_e^m$

ここで、軸方向変位 u_x は次のように要素内でアイソパラメトリック補間されるのであった。

$u_x = (1-\xi)u^1_x + \xi u^2_x$

これを用いると軸方向の伸縮による歪みエネルギーは節点変位を使って次のようになる。

$W_e^m = \frac{1}{2}\int_{a_i}^{b_i} AE \left(\frac{\partial u}{\partial x}\right)^2 dx = \frac{1}{2}\int_{0}^{1} AE \left(\frac{u^2_x - u^1_x}{L}\right)^2 Ld\xi = \frac{AE}{2L}(u^1_x - u^2_x)^2\\\;\; = \frac{1}{2}\{\begin{array}u^1_x\\ u^2_x\end{array}\}^T \left[ \frac{AE}{L}$\begin{array}1 & -1\\ -1 & 1 \end{array}$\right] \{\begin{array}u^1_x\\ u^2_x\end{array}\} = \frac{1}{2}\{\begin{array}u^1_x\\ u^2_x\end{array}\}^T \b{K}^b_e \{\begin{array}u^1_x\\ u^2_x\end{array}\}$

但し、 $\b{K}^b_e$ は剛性行列である。

$\b{K}^b_e = \frac{AE}{L}$\begin{array}1 & -1\\ -1 & 1 \end{array}$$

質量行列

軸方向の変位速度 $\dot{u}_x$ による力学的エネルギー T^m は次のようになる。

$T^m = \frac{1}{2}\int_0^L A\rho (\dot{u}_x)^2 dx = \sum_e\frac{1}{2}\int_{a_i}^{b_i} A\rho (\dot{u}_x)^2 dx = \sum_e T_e^m$

ここで T_e^m は要素内での力学的エネルギーで、

$T_e^m = \frac{1}{2}\int_{a_i}^{b_i} A\rho (\dot{u}_x)^2 dx = \frac{1}{2}A\rho L\int^1_0 (\dot{u}_x)^2d\xi\\\;\;\; = \frac{1}{2}A\rho L\int^1_0 \{(1-\xi)\dot{u}^1_x + \xi \dot{u}^2_x\}^2d\xi\\\;\;\; = \frac{1}{2}\{\begin{array}u^1_x\\ u^2_x\end{array}\}^T \left[ A\rho L$\begin{array}\frac{1}{3} & \frac{1}{6}\\ \frac{1}{3} & \frac{1}{3} \end{array}$\right] \{\begin{array}u^1_x\\ u^2_x\end{array}\} = \frac{1}{2}\{\begin{array}u^1_x\\ u^2_x\end{array}\}^T \b{M}^b_e \{\begin{array}u^1_x\\ u^2_x\end{array}\}$

但し、 M^b_e は次のような軸方向変位速度に対する要素質量行列である。

$M^b_e = A\rho L$\begin{array}\frac{1}{3} & \frac{1}{6}\\ \frac{1}{3} & \frac{1}{3} \end{array}$$

非線形解析

節点辺りに法線の回転角 $\theta^i$ と変位 $\b{u}^i$ を持つ。

$\b{n}^i = R(\theta^i) \b{N}^i$

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

ここで、

$\b{V} = \b{X}^1 - \b{X}^2\;\;\;\b{v} = \b{x}^1 - \b{x}^2 = V^{12} + \b{u}^1 -\b{u}^2$

とする。要素の変形前の長さと変形後の長さは次のとおり

$l = |\b{x}^1-\b{x}^2| = |\b{v}|$

$L = |\b{X}^1-\b{X}^2| = |\b{V}|$

回転角

要素内の回転角は、微小かつ要素の長さの変化が小さいと仮定して次のように決める。

$\psi^i = \frac{(\b{x}^1-\b{x}^2)\cdot \b{n}^i}{|\b{X}^1-\b{X}^2|} = \frac{1}{L}(\b{v}\cdot \b{n}^i)$

エネルギー

さて、ここで微小歪みを仮定して、要素内の歪みエネルギー W_e が曲げによるエネルギー W^b_e と、軸方向歪みによるエネルギー W^m_e の次のような和であると仮定する。

W_e = W^b_e+W^m_e

$W^b_e =\frac{1}{2} \{\psi^1, \psi^2\}\frac{2EI}{L^3}\[\begin{array}2 & 1\\ 1 & 2 \end{array}\]\{\begin{array}\psi^1\\ \psi^2\end{array}\}$

$W^m_e = \frac{1}{2}\frac{A\rho}{L} (l-L)^2$

回転角の微分

法線ベクトルの微分

$R(\theta) = \left(\begin{array}\cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{array}\right)$

$R(\theta)' = \left(\begin{array}-\sin\theta & -\cos\theta \\ \cos\theta & -\sin\theta \end{array}\right) \;\;\Right\;\; (\b{n}^i)' = R(\theta)'N^i = \left(\begin{array}-n^i_y\\n^i_x\end{array}\right) = \b{n}^i_\bot$

$R(\theta)'' = \left(\begin{array}-\cos\theta & \sin\theta \\ -\sin\theta & -\cos\theta \end{array}\right) \;\;\Right\;\; (\b{n}^i)'' = R(\theta)''N^i = -R(\theta)N^i = -\b{n}^i$

頂点間のベクトルの微分

$\frac{\partial v^{12}}{\partial u^1} = \b{I},\;\;\frac{\partial v^{12}}{\partial u^2} = -\b{I}$

$\frac{\partial^2 v^{12}}{\partial (u^i)^2} = 0,\;\; \frac{\partial^2 v^{12}}{\partial (\theta^i)^2} = 0,\;\; \frac{\partial^2 v^{12}}{\partial \theta^i \partial u^j} = 0$

要素内の回転角の微分

$\frac{\partial\psi^i}{\partial\b{U}_e} = \frac{1}{L} (\frac{\partial\b{v}}{\partial\b{U}_e}\cdot \b{n}^i + \b{v}\cdot \frac{\partial\b{n}^i}{\partial\b{U}_e})$

$\begin{array}{lllll} \frac{\partial\psi^1}{\partial\b{U}_e} = \frac{1}{L}\{\b{n}^1,\;\;& (\b{v}\cdot\b{n}^1_\bot),\;\; & -\b{n}^1,& 0 & \}^T\\ \frac{\partial\psi^2}{\partial\b{U_e}} =\frac{1}{L}\{\b{n}^2,\;\;& 0,\;\;& -\b{n}^2,\;\;&(\b{v}\cdot\b{n}^2_\bot) & \}^T\end{array}$

$\frac{\partial^2\psi^1}{\partial\b{U}_e\partial\b{U}_e} = \frac{1}{L}$\begin{array} 0 & \b{n}^1_\bot & 0 & 0 \\ {\b{n}^1_\bot}^T & -(\b{v}\cdot\b{n}^1) & -{\b{n}^1_\bot}^T & 0 \\ 0 & -\b{n}^1_\bot & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array}$\;\;\; \frac{\partial^2\psi^1}{\partial\b{U}_e\partial\b{U}_e} = \frac{1}{L}$\begin{array} 0 & 0 & 0 & \b{n}^2_\bot \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\b{n}^2_\bot \\ {\b{n}^2_\bot}^T & 0 & -{\b{n}^2_\bot}^T & -(\b{v}\cdot\b{n}^2) \end{array}$$

曲げエネルギーの微分

微分の連鎖則から次のように曲げエネルギーの１回微分（＝曲げに対する内力）が計算できる。

$\b{Q}^b_e = \frac{\partial W^b_e}{\partial U_e} = \frac{\partial W^b_e}{\partial\psi^1}\frac{\partial \psi^1}{\partial U_e} + \frac{\partial W^b_e}{\partial\psi^2}\frac{\partial \psi^2}{\partial U_e} = \frac{2EI}{L^3}\left{$2\psi^1+\psi^2$\frac{\partial \psi^1}{\partial U_e} + $\psi^1+2\psi^2$\frac{\partial \psi^2}{\partial U_e}\right}$

$\b{K}^b_e = \frac{\partial^2 W^b_e}{\partial U_e\partial U_e}\\\;\;\; = \frac{\partial^2 W^b_e}{\partial\psi^1\partial\psi^1}\frac{\partial \psi^1}{\partial U_e}$\frac{\partial \psi^1}{\partial U_e}$^T + \frac{\partial^2 W^b_e}{\partial\psi^1\partial\psi^2}\frac{\partial \psi^1}{\partial U_e}$\frac{\partial \psi^2}{\partial U_e}$^T + \frac{\partial^2 W^b_e}{\partial\psi^1\partial\psi^2}\frac{\partial \psi^2}{\partial U_e}$\frac{\partial \psi^1}{\partial U_e}$^T + \frac{\partial^2 W^b_e}{\partial\psi^2\partial\psi^2}\frac{\partial \psi^2}{\partial U_e}$\frac{\partial \psi^2}{\partial U_e}$^T\;\cdots\;initial\;strain\;stiffness\\\;\;\;\;\;\; + \frac{\partial W^b_e}{\partial \psi^1}\frac{\partial\psi^1}{\partial\b{U}_e\partial\b{U}_e}+ \frac{\partial W^b_e}{\partial \psi^2}\frac{\partial\psi^2}{\partial\b{U}_e\partial\b{U}_e}\;\cdots\;initial\;stress\;stiffness \\\;\;\; = \frac{2EI}{L^3}\left{ 2\frac{\partial \psi^1}{\partial U_e}$\frac{\partial \psi^1}{\partial U_e}$^T + \frac{\partial \psi^1}{\partial U_e}$\frac{\partial \psi^2}{\partial U_e}$^T + \frac{\partial \psi^2}{\partial U_e}$\frac{\partial \psi^1}{\partial U_e}$^T + 2\frac{\partial \psi^2}{\partial U_e}$\frac{\partial \psi^2}{\partial U_e}$^T\right}\;\cdots\;initial\;strain\;stiffness\\\;\;\;\;\;\; + \frac{2EI}{L^3}\left{(2\psi^1+\psi^2)\frac{\partial\psi^1}{\partial\b{U}_e\partial\b{U}_e}+ (\psi^1+2\psi^2)\frac{\partial\psi^2}{\partial\b{U}_e\partial\b{U}_e}\right}\;\cdots\;initial\;stress\;stiffness$

軸方向歪みエネルギーの微分

内力

軸方向の歪みエネルギーは次のように、梁の長さの変化で定義した。

$W^m_e = \frac{1}{2}\frac{AE}{L} (l-L)^2$

このエネルギーの微分を計算するために、梁の長さを両端のベクトルで微分してみよう。

$\frac{\partial l}{\partial \b{v}} = \frac{\partial|\b{v}|}{\partial \b{v}} = \frac{\b{v}}{|\b{v}|}$

これを用いるとエネルギーの梁の両端を結ぶベクトルでのでの微分は次のとおり。

$\frac{\partial W^m_e}{\partial \b{v}} = \frac{AE}{L} (l-L)\frac{\partial l}{\partial \b{v}} = +\frac{AE}{L} (l-L)\frac{\b{v}}{|\b{v}|}$

これを用いるとエネルギーの梁の節点変位における微分は次のとおり。

$\{\begin{array}\frac{\partial W^m_e}{\partial \b{u}^1} = \frac{\partial W^m_e}{\partial \b{v}}\frac{\partial\b{v}}{\partial \b{u}^1} = +\frac{\partial W^m_e}{\partial \b{v}}\\ \frac{\partial W^m_e}{\partial \b{u}^2} = \frac{\partial W^m_e}{\partial \b{v}}\frac{\partial\b{v}}{\partial \b{u}^2} = -\frac{\partial W^m_e}{\partial \b{v}}\end{array}$

最終的に軸方向の歪みに起因する内力は次のようになる。

$\b{Q}_e^m = \{\begin{array}+\frac{AE}{L} (l-L)\frac{\b{v}}{|\b{v}|}\\ 0\\ -\frac{AE}{L} (l-L)\frac{\b{v}}{|\b{v}|}\\ 0 \end{array}\}$

接線剛性行列

接線剛性行列は歪みエネルギーの２回微分として定義される。

２頂点の長さの２頂点を結ぶベクトルで２回微分してみよう。

$\frac{\partial^2 l}{\partial\b{v}\partial\b{v}} = \frac{\b{I}}{|\b{v}|} - \frac{\b{v}\b{v}^T}{|\b{v}|^3}$

これを用いると軸方向エネルギーの２回微分は次のように書ける。

$\frac{\partial^2 W^m_e}{\partial \b{v}\partial \b{v}} = \frac{AE}{L}\{\frac{\b{v}\b{v}^T}{|\b{v}|^2} + (l-L)$\frac{I}{|\b{v}|}-\frac{\b{v}\b{v}^T}{|\b{v}|^3}$\} = \frac{AEL}{l^3}\b{v}\b{v}^T + \frac{AE(l-L)}{Ll}\b{I}$

上の第一項が初期変位剛性、第二項は初期応力剛性に対応している。これらを用いて節点変位での微分を計算すると次のようになる。

$\{\begin{array}{l}\frac{\partial^2 W^m_e}{\partial \b{u}^1\partial \b{u}^1} = \frac{\partial^2 W^m_e}{\partial \b{v}\partial \b{v}}\frac{\partial\b{v}}{\partial \b{u}^1}\frac{\partial\b{v}}{\partial \b{u}^1} = \frac{\partial^2 W^m_e}{\partial \b{v}\partial \b{v}}\\ \frac{\partial^2 W^m_e}{\partial \b{u}^1\partial \b{u}^2} = \frac{\partial^2 W^m_e}{\partial \b{u}^2\partial \b{u}^1} = \frac{\partial^2 W^m_e}{\partial \b{v}\partial \b{v}}\frac{\partial\b{v}}{\partial \b{u}^1}\frac{\partial\b{v}}{\partial \b{u}^2} = -\frac{\partial^2 W^m_e}{\partial \b{v}\partial \b{v}}\\ \frac{\partial^2 W^m_e}{\partial \b{u}^2\partial \b{u}^2} = \frac{\partial^2 W^m_e}{\partial \b{v}\partial \b{v}}\frac{\partial\b{v}}{\partial \b{u}^2}\frac{\partial\b{v}}{\partial \b{u}^2} = \frac{\partial^2 W^m_e}{\partial \b{v}\partial \b{v}}\end{array}$

これらを纏めると、軸方向の歪みに起因する接線剛性行列はは以下のとおり。

$\b{K}^m_e = $\begin{array}\frac{AEL}{l^3}\b{v}\b{v}^T + \frac{AE(l-L)}{Ll}\b{I} & 0 & -\frac{AEL}{l^3}\b{v}\b{v}^T - \frac{AE(l-L)}{Ll}\b{I} & 0\\ 0 & 0 & 0 & 0\\ -\frac{AEL}{l^3}\b{v}\b{v}^T - \frac{AE(l-L)}{Ll}\b{I} & 0 & \frac{AEL}{l^3}\b{v}\b{v}^T + \frac{AE(l-L)}{Ll}\b{I} & 0 \\ 0 & 0 & 0 & 0 \end{array}$$

離散化された方程式とその解法

要素内で次のようにエネルギーが定義されていた。