动力学程序讲解

动力学程序（附程序的详细翻译）

第一章

第8页，有阻尼单自由度问题

%第8页，有阻尼单自由度问题

function VTB1(m,c,k,x0,v0,tf) %m是质量，c,k刚度，x0质量块位移，v0质量块速度clc %清屏

wn=sqrt(k/m);

z=c/2/m/wn;%ζ=n/wn,n=c/2m

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

fprintf('固有频率为%.3g.rad/s.\n',wn);%输出wn

fprintf('阻尼系数为%.3g.\n',z);%输出ξ

fprintf('有阻尼的固有频率为%.3g.rad/s.\n',wd);%输出wd

t=0:tf/10000:tf;%时域的长短，决定很坐标的长短

if z<1 %判断ζ的取值，如果ζ<1,弱阻尼情况，按以下1-10到1-12公式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);%输出A

fprintf('phi=%.3g\n',phi);%输出ψ

elseif z==1 %临界阻尼情况，按照1-14到1-15公式计算

a1=x0;% C1

a2=v0+wn*x0;%C2

fprintf('a1=%.3g\n',a1);%输出a1

fprintf('a2=%.3g\n',a2);%输出a2

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

else %过阻尼，按照1-13公式计算

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

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

fprintf('a1=%.3g\n',a1);%输出a1

fprintf('a2=%.3g\n',a2);%输出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 on

xlabel('时间(s)')

ylabel('位移(m)')

title('位移相对时间的关系')

单自由度系统的谐迫振动（P11业）

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

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

clc %清屏

wn=sqrt(k/m);

z=c/2/m/wn;%ζ=n/wn,n=c/2m

lan=w/wn %λ的求法1-18公式

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

fprintf('固有频率为%.3g.rad/s.\n',wn); %输出Wn

fprintf('阻尼系数为%.3g.\n',z);%输出ζ

fprintf('有阻尼的固有频率为%.3g.rad/s.\n',wd); %输出Wd

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

t=0:tf/1000:tf;

phi1=atan(x0*wd/(v0+z*wn*x0));%按1-12求ψ

phi2=atan(2*z*lan/(1-lan^2));%求Φ

b=wn^2*f0/k/sqrt((wn^2-w^2)+(2*z*wn*w)^2);%求B的稳态响应的振幅

x=a*exp(-z*wn*t).*sin(sqrt(1-z^2)*wn*t+phi1)+b*sin(w*t-phi2);%响应方程

plot(real(t),real(x)),grid

xlabel('时间(s)')

ylabel('位移(cm)')

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

第二章

一、矩阵迭代法（P38业，例2-3）

function jzdd %矩阵迭代法

clear all %清空当前所有数据

close all %关闭当前所有的绘图窗口

fid1=fopen('A-1','wt'); %建立第一个名为“A-1”的可写文档

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; %inv表示对矩阵取逆，公式2-65

A=ones(3,1);%定义一个初始迭代阵型，ones()函数表示是3个位为1的单列矩阵，ones(i,j)

%则是i行两列都是j都是1的数组！在这方法中一般取ones(i,1),i=质量各数

for i=1:3 %进入for循环，将要执行书本上37业的迭代，循环次数为1到i，没执行一次循环代表一个质量块的主阵型

pp0=0; %令初始的P（0）=0，并且输出Ｐ（０）

i %输出i

B=D*A; %代入2-60的公式

pp=1.0/B(3); %同上，B(3),是你当时定义主阵型时，第i个元素全为1，当然也可以是第一个

A=B/B(3); %同上

%以上三段代码，是初始赋值

while abs((pp-pp0)/pp)>1e-6 % 判断是否满足迭代的条件|P（k）^2-P(k-1)^2|<δ

%满足后退出while循环

pp0=pp; %P(k)^2=上诉最终的pp值

B=D*A;

pp=1.0/B(3);

A=B/B(3);

%真正地迭代代码

end %结束while循环

f=sqrt(pp)/2/pi %公式2-61

fprintf(fid1,'%20.5f',A); %fid为文件句柄，指定要写入数据的文件，format是用来控制所写数据格式的格式符，与fscanf函数相同，A是用来存放数据的矩阵。

fprintf(fid2,'%20.5f',f); %储存f

D=D-A*A'*M/(A'*M*A*pp); %下一阶的动力矩阵【D*】，A’代表A的逆矩阵

end%结束for循环

fid1=fopen('A-1','rt');%读取A-1文档

A=fscanf(fid1,'%f',[3,3]);%输入3x3的A矩阵，来源为fid1

fid2=fopen('B-1','rt');%读取B-1文档

f=fscanf(fid2,'%f',[3,1]);%输入（拾取）3x1的f矩阵,来源是fid2

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 %bar代表绘制长条图t代表起点横坐标，f(t,1)纵坐标,输出各阶频率阶次的的固有频率

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

bar(t,A(t,1)),xlabel('自由度（质量块）'),ylabel('振型向量'),%输出对应的X阶主振型

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

pause(0.1)%暂停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)%暂停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