动力学第三章

第2章

function VTB2(m,c,k,x0,v0,tf,w,f0)

%单自由度系统得谐迫振动

clc

wn=sqrt(k/m);

z=c/2/m/wn;

lan=w/wn

wd=wn*sqrt(1-z^2);

A=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);

t=0:tf/1000:tf;

phi1=atan2(x0*wd,v0+z*wn*x0);

phi2=atan2(2*z*lan,1-lan^2);

B=wn^2*f0/k/sqrt((wn^2-w^2)^2+(2*z*wn*w)^2);

x=A*exp(-z*wn*t)、*sin(sqrt(1-z^2)*wn*t+phi1)+B*sin(w*t-phi2); plot(t,x),grid

xlabel('时间(s)')

ylabel('位移')

title('位移与时间得关系')

function VTB1(m,c,k,x0,v0,tf)

%VTB1用来计算单自由度有阻尼自由振动系统得响应

%VTB1绘出单自由度有阻尼自由振动系统得响应图

%m为质量;c为阻尼;k为刚度;x0为初始位移;v0为初始速度;tf为仿真时间

%VTB1(zeta,w,x0,v0,tf)绘出单自由度有阻尼自由振动系统得响应图

%程序中z为阻尼系数ξ;wn为固有频率ωn;A为振动幅度;phi为初相位θ

clc

wn=sqrt(k/m);

z=c/2/m/wn;

wd=wn*sqrt(1-z^2);

fprintf('固有频率为%、3g、rad/s、\n',wd);

fprintf('阻尼系数为%、3g、\n',z);

fprintf('有阻尼得固有频率为%、3g、rad/s、\n',wd);

t=0:tf/1000:tf;

if z<1

A=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);

phi=atan2(x0*wd,v0+z*wn*x0)

x=A*exp(-z*wn*t)、*sin(wd*t+phi);

fprintf('A=%、3g\n',A);

elseif z==1

a1=x0;

a2=v0+wn*x0;

fprintf('a1=%、3g\n',a1);

fprintf('a2=%、3g\n',a2);

x=(a1+a2*t)、*exp(-wn*t);

else

a1=(-v0+(-z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);

a2=(v0+(z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);

fprintf('a1=%、3g\n',a1);

fprintf('a2=%、3g\n',a2);

x=exp(-z*wn*t)、*(a1*exp(-wn*sqrt(z^2-1)*t)+a2*exp(wn*sqrt(z^2-1)*t)); end

plot(t,x),grid

xlabel('时间(s)')

ylabel('位移')

title('位移与时间得关系')

function jzdd

%矩阵迭代法求系统得三阶固有频率与主阵型

clear all

close all

fid1=fopen('A-1','wt'); %建立主振型文件

fid2=fopen('B-1','wt'); %建立固有频率文件

%输入质量矩阵

M(1,1)=2;

M(2,2)=1、5;

M(3,3)=1;

%输入刚度矩阵

K(1,1)=5;K(1,2)=-2;K(2,1)=-2;K(2,2)=3;K(2,3)=-1;K(3,2)=-1;K(3,3)=1

%计算特征值与特征向量

D=inv(K)*M; %原始动力矩阵

A=ones(3,1); %初始振型

for i=1:3 %计算三阶固有频率与主振型

pp0=0;

i

B=D*A;

pp=1、0/B(3); %B(3)为B中得最后一个元素

A=B/B(3);

while abs((pp-pp0)/pp)>1、e-6

pp0=pp;

B=D*A;

pp=1、0/B(3);

A=B/B(3);

end

f=sqrt(pp)/2/pi %固有频率单位转换为Hz

fprintf(fid1,'%20、5f',A); %输入主振型数据

fprintf(fid2,'%20、5f',f); %输入固有频率数据

D=D-A*A'*M/(A'*M*A*pp);

end

fid1=fopen('A-1','rt'); %打开主振型文件

A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵

fid2=fopen('B-1','rt'); %打开固有频率文件

f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵

t=1:3;

h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);

bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),

title('固有频率'),hold on,grid

h1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);

bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),

title('1阶主振型'),hold on,grid

pause(0、1)

h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);

bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),

title('2阶主振型'),hold on,grid

pause(0、1)

h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);

bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),

title('3阶主振型'),hold on,grid

pause(0、1)

end

%chuandijuzhen、m; %传递矩阵方法求固有频率

clear all,clear close

J1=1;J2=1;J3=2;

k2=1100000;k3=1200000;k4=100000;

fid=fopen('chuandi','wt'); %建立(打开)速度文件

M1L=0;

for WN=0:0、01:2000

shita1R=1;M1R=-WN^2*J1;

shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;

shita3R=shita2R+1/k3*M2R;M3R=shita2R*(-WN^2*J3)+(1+(-WN^2*J3)/k3)*M2R; shita4R=shita3R+1/k4*M3R;

if abs(shita4R)<0、005

WN %搜索到得固有频率(rad/s),并显示