《有限元分析基础教程》（曾攀）笔记二-梁单元有限元方程推导
不得不说，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值即可。
