关于梁挠曲线方程解析解与有限元解的简单对比及其他

一、对比

对于受均布荷载的简支梁，假设梁的长度L=1，均布荷载大小为1，弹模E为1，惯性矩I为1，那么梁的挠曲线方程的解析解为

$$v(x)=-\left(\frac{x^4}{24}-\frac{x^3}{12}+\frac{x}{24}\right)$$

梁单元的形函数如下

根据形函数的定义，只要知道了单元两端节点的位移，即可用形函数求出单元上任意位置的挠度。即

$$v(x)=N(x)\cdot q^{e}$$

其中$q^{e}$为节点的位移列阵。

下面对比一下，由形函数得到的挠度曲线与挠曲线的解析解之间的差别，为了对比方便，这里采用形函数时，节点的位移均由挠曲线解析解求出，换言之，在节点处的位移是精确的。

1.假设将上图结构中梁看作一个单元

$$\begin{split} v(0)=0\quad & v'(0)=-\dfrac{1}{24}\\ v(1)=0\quad & v'(1)=\dfrac{1}{24} \end{split}$$

把形函数行矩阵与上面的节点位移列阵相乘，可以得到挠曲线方程

$$v1(x)=\frac{x^2}{24}-\frac{x}{24}$$

2.假设将上图中梁看着2个单元

每个单元的长度L=1/2
利用MM求得挠曲线方程的分段函数为

分别把解析解，1个单元，2个单元求得的挠曲线方程画出来

可以看到，采用2个单元时，挠曲线方程与解析解十分接近，几乎重合；而采用1个单元时，较解析解差别较大。考察在1/4处的数值

类型	1/4处挠度值
解析解	-0.00927734
1个单元	-0.0078125
2个单元	-0.00911458

差别很明显。

 

二、假设采用有限元法计算简支梁两端节点的位移，同样测试采用1个单元和采用2个单元所得的结果

1、采用1个单元

单刚矩阵为

$$\left( \begin{array}{cccc} 12 & 6 & -12 & 6 \\ 6 & 4 & -6 & 2 \\ -12 & -6 & 12 & -6 \\ 6 & 2 & -6 & 4 \\ \end{array} \right)$$

由于$v_{1}$和$v_{2}$均为0，所以最终的求解矩阵为

$$\begin{bmatrix}4 & 2\\ 2 & 4 \end{bmatrix}\begin{bmatrix}\theta_{1}\\ \theta_{2} \end{bmatrix}=\begin{bmatrix}-\frac{1}{12}\\ \frac{1}{12} \end{bmatrix}$$

求得$\theta_{1}=-\frac{1}{24}$，$\theta_{2}=\frac{1}{24}$

可以看到，和解析解方程求得该处的转角是相等的。

2、采用2个单元

由于长度L为1/2，所以每个单元的单刚矩阵均为

$$\left( \begin{array}{cccc} 96 & 24 & -96 & 24 \\ 24 & 8 & -24 & 4 \\ -96 & -24 & 96 & -24 \\ 24 & 4 & -24 & 8 \\ \end{array} \right)$$

对于3个节点的位移列阵，$v_{1}$，$v_{3}$均为零，删除相对应的矩阵行和列，最终的求解矩阵为

$$\begin{bmatrix}8 & -24 & 4 & 0\\ -24 & 192 & 0 & 24\\ 4 & 0 & 16 & 4\\ 0 & 24 & 4 & 8 \end{bmatrix}\begin{bmatrix}\theta_{1}\\ v_{2}\\ \theta_{2}\\ \theta_{3} \end{bmatrix}=\begin{bmatrix}-\frac{1}{48}\\ -\frac{1}{2}\\ 0\\ \frac{1}{48} \end{bmatrix}$$

求得最终结果为

$$\begin{bmatrix}\theta_{1}\\ v_{2}\\ \theta_{2}\\ \theta_{3} \end{bmatrix}=\begin{bmatrix}-\frac{1}{24}\\ -\frac{5}{384}\\ 0\\ \frac{1}{24} \end{bmatrix}$$

可以看到求得的结果与解析解求得结果是一样的。

3、总结

有限元法也是采用了试函数，属于近似方法，为什么在节点处求得的结果与解析解是一样的？

具体原因我也回答不上来，我个人的理解是因为试函数是采用了4次多项式，而解析解同样也是4次多项式，但这也无法解释上述问题。期待解答。

假如有限元方法在所有单元节点处的位移是精确的，那么把单元划分的越多，结果就与真实解约接近。