基于Rayleigh-Ritz法求解均布荷载下弹性梁的挠度值（Mathematica推导）

本篇推文分享的是夸克工作室的《有限元分析基础篇ANSYS与Mathematica》中的Rayleigh-Ritz法介绍，觉得蛮有意思，就分享出来供大家茶余饭后小看几眼~本文以均布荷载下弹性梁的挠度值为例，分别将三角函数、幂级数、形函数作为一系列的基函数，应用Rayleigh-Ritz法对其求解，并附上Mathematica代码供读者验证推导。

原理

在Rayleigh-Ritz方法中，首先假设符合于边界条件的试探函数（Trial Solution Function），并将其函数代入势能函数中，再对试函数的各项未知系数作微分并令之为0，最终求得试函数的各项系数。

求解均布荷载下弹性梁的挠度值

求解步骤：

1. 1. 构造试函数      ;

2. 2. 挠度值表示为试函数的线性组合，其中$c_i$为$i$项的系数，挠度 

3.      
   

4. 3. 将上式带入到梁单元势能函数中：

5.      

6.      
   

7. 4. 最后将势能      取最小值，即对函数的系数      取偏微分，并令之为0，得到未知系数      的值，带回到上式得出梁元挠度表达式。

以三角函数为试探函数：

当    为偶数值时，若令    则函数值为0，不满足悬臂梁的边界条件要求，因此，假设试探函数中的    为奇数。

程序绘制出位置    从0到    的挠度值，如下图所示。图中显示位移值约为1.2，而理论值应该为0.21368，虽然曲线趋势是正确的，然而数值仍然有相当的差距。若将梁体划分为两个“单元块”，则误差将会大大降低，读者可自行验证。

图1 三角函数为基函数的单元素弹性梁挠度曲线

以幂级数为试探函数：

由于悬臂梁要求固定的边界条件，约束点上的位移与斜率都是0，故幂级数中的  为0或者1时，由原点位移值  可知，  的零次项与1次项系数为0，可假设试探函数从二次项开始。

程序绘制出位置    从0到    的挠度值，如下图所示。图中显示位移值约为0.8，而理论值应该为0.21368，虽然曲线趋势是正确的，然而数值仍然有相当的差距。若将梁体划分为两个“单元块”，则误差将会大大降低，读者可自行验证。

图2 幂函数为基函数的单元素弹性梁挠度曲线

以形函数为试探函数：

由于悬臂梁要求固定的边界条件，约束点上的位移与斜率都是0，故幂级数中的 为0或者1时，由原点位移值 与斜率值可知可知，第一项与第三项系数为0，可知试探函数只有第二项和第四项。

程序绘制出位置    从0到    的挠度值，如下图所示。图中显示位移值约为0.8，而理论值应该为0.21368，虽然曲线趋势是正确的，然而数值仍然有相当的差距。若将梁体划分为两个“单元块”，则误差将会大大降低，读者可自行验证。

图3形函数为基函数的单元素弹性梁挠度曲线