基于MIMO的悬臂梁振动响应有限元计算原理及应用

“本文介绍了梁的有限元动力学分析基本原理，并基于梁有限元模型，运用MIMO（多输入多输出）算法，计算梁在多个输入力下的振动响应。单自由度质量-弹簧-阻尼系统的振动动力学方程的计算和求解是深入理解本文的基础。”

01

—

目的：计算悬臂梁在多个输入下的振动响应

图1是一个悬臂梁及其有限元参数，在中间位置和悬臂远端均受动态力，求整个悬臂梁在两个（或多个）力下的振动响应（位移，速度，或加速度）。

图1

要进行振动响应计算，首先要建立动力学方程，这就要借助有限元的方法。

02

—

自由状态梁的动力学方程有限元计算原理

本文算法是将图1中的梁划分为36个梁单元。本节仅针对3个梁单元的有限元计算进行举例。

图2是Euler-Bernoulli梁单元的简化模型，梁单元左侧点(Node1)受到f1的力和T1的扭矩，产生v1的位移及theta1的扭转；同样右侧点(Node2)受到f2的力和T2的扭矩，产生v2的位移及theta2的扭转。

图2

构造整个梁模型的动力学方程，最重要的是组装M(质量),C(阻尼),K(刚度)矩阵（学过单自由度质量-弹簧-阻尼系统振动的同学一定不会感到陌生）。整个过程如同将大象装冰箱里一样，分为流畅的三步：

一，翻书查找梁单元理论公式：

理论上（太多公式，不作介绍

），该梁单元的质量矩阵Me，刚度矩阵Ke （上标e代表element）如图3，为了简单方便，将阻尼矩阵Ce设置成比例阻尼（如图3）。那么梁单元的动力学方程就如图3红字部分。

图3

二，将多个单元矩阵组装成整个梁矩阵：

将3个梁单元组装起来，也就是将Me, Ke, Ce 组装成M123（图4），K123（图5），C123（比例阻尼，组装好M123, K123后直接算）。

图4

图5

三，基于组装好的矩阵建立动力学方程：

组装质量矩阵M，刚度矩阵K，是建立动力学方程的关键。图6中红框内的矩阵取阻尼矩阵C=0时，求取矩阵的特征值和特征向量，即该梁自由状态下的共振频率和振型(包括位移和转角)。

图6

对图6中红框内的矩阵求逆，即外力和响应的传递函数矩阵，如图7红框内矩阵。

图7

03

—

悬臂梁的动力学方程有限元计算原理

本文要介绍的是实际应用中比较常见的约束方式：悬臂梁。怎样在模型中加入约束条件呢？方法就在矩阵的构造内：

因为Node0处，位移和转角都是0，所以原先8*8的矩阵（图8），可以不用考虑前2行前2列，变成了6*6的矩阵（图9）。

图8

图9

同样的，对图9中该6*6的矩阵取阻尼矩阵C=0时，求取矩阵的特征值和特征向量，即该悬臂梁的共振频率和振型(包括位移和转角)。

对该矩阵求逆，得到悬臂状态下，外力和响应的传递函数矩阵，如图10右侧6*6矩阵。

图10

更进一步，我们可以忽略各点扭矩的影响，那么该矩阵可以再次从6*6缩减到3*3矩阵，如图11。

图11

04

—

悬臂梁的动力学方程有限元计算算例

用有限元法算出了36*36的传递函数矩阵后，可求出在多个力下的悬臂梁振动响应，如图12。

图12

请注意，该36*36矩阵是包含频率的复数矩阵，形如图13。举例：蓝色表示f=1Hz时的36*36的复数矩阵；红色表示f=2Hz时的36*36的复数矩阵；绿色表示f=3Hz时的36*36的复数矩阵……

图13

将所有的第1行第1列复数矩阵沿着频率方向连起来，即图13中的蓝色线，即是h1,1

将所有的第1行第6列复数矩阵沿着频率方向连起来，即图13中的黑色线，即是h1,6

……

这些连线即是传递函数（复数），包含幅值比和相位差（如图14，15）。

图14

图15

如果将36*36矩阵某一行或者某一列传递函数复数的虚部依次画出来（例如：h1,1; h1,2; h1,3; …… h1,36)，则可以从中看到振型，如图16。

图16

图17是根据图9中的矩阵算出的特征向量（归一化后），特征向量即振型。可以看出，该振型和图16一致。

图17

图18

在f18,f36力下的整个悬臂梁振动响应计算步骤，依然是帅气的三步（参见图18）：

一，对f18, f36进行傅立叶变换，得到复数输入力f18,f36 (包含幅值和相位信息）；

二，将传递函数矩阵（复数）乘以输入力f18,f36（复数），得到梁上各点响应频谱（复数）；

三，将梁上各点响应频谱（复数）反傅立叶变换，得到梁上各点时域振动响应。

05

—

悬臂梁振动响应计算结果及总结

以下直接给出计算结果：

图19

1）图19：

f18=1*sin(2*pi*30*t); f36=1*sin(2*pi*30*t)。这两个力在梁的不同位置，虽然同相位，但是起到互相抵消的作用。合成振动仍为第2阶振型。

图20

2）图20:

f18=1*sin(2*pi*30*t); f36=1*sin(2*pi*30*t+pi)。这两个力在梁的不同位置，相位相反，合力起到放大作用。合成振动仍为第2阶振型。

图21

3）图21：

f18=1*sin(2*pi*30*t); f36=1*sin(2*pi*166*t)。这两个力激发起来悬臂梁的两个固有频率，所以振动特性体现了两个振型的叠加。

总结：结构在多个动态力下的振动响应受到两大因素影响：

1）激励力：包括力的位置，频率，幅值，相位的综合影响。

2）结构本身的动力学特性：特别是共振频率，及共振频率下结构振型。

06

—

悬臂梁案例的引申应用

前面简单介绍了梁的振动响应计算方法。可能大家会问这样的问题：为什么长度取 0.628m，为什么将其划分为36个梁单元？

我们不妨看如下的动画（图22）：

图22

将该梁卷成一圈，正好是半径为0.1m的圆。同时36个单元，每个单元对应的是 10度。

图23

图23是圆内的多个激励力。图24，图25分别是圆的某阶共振频下的振型。

图24

图25

至此，点到为止～