基于MATLAB多自由度单向串联振动系统的计算分析

于 翔，周 松

(沈阳航空航天大学 机电工程学院，辽宁 沈阳 110136)

基于MATLAB多自由度单向串联振动系统的计算分析

于 翔，周 松

(沈阳航空航天大学 机电工程学院，辽宁 沈阳 110136)

在工程振动问题中，固有特性(固有频率和主振型)是描述振动系统的主要参数。根据矩阵迭代法，利用MATLAB/GUI界面设计平台，提出了针对多自由度单向串联振动系统固有特性的计算方法，并通过算例求解和误差分析，完成了合理性验证，为复杂振动系统响应的初步解耦奠定了基础。

MATLAB；多自由度；振动系统；矩阵迭代法；固有特性

振动力学不仅是近代应用力学的一个重要分支，而且是现今工业生产中分析和解决机械产品运转稳定和寿命持久问题的重要手段。随着科学技术的高速发展，振动力学在机械加工、航空航天以及交通运输等工业技术领域中占有愈来愈重要的地位[1]。其中，固有特性(固有频率和主振型)作为描述振动系统的主要参数[2-5]，在响应模拟和仿真预测中发挥着重要的作用。因此，通过研究振动系统固有特性的求解方法，可以为防止系统共振提供理论依据[6-8]，也为进行系统响应的初步解耦和机械结构的强度分析奠定基础[9]。然而，对于工程振动问题，无论是复杂的运动体还是离散的弹性体都可以通过有限元法化解为理想节点系统的振动模型，但是由于系统自由度不唯一，并且以多自由度单向串联振动系统最为常见，因此如何准确地求得振动系统任意阶固有特性的有效解，一直都是研究振动力学问题的关键。

1 固有特性的求解方法

多自由度振动系统固有特性的算法主要有直接计算和近似计算2种。其中，对低自由度振动系统，直接计算作为首选方法，既简便快捷又准确有效，但是在实际工程中，大多数系统的振动形式非常复杂且自由度很高，用频率方程(特征方程)很难进行直接求解。然而，通过相关研究发现，近似计算中的矩阵迭代法因为其精确度不依赖于假设振型，假设初值的好坏只影响迭代次数，与所求固有特性的精度无关，即迭代运算总是收敛于振动系统的最低阶固有特性[1,10]。因此，对于实际振动问题，往往采用矩阵迭代的方法求解系统的固有特性。

在多自由度正定系统的自由振动中，系统的主振型方程可分别表示为：

ω2{A}=[M]-1[K]{A}

(1)

(2)

式中：ω为固有频率；{A}为主振型；[M]为质量矩阵；[K]为刚度矩阵；[δ]为柔度矩阵。矩阵迭代法是从假设主振型出发，对以上两式进行矩阵迭代运算。因为工程上对系统的最低阶或较低阶固有特性比较重视[2-3]，所以通常只对式(2)进行迭代运算。引进动力矩阵：

[D]=[δ][M]

(3)

通过矩阵迭代，在规定的有效位数内，当发现{A}k≈{A}k-1时，运算结束。此时，{A}k或{A}k-1即为系统第一阶主振型{A(1)}的近似值；系数ak即为系统第一阶固有频率平方倒数的近似值，即：

{A(1)}≈{A}k-1={A}k

(4)

(5)

利用式(2)不仅可以求出最低阶固有特性，而且还能依次求得较低阶固有特性。其方法是在求系统第(s+1)阶固有特性时，每次迭代前都乘以清型变换后的动力矩阵[D]s。根据式(2)，则有：

[D]s=[D][Q]s=

(6)

式中：ωj为第j阶固有频率；{A(j)}为第j阶主振型；Mj为系统中第j个物体的质量。通过上述方法得到清型变换后的动力矩阵[D]s。然后，根据第一阶固有特性的算法求出第(s+1)阶固有特性，并以此类推求解系统的任意阶固有特性，其对应的程序流程图如图1所示。

图1 矩阵迭代法求解固有特性的程序流程图

2 程序设计

