三自由系统的动力学分析

石家庄铁道大学SHIJIAZHUANG TIEDAO UNIVERSITY

《振动理论》课程论文

培养单位_ 机械工程学院

学科专业_ 机械电子工程

课程名称振动理论

任课教师李韶华

学生姓名赵

学号

提交日期 2010.01.17

三自由系统的动力学分析

摘要

工程上较复杂的振动问题多数需要用多自由度系统的振动理论来解决。我们熟悉的教材上给出的都是理论求解的方法，本文旨在进行三自由系统的动力学分析。本文将先分析三自由系统的固有振动，其中采用大家熟悉的振型叠加法研究系统的响应，关键是利用Matlab软件求解三自由系统的理论解与数值解，绘图并分析两者的差异和规律。

关键词：三自由系统 Matlab 理论解数值解

Abstract

On the engineering ，more complicated vibration problem need to use multi-freedom degree system to solve. The teaching material that we acquaint with offer the theory method. This text aims at carrying on the dynamics analysis of three-free systems. This text will analyze the proper vibration of three free systems first and adopt fold responding to research system, the key is the theory solution and number-solution that makes use of Matlab software to solve three free systems, paint and analyze the difference and regulation.

Key words:three-freedom degree system Matlab

number-solution theory solution

1

2 引言

一个具有n个自由度的系统，它在任一瞬时的运动形态需要用n个相互耦合的二阶常微分方程组成的方程组。对本文的三自由度的阻尼系统而言，它具有3个固有频率（有可能出现重值），当系统按任意一个固有频率作自由振动时，系统的运动是一种同步运动，称为主振动。系统运动时具有的振动形态称为主振型，或称为模态。利用主坐标，n自由度系统的振动可以当作n个单自由度系统的振动来考虑，然后通过叠加得到原来系统的振动，这种分析方法即为振型叠加法。运用MATLAB 强大、完善的通用数值分析、矩阵运算和图形功能 ,以及MATLAB 控制系统工具箱函数 ,对机械振动系统进行仿真分析 ,可为机械振动系统的分析与设计提供一个强大而方便的工具。

1. 运动微分方程的建立

本文主要采用影响系数方法。无阻尼作用力方程：

M+KX=P(t)

刚度矩阵 K中的元素kij是使系统仅在第 j 个坐标上产生单位位移而相应于第i 个坐标上所需施加的力；质量矩阵 M 中的元素mij是使系统仅在第j 个坐标上产生单位加速度而相应于第i个坐标上所需施加的力，mij，kij分别称为质量影响系数及刚度影响系数，根据它们的物理意义可以直接写出矩阵 M 和 K，从而建立作用力方程，这种方法称为影响系数方法。

引例：如图1-1所示为三自由度的弹簧质量系统，求它的作用力方程

3

图 1-1

由影响系数法不难得到系统的微分方程为

:

2. 求固有频率和振型

令主振动：

或直接用，得

令，得

4

特征方程：（3-α）（α-5α+4）=0得α1=1；α2=3；α3=4

求振型：分别代入α1=1；α2=3；α3=4 的值得

3. 求系统响应

假设图1-1系统中左边第一个质量上作用有激振力P 1（t ）=P

0，其中w=1.7,试求系统的稳态响应。

用振型叠加法，变换成正则坐标下的方程，求解正则坐标下的稳态响应，再变换成物理坐标下的稳态响应，即为系统稳态解，具体过程参看下面理论解的求解方法。

4. 用Matlab 求解数值解和理论解

4.1 理论解

假设k=3，m=3，p 0=2，用Matlab 求解得稳态响应为

X=(2195800*sin((17*t)/10))/6923

-(20000*sin((17*t)/10))/62937

-(2000000*sin((17*t)/10))/692307 ）T

稳态响应曲线如图：

图4-1

4.2 数值解

用Matlab的simulink仿真得到的图形为：

图4-2

4.3 分析比较

在对图4-1和图4-2进行的比较中可以清楚地看出：图4-1理论解得到的是系统的稳态响应曲线，图4-2数值解中不仅可以看到暂态响应的变化还有稳态响应曲线，从图4-2可以清晰地看出在接近90时系统趋于稳定，而且稳定以后和理论解的稳态响应曲线非常一致，这说明理论解求解稳态响应的方法很可靠。

5

附录

求理论解的具体程序

syms p，t

m=3;k=3;p0=2;

M=[m 0 0;0 m 0;0 0 m];K=[3*k -k 0;-k 2*k -k;0 -k 3*k]; P=[p0; 0 ;0]

B=K-p.*M

eval(solve(det(B)))

[F,D]=eig(B)

F=F*[0 1 0;0 0 1;1 0 0] %F为振型

KP=F'*K*F; MP=F'*M*F %主刚度，主质量矩阵

H= F* inv( sqrt(MP))%正则振型

V=MP\KP %谱矩阵

R0=H'*P %正则坐标下激励的幅值

e1=R0(1,1)/(1-1.7^2)*sin(1.7*t) % 正则振型下的稳态解

e2=R0(2,1)/(3-1.7^2)*sin(1.7*t)

e3=R0(3,1)/(4-1.7^2)*sin(1.7*t)

X=H*[e1;e2;e3] %理论解

clear

t=0:0.05:10;%绘图

plot(t,(2195800*sin((17*t)/10))/692307)

hold on

plot(t, -(20000*sin((17*t)/10))/62937,'r') hold on

plot(t, -(2000000*sin((17*t)/10))/692307,'g')

求数值解的具体程序（用simulink仿真的方法（也可以用sim命令做））m=3;k=3;

M=[m 0 0;0 m 0;0 0 m];K=[3*k -k 0;-k 2*k -k;0 -k 3*k]; A=cat(1,cat(2,zeros(3,3),eye(3)),cat(2,-inv(M)*K,ze ros(3,3)));

G=eye(3);B=cat(1, zeros(3,3),inv(M)*G);

C=cat(2,eye(3), zeros(3,3));

D= zeros(3,3);%以上编辑为M文件，仿真时在命令行输入M文件名即可

6