动力学程序.docx

1、动力学程序动力学程序（附程序的详细翻译）第一章第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/2mwd=wn*sqrt(1-z2);fprintf(固有频率为%.3g.rad/s.n,wn);%输出wnfprintf(阻尼系数为%.3g.n,z);%输出fprintf(有阻尼的固有频率为%.3g.rad/s.n,wd);%输出wdt=0:tf/10000:tf;%时域的长短，决定很坐标的长短i

2、f z1 %判断的取值，如果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-61fprintf(fid1,%20.5f,A); %fid为文件句柄，指定要写入数据的文件，format是用来控制所写数据格式的格式符，与fscanf函数相同，A是用来存放数据的矩阵。fprintf(fid2,%20.5f,f); %储存fD=D-A*A*M/(A*M*A*pp

3、); %下一阶的动力矩阵【D*】，A代表A的逆矩阵end%结束for循环fid1=fopen(A-1,rt);%读取A-1文档A=fscanf(fid1,%f,3,3);%输入3x3的A矩阵，来源为fid1fid2=fopen(B-1,rt);%读取B-1文档f=fscanf(fid2,%f,3,1);%输入（拾取）3x1的f矩阵,来源是fid2t=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

4、代表绘制长条图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,gridpause(0.1)%暂停0.1秒h1=figure(numbertitle,off,name,2,pos,50 200 420 420);bar(t,A(t,2),xlabel(自由度（质量块）),ylabel(振型向量),title(2阶主振型),ho

5、ld on,gridpause(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二、传递矩阵法（P45页，例2-5）function cdjzclear all %清空当前所有数据close all %关闭当前所有的绘图窗口J1=1;J2=1;J3=2; %定义各个转动惯量Jk2=1100000;k3=1200000;k4=100000; %定义各个刚度的具体数值fid=fopen

6、(chuandi,wt); %建立打开速度文件M1L=0; %定义变量M1L的初始值为0，M1L为第一个质量块左边的转矩大小for WN=0:.01:2000 %for 循环 shita1R=1;M1R=-WN2*J1; %定义1，同时利用公式2-73求出质量块右边的转矩大小M1R shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN2*J2)+(1+(-WN2*J2)/k2)*M1R;%求2，M2R，应用到2-78公式，问题是这个公式是L的 shita3R=shita2R+1/k2*M2R;M3R=shita2R*(-WN2*J3)+(1+(-WN2*J3)/k

7、3)*M2R; %依次类推 shita4R=shita3R+1/k4*M3R; %（待分析，不知神马东东） if abs(shita4R)1e-6 % 判断是否满足迭代的条件|P（k）2-P(k-1)2|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-61fprintf(fid1,%20.5f,A); %fid为文件句柄，指定要写入数据的文件，format是用来控

8、制所写数据格式的格式符，与fscanf函数相同，A是用来存放数据的矩阵。fprintf(fid2,%20.5f,f);%储存fD=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矩阵，来源为fid1fid2=fopen(B-3,rt);%读取B-1文档f=fscanf(fid2,%f,3,1);%输入（拾取）3x1的f矩阵,来源是fid2t=1:3;h1=figure(numbertitle,off,name,0,pos,

9、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,gridpause(0.1)%暂停0.1秒h1=figure(numbertitle,off,na

10、me,2,pos,50 200 420 420);bar(t,A(t,2),xlabel(自由度（质量块）),ylabel(振型向量),title(2阶主振型),hold on,gridpause(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,grid2-14作业题2-14作业题function cdjz2 %clear all, close all%清空当前所有数据，关闭当前所有的绘图

11、窗口J1=1/2;J2=1;%定义各个转动惯量J，I=1k1=100000;k2=100000;%定义各个刚度的具体数值fid=fopen(chuandi,wt) %建立打开速度文件M1L=0 %定义变量M1L的初始值为0，M1L为第一个质量块左边的转矩大小for WN=0:.01:2000 %for 循环 shita1R=1;M1R=-WN2*J1;%定义1=1，同时利用公式2-73求出质量块右边的转矩大小M1R shita2R=shita1R+1/k1*M1R;M2R=shita1R*(-WN2*J2)+(1+(-WN2*J2)/k1)*M1R;%求2，M2R，应用到2-78公式，问题是这

12、个公式是L的 shita3R=shita2R+1/k2*M2R; %依次类推 if abs(shita3R)0.005 %判断精度是否满足，满足后执行 WN %搜索到的固有频率（rad/s），并显示 shita=shita1R;shita2R;shita3R %搜索到振型，并显示 bar(shita),xlabel(对应的质量块),ylabel(振型向量)%画矩阵图形，条形图 pause(0.1) %1秒一暂停 endfprintf(fid,%30.15f,shita3R);%储存shita3R endfid=fopen(chuandi,rt); %打开可读写chuandi文件x=fscanf

13、(fid,%f,1,200001); %拾取数据t=1:200001;plot(.01*t,x),grid,xlabel(频率rad/s),ylabel(第三个质量块的转角(rad),title(用传递矩阵法求固有频率)第三章一、欧拉法（P52业，例3-1）function vtb3(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 for t

14、=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页）function vtb4(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 for t=0:delt:tf %delt为时