不得不说，Mathematica真是个好东西，以前学习有限元的时候，对于书中的方程推导，看到了就看过去了，从没有想过要自己推导一遍，原因是手工推导太复杂。有了MM，原来很复杂的东西突然变得简单了。
1.单元几何描述
上图是纯弯梁单元，长度l，弹模E，面积A，惯性矩I。两个节点1和2的位移列阵为
[
q^{e}=[v_{1}, heta_{1},v_{2}, heta_{2}]^{T}
]
$v$是挠度(defection)，或者叫位移；$ heta$是转角(slope)。需注意的是$v$和$ heta$的方向，一个是向上，一个是逆时针。
两个节点的节点力矩阵为
[
P^{e}=[P_{v1},M_{1},P_{v2},M_{2}]^{T}
]
当然实际情况往往是在梁的长度方向上作用有荷载，而不是只在节点处有，这时就要进行荷载等效，后面会有说明。注意这两个矩阵都是列矩阵。
需要注意的是，节点力矩阵表示的的是节点上的所有的力，不仅包括荷载引起的等效节点力，还包括节点的反力，反力矩等。
2.单元位移场表达
由于有4个位移节点的已知条件，那么假设纯弯曲梁单元的位移挠度函数具有四个待定系数，如下形式
egin{equation}
v(x)=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}
end{equation}
对于两端节点，位移和转角分别为$v_{1}, heta_{1},v_{2}, heta_{2}$，注意挠曲线方程在一点出的导数值即为改点的转角，所以四个边界条件为
$$
egin{cases}
v(0)=v_{1} & v'(0)= heta_{1}\
v(L)=v_{2} & v'(L)= heta_{2}
end{cases}
$$
使用MM求解方程组
将求得的待定系数带入原方程，可得
将四个位移合并同类项，可以得到
即最终的挠曲线方程vfea为
$$
vfea= ext{$ heta $1} left(frac{x^3}{L^2}-frac{2 x^2}{L}+x
ight)+ ext{$ heta $2} left(frac{x^3}{L^2}-frac{x^2}{L}
ight)+ ext{v1} left(frac{2 x^3}{L^3}-frac{3 x^2}{L^2}+1
ight)+ ext{v2} left(frac{3 x^2}{L^2}-frac{2 x^3}{L^3}
ight)
$$
如果令$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)=Ecdot B(x)q^{e}=S(x)q^{e}
$$
其中$S(x)=Ecdot B(x)$，叫做单元的应力矩阵，表示单元的应力与位移的关系。
4.单元势能表达式
单元的势能为
egin{equation} label{势能方程}
varPi^{e}=U^{e}-W^{e}
end{equation}
其中应变能$U^{e}$的表达式为
$$
egin{split}U^{e} & =dfrac{1}{2}int_{0}^{L}int_{A}sigmacdotvarepsilon dAcdot dx\
& =dfrac{1}{2}int_{0}^{L}int_{A}Ecdot B(x)q^{e}cdot B(x)q^{e}dAcdot dx\
& =dfrac{1}{2}int_{0}^{L}int_{A}Ecdot(q^{e})^{T}cdot B(x)^{T}cdot B(x)q^{e}dAcdot dx\
& =dfrac{1}{2}q^{eT}[int_{0}^{L}int_{A}Ecdot B(x)^{T}cdot B(x)dAcdot 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}},Lleft(frac{6x}{L^{3}}-frac{4}{L^{2}}
ight),frac{6}{L^{2}}-frac{12x}{L^{3}},Lleft(frac{6x}{L^{3}}-frac{2}{L^{2}}
ight)]
$$
$$
K^{e}=EIcdotint_{0}^{L}B(x)^{T}cdot B(x)dx=EIcdotleft(egin{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}
ight)=dfrac{EI}{L^{3}}left(egin{array}{cccc}
12 & 6L & -12 & 6L\
6L & 4L^{2} & -6L & 2L^{2}\
-12 & -6L & 12 & -6L\
6L & 2L^{2} & -6L & 4L^{2}
end{array}
ight)
$$
由于所有的力都简化成了节点力，所以外力功为
$$
W^{e}=P^{eT}cdot q^{e}
$$
5.单元的刚度方程
由最小势能原理，对式 ef{势能方程}中的$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值即可。