《有限元分析基础教程》（曾攀）笔记二-梁单元方程推导（二）：简支梁挠曲线近似解
一、“近似”的两种分类
一个复杂的函数，可以通过一系列的“基底函数”（base function）的组合来近似，也就是函数逼近，有两种典型的方法：
- 基于全域的逼近，如傅立叶级数展开；
- 基于子域的分段函数组合，如有限元方法。
第一种函数逼近方式，就是力学分析中经典的瑞利-里兹方法（Rayleigh-Ritz），这种方法的特点是基底函数比较复杂，一般是高阶连续函数，通常仅需采用前面几阶函数组合即可得到较高的逼近精度，比如展开为傅立叶级数。
第二种函数逼近方式就是现代力学分析中的有限元思想，即“分段逼近”，每一个分段函数一般比较简单，采用线性函数或者二次函数即可，但是需要较多的分段才能得到逼近效果，工作量比较大。
关于两种逼近方法更加形象的说明，可以参考曾攀老师书中的一个图片，也就是下面这张图
1.1 利用基于全域的逼近方法求近似解
瑞利-里兹法的核心观点就是选取“试函数”，这个试函数必须首先满足位移边界条件，当然还会带有一些待定系数，然后将其他的变量都用这个“试函数”来表达，通过其他的边界条件或者能量方法（虚功原理或极小势能原理）来求解待定系数。
1.1.1 利用虚功原理求解近似解
简单来说，虚功原理的含义就是如果系统有一个虚位移，那么外力在虚位移上所做的外虚功应该内力所做的内虚功。
内虚功
\begin{equation}
\delta U=\intop_{\varOmega}\sigma_{x}\delta\varepsilon_{x}d\varOmega
\end{equation}
外力虚功
\begin{equation}
\delta W=\intop_{l}p\delta vdx
\end{equation}
由于选取的“试函数”只需满足位移边界条件即可，所以“试函数”的选取条件非常宽泛。但是不同的试函数求得最终结果的精度确实相差很大的，所以选取一个合理的试函数非常关键。
比如，对于纯弯曲的简支梁而言，它的位移边界条件就是在0和L处的位移为零，而且位移的形状应该是中间大，然后逐渐向两端减小为0。很自然的首先想到的就是二次函数抛物线，另外一种就是三角函数正弦曲线。为简化，这里将梁的长度L均取为1
- 试函数$v1(x)=-x^{2}+x$
- 试函数$v2(x)=\frac{1}{4}\sin(\pi x)$
画出上述两个试函数在0~1之间的图形
可以看到二者都符合位移边界条件，而且形状的大概模样和预想的梁挠度曲线是一致的。
当试函数采用$c_{1}\sin\left(\frac{\pi x}{L}\right)$时，根据虚功原理，可以求得系数c1，得到最终的挠曲线为
\begin{equation}
\frac{4L^{4}p_{0}}{\pi^{5}\text{EI}}\sin\left(\frac{\pi x}{L}\right)
\end{equation}
当试函数采用$c_1 \left(L x-x^2\right)$时，得到的最终挠曲线为
\begin{equation}\frac{L^{2}p_{0}}{24\text{EI}}(Lx-x^{2})
\end{equation}
令$\frac{p_{0}}{EI}=1$和$L=1$，分别画出0~L直接解析解、采用两种试函数得到的近似解的图形
可以看出采用$c_{1}\sin\left(\frac{\pi x}{L}\right)$作为试函数，最终的到曲线v1(x)与解析解v0(x)是非常接近的，在上图中几乎看不出差别，但是采用$c_1 \left(L x-x^2\right)$作为试函数求解出的近似解却和解析解相差比较大，进一步放大图形
可以看到，当图形放大之后，才能看出v0(x)和v1(x)直接的细微差别。
1.1.2 利用极小势能原理求近似解
通过上面采用虚功原理求近似解，可以看出试函数采用三角函数得到的结果更加准确，事实上这就是是傅立叶级数的展开式。
极小势能原理，也就是势能驻值原理，简而言之就是系统的势能极小时，系统是稳定的。
我们一般求系统的势能是用应变能减去外力功，即$\varPi=U-W$。
其中U表示应变能，它的表达式为
\begin{equation}
U=\dfrac{1}{2}\intop_{\varOmega}\sigma_{x}\varepsilon_{x}d\varOmega
\end{equation}
W表示外力功，表达式为
\begin{equation}
W=\intop_{l}pvdx
\end{equation}
基于上面虚功原理，可以看出试函数选取三角函数更加合理，事实上，由于解析解是采用幂级数形式的方程，采用傅立叶级数（三角函数）进行近似更加合理，后面我们也会看到，实际上试函数的待定系数就是解析解的傅立叶系数。
选取试函数为
\begin{equation}
v(x)=c_{1}sin\dfrac{\pi x}{l}+c_{2}sin\dfrac{2\pi x}{l}+c_{3}sin\dfrac{3\pi x}{l}
\end{equation}
采用Mathematica进行求解过程如下
\[
\text{v3}(\text{x$\_$})\text{:=}c_1 \sin \left(\frac{\pi x}{L}\right)+c_2 \sin \left(\frac{2 \pi x}{L}\right)+c_3 \sin \left(\frac{3 \pi x}{L}\right)
\]
\[
\text{U3}=\frac{1}{2} \text{EI} \int_0^L \text{v3}''(x)^2 \, dx=\frac{\pi ^4 \left(c_1^2+16 c_2^2+81 c_3^2\right) \text{EI}}{4 L^3}
\]
\[
\text{W3}=p_0 \int_0^L \text{v3}(x) \, dx=\frac{2 \left(3 c_1+c_3\right) L p_0}{3 \pi }
\]
\begin{equation}
\Pi =\text{U3}-\text{W3}=\frac{\pi ^4 \left(c_1^2+16 c_2^2+81 c_3^2\right)
\text{EI}}{4 L^3}-\frac{2 \left(3 c_1+c_3\right) L p_0}{3 \pi }
\end{equation}
将势能表达式对c1,c2,c3求驻值，即对各自的偏导为0，可以求出c1，c2，c3。
\[
\text{Solve}\left[\left\{\frac{\partial \Pi }{\partial c_1}=0,\frac{\partial \Pi }{\partial c_2}=0,\frac{\partial \Pi }{\partial c_3}=0\right\},\left\{c_1,c_2,c_3\right\}\right]
\]
\[
\left\{\left\{c_1\to \frac{4 L^4 p_0}{\pi ^5 \text{EI}},c_2\to 0,c_3\to \frac{4 L^4 p_0}{243 \pi ^5 \text{EI}}\right\}\right\}
\]
可以看到c1的值与使用虚功原理求得的值是一样的。由于采用了两项三角函数组成的试函数，最终求得的挠曲线相比虚功原理一项三角函数得到的挠曲线更加精确。根本的原因是这样的，上述c1，c2，c3就是解析解的按照傅立叶级数展开之后的傅立叶系数。
根据周期为$2l$的周期函数的傅立叶级数公式（《高等数学》同济六版下册P316）
由于解析解为
\begin{equation}v(x)=\dfrac{p_{0}}{EI}(\dfrac{1}{24}x^{4}-\dfrac{1}{12}lx^{3}+\dfrac{1}{24}l^{3}x)
\end{equation}
为了简便，将$\dfrac{p_{0}}{EI}$都看作1，得到简化后的挠曲线解析解为
\begin{equation}
v(x)=\dfrac{1}{24}x^{4}-\dfrac{1}{12}lx^{3}+\dfrac{1}{24}l^{3}x\end{equation}
这个函数在[0,1]的图形可以做奇函数的周期拓展，所以傅立叶级数展开式只有正弦项，没有余弦项。其傅立叶系数为
\[
v(\text{x$\_$})\text{:=}\frac{L^3 x}{24}-\frac{L x^3}{12}+\frac{x^4}{24}
\]
\[
\text{b1}=\frac{2 \int_0^L v(x) \sin \left(\frac{\pi x}{L}\right) \, dx}{L}=\frac{4 L^4}{\pi ^5}
\]
\[
\text{b2}=\frac{2 \int_0^L v(x) \sin \left(\frac{2 \pi x}{L}\right) \, dx}{L}=0
\]
\[
\text{b3}=\frac{2 \int_0^L v(x) \sin \left(\frac{3 \pi x}{L}\right) \, dx}{L}=\frac{4 L^4}{243 \pi ^5}
\]
可以看出根据傅立叶公式求出的傅立叶系数与根据势能原理求出的试函数的待定系数是一样一样的。
下面一篇会从有限元方法推导梁单元。To be continue………….
