动力学程序讲解Word格式.docx
《动力学程序讲解Word格式.docx》由会员分享,可在线阅读,更多相关《动力学程序讲解Word格式.docx(21页珍藏版)》请在冰豆网上搜索。
![动力学程序讲解Word格式.docx](https://file1.bdocx.com/fileroot1/2022-11/23/68ea95f7-f60e-46e0-ad94-63e686d80096/68ea95f7-f60e-46e0-ad94-63e686d800961.gif)
%输出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,并且输出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值
%真正地迭代代码
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