Reigh-Ritz法求解欧拉伯努利梁问题

1. 问题描述

一端固定的悬臂梁,同时承受:

• 均布载荷 q = −45 KN/m(全长)
• 集中载荷 P = −100 KN(作用在 x₁ = 4 m 处)

参数	符号	数值
杨氏模量	E	20·10⁶ KPa
梁长	L	6 m
截面宽	b	0.3 m
截面高	h	0.5 m
截面惯性矩	I = bh³/12	0.003125 m⁴

2. Rayleigh-Ritz 法思路

直接解 Euler-Bernoulli 微分方程

$$EI\frac{d^{4}y}{d{x}_1^4}=q(x_1)$$

在有集中载荷、分布载荷混合时,需要奇异函数逐段积分。
Rayleigh-Ritz 绕开这一点,把它转为一个能量最小化问题:

1. 假设挠度形状 $y_{\text{approx}}(x_1)={\sum }_{i=0}^n{a}_i{x}_1^i$(这里取 n = 6)。
2. 施加几何边界条件(固定端 $y(0)=0,y^{\prime }(0)=0$),消去 a₀、a₁。
3. 计算总势能 $\Pi =U-W$(应变能 − 外力做功)。
4. 对剩余系数取变分 $\partial \Pi /\partial a_i=0$,得到线性方程组,解出所有 aᵢ。

3. 外力做功

分布载荷沿梁长积分,集中载荷只在作用点取值:

$$
W = \underbrace{\int_0^L q, y(x_1), dx_1}_{W_1}

• \underbrace{P, y(d)}_{W_2}
  $$
  $W_1$是均布载载荷做功，$W_2$是集中力做功载荷。

4. 内部应变能

$x_1$是中性轴方向，$x_2$是垂直中性轴方向。
$$U={\int }_0^L\frac{EI}{2}{\left(\frac{d^{2}y}{d{x}_1^2}\right)}^{2}d{x}_1$$

为什么应变能写成上述形式:

4.1 应变能密度(三维弹性体)

线弹性材料单位体积的储能:

$$u_0=\frac{1}{2}{\sigma }_{ij}{\epsilon }_{ij}$$

对纯单轴弯曲,只有轴向应力 ${\sigma }_{11}$ 与应变 ${\epsilon }_{11}$ 起作用:

$$u_0=\frac{1}{2}{\sigma }_{11}{\epsilon }_{11}=\frac{1}{2}E{\epsilon }_{11}^2$$

4.2 Euler-Bernoulli 梁假设

“平面截面保持平面、且始终垂直于中性轴” ⇒ 离中性轴距离为 $x_2$ 处的轴向应变
等于曲率 κ 乘以 $-x_2$:

$${\epsilon }_{11}(x_1,x_2)=-x_2\kappa (x_1),\kappa (x_1)\approx \frac{d^{2}y}{d{x}_1^2}$$

(小变形下 $\kappa =y^{\prime \prime }/(1+y^{\prime 2}{)}^{3/2}\approx y^{\prime \prime }$。)

4.3 对截面积分 → 线密度

在 x₁ 处、厚度为 dx₁ 的一小段梁,其储能为

$$dU=\left({\iint }_A{u}_{0}dA\right)d{x}_1=\frac{E}{2}{\kappa }^2{\iint }_A{x}_2^{2}dAd{x}_1$$

其中面积积分正是截面的惯性矩定义:

$$I={\iint }_A{x}_2^{2}dA$$

于是每段弯曲应变能的线密度为

$$\frac{dU}{d{x}_1}=\frac{1}{2}EI{\kappa }^2=\frac{1}{2}EI{\left(\frac{d^{2}y}{d{x}_1^2}\right)}^2$$

4.4 沿梁长积分

最终整根梁的弯曲应变能:

$$\boxed{U={\int }_0^L\frac{EI}{2}{\left(\frac{d^{2}y}{d{x}_1^2}\right)}^{2}d{x}_1}$$

5. 变分与求解

总势能:

$$\Pi (a_2,a_3,a_4,a_5,a_6)=U-W$$

稳定性(能量极小)要求:

$$\frac{\partial \Pi }{\partial a_i}=0,i=2,\dots ,6$$

6. 欧拉伯努利方程-解析解

对 Euler-Bernoulli 方程四次积分,加上端点反力 $R_1=-qL-P$、固端弯矩
$M_1=Pd+q{L}^2/2$,并用奇异函数处理集中载荷:

$$
y_{\text{act}}(x_1) = \frac{1}{E I}!\left[
\frac{M_1 x_1^{2}}{2} + \frac{R_1 x_1^{3}}{6}

• \frac{q x_1^{4}}{24} + \frac{P,\langle x_1 - d\rangle^{3}}{6}
  \right]
  $$

其中 $⟨x_1-d{⟩}^n$ 是 Macaulay 奇异函数（括号函数），当 $x_1<d$ 时取值为 0。

7. 弯矩、剪力、应力

Euler-Bernoulli 梁的内力关系:

$$M(x_1)=EI{y}^{\prime \prime }(x_1),V(x_1)=EI{y}^{\prime \prime \prime }(x_1)$$

弯曲应力、剪应力(矩形截面):

$${\sigma }_{11}=-\frac{M{x}_2}{I},{\sigma }_{12}=\frac{VQ(x_2)}{Ib},Q(x_2)=b\left(\frac{h}{2}-x_2\right)\left(\frac{h}{4}+\frac{x_2}{2}\right)$$

8. 结果分析

• 红色是解析解，黑色为近似解

六次多项式已足够逼近这个问题的精确解。但是剪力图作为更高次的求导结果在 x₁=4 m 处,解析解有阶跃,这是因为用的 Ritz 试函数是全局光滑多项式，而真实解在点载荷处本来就不光滑。提高阶数只能在这个前提下“缓解”，但解决不了本质问题。

• 左边列是解析解，右边列为近似解

通过上图可以看到 Ritz 方法用全局光滑函数逼近高阶导数时，在自由端会产生误差，这是因为 Ritz 方法保证的是“能量最优”，不是“导数局部精确”，而剪力：是三阶导，对局部误差极其敏感，在边界最容易由于积累暴露问题。

感谢阅读​！如果觉得本文对你有帮助，欢迎关注我的csdn账号@「Ahicert」和vx_gzh「工业界的程序猿」。在vx_gzh后台回复 ​“Ritz”关键词​即可免费获取完整的代码示例