本ポストはこちらの続きです。
梁の曲げひずみエネルギー
梁の曲げひずみエネルギー \(U\) は、たわみ \(w\) の2階微分(曲率)を使って次のように表せます。
\[U = \frac{EI}{2}\int_0^{L_e} \left(\frac{d^2w}{dx^2}\right)^2 dx\]
梁要素の3次エルミート補間の形状関数
梁要素の3次エルミート補間の形状関数(要素内座標 \(\xi\) は \(0\)〜\(1\)、\(\xi=0\) が節点1、\(\xi=1\) が節点2に対応、\(L_e\) は要素の長さ)は次のとおりです。
\[\begin{eqnarray}
N_1(\xi) &=& 1 - 3\xi^2 + 2\xi^3\\
N_2(\xi) &=& L_e\left(\xi - 2\xi^2 + \xi^3\right)\\
N_3(\xi) &=& 3\xi^2 - 2\xi^3\\
N_4(\xi) &=& L_e\left(-\xi^2 + \xi^3\right)
\end{eqnarray}\]
節点自由度ベクトル
1つの梁要素には両端に2つの節点があり、節点1には「たわみ \(w_1\)」と「傾き \(\theta_1\)」の自由度が、節点2には「たわみ \(w_2\)」と「傾き \(\theta_2\)」の自由度があるため、1つの梁要素の自由度は計4個となり、節点自由度ベクトル \(d\) として表すと以下のとおりになります。
\[\mathbf{d} = \begin{bmatrix}w_1\\ \theta_1\\ w_2\\ \theta_2\end{bmatrix}\]
前述の4つの形状関数それぞれに対応する節点自由度ベクトルの各要素は以下のとおりなります、
\[N_1:w_1=1,\theta_1=0,w_2=0,\theta_2=0\\
N_2:w_1=0,\theta_1=1,w_2=0,\theta_2=0\\
N_3:w_1=0,\theta_1=0,w_2=1,\theta_2=0\\
N_4:w_1=0,\theta_1=0,w_2=0,\theta_2=1\]
たわみ \(w\) を形状関数と節点自由度で表す
要素内のたわみ \(w(\xi)\) は形状関数の重ね合わせになりますので、節点自由度ベクトル \(\mathbf{d} = \begin{bmatrix}w_1\\ \theta_1\\ w_2\\ \theta_2\end{bmatrix}\) と形状関数の行ベクトル \(\mathbf{N}(\xi) = \begin{bmatrix}N_1 & N_2 & N_3 & N_4\end{bmatrix}\) から、
\[w(\xi) = \mathbf{N}(\xi)\,\mathbf{d}\] と表せます。
形状関数の2階微分
べき乗の微分公式 \(\dfrac{d}{d\xi}\xi^n = n\xi^{n-1}\) を利用します。
\(N_1(\xi) = 1 - 3\xi^2 + 2\xi^3\)
1階微分(各項を個別に微分)
\[\frac{dN_1}{d\xi} = \underbrace{0}_{1の微分} - \underbrace{6\xi}_{3\xi^2の微分} + \underbrace{6\xi^2}_{2\xi^3の微分} = -6\xi + 6\xi^2\]
2階微分(もう一度微分)
\[\frac{d^2N_1}{d\xi^2} = -6 + 12\xi\]
\(N_2(\xi) = L_e(\xi - 2\xi^2 + \xi^3)\)
\(L_e\) は \(\xi\) に無関係な定数ですので、微分の外に出したまま、中身だけ微分します。
1階微分
\[\frac{dN_2}{d\xi} = L_e\left(\underbrace{1}_{\xiの微分} - \underbrace{4\xi}_{2\xi^2の微分} + \underbrace{3\xi^2}_{\xi^3の微分}\right) = L_e(1 - 4\xi + 3\xi^2)\]
2階微分
\[\frac{d^2N_2}{d\xi^2} = L_e(-4 + 6\xi)\]
\(N_3(\xi) = 3\xi^2 - 2\xi^3\)
1階微分
\[\frac{dN_3}{d\xi} = \underbrace{6\xi}_{3\xi^2の微分} - \underbrace{6\xi^2}_{2\xi^3の微分} = 6\xi - 6\xi^2\]
2階微分
\[\frac{d^2N_3}{d\xi^2} = 6 - 12\xi\]
\(N_4(\xi) = L_e(-\xi^2 + \xi^3)\)
1階微分
\[\frac{dN_4}{d\xi} = L_e\left(\underbrace{-2\xi}_{-\xi^2の微分} + \underbrace{3\xi^2}_{\xi^3の微分}\right) = L_e(-2\xi + 3\xi^2)\]
2階微分
\[\frac{d^2N_4}{d\xi^2} = L_e(-2 + 6\xi)\]
まとめ
\[\begin{eqnarray}
N_1'' &=& -6 + 12\xi\\
N_2'' &=& L_e(-4 + 6\xi)\\
N_3'' &=& 6 - 12\xi\\
N_4'' &=& L_e(-2 + 6\xi)
\end{eqnarray}\]
座標を \(x\) から \(\xi\) に変換する(\(d^2w/dx^2\) の書き換え)
\(\xi = x/L_e\)(つまり \(x = L_e\xi\))という関係がありますので、微分の連鎖律により、
\[\frac{d}{dx} = \frac{1}{L_e}\frac{d}{d\xi} \quad\Longrightarrow\quad \frac{d^2}{dx^2} = \frac{1}{L_e^2}\frac{d^2}{d\xi^2}\]
したがって、
\[\frac{d^2w}{dx^2} = \frac{1}{L_e^2}\frac{d^2\mathbf{N}}{d\xi^2}\,\mathbf{d} = \frac{1}{L_e^2}\begin{bmatrix}N_1'' & N_2'' & N_3'' & N_4''\end{bmatrix}\mathbf{d}\]
積分変数を \(x\) から \(\xi\) に変換する
\(x = L_e\xi\) より \(dx = L_e\,d\xi\)、積分区間は \(x:0\to L_e\) が \(\xi:0\to1\) に対応します。
梁の曲げひずみエネルギー \(U\) に代入
\[U = \frac{EI}{2}\int_0^{L_e}\left(\frac{d^2w}{dx^2}\right)^2 dx = \frac{EI}{2}\int_0^{1}\left(\frac{1}{L_e^2}\begin{bmatrix}N_1'' & N_2'' & N_3'' & N_4''\end{bmatrix}\mathbf{d}\right)^2 L_e\,d\xi\]
\(\left(\cdot\right)^2\) の部分は「スカラー」ですので、
\[(\text{スカラー})^2 = (\text{スカラー})^{\mathrm T}(\text{スカラー})\]
と書き直せます。
\(\begin{bmatrix}N_1''&N_2''&N_3''&N_4''\end{bmatrix}\mathbf{d}\) の転置は \[\mathbf{d}^{\mathrm T}\begin{bmatrix}N_1''\\N_2''\\N_3''\\N_4''\end{bmatrix}\] ですので、
\[U = \frac{EI}{2}\cdot\frac{L_e}{L_e^4}\int_0^1 \mathbf{d}^{\mathrm T}\begin{bmatrix}N_1''\\N_2''\\N_3''\\N_4''\end{bmatrix}\begin{bmatrix}N_1'' & N_2'' & N_3'' & N_4''\end{bmatrix}\mathbf{d}\;d\xi\]
\(\mathbf{d}\) は \(\xi\) に無関係(節点の値は積分中は定数)ですので、積分の外に出せます。
\[U = \frac{1}{2}\,\mathbf{d}^{\mathrm T}\left[\frac{EI}{L_e^3}\int_0^1 \begin{bmatrix}N_1''\\N_2''\\N_3''\\N_4''\end{bmatrix}\begin{bmatrix}N_1'' & N_2'' & N_3'' & N_4''\end{bmatrix}d\xi\right]\mathbf{d}\]
剛性行列 \(K_e\) を定義する
\[K_e = \frac{EI}{L_e^3}\int_0^1 \begin{bmatrix}N_1''\\N_2''\\N_3''\\N_4''\end{bmatrix}\begin{bmatrix}N_1'' & N_2'' & N_3'' & N_4''\end{bmatrix} d\xi\]
剛性行列 \(K_e\) の積分部分を求める
\(K_e\) の \((i,j)\) 成分は \(\dfrac{EI}{L_e^3}\displaystyle\int_0^1 N_i''(\xi)N_j''(\xi)\,d\xi\) で計算されます。
\[N_1''=-6+12\xi,\quad N_2''=L_e(-4+6\xi),\quad N_3''=6-12\xi,\quad N_4''=L_e(-2+6\xi)\]
かつ
\[N_3''=-N_1''\]
から各成分は以下のとおりになります。
(1,1)成分: \(\int_0^1 (N_1'')^2\,d\xi\)
\[(N_1'')^2=(-6+12\xi)^2=36-144\xi+144\xi^2\] \[\int_0^1(36-144\xi+144\xi^2)d\xi=36-72+48=12\]
(1,2)成分: \(\int_0^1 N_1''N_2''\,d\xi\)
\[N_1''N_2''=L_e(-6+12\xi)(-4+6\xi)=L_e(24-84\xi+72\xi^2)\] \[\int_0^1 L_e(24-84\xi+72\xi^2)d\xi=L_e(24-42+24)=6L_e\]
(1,3)成分: \(\int_0^1 N_1''N_3''\,d\xi\)
\(N_3''=-N_1''\) ですので、
\[\int_0^1 N_1''\cdot(-N_1'')\,d\xi=-\int_0^1(N_1'')^2 d\xi=-12\]
(1,4)成分: \(\int_0^1 N_1''N_4''\,d\xi\)
\[N_1''N_4''=L_e(-6+12\xi)(-2+6\xi)=L_e(12-60\xi+72\xi^2)\] \[\int_0^1 L_e(12-60\xi+72\xi^2)d\xi=L_e(12-30+24)=6L_e\]
(2,2)成分: \(\int_0^1 (N_2'')^2\,d\xi\)
\[(N_2'')^2=L_e^2(-4+6\xi)^2=L_e^2(16-48\xi+36\xi^2)\] \[\int_0^1 L_e^2(16-48\xi+36\xi^2)d\xi=L_e^2(16-24+12)=4L_e^2\]
(2,3)成分: \(\int_0^1 N_2''N_3''\,d\xi\)
\(N_3''=-N_1''\) ですので、(1,2)成分の結果を使って、
\[\int_0^1 N_2''\cdot(-N_1'')\,d\xi=-6L_e\]
(2,4)成分: \(\int_0^1 N_2''N_4''\,d\xi\)
\[N_2''N_4''=L_e^2(-4+6\xi)(-2+6\xi)=L_e^2(8-36\xi+36\xi^2)\] \[\int_0^1 L_e^2(8-36\xi+36\xi^2)d\xi=L_e^2(8-18+12)=2L_e^2\]
(3,3)成分: \(\int_0^1 (N_3'')^2\,d\xi\)
\(N_3''=-N_1''\) ですので \((N_3'')^2=(N_1'')^2\)。よって(1,1)成分と同じ:
\[\int_0^1(N_3'')^2 d\xi=12\]
(3,4)成分: \(\int_0^1 N_3''N_4''\,d\xi\)
\(N_3''=-N_1''\) ですので、(1,4)成分の結果を使って、
\[\int_0^1(-N_1'')\cdot N_4''\,d\xi=-6L_e\]
(4,4)成分: \(\int_0^1 (N_4'')^2\,d\xi\)
\[(N_4'')^2=L_e^2(-2+6\xi)^2=L_e^2(4-24\xi+36\xi^2)\] \[\int_0^1 L_e^2(4-24\xi+36\xi^2)d\xi=L_e^2(4-12+12)=4L_e^2\]
積分結果のまとめ(対称性より残り6成分も決定)
\[\int_0^1 \begin{bmatrix}N_1''\\N_2''\\N_3''\\N_4''\end{bmatrix}\begin{bmatrix}N_1'' & N_2'' & N_3'' & N_4''\end{bmatrix}d\xi=\begin{bmatrix}12 & 6L_e & -12 & 6L_e\\6L_e & 4L_e^2 & -6L_e & 2L_e^2\\-12 & -6L_e & 12 & -6L_e\\6L_e & 2L_e^2 & -6L_e & 4L_e^2\end{bmatrix}\]
以上です。