MATLAB 作为一款数据分析软件，对矩阵迭代运算尤其合适[3，11]。通过GUI界面设计，能为用户提供人性化的操作界面[12]，便于针对多自由度振动系统的计算分析。

2.1GUI界面操作

根据矩阵迭代法，基于MATLAB软件设计的计算多自由度单向串联振动系统固有特性的GUI操作界面如图2所示。通过输入各单元的质量和刚度以及振动系统的初始假设振型，就可以自动生成系统本身固有的质量矩阵和刚度矩阵。同时，可以根据用户需要，求出系统的任意阶固有特性。

图2 软件操作界面

2.2MATLAB程序

根据矩阵迭代法计算系统固有特性的关键在于能否求出动力矩阵，而动力矩阵是由质量矩阵和柔度矩阵进行矩阵乘法运算得到的[2]。对于多自由度单向串联振动系统，因为柔度矩阵是拐角矩阵[12]，而拐角矩阵具有沿主对角线分层的特点，如图3所示，所以在程序语言的编写过程中，根据单层循环结构很难完成矩阵建立。

图3 柔度矩阵的特有形式拐角矩阵

通过对多自由度单向串联振动系统的研究发现，柔度矩阵主对角线上元素为各自由度上刚度倒数值的累加和，因此通过双层循环嵌套和矩阵元素的行列转换，可以准确建立任意阶数的柔度矩阵。其程序指令如下：

i=1;w=0;

while i<=n %层间循环

j=i;w=w+1/k(i,1); %刚度倒数求和

while j<=n %层内循环

K0(j,i)=w; K0(i,j)=w; %主对角线对称元素赋值

j=j+1;

end

i=i+1;

end

由于矩阵迭代法的主要计算方式就是归一化后根据程序设计精度反复求积和不断赋值，所以在MATLAB程序语言的编写过程中，通过矩阵中列与列的对比、运算和变换，可以在提高运算效率的同时减少计算空间，避免程序语句过于复杂而引起的变量重复调用。其对应程序段如下：

%初始假设振型归一化

a=A0(end,1);A0(:,1)=A0(:,1)./a;

D=K0*M0; %初始动力矩阵

t=1;

while t<=s %求解第s阶固有频率和主振型

A(:,1)=A0(:,1); %初始假设振型赋值

A(:,2)=D*A(:,1);a=A(end,2);A(:,2)=A(:,2)./a;

%矩阵迭代法

while abs(sum(A(:,1)-A(:,2)))>=0.000001

A(:,1)=A(:,2);A(:,2)=D*A(:,1);

a=A(end,2);A(:,2)=A(:,2)./a;

end

M=A(:,2)'*M0*A(:,2); %第s阶主质量

%清型变换后的动力矩阵

