基于Python的Rayleigh梁系统控制方程仿真异常求助

问题背景 #

针对配备弹性Pasternak内层、两端弹性约束的含移动质量Rayleigh梁系统，控制方程、边界条件及初始条件如下：

控制方程 #

$$
EI \frac{\partial^4 W_1(x,t)}{\partial x^4} + \mu \frac{\partial^2 W_1(x,t)}{\partial t^2} - T_{c1}^r \mu R_1 \frac{\partial^4 W_1(x,t)}{\partial x^2 \partial t^2} + T_{c2}^r \left[ K_1 (W_1(x,t)-W_2(x,t)) - G_1 \left( \frac{\partial^2 W_1(x,t)}{\partial x^2} - \frac{\partial^2 W_2(x,t)}{\partial x^2} \right) \right] = F_1(x,t)
$$

边界条件 #

$$
\frac{\partial^2 W_1(0,t)}{\partial x^2} = 0, \quad \frac{\partial^2 W_1(L,t)}{\partial x^2} = 0
$$
$$
EI \frac{\partial^3 W_i(0,t)}{\partial x^3} = K_1 W_i(0,t), \quad EI \frac{\partial^3 W_i(L,t)}{\partial x^3} = K_1 W_i(L,t) \quad (i=1,2)
$$

初始条件 #

$$
W_1(x,0) = 0, \quad \left. \frac{\partial^2 W_1(x,t)}{\partial t^2} \right|{t=0} = 0
$$
$$
W_2(x,0) = 0, \quad \left. \frac{\partial^2 W_2(x,t)}{\partial t^2} \right|{t=0} = 0
$$

采用有限元法(FEM)、Crank-Nicolson法、有限差分法(FDM)仿真后结果均为直线，说明系统未产生预期动态响应，以下是针对性排查与解决步骤：

排查与解决步骤 #

1. 载荷与移动质量的有效性验证 #

• 确认移动质量的载荷项是否正确引入：移动质量对梁的作用为随时间变化的集中力，形式为 $F_1(x,t) = m g \delta(x - x_m(t)) + m \frac{d^2 W_1(x_m(t),t)}{dt^2}$（$\delta$ 为狄拉克函数）。检查离散化实现：
  • FEM中是否将集中载荷等效为对应节点的载荷向量；
  • FDM中是否在质量当前位置的网格点施加了对应力项；
  • 若未施加任何载荷（包括移动质量），系统在齐次初始条件下会保持零位移，需确认载荷代码是否生效。

2. 初始条件的完整性检查 #

• 动态仿真需明确初始速度（$\frac{\partial W_1(x,0)}{\partial t}$），当前仅给出初始位移和初始加速度为0，若初始速度也设为0且无外部激励，系统无振动响应。可尝试施加微小初始位移或初始速度作为扰动，验证是否产生振动。
• 检查代码中初始条件的赋值逻辑，避免所有初始状态（位移、速度、加速度）全为0。

3. 边界条件的离散化准确性排查 #

• 弯矩边界条件（$\frac{\partial^2 W_1}{\partial x^2}=0$）：
  • FDM中需确认边界点的二阶导数差分格式是否正确实现，比如采用单边差分或调整边界点的方程；
  • FEM中需选用满足弯矩连续的单元（如Hermite三次单元），确保形函数的二阶导数在边界处满足零弯矩条件。
• 弹性约束剪力条件（$EI \frac{\partial^3 W_i}{\partial x^3} = K_1 W_i$）：
  • FDM中需调整边界点的三阶导数差分方程，将弹性约束转化为边界节点的力平衡；
  • FEM中需在边界节点添加对应弹簧刚度的贡献到整体刚度矩阵中，避免误将边界设为固支或自由端。

4. 控制方程的耦合项与特殊项处理 #

• Rayleigh梁转动惯量项（$-T_{c1}^r \mu R_1 \frac{\partial^4 W_1}{\partial x^2 \partial t^2}$）：检查时间-空间混合导数的离散方式，比如Crank-Nicolson法中是否正确处理了二阶时间导数与二阶空间导数的耦合，FEM中是否将该项纳入单元质量矩阵或附加刚度矩阵。
• Pasternak内层耦合项：需确认是否建立了$W_1$与$W_2$的耦合方程组并同步求解，若仅单独求解$W_1$的方程，忽略$W_2$的约束作用，会导致系统响应异常。

5. 物理参数的合理性检查 #

• 确认所有参数（$EI$、$\mu$、$T_{c1}^{r$、$R_1$、$T_{c2}}r$、$K_1$、$G_1$、$L$、移动质量$m$等）的数值是否合理：
  • 若$EI$（弯曲刚度）过大，梁刚性极强无变形；若过小，可能因数值精度问题导致响应被掩盖；
  • 移动质量$m$若设为0，自然无动态载荷输入；
  • 建议对参数进行无量纲化处理，避免因数量级差异过大导致数值舍入误差。

6. 数值方法的稳定性与收敛性验证 #

• FDM/Crank-Nicolson法：检查时间步长$\Delta t$与空间步长$\Delta x$是否满足稳定性条件，显式FDM需符合CFL条件，Crank-Nicolson法虽无条件稳定，但步长过大可能导致数值阻尼过大，抑制振动响应。
• FEM：确认单元数量是否足够，单元过少会导致离散精度不足，无法捕捉梁的振动模态，结果趋近于直线。

7. 代码调试技巧 #

• 简化问题验证：先移除移动质量和Pasternak内层，求解简单Rayleigh梁的自由振动，确认能得到预期振动曲线后，再逐步添加复杂项排查问题。
• 输出中间结果：打印初始时刻、第一个时间步的位移数组，检查是否有非零值；输出离散后的系统矩阵（刚度矩阵、质量矩阵），排查是否存在奇异或全零的行/列。

内容的提问来源于stack exchange，提问作者Racheal Babatunde