×

[PR]この広告は3ヶ月以上更新がないため表示されています。
ホームページを更新後24時間以内に表示されなくなります。

|TOP Page|


FEM for 2D thin-walled beam

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




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




Last Update:2009年10月25日


目次

概要

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

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

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

静的線形曲げ問題

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

両端固定梁

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

w(0)=w(D)=w'(0)=w'(D)=0

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

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

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

弱形式化

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

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

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

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

左辺第1,2項はv(0)=v(D)=v'(0)=v'(D)=0から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]となっていることにしよう。つまり、要素の中で、vC^2級関数であるとすると、vの2階微分は区間連続である。さて、vが連続であるためには、節点を区間要素の両端に置けば、隣合う要素同士で値を共有するので良かった。同様にv'が連続であるためには勾配を表す節点を区間要素の両端に置けばよい

要素I_iの中で補間関数の2階微分が定数であるとしよう。つまり、要素間で補間関数は3次関数となる。要素I_iの4つの節点(両端における変位と変位勾配)に対応するアイソパラメトリックな補間関数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

ここで曲げ剛性EIを要素内一定であるとして積分の外にくくりだす。さらに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

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

\{\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}は頂点での変位とその1次微分と一般座標\xiにおける変位の2回微分を結びつける行列である。

\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は次のように書くことができる。

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

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

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

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}

2次元線形梁要素

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

剛性行列

軸方向歪みによるエネルギー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}\)

非線形解析

節点i辺りに法線の回転角\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と変形後の長さlは次のとおり

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}\)

曲げエネルギーの微分

微分の連鎖則から次のように曲げエネルギーの1回微分(=曲げに対する内力)が計算できる。

\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}\}

接線剛性行列

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

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

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

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

\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}\)

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

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

W_e = W^b_e+W^m_e

内力と接線剛性行列は、全体のエネルギーの微分であり、エネルギーは要素ごとに曲げと軸方向の歪みを足し合わせて得られた、次のように定義される。

\b{Q} = \frac{\partial W}{\partial\b{U}} = \sum_e \b{Q}^b_e +\b{Q}^m_e

\b{K} = \frac{\partial^2 W}{\partial\b{U}\partial\b{U}} = \sum_e \b{K}^b_e +\b{K}^m_e

参考にしたもの

Links

東京大学新領域創成科学研究科、久田・杉浦研究室、渡邊先生の講義資料
http://www.sml.k.u-tokyo.ac.jp/members/nabe/
『 固 体 力 学 の 達 人 (基 礎 編) 』中 谷 彰 宏 著
http://www.md.ams.eng.osaka-u.ac.jp/~nakatani/Lectures/Fundamentals_of_Solid_Mechanics/TEXT/HTML/

Books

        

非線形有限要素法のためのテンソル解析の基礎 久田俊明 著
連続体力学―簡明な理論と例題 (1979年) (理工学海外名著シリーズ〈31〉) P.チャドウィック 著
テンソル解析と連続体力学 (1979年) (理工学海外名著シリーズ〈30〉) W.フリューゲ 著


Last Update:2009年10月25日
Made by Nobuyuki UMETANI  梅谷 信行
n.umetani@gmail.com