动力学程序.docx
《动力学程序.docx》由会员分享，可在线阅读，更多相关《动力学程序.docx（21页珍藏版）》请在冰豆网上搜索。
动力学程序
动力学程序（附程序的详细翻译）
第一章
第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，并且输出Ｐ（０）
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，并且输出Ｐ（０）
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，并且输出Ｐ（０）
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为时
升级会员 