动力学程序.docx
《动力学程序.docx》由会员分享,可在线阅读,更多相关《动力学程序.docx(21页珍藏版)》请在冰豆网上搜索。
![动力学程序.docx](https://file1.bdocx.com/fileroot1/2023-1/29/17e084f1-a56f-4ca3-b7cc-736f4075b488/17e084f1-a56f-4ca3-b7cc-736f4075b4881.gif)
动力学程序
动力学程序(附程序的详细翻译)
第一章
第8页,有阻尼单自由度问题
%第8页,有阻尼单自由度问题
functionVTB1(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;%时域的长短,决定很坐标的长短
ifz<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);%输出ψ
elseifz==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),gridon
xlabel('时间(s)')
ylabel('位移(m)')
title('位移相对时间的关系')
单自由度系统的谐迫振动(P11业)
functionvtb2(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)
functionjzdd%矩阵迭代法
clearall%清空当前所有数据
closeall%关闭当前所有的绘图窗口
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=质量各数
fori=1:
3%进入for循环,将要执行书本上37业的迭代,循环次数为1到i,没执行一次循环代表一个质量块的主阵型
pp0=0;%令初始的P(0)=0,并且输出P(0)
i%输出i
B=D*A;%代入2-60的公式
pp=1.0/B(3);%同上,B(3),是你当时定义主阵型时,第i个元素全为1,当然也可以是第一个
A=B/B(3);%同上
%以上三段代码,是初始赋值
whileabs((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',[50200420420]);
bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),
title('固有频率'),holdon,grid%bar代表绘制长条图t代表起点横坐标,f(t,1)纵坐标,输出各阶频率阶次的的固有频率
h1=figure('numbertitle','off','name','1','pos',[50200420420]);
bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),%输出对应的X阶主振型
title('1阶主振型'),holdon,grid
pause(0.1)%暂停0.1秒
h1=figure('numbertitle','off','name','2','pos',[50200420420]);
bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('2阶主振型'),holdon,grid
pause(0.1)%暂停0.1秒
h1=figure('numbertitle','off','name','3','pos',[50200420420]);
bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('3阶主振型'),holdon,grid
二、传递矩阵法(P45页,例2-5)
functioncdjz
clearall%清空当前所有数据
closeall%关闭当前所有的绘图窗口
J1=1;J2=1;J3=2;%定义各个转动惯量J
k2=1100000;k3=1200000;k4=100000;%定义各个刚度的具体数值
fid=fopen('chuandi','wt');%建立打开速度文件
M1L=0;%定义变量M1L的初始值为0,M1L为第一个质量块左边的转矩大小
forWN=0:
.01:
2000%for循环
shita1R=1;M1R=-WN^2*J1;%定义θ1,同时利用公式2-73求出质量块右边的转矩大小M1R
shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;%求θ2,M2R,应用到2-78公式,问题是这个公式是θL的
shita3R=shita2R+1/k2*M2R;M3R=shita2R*(-WN^2*J3)+(1+(-WN^2*J3)/k3)*M2R;%依次类推
shita4R=shita3R+1/k4*M3R;%(待分析,不知神马东东)
ifabs(shita4R)<0.005%abs()函数,对括号内的所有元素平方后取绝对值;判断精度是否满足,满足后执行
WN%搜索到的固有频率(rad/s),并显示
shita=[shita1R;shita2R;shita3R;shita4R]%搜索到阵型,并显示
bar(shita),xlabel('对应的质量块'),ylabel('振型向量')%画θ矩阵图形
pause(1.0)%1秒一暂停
end%结束if语句
fprintf(fid,'%30.15f',shita4R);%将θ4输到文档中,数据类型是30位数,15位小数的浮点数
end%结束for循环
fid=fopen('chuandi','rt');%此时fid有返回值,当是正数时代表打开文件成功,-1代表失败,rt表示读写
x=fscanf(fid,'%f',[1,200001]);%要求输入1到2000001的数
t=0:
.01:
2000;
plot(t,x),grid,xlabel('频率rad/s'),ylabel('第四个质量块的转角(rad)'),title('用传递矩阵法求固有频率')
作业(P48-P49)
2-8作业题
functionjzdd2%第2-8作业题
clearall%清空当前所有数据
closeall%关闭当前所有的绘图窗口
fid1=fopen('A-2','wt');%建立第一个名为“A-1”的可写文档
fid2=fopen('B-2','wt');%同上
M(1,1)=1;
M(2,2)=1;
M(3,3)=1;%以上三段代码是为了输入质量矩阵,依题意均为1
K(1,1)=2;
K(1,2)=-1;K(2,1)=-1;
K(2,2)=2;
K(2,3)=-1;K(3,2)=-1;
K(3,3)=1;%输入刚度矩阵,依题意均为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=质量各数
fori=1:
3%进入for循环,将要执行书本上37业的迭代,循环次数为1到i,每执行一次循环代表一个质量块的主阵型
pp0=0;%令初始的P(0)=0,并且输出P(0)
i%输出i
B=D*A;%代入2-60的公式
pp=1.0/B(3);%同上,B(3),是你当时定义主阵型时,第i个元素全为1,当然也可以是第1个
A=B/B(3);%同上
%以上三段代码,是初始赋值
whileabs((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-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-2','rt');%读取A-2文档
A=fscanf(fid1,'%f',[3,3]);%输入3x3的A矩阵,来源为fid1
fid2=fopen('B-2','rt');%读取B-2文档
f=fscanf(fid2,'%f',[3,1]);%输入(拾取)3x1的f矩阵,来源是fid2
t=1:
3;
h1=figure('numbertitle','off','name','0','pos',[50200420420]);
bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),
title('固有频率'),holdon,grid%bar代表绘制长条图t代表起点横坐标,f(t,1)纵坐标,输出各阶频率阶次的的固有频率
h1=figure('numbertitle','off','name','1','pos',[50200420420]);
bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),%输出对应的X阶主振型
title('1阶主振型'),holdon,grid
pause(0.1)%暂停0.1秒
h1=figure('numbertitle','off','name','2','pos',[50200420420]);
bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('2阶主振型'),holdon,grid
pause(0.1)%暂停0.1秒
h1=figure('numbertitle','off','name','3','pos',[50200420420]);
bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('3阶主振型'),holdon,grid
2-9作业题
functionjzdd3
clearall%清空当前所有数据
closeall%关闭当前所有的绘图窗口
fid1=fopen('A-3','wt');%建立第一个名为“A-3”的可写文档
fid2=fopen('B-3','wt');%同上
M(1,1)=4;
M(2,2)=2;
M(3,3)=1;%以上三段代码是为了输入质量矩阵
K(1,1)=4;
K(1,2)=-1;K(2,1)=-1;
K(2,2)=2;
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=质量各数
fori=1:
3%进入for循环,将要执行书本上37业的迭代,循环次数为1到i,每执行一次循环代表一个质量块的主阵型
pp0=0;%令初始的P(0)=0,并且输出P(0)
i%输出i
B=D*A;%代入2-60的公式
pp=1.0/B
(1);%同上,B
(1),是你当时定义主阵型时,第i个元素全为1,当然也可以是第一个
A=B/B
(1);%同上
%以上三段代码,是初始赋值
whileabs((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
(1);
A=B/B
(1);
%真正地迭代代码
end%结束while循环
f=sqrt(pp)%公式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-3','rt');%读取A-3文档
A=fscanf(fid1,'%f',[3,3]);%输入3x3的A矩阵,来源为fid1
fid2=fopen('B-3','rt');%读取B-1文档
f=fscanf(fid2,'%f',[3,1]);%输入(拾取)3x1的f矩阵,来源是fid2
t=1:
3;
h1=figure('numbertitle','off','name','0','pos',[50200420420]);
bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),
title('固有频率'),holdon,grid%bar代表绘制长条图t代表起点横坐标,f(t,1)纵坐标,输出各阶频率阶次的的固有频率
h1=figure('numbertitle','off','name','1','pos',[50200420420]);
bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),%输出对应的X阶主振型
title('1阶主振型'),holdon,grid
pause(0.1)%暂停0.1秒
h1=figure('numbertitle','off','name','2','pos',[50200420420]);
bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('2阶主振型'),holdon,grid
pause(0.1)%暂停0.1秒
h1=figure('numbertitle','off','name','3','pos',[50200420420]);
bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('3阶主振型'),holdon,grid
2-14作业题2-14作业题
functioncdjz2%
clearall,closeall%清空当前所有数据,关闭当前所有的绘图窗口
J1=1/2;J2=1;%%定义各个转动惯量J,I=1
k1=100000;k2=100000;%定义各个刚度的具体数值
fid=fopen('chuandi','wt')%建立打开速度文件
M1L=0%定义变量M1L的初始值为0,M1L为第一个质量块左边的转矩大小
forWN=0:
.01:
2000%for循环
shita1R=1;M1R=-WN^2*J1;%定义θ1=1,同时利用公式2-73求出质量块右边的转矩大小M1R
shita2R=shita1R+1/k1*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k1)*M1R;%求θ2,M2R,应用到2-78公式,问题是这个公式是θL的
shita3R=shita2R+1/k2*M2R;%依次类推
ifabs(shita3R)<0.005%判断精度是否满足,满足后执行
WN%搜索到的固有频率(rad/s),并显示
shita=[shita1R;shita2R;shita3R]%搜索到振型,并显示
bar(shita),xlabel('对应的质量块'),ylabel('振型向量')%画θ矩阵图形,条形图
pause(0.1)%1秒一暂停
end
fprintf(fid,'%30.15f',shita3R);%储存shita3R
end
fid=fopen('chuandi','rt');%打开可读写chuandi文件
x=fscanf(fid,'%f',[1,200001]);%拾取数据
t=1:
200001;
plot(.01*t,x),grid,xlabel('频率rad/s'),ylabel('第三个质量块的转角(rad)'),title('用传递矩阵法求固有频率')
第三章
一、欧拉法(P52业,例3-1)
functionvtb3(m,c,k,x0,v0,tf,w,f0,delt)%请运行vtb3(18.2,1.49,43.8,1,1,100,15,44.5,1)
%用欧拉法计算单自由度系统谐迫振动响应
wn=sqrt(k/m);%计算固有频率wn
fid1=fopen('disp','wt');%建立一个位移文件disp.dat
fort=0:
delt:
tf%delt为时间步长
xdd=(f0*sin(w*t)-k*x0-c*v0)/m;%计算加速度,代入公式3-1第三条公式,自下而上
xd=v0+xdd*delt;%计算速度,
x=x0+xd*delt;%计算位移x
fprintf(fid1,'%10.4f',x0);%向文件中写入数据
x0=x;v0=xd;
t
end
fid2=fopen('disp','rt');%打开disp.dat文件
n=tf/delt;%disp.dat文件中位移的个数
x=fscanf(fid2,'%f',[1,n]);%将disp.dat文件中位移写成矩阵
t=1:
n;
plot(t,x),grid
xlabel('时间(s)')
ylabel('位移(m)')
title('位移与时间的关系')
二、欧拉法的改进(P52页)
functionvtb4(m,c,k,x0,v0,tf,w,f0,delt)%请运行vtb4(282,249,43.8,1,1,100,15,4.5,1)
%用改进的欧拉法计算单自由度系统谐迫振动响应
wn=sqrt(k/m);%计算固有频率wn
fid1=fopen('disp','wt');%建立一个位移文件disp.dat
fort=0:
delt:
tf%delt为时