动力学程序讲解Word格式.docx
《动力学程序讲解Word格式.docx》由会员分享，可在线阅读，更多相关《动力学程序讲解Word格式.docx（21页珍藏版）》请在冰豆网上搜索。
%输出A
phi=%.3g\n'
phi）;
%输出ψ
elseifz==1%临界阻尼情况，按照1-14到1-15公式计算
a1=x0;
%C1
a2=v0+wn*x0;
%C2
a1=%.3g\n'
a1）;
%输出a1
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）;
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）
%单自由度系统的谐迫振动
lan=w/wn%λ的求法1-18公式
%输出Wn
%输出ζ
%输出Wd
a=sqrt（（（v0+z*wn*x0）^2+（x0*wd）^2）/wd^2）;
tf/1000:
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
位移（cm）'
位移与时间的关系'
第二章
一、矩阵迭代法（P38业，例2-3）
functionjzdd%矩阵迭代法
clearall%清空当前所有数据
closeall%关闭当前所有的绘图窗口
fid1=fopen（'
A-1'
'
wt'
）;
%建立第一个名为“A-1”的可写文档
fid2=fopen（'
B-1'
%同上
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值
%真正地迭代代码
end%结束while循环
f=sqrt（pp）/2/pi%公式2-61
fprintf（fid1,'
%20.5f'
%fid为文件句柄，指定要写入数据的文件，format是用来控制所写数据格式的格式符，与fscanf函数相同，A是用来存放数据的矩阵。
fprintf（fid2,'
f）;
%储存f
D=D-A*A'
*M/（A'
*M*A*pp）;
%下一阶的动力矩阵【D*】，A’代表A的逆矩阵
end%结束for循环
rt'
%读取A-1文档
A=fscanf（fid1,'
%f'
[3,3]）;
%输入3x3的A矩阵，来源为fid1
%读取B-1文档
f=fscanf（fid2,'
[3,1]）;
%输入（拾取）3x1的f矩阵,来源是fid2
t=1:
3;
h1=figure（'
numbertitle'
off'
name'
0'
pos'
[50200420420]）;
bar（t,f（t,1））,xlabel（'
频率阶级次'
）,ylabel（'
Hz'
）,
固有频率'
）,holdon,grid%bar代表绘制长条图t代表起点横坐标，f（t,1）纵坐标,输出各阶频率阶次的的固有频率
1'
bar（t,A（t,1））,xlabel（'
自由度（质量块）'
振型向量'
）,%输出对应的X阶主振型
1阶主振型'
）,holdon,grid
pause（0.1）%暂停0.1秒
2'
bar（t,A（t,2））,xlabel（'
2阶主振型'
3'
bar（t,A（t,3））,xlabel（'
3阶主振型'
二、传递矩阵法（P45页，例2-5）
functioncdjz
J1=1;
J2=1;
J3=2;
%定义各个转动惯量J
k2=1100000;
k3=1200000;
k4=100000;
%定义各个刚度的具体数值
fid=fopen（'
chuandi'
%建立打开速度文件
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（'
对应的质量块'
）%画θ矩阵图形
pause（1.0）%1秒一暂停
end%结束if语句
fprintf（fid,'
%30.15f'
shita4R）;
%将θ4输到文档中，数据类型是30位数，15位小数的浮点数
end%结束for循环
%此时fid有返回值，当是正数时代表打开文件成功，-1代表失败，rt表示读写
x=fscanf（fid,'
[1,200001]）;
%要求输入1到2000001的数
2000;
plot（t,x）,grid,xlabel（'
频率rad/s'
第四个质量块的转角（rad）'
）,title（'
用传递矩阵法求固有频率'
作业（P48-P49）
2-8作业题
functionjzdd2%第2-8作业题
A-2'
B-2'
M（1,1）=1;
M（2,2）=1;
%以上三段代码是为了输入质量矩阵，依题意均为1
K（1,1）=2;
K（1,2）=-1;
K（2,1）=-1;
K（2,2）=2;
%输入刚度矩阵，依题意均为1
%定义一个初始迭代阵型，ones（）函数表示是3个位为1的单列矩阵，ones（i,j）
%则是i行两列都是j都是1的数组！
3%进入for循环，将要执行书本上37业的迭代，循环次数为1到i，每执行一次循环代表一个质量块的主阵型
%同上，B（3）,是你当时定义主阵型时，第i个元素全为1，当然也可以是第1个
f=sqrt（pp）%公式2-61
%下一阶的动力矩阵【D*】，A’代表A的逆矩阵
end%结束for循环
%读取A-2文档
%输入3x3的A矩阵，来源为fid1
%读取B-2文档
）,%输出对应的X阶主振型
），
2-9作业题
functionjzdd3
A-3'
%建立第一个名为“A-3”的可写文档
B-3'
M（1,1）=4;
M（2,2）=2;
K（1,1）=4;
%输入刚度矩阵
pp=1.0/B
（1）;
%同上，B
（1）,是你当时定义主阵型时，第i个元素全为1，当然也可以是第一个
A=B/B
（1）;
δ，满足后退出while循环
end%结束while循环
%储存f
%读取A-3文档
2-14作业题2-14作业题
functioncdjz2%
clearall,closeall%清空当前所有数据，关闭当前所有的绘图窗口
J1=1/2;
%%定义各个转动惯量J，I=1
k1=100000;
k2=100000;
%定义各个刚度的具体数值
）%建立打开速度文件
M1L=0%定义变量M1L的初始值为0，M1L为第一个质量块左边的转矩大小
%定义θ1=1，同时利用公式2-73求出质量块右边的转矩大小M1R
shita2R=shita1R+1/k1*M1R;
M2R=shita1R*（-WN^2*J2）+（1+（-WN^2*J2）/k1）*M1R;
ifabs（shita3R）<
0.005%判断精度是否满足，满足后执行
WN%搜索到的固有频率（rad/s），并显示
shita3R]%搜索到振型，并显示
）%画θ矩阵图形，条形图
pause（0.1）%1秒一暂停
end
fprintf（fid,'
shita3R）;
%储存shita3R
%打开可读写chuandi文件
%拾取数据
200001;
plot（.01*t,x）,grid,xlabel（'
第三个质量块的转角（rad）'
第三章
一、欧拉法（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'
%建立一个位移文件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
fid2=fopen（'
%打开disp.dat文件
n=tf/delt;
%disp.dat文件中位移的个数
x=fscanf（fid2,'
[1,n]）;
%将disp.dat文件中位移写成矩阵
t=1:
n;
plot（t,x）,grid
xlabel（'
ylabel（'
title（'
二、欧拉法的改进（P52页）
functionvtb4（m,c,k,x0,v0,tf,w,f0,delt）%请运行vtb4（282,249,43.8,1,1,100,15,4.5,1）
%用改进的欧拉法计算单自由度系统谐迫振动响应
tf%delt
升级会员 