D=D-(A(:,2)*A(:,2)'*M0).*a./M;

t=t+1;

end

a=sqrt(1/a);Ai=A(:,2); %第s阶固有特性

3 算例分析

在如图4所示的三自由度简化振动系统中，假设不考虑单位之间的相互换算，即各振动单元的刚度k和质量m均为1，并且各单元在起始时刻都处于零位[3]。根据MATLAB程序获得的系统第三阶固有特性见表1～3。

图4 三自由度简化振动系统

表1 设计精度为10-3时的计算值与迭代值对比表

表2 设计精度为10-4时的计算值与迭代值对比表

表3 设计精度为10-5时的计算值与迭代值对比表

通过以上数据表的对比可以发现，矩阵迭代法对于计算多自由度单向串联振动系统固有频率的影响很小，但是对于主振型的影响很大。随着系统阶数的不断提高，误差呈现出逐级放大的趋势。这可能是因为在矩阵迭代法中，无论求解第几阶主振型都需要其对应的上一阶主振型参与计算，从而造成误差累积。因此为了提高计算精度，需将程序的设计精度至少提高两个数量级，例如，三自由度振动系统的设计精度应为10-5。但是，这在一定程度上增加了计算量，不过现有的高性能计算机可以很好地弥补矩阵迭代法的不足并发挥其长处。

根据MATLAB/GUI软件程序获得的三自由度振动系统的各单元振动图像如图5所示。通过观察可以发现，由于各单元之间的相互影响，运动形式不再是正弦函数，并且单元排列阶数越高，其波形越不稳定，如图5中第三阶主振型(单元M3密集虚线)所示，从而验证了随自由度增加系统复杂性提高的结论[7，9]。

图5 各单元振动模拟图像

4 结束语

本文采用MATLAB软件，根据矩阵迭代法，计算出多自由度单向串联振动系统的固有特性，提出了针对特殊矩阵以及系统任意阶固有特性的编程方法。通过对三自由度振动系统的算例分析，对比了直接计算法和矩阵迭代法的区别以及各单元的模拟振动响应曲线。计算结果表明，随着系统阶数的提高，矩阵迭代法的误差呈现出逐级放大的趋势，但是通过提高设计精度，可以有效降低误差。同时，根据图像分析可以发现，随着单元阶数等级的提高，其对应单元的波形越不稳定，从而验证了随自由度增加系统复杂性提高的结论，为设计人员分析多自由度系统自振特性及受迫振动提供了参考依据，具有一定的实际应用价值。

[1] 刘习军. 工程振动与测试技术[M]. 天津：天津大学出版社，1999：1-57.

[2] 文涛，胡青春. 基于MATLAB语言的多自由度振动系统的固有频率及主振型计算分析[J]. 中国制造业信息化，2007，36(1)：2-19.

[3] 潘延亮，董大伟，闫兵，等. 发动机隔振系统振动固有特性的优化计算及影响因素分析[J]. 车用发动机，2012(2)：17-25.

[4] 周斌. 内燃机弹性基础隔振系统固有频率偏差分析[J]. 西南交通大学学报，1988，33(5)：503-507.

[5] 王秀，刘初升，彭利平，等. 含横向导支端Euler-Bernoulli梁弯曲振动的固有特性分析[J]. 机械设计，2014，31(2)：81-85.

[6] 张鹏，安子军. 摆线钢球行星传动动力学建模与固有特性分析[J]. 中国机械工程，2014，25(2)：157-162.

[7] 赵克刚，李健，姚伟浩. 车辆传动系统固有特性的识别[J]. 重庆理工大学学报：自然科学，2014，28(4)：5-9.

[8] 唐贵基，陈秀娟. 旋转机械轴系扭振固有特性分析[J]. 汽轮机技术，2010，52(5)：339-344.

[9] 张宪，何洋，钟江，等. 疲劳振动试验台的模态与谐响应分析[J]. 机械设计与制造，2008(4)：12-14.

[10] 黄廷祝. 特殊矩阵分析及应用[M]. 北京：科学出版社，2007：89-95.

[11] Balakumar B，Jame H D. MATLAB原理与工程应用[M].2版. 北京：电子工业出版社，2006：292-341.

[12] 邱金蕙，王矞辉，李振全. 基于MATLAB/GUI的新型界面开发方式[J]. 河北工业科技，2008，25(4)：16-23.

InherentCharacteristicsoftheOne-wayTandemVibrationSystemwithMultipleDegreesofFreedomBasedonMATLAB/GUIInterface

YU Xiang, ZHOU Song

(Shenyang Aerospace University, Liaoning Shenyang, 110136, China)

Inherent characteristics such as natural frequencies and mode shapes are the main parameters of vibration system in engineering vibration field. In this paper, it utilizes the MATLAB/GUI interface design platform and proposes a calculation method for one-way tandem vibration system with multiple degrees of freedom based on matrix iteration. In addition, it validates the method's rationality by solving practical examples and analyzing errors. The method lays a foundation for preliminary decoupling the response of complex vibration systems.

MATLAB; Multiple Segrees of Freedom; Vibration System; Matrix Iteration; Inherent Characteristics

10.3969/j.issn.2095-509X.2014.09.015

2013-11-07

于翔(1988—)，男，山东烟台人，沈阳航空航天大学硕士研究生，主要研究方向为振动载荷下材料疲劳寿命预测。

O327；TB115

A

2095-509X(2014)09-0062-04