不得不說，Mathematica真是個好東西，以前學習有限元的時候，對於書中的方程推導，看到了就看過去了，從沒有想過要自己推導一遍，原因是手工推導太復雜。有了MM，原來很復雜的東西突然變得簡單了。
1.單元幾何描述
上圖是純彎梁單元，長度l，彈模E，面積A，慣性矩I。兩個節點1和2的位移列陣為
\[
q^{e}=[v_{1},\theta_{1},v_{2},\theta_{2}]^{T}
\]
$v$是撓度(defection)，或者叫位移；$\theta$是轉角(slope)。需注意的是$v$和$\theta$的方向，一個是向上，一個是逆時針。
兩個節點的節點力矩陣為
\[
P^{e}=[P_{v1},M_{1},P_{v2},M_{2}]^{T}
\]
當然實際情況往往是在梁的長度方向上作用有荷載，而不是只在節點處有，這時就要進行荷載等效，后面會有說明。注意這兩個矩陣都是列矩陣。
需要注意的是，節點力矩陣表示的的是節點上的所有的力，不僅包括荷載引起的等效節點力，還包括節點的反力，反力矩等。
2.單元位移場表達
由於有4個位移節點的已知條件，那么假設純彎曲梁單元的位移撓度函數具有四個待定系數，如下形式
\begin{equation}
v(x)=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}
\end{equation}
對於兩端節點，位移和轉角分別為$v_{1},\theta_{1},v_{2},\theta_{2}$，注意撓曲線方程在一點出的導數值即為改點的轉角，所以四個邊界條件為
$$
\begin{cases}
v(0)=v_{1} & v'(0)=\theta_{1}\\
v(L)=v_{2} & v'(L)=\theta_{2}
\end{cases}
$$
使用MM求解方程組
將求得的待定系數帶入原方程，可得
將四個位移合並同類項，可以得到
即最終的撓曲線方程vfea為
$$
vfea=\text{$\theta $1} \left(\frac{x^3}{L^2}-\frac{2 x^2}{L}+x\right)+\text{$\theta $2} \left(\frac{x^3}{L^2}-\frac{x^2}{L}\right)+\text{v1} \left(\frac{2 x^3}{L^3}-\frac{3 x^2}{L^2}+1\right)+\text{v2} \left(\frac{3 x^2}{L^2}-\frac{2 x^3}{L^3}\right)
$$
如果令$\zeta=\frac{x}{L}$，上式中位移前的系數組成的矩陣稱之為形函數矩陣，也就是常說的形函數。
即
$$
v(x)=N(x)q^{e}
$$
3.單元應變場，應力場的表達
應變的表達式為
$$
\varepsilon=-yv''(x)=B(x)q^{e}
$$
其中$B(x)=-yN''(x)$，$B(x)$叫做單元的幾何矩陣，表示應變與位移的幾何關系。
應力的表達式為
$$
\sigma(x)=E\cdot B(x)q^{e}=S(x)q^{e}
$$
其中$S(x)=E\cdot B(x)$，叫做單元的應力矩陣，表示單元的應力與位移的關系。
4.單元勢能表達式
單元的勢能為
\begin{equation} \label{勢能方程}
\varPi^{e}=U^{e}-W^{e}
\end{equation}
其中應變能$U^{e}$的表達式為
$$
\begin{split}U^{e} & =\dfrac{1}{2}\int_{0}^{L}\int_{A}\sigma\cdot\varepsilon dA\cdot dx\\
& =\dfrac{1}{2}\int_{0}^{L}\int_{A}E\cdot B(x)q^{e}\cdot B(x)q^{e}dA\cdot dx\\
& =\dfrac{1}{2}\int_{0}^{L}\int_{A}E\cdot(q^{e})^{T}\cdot B(x)^{T}\cdot B(x)q^{e}dA\cdot dx\\
& =\dfrac{1}{2}q^{eT}[\int_{0}^{L}\int_{A}E\cdot B(x)^{T}\cdot B(x)dA\cdot dx]q^{e}\\
& =\dfrac{1}{2}q^{eT}K^{e}q^{e}
\end{split}
$$
注意由於$B(x)$和$q^{e}$是單行單列矩陣，所以相乘結果與轉置后調換點乘順序后的結果一樣。
對於中間的一大坨，用$K^{e}$代替，$K^{e}$就是單元剛度矩陣。利用MM求它的具體結果，如下
$$
B(x)=N’’(x)=[\frac{12x}{L^{3}}-\frac{6}{L^{2}},L\left(\frac{6x}{L^{3}}-\frac{4}{L^{2}}\right),\frac{6}{L^{2}}-\frac{12x}{L^{3}},L\left(\frac{6x}{L^{3}}-\frac{2}{L^{2}}\right)]
$$
$$
K^{e}=EI\cdot\int_{0}^{L}B(x)^{T}\cdot B(x)dx=EI\cdot\left(\begin{array}{cccc}
\frac{12}{L^{3}} & \frac{6}{L^{2}} & -\frac{12}{L^{3}} & \frac{6}{L^{2}}\\
\frac{6}{L^{2}} & \frac{4}{L} & -\frac{6}{L^{2}} & \frac{2}{L}\\
-\frac{12}{L^{3}} & -\frac{6}{L^{2}} & \frac{12}{L^{3}} & -\frac{6}{L^{2}}\\
\frac{6}{L^{2}} & \frac{2}{L} & -\frac{6}{L^{2}} & \frac{4}{L}
\end{array}\right)=\dfrac{EI}{L^{3}}\left(\begin{array}{cccc}
12 & 6L & -12 & 6L\\
6L & 4L^{2} & -6L & 2L^{2}\\
-12 & -6L & 12 & -6L\\
6L & 2L^{2} & -6L & 4L^{2}
\end{array}\right)
$$
由於所有的力都簡化成了節點力，所以外力功為
$$
W^{e}=P^{eT}\cdot q^{e}
$$
5.單元的剛度方程
由最小勢能原理，對式\ref{勢能方程}中的$\varPi^{e}$對$q^{e}$取極小值，可以得到單元的剛度方程
$$
K^{e}\cdot q^{e}=P^{e}
$$
由於上面的這個方程式基於能量的原理建立的，所以是具有普遍性的。注意$P^{e}$是節點所有的力矩陣，包括外荷載的等效節點力和節點的反力。
其實對於怎么把$\varPi^{e}$對$q^{e}$取極小值我現在也不得其解，這個應該是用到了變分原理的東西，有點類似求導，這個問題待以后學習變分原理的時候再說。
6.等效節點荷載
由於有限元方法只能在單元的節點上存在荷載，所以如果在單元上有荷載，要將其轉化為等效的節點荷載，轉換的原則是做功相等。比如對於一個受到均布荷載的梁單元
可以將其等效成節點的力和彎矩
這時，所有等效之前外荷載在單元上所做的功，應該等於等效后節點力在節點位移上所做的功。
單元的撓度位移場$v(x)=N(x)q^{e}$。
等效之前外荷載所做的功
$$
W=\int_{0}^{L}p(x)v(x)dx=[\int_{0}^{L}(-p_{0})N(x)dx]\cdot q^{e}
$$
等效之后外荷載所做的功
$$
W=P^{eT}\cdot q^{e}
$$
由於二者相等，所以
$$
P^{eT}=\int_{0}^{L}(-p_{0})N(x)dx==[-\dfrac{p_{0}L}{2},-\dfrac{p_{0}L^{2}}{12},-\dfrac{p_{0}L}{2},\dfrac{p_{0}L^{2}}{12}]
$$
注意上面的這個是均布荷載$p_{0}$在梁單元上的等效節點荷載，方程建立的過程與兩端的約束是無關的，所以等效節點力也與兩端的約束無關。不可將“等效節點荷載”與“荷載引起的節點反力”混淆，比如對於簡支梁，在兩端的鉸支座上是沒有彎矩的，但是等效的節點彎矩在節點上是存在的，相當於等效的節點彎矩引起了端點出的轉角。
如果是在某一點上的集中力進行等效，連積分都不用，直接用力F乘以該點的位移，該點的位移由形函數給出，並給出形函數中的x值即可。