哈工大机械动力学作业

方法1、用模态分析的方法求解： 1、 频响函数矩阵

二自由度系统频响函数矩阵为一2×2方阵，其表达式为：21()()H K M j C ωωω-=-+，其

中9001M ??=????； 2.70.30.30.3C -??=??-??；27333K -??

=??-??

写成矩阵形式：

1

211122

2122()()279 2.730.3()()()30.330.3H H j j H H H j j ωωωω

ωωωωωωω-??-+--??==????---+????

方法2、用模态分析的方法求解：

1）求该二阶系统的一、二阶固有频率及其阻尼比

利用求解该系统振动微分方程()M x C x Kx f t ??

?

++=的特征值i λ来确定固有频率及其阻尼比，下面利用Matlab 来计算固有频率0i ω和阻尼比i ξ：

syms x ;

m1=9; m2=1; k1=24; k2=3;k3=0; c1=2.4;c2=0.3;c3=0; M=[m1 0;0 m2];

C=[c1+c2,-c2;-c2,c2+c3]; K=[k1+k2,-k2;-k2,k2+k3]; y=det(M*x^2+C*x+K);

得到解：

由特征解可以知道：011 1.4142w λ==

1

101

0.1

0.07091.4107

w σζ=

=

=

0222w λ=== 2

202

0.2

0.10051.9900

w σζ=

=

=

1）用解耦的方法

M=[9,0;0,1],C=[2.7,-0.3;-0.3,0.3],K=[27,-3;-3,3];x0=[1;0];v0=[0;0];Ft =[0;1]

Minv2=inv(sqrt(M)); Kt=Minv2*K*Minv2; Ct=Minv2*C*Minv2; [P ,D ]=eig(Kt);

[lambda,I]=sort(diag(D)); P=P(:,I); C1=P'*Ct*P; K1=P'*Kt*P;

ft=P'*inv （M ）*P*Ft;

得到结果：

由振动()Mx Cx Kx F t ++=，令1/2

*x M q -=得：

11()x C x K x f t ++=

由于C1和K1都是对角矩阵，两个振动得到解耦，可以利用单自由度有阻尼有外作用力下的质量弹簧系统分别解得：

0.220.7071sin 30.440.7071sin 3q q q t

q q q t

++=-?++=?

120.07070.40.1

22

ζζ=

===?

12 1.411.99d d ωωωω===

初始条件：1/21/2

-2.12130(0)'(0),(0)'(0)-2.12130q P M x q P M x ????====?

???

????

得解：0.11()(0.06234sin1.41 2.1299cos1.41)0.1003sin30.0086cos3t q t e t t t t -=--++

0.2

2()(0.0084s i n 1.99

2.0892c o s 1.99)0.1337i n 3

0.0321c o s 3t q t e t t t t -=----

1/2()()x t M Pq t -=

0.11()(0.0167sin1.410.9945cos1.41)0.0079sin30.0055cos3t x t e t t t t

-=+++

0.22()(0.0381sin1.99 +0.0288cos1.99)-0.1655sin3-0.0288cos3t x t e t t t t -=

由MATLAB ：

t=0:0.01:40;

x1=-0.2357*(exp(-0.1*t).*(-0.06234*sin(1.41*t)+(-2.1299)*cos(1.41*t))+0.0086*sin(3*t)+0.0086*cos(3*t))-0.2357*(exp(-0.2*t).*(-0.0084*sin(1.99*t)+(-2.0892)*cos(1.99*t))-0.1337*sin(3*t)-0.0321*cos(3*t)); x2=-0.7071*(exp(-0.1*t).*(-0.06234*sin(1.41*t)+(-2.1299)*cos(1.41*t))+0.0086*sin(3*t)+0.0086*cos(3*t))+0.7071*(exp(-0.2*t).*(-0.0084*sin(1.99*t)+(-2.0892)*cos(1.99*t))-0.1337*sin(3*t)-0.0321*cos(3*t)); plot(t,x1,'r',t,x2,'b'); xlabel('time(sec)'); legend('x1','x2'); grid on;

2用振型叠加的方法求解： 由()M x C x Kx F t ??

?

++=，

求固有频率和主振型：{}2

()*0n N K M x ω-=得：

2212121/31/32;4;X ;X 11n n ωω-????

====

? ?????

{}{}{}{}112212X 'X 2,X 'X 2M M M M ====

而X i i

N =可知121/31/3X ,X 11N N -??==????即1/31/3X 11N -?=??

令()*()N x t X q t =，得到q(t)(t)Kq(t)f(t)Cq ++=

0.220.7071sin 30.440.7071sin 3q q q t

q q q t

++=?++=?

121/31/3,11-??==??

?

?v v 120.4

0.0707;0.1

22ζζ=

===?

121.41; 1.99d d ωωωω====,

11

2.12130(0)(0), and (0)(0)2.12130N N X X --????==== ? ?-????

q x q x ，

通解

12()0.1003sin3t 0.0086cos3t ()0.1337sin3t 0.0321cos3t

q t q t =--=--

0.110.22()(0.3645sin1.41 2.1299cos1.41)0.1003sin 3t 0.0086cos3t ()(0.4148sin1.99 2.0892cos1.99t)0.1337sin 30.0321cos3t t t q t e t t q t e t t --=+--=---()()()()()()0.10.210.12x(t)X '*()

x (t)0.2357*(0.3645sin1.41 2.1299cos1.41)0.2357*(0.2074sin1.99-1.9793co 0.2357*0.1003*3*0.0086s1.99t)x (t)0.7071**3*0.2357*0.1337*3*0.031*3*(2N t t t sin t cos t sin t cos t q t e t t e t e ---+---=--==+-0.20.3645sin1.41 2.1299cos1.41)0.7071*(0.2074sin1.99-1.9793cos1.99t) +0.7071*exp(-0.2*t).*(0.274*sin(1.99*t)+(-1.9793)*cos(1.99*t))-0.0709*sin(3*t)+0.0247*cos(3*t)

t t t e t -++

t=0:0.01:40;

x1=0.2357*exp(-0.1*t).*(0.3645*sin(1.41*t)+(2.1299)*cos(1.41*t))-0.2357*exp(-0.2*t).*(0.274*si n(1.99*t)+(-1.9793)*cos(1.99*t))+0.2357*(-0.1003*sin(3*t)-0.0086*cos(3*t))-0.2357*(-0.1337*si n(3*t)-0.0321*cos(3*t));

x2=0.7071*exp(-0.1*t).*(0.3645*sin(1.41*t)+(2.1299)*cos(1.41*t))+0.7071*exp(-0.2*t).*(0.274*s in(1.99*t)+(-1.9793)*cos(1.99*t))+0.7071*(-0.1003*sin(3*t)-0.0086*cos(3*t))+0.7071*(-0.1337*s in(3*t)-0.0321*cos(3*t)); plot(t,x1,'r',t,x2,'b'); xlabel('time(sec)'); legend('x1','x2');

grid on;

方法3：基于状态空间理论的matlab 动力响应分析方法：

可以看成是一个以下的

M1=9kg M2=1Kg C1=2.4，C2=0.3，C3=0

K1=24N/m ，K2=3N/m ，K3=0

理论依据：多自由度振动系统的运动方程为

)(t f kq q C q

M =++ （1） 引入一个辅助方程

0q M q

M =- （2）

把方程（1）、（2）组合起来，有

?

?????=????????????-+????????????0f(t)q q M 00K q q 0M M C （3） 可表示为 )(t E f Hx x

G =+

（4）

式中

??????=0M M C G ，??????-=M 00K H ，??????=0I E ，?

??

???=q q x 式（4）又可以表示为

)(t Bu Ax x

+= （5） 式中

H G A 1--=，E G B 1-=，)()(t t f u =

式（5）即为n 自由度振动系统的状态方程。我们可以根据具体的要求定义相应的输出方程

)(t Du CX y += (6)

注意这里的矩阵C 并非是系统的阻尼矩阵。式（5）、式（6）一起称为系统的状态空间

模型。我们可以利用Matlab的控制工具箱和状态空间理论方便的研究机械振动系统的固有特性与振动响应。

由MATLAB：

k1=24;k2=3;k3=0;

m1=9;m2=1;

c1=2.4;c2=0.3;c3=0;

M=[m1,0;0,m2];

K=[k1+k2,-k2;-k2,k2+k3];

C=[c1+c2,-c2;-c2,c2+c3];

n=size(M,1);

Z=zeros(n);

G=[C M;M Z];

H=[K,Z;Z,-M];

E=[eye(n);Z];

ssA=-inv(G)*H;

ssB=G\E;

ssC=in[1,0,0,0;0,1,0,0];

ssD=0*inv(M);

t=0:0.1:40;

f1=0*t;

f2=1*sin(3*t);

u=[f1;f2];

sys=ss(ssA,ssB,ssC,ssD);

x0=[1,0,0,0];

x=lsim(sys,u,t,x0);

plot(t,x(:,1),'r--',t,x(:,2));

legend('x1','x2');

grid on

验证分析1：直接法认证：

t=0:0.1:40;

s11=-0.1+1.41i;s12=-0.1-1.41i;s21=-0.2+1.99i;s22=-0.2-1.99i;

c11=2.7;c12=-0.3;c21=-0.3;c22=0.3;m11=9;m22=1;k11=27;k12=-3;k21=-3;k2 2=3;

u11=(-c12*s11-k12)/(m11*s11^2+c11*s11+k11);

u12=(-c12*s12-k12)/(m11*s12^2+c11*s12+k11);

u21=(-c12*s21-k12)/(m11*s21^2+c11*s21+k11);

u22=(-c12*s22-k12)/(m11*s22^2+c11*s22+k11);

u11=0.3327+0.0001i;u12=0.3327-0.0001i;u21=-0.3333-0.0000i;u22=-0.3333 +0.0000i;

M=[u11,u12,u21,u22;1,1,1,1;u11*s11,u12*s12,u21*s21,u22*s22;s11,s12,s2 1,s22];

B=[0.9931;0.1662;(-3*0.0083);(3*0.0277)];

A1= 0.7871 - 0.0574i;A2= 0.7871 + 0.0574i;A3= -0.7040 + 0.0510i;A4=-0.7040 - 0.0510i;

x1=u11*A1*exp(s11*t)+u12*A2*exp(s12*t)+u21*A3*exp(s21*t)+u22*A4*exp(s 22*t)+ 0.0031*cos(3*t)+0.0083*sin(3*t);

x2=A1*exp(s11*t)+A2*exp(s12*t)+A3*exp(s21*t)+A4*exp(s22*t)

-0.1662*cos(3*t)-0.0277 *sin(3*t);

plot(t,x1,'r--',t,x2,'b-');

legend('x1','x2');

grid on

2用simscape建立实际的物理模型来进行仿真分析：

得到的结果为：