动力学第三章
- 格式:doc
- 大小:53.00 KB
- 文档页数:12
第三章化学动力学(总15页)--本页仅作为文档封面,使用时请直接删除即可----内页可以根据需求调整合适字体及大小--第三章 化学动力学3-1.在1 100 K 时,3NH (g)在金属钨丝上发生分解。
实验测定,在不同的3NH (g)的初始压力0p 下所对应的半衰期12t ,获得下列数据 0/Pa p ×104 ×104 ×104 12/min t 试用计算的方法,计算该反应的级数和速率系数。
解: 根据实验数据,反应物3NH (g)的初始压力不断下降,相应的半衰期也不断下降,说明半衰期与反应物的起始浓度(或压力)成正比,这是零级反应的特征,所以基本可以确定是零级反应。
用半衰期法来求反应的级数,根据半衰期法的计算公式12121,121,2n t a t a -⎛⎫= ⎪⎝⎭即 ()12,112,221ln /1ln(/)t t n a a =+把实验数据分别代入,计算得()()12,112,2440,20,1ln /ln 7.6/3.7110ln(/)ln(1.710/3.510)t t n p p --=+=+≈⨯⨯ 同理,用后面两个实验数据计算,得 ()ln 3.7/1.710ln(0.75/1.7)n =+≈所以,该反应为零级反应。
利用零级反应的积分式,计算速率系数。
正规的计算方法应该是分别用3组实验数据,计算得3个速率系数,然后取平均值。
这里只列出用第一组实验数据计算的结果,即010022p at k k == 4310012 3.510Pa 2.310 Pa min 227.6 minp k t -⨯===⨯⋅⨯3-2.某人工放射性元素,能放出α粒子,其半衰期为15 min 。
若该试样有80%被分解,计算所需的时间解:放射性元素的蜕变,符合一级反应的特征。
对于一级反应,已知半衰期的数值,就能得到速率系数的值,因为一级反应的半衰期是与反应物浓度无关的常数。
然后再根据一级反应的定积分式,计算分解80%所需的时间。
第三章充气轮胎动力学§3-1 概述轮胎是车辆重要的组成部分,直接与地面接触。
其作用是支承整车的重量,与悬架共同缓冲来自路面的不平度激励,以保证车辆具有良好的乘坐舒适性和行驶平顺性;保证车轮和路面具有良好的附着性,以提高车辆驱动性、制动性和通过性,并为车辆提供充分的转向力。
一、轮胎运动坐标系二、车轮运动参数1.滑动率2.轮胎侧偏角a3.轮胎径向变形§3-2 轮胎的功能、结构及发展轮胎的基本功能包括:1)支撑整车重量;2)与悬架元件共同作用,衰减由路面不平引起的振动与冲击;3)传递纵向力,以实现驱动和制动;4)传递侧向力,以使车辆转向并保证行驶稳定性。
为实现以上功能,任何一个充气轮胎都必须具备以下基本结构:(1)胎体(2)胎圈(3)胎面常用的车用充气轮胎有两种,即斜交轮胎和子午线轮胎。
二者在结构上有明显不同,主要区别在于胎体帘线角度的不同。
所谓“帘线角”即为胎体帘布层单线与车轮中心线形成的夹角。
根据车辆动力学研究内容的不同,轮胎模型可分为:(1)轮胎纵滑模型主要用于预测车辆在驱动和制动工况时的纵向力。
(2)轮胎侧偏模型和侧倾模型主要用于预测轮胎的侧向力和回正力矩,评价转向工况下低频转角输入响应。
(3)轮胎垂向振动模型主要用于高频垂向振动的评价,并考虑轮胎的包容特性(包含刚性滤波和弹性滤波特性)。
这里仅对几种常用的轮胎模型给予介绍。
(1)幂指数统一轮胎模型幂指数统一轮胎模型的特点是:。
1)采用了无量纲表达式,其优点在于由纯工况下的一次台架试验得到的试验数据可应用于各种不同的路面。
当路面条件改变时,只要改变路面的附着特性参数,代人无量纲表达式即可得该路面下的轮胎特性。
2)无论是纯工况还是联合工况,其表达式是统一的。
3)可表达各种垂向载荷下的轮胎特性。
4)保证了可用较少的模型参数实现全域范围内的计算精度,参数拟合方便,计算量小。
在联合工况下,其优势更加明显。
5)能拟合原点刚度。
(2)“魔术公式”轮胎模型“魔术公式”轮胎模型的特点是:1)用一套公式可以表达出轮胎的各向力学特性,统一性强,编程方便,需拟合参数较少,且各个参数都有明确的物理意义,容易确定其初值。
第2章function VTB2(m,c,k,x0,v0,tf,w,f0)%单自由度系统的谐迫振动clcwn=sqrt(k/m);z=c/2/m/wn;lan=w/wnwd=wn*sqrt(1-z^2);A=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);t=0:tf/1000:tf;phi1=atan2(x0*wd,v0+z*wn*x0);phi2=atan2(2*z*lan,1-lan^2);B=wn^2*f0/k/sqrt((wn^2-w^2)^2+(2*z*wn*w)^2);x=A*exp(-z*wn*t).*sin(sqrt(1-z^2)*wn*t+phi1)+B*sin(w*t-phi2);plot(t,x),gridxlabel('时间(s)')ylabel('位移')title('位移与时间的关系')function VTB1(m,c,k,x0,v0,tf)%VTB1用来计算单自由度有阻尼自由振动系统的响应%VTB1绘出单自由度有阻尼自由振动系统的响应图%m为质量;c为阻尼;k为刚度;x0为初始位移;v0为初始速度;tf为仿真时间%VTB1(zeta,w,x0,v0,tf)绘出单自由度有阻尼自由振动系统的响应图%程序中z为阻尼系数ξ;wn为固有频率ωn;A为振动幅度;phi为初相位θclcwn=sqrt(k/m);z=c/2/m/wn;wd=wn*sqrt(1-z^2);fprintf('固有频率为%.s.\n',wd);fprintf('阻尼系数为%.3g.\n',z);fprintf('有阻尼的固有频率为%.s.\n',wd);t=0:tf/1000:tf;if z<1A=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);elseif z==1a1=x0;a2=v0+wn*x0;fprintf('a1=%.3g\n',a1);fprintf('a2=%.3g\n',a2);x=(a1+a2*t).*exp(-wn*t);elsea1=(-v0+(-z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);a2=(v0+(z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);fprintf('a1=%.3g\n',a1);fprintf('a2=%.3g\n',a2);x=exp(-z*wn*t).*(a1*exp(-wn*sqrt(z^2-1)*t)+a2*exp(wn*sqrt(z^2-1)*t)); endplot(t,x),gridxlabel('时间(s)')ylabel('位移')title('位移与时间的关系')function jzdd%矩阵迭代法求系统的三阶固有频率和主阵型clear allclose allfid1=fopen('A-1','wt'); %建立主振型文件fid2=fopen('B-1','wt'); %建立固有频率文件%输入质量矩阵M(1,1)=2;M(2,2)=;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; %原始动力矩阵A=ones(3,1); %初始振型for i=1:3 %计算三阶固有频率和主振型pp0=0;iB=D*A;pp=B(3); %B(3)为B中的最后一个元素A=B/B(3);while abs((pp-pp0)/pp)>pp0=pp;B=D*A;pp=B(3);A=B/B(3);endf=sqrt(pp)/2/pi %固有频率单位转换为Hzfprintf(fid1,'%',A); %输入主振型数据fprintf(fid2,'%',f); %输入固有频率数据D=D-A*A'*M/(A'*M*A*pp);endfid1=fopen('A-1','rt'); %打开主振型文件A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵fid2=fopen('B-1','rt'); %打开固有频率文件f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵t=1:3;h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]); bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),title('固有频率'),hold on,gridh1=figure('numbertitle','off','name','1','pos',[50 200 420 420]); bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'), title('1阶主振型'),hold on,gridpauseh1=figure('numbertitle','off','name','2','pos',[50 200 420 420]); bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'), title('2阶主振型'),hold on,gridpauseh1=figure('numbertitle','off','name','3','pos',[50 200 420 420]); bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'), title('3阶主振型'),hold on,gridpauseend%; %传递矩阵方法求固有频率clear all,clear closeJ1=1;J2=1;J3=2;k2=1100000;k3=1200000;k4=100000;fid=fopen('chuandi','wt'); %建立(打开)速度文件M1L=0;for WN=0::2000shita1R=1;M1R=-WN^2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;shita3R=shita2R+1/k3*M2R;M3R=shita2R*(-WN^2*J3)+(1+(-WN^2*J3)/k3)*M2R;shita4R=shita3R+1/k4*M3R;if abs(shita4R)<WN %搜索到的固有频率(rad/s),并显示shita=[shita1R;shita2R;shita3R;shita4R];%搜索到振型,并显示bar(shita),xlabel('对应的质量块'),ylabel('振型向量')pauseendfprintf(fid,'%',shita4R);endfid=fopen('chuandi','rt');x=fscanf(fid,'%f',[1,200001]);t=1:200001;plot*t,x);grid,xlabel('频率rad/s'),ylabel('第四个质量块的转角(rad/s)'),title('用传递矩阵法求固有频率')function cdjz2%; %传递矩阵方法求固有频率clear all,clear closeJ1=;J2=1;k2=100000;k3=100000;fid=fopen('chuandi3','wt'); %建立(打开)速度文件M1L=0;for WN=0::2000shita1R=1;M1R=-WN^2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;shita3R=shita2R+1/k3*M2R;if abs(shita3R)<WN %搜索到的固有频率(rad/s),并显示shita=[shita1R;shita2R;shita3R;] %搜索到振型,并显示bar(shita),xlabel('对应的质量块'),ylabel('振型向量')pauseendfprintf(fid,'%',shita3R);endfid=fopen('chuandi3','rt');x=fscanf(fid,'%f',[1,200001]);t=1:200001;plot*t,x);grid,xlabel('频率rad/s'),ylabel('第四个质量块的转角(rad/s)'),title('用传递矩阵法求固有频率')function zuoye8%矩阵迭代法求系统的三阶固有频率和主阵型clear allclose allfid1=fopen('A-2','wt'); %建立主振型文件fid2=fopen('B-2','wt'); %建立固有频率文件%输入质量矩阵M(1,1)=1;M(2,2)=1;M(3,3)=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%计算特征值与特征向量D=inv(K)*M; %原始动力矩阵A=ones(3,1); %初始振型for i=1:3 %计算三阶固有频率和主振型pp0=0;iB=D*A;pp=B(3); %B(3)为B中的最后一个元素A=B/B(3);while abs((pp-pp0)/pp)>pp0=pp;B=D*A;pp=B(3);A=B/B(3);endf=sqrt(pp)fprintf(fid1,'%',A); %输入主振型数据fprintf(fid2,'%',f); %输入固有频率数据D=D-A*A'*M/(A'*M*A*pp);endfid1=fopen('A-2','rt'); %打开主振型文件A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵fid2=fopen('B-2','rt'); %打开固有频率文件f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵t=1:3;h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),title('固有频率'),hold on,gridh1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('1阶主振型'),hold on,gridpauseh1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('2阶主振型'),hold on,gridpauseh1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('3阶主振型'),hold on,gridpauseendfunction zuoye9%矩阵迭代法求系统的三阶固有频率和主阵型clear allclose allfid1=fopen('A-3','wt'); %建立主振型文件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; %原始动力矩阵U=inv(K)A=ones(3,1); %初始振型for i=1:3 %计算三阶固有频率和主振型pp0=0;iB=D*A;pp=B(1); %B(1)为B中的第一个元素A=B/B(1);while abs((pp-pp0)/pp)>pp0=pp;B=D*A;pp=B(1);A=B/B(1);endf=sqrt(pp)fprintf(fid1,'%',A); %输入主振型数据fprintf(fid2,'%',f); %输入固有频率数据D=D-A*A'*M/(A'*M*A*pp);endfid1=fopen('A-3','rt'); %打开主振型文件A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵fid2=fopen('B-3','rt'); %打开固有频率文件f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵t=1:3;h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]); bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),title('固有频率'),hold on,gridh1=figure('numbertitle','off','name','1','pos',[50 200 420 420]); bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'), title('1阶主振型'),hold on,gridpauseh1=figure('numbertitle','off','name','2','pos',[50 200 420 420]); bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'), title('2阶主振型'),hold on,gridpauseh1=figure('numbertitle','off','name','3','pos',[50 200 420 420]); bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'), title('3阶主振型'),hold on,gridpauseendfunction cdjz2%; %传递矩阵方法求固有频率clear all,clear closeJ1=;J2=1;k2=10000;k3=10000;fid=fopen('chuandi3','wt'); %建立(打开)速度文件M1L=0;for WN=0::2000shita1R=1;M1R=-WN^2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;shita3R=shita2R+1/k3*M2R;if abs(shita3R)<WN %搜索到的固有频率(rad/s),并显示shita=[shita1R;shita2R;shita3R;] %搜索到振型,并显示bar(shita),xlabel('对应的质量块'),ylabel('振型向量')pauseendfprintf(fid,'%',shita3R);endfid=fopen('chuandi3','rt');x=fscanf(fid,'%f',[1,200001]);t=1:200001;plot*t,x);grid,xlabel('频率rad/s'),ylabel('第三个质量块的转角(rad/s)'),title('用传递矩阵法求固有频率')第3章function vtb3(m,c,k,x0,v0,tf,w,f0,delt)%用欧拉法计算单自由度系统谐迫振动响应wn=sqrt(k/m); %计算固有频率fid1=fopen('disp','wt'); %建立一个位移文件for t=0:delt:tf; %delt为时间步长xdd=(f0*sin(w*t)-k*x0-c*v0)/m; %计算加速度xd=v0+xdd*delt; %计算速度x=x0+xd*delt; %计算位移xfprintf(fid1,'%',x0); %向文件中写数据x0=x;v0=xd;tendfid2=fopen('disp','rt'); %打开文件n=tf/delt; %文件中位移的个数x=fscanf(fid2,'%f',[1,n]); %将文件中文艺写成矩阵t=1:n;plot(t,x),gridxlabel('时间(s)')ylabel('位移(s)')title('位移与时间的关系')function vtb4(m,c,k,x0,v0,tf,w,f0,delt)%用改进的欧拉法计算单自由度系统谐迫振动响应wn=sqrt(k/m); %计算固有频率fid1=fopen('disp','wt'); %建立一个位移文件for t=0:delt:tf; %delt为时间步长xdd=(f0*sin(w*t)-k*x0-c*v0)/m; %计算加速度x3d=(f0*w*cos(w*t)-k*v0-c*xdd)/m;xd=v0+xdd*delt+x3d*delt^2/2; %计算速度x=x0+xd*delt+xdd*delt^2/2; %计算位移xfprintf(fid1,'%',x0); %向文件中写数据x0=x;v0=xd;tendfid2=fopen('disp','rt'); %打开文件n=tf/delt; %文件中位移的个数x=fscanf(fid2,'%f',[1,n]); %将文件中文艺写成矩阵t=1:n;plot(t,x),gridxlabel('时间(s)')ylabel('位移(s)')title('位移与时间的关系')function vtb5(tf,delt) %用线性加速度法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen('disp5','wt'); %建立一个位移文件m=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵c=*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵x0=[1 1 1]'; %初始位移v0=[1 1 1]'; %初始速度md=inv(m+delt/2*c+1/6*delt^2*k);m1=inv(m);[E,F]=eig(m1*k);diag(sqrt(F)); %计算固有频率(rad/s)for t=0:delt:tf; %delt为时间步长f=[*sin*t) *cos*t) *sin*t)]';if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度elsex=md*(m*(x0+delt*v0+delt^2/3*xdd0)+c*(delt/2*x0+delt^2/3*v0+delt^3/12*xdd0)+delt^2/6*f);%计算位移xdd=6/delt^2*(x-(x0+delt*v0+delt^2/3*xdd0)); %计算加速度xd=v0+delt/2*(xdd0+xdd);%计算速度xdd0=xdd;v0=xd;x0=x;fprintf(fid1,'%',x0);%向文件中写数据t %显示计算时间步长endendfid2=fopen('disp5','rt'); %打开文件n=tf/delt; %文件中位移的个数x=fscanf(fid2,'%f',[3,n]); %将文件中的位移写成矩阵t=1:n;figure('numbertitle','off','name','自由度-1的位移','pos',[450 180 400 420]);plot(t,x(1,t)),grid,xlabel('时间*秒'),title('自由度-1的位移与时间的关系')figure('numbertitle','off','name','自由度-2的位移','pos',[350 160 400 420]);plot(t,x(2,t)),grid,xlabel('时间*秒'),title('自由度-3的位移与时间的关系')figure('numbertitle','off','name','自由度-3的位移','pos',[250 140 400 420]);plot(t,x(3,t)),grid,xlabel('时间*秒'),title('自由度-3的位移与时间的关系')function vtb6(tf,delt) %用线纽马克-β法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen('disp6','wt'); %建立一个位移文件m=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵c=*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵x0=[1 1 1]'; %初始位移v0=[1 1 1]'; %初始速度bita=1/6;md=inv(m+delt/2*c+bita*delt^2*k);m1=inv(m);[E,F]=eig(m1*k);diag(sqrt(F)); %计算固有频率(rad/s)for t=0:delt:tf; %delt为时间步长f=[*sin*t) *cos*t) *sin*t)]';if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度elsexdd=md*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-bita)*delt^2*xdd0)); %计算加速度xd=v0+delt/2*(xdd0+xdd);%计算速度x=x0+delt*v0+delt^2/2*xdd0+bita*delt^3*(xdd-xdd0)/delt; %计算位移v0=xd;x0=x;fprintf(fid1,'%',x0);%向文件中写数据t %显示计算时间步长endendfid2=fopen('disp6','rt'); %打开文件n=tf/delt; %文件中位移的个数x=fscanf(fid2,'%f',[3,n]); %将文件中的位移写成矩阵t=1:n;figure('numbertitle','off','name','自由度-1的位移','pos',[450 180 400 420]);plot(t,x(1,t)),grid,xlabel('时间*秒'),title('自由度-1的位移与时间的关系')figure('numbertitle','off','name','自由度-2的位移','pos',[350 160 400 420]);plot(t,x(2,t)),grid,xlabel('时间*秒'),title('自由度-2的位移与时间的关系')figure('numbertitle','off','name','自由度-3的位移','pos',[250 140 400 420]);plot(t,x(3,t)),grid,xlabel('时间*秒'),title('自由度-3的位移与时间的关系')function vtb7(tf,delt) %用威尔逊θ法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen('disp7','wt'); %建立一个位移文件m=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵c=*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵x0=[1 1 1]'; %初始位移v0=[1 1 1]'; %初始速度theta=;md=inv(k+3*c/theta/delt+6/(theta*delt^2)*m);m1=inv(m);[E,F]=eig(m1*k);diag(sqrt(F)); %计算固有频率(rad/s)for t=0:delt:tf; %delt为时间步长f=[*sin*t) *cos*t) *sin*t)]';if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度elsextheta=md*(m*(2*xdd0+6/theta/delt*v0+6/(theta*delt)^2*x0)+c*(theta*delt/2*xdd0+2*v0+3/theta/delt* x0)+f); %计算(t+θdelt)时刻的速度xddtheta=6/(theta*delt)^2*(xtheta-x0)-6/theta/delt*v0-2*xdd0; %计算(t+θdelt)时刻的加速度xdd=(1-1/theta)*xdd0+1/theta*xddtheta; %计算(t+delt)时刻的加速度xd=v0+delt/2*(xdd0+xdd); %计算(t+delt)速度x=x0+delt*v0+delt^2*(2*xdd0+xdd)/6; %计算(t+delt)位移v0=xd;x0=x;xdd0=xdd;fprintf(fid1,'%',x0);%向文件中写数据t %显示计算时间步长endendfid2=fopen('disp7','rt'); %打开文件n=tf/delt; %文件中位移的个数x=fscanf(fid2,'%f',[3,n]); %将文件中的位移写成矩阵t=1:n;figure('numbertitle','off','name','自由度-1的位移','pos',[450 180 400 420]);plot(t,x(1,t)),grid,xlabel('时间*秒'),title('自由度-1的位移与时间的关系')figure('numbertitle','off','name','自由度-2的位移','pos',[350 160 400 420]);plot(t,x(2,t)),grid,xlabel('时间*秒'),title('自由度-2的位移与时间的关系')figure('numbertitle','off','name','自由度-3的位移','pos',[250 140 400 420]);plot(t,x(3,t)),grid,xlabel('时间*秒'),title('自由度-3的位移与时间的关系')。
第2章function VTB2(m,c,k,x0,v0,tf,w,f0)%单自由度系统得谐迫振动clcwn=sqrt(k/m);z=c/2/m/wn;lan=w/wnwd=wn*sqrt(1-z^2);A=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);t=0:tf/1000:tf;phi1=atan2(x0*wd,v0+z*wn*x0);phi2=atan2(2*z*lan,1-lan^2);B=wn^2*f0/k/sqrt((wn^2-w^2)^2+(2*z*wn*w)^2);x=A*exp(-z*wn*t)、*sin(sqrt(1-z^2)*wn*t+phi1)+B*sin(w*t-phi2); plot(t,x),gridxlabel('时间(s)')ylabel('位移')title('位移与时间得关系')function VTB1(m,c,k,x0,v0,tf)%VTB1用来计算单自由度有阻尼自由振动系统得响应%VTB1绘出单自由度有阻尼自由振动系统得响应图%m为质量;c为阻尼;k为刚度;x0为初始位移;v0为初始速度;tf为仿真时间%VTB1(zeta,w,x0,v0,tf)绘出单自由度有阻尼自由振动系统得响应图%程序中z为阻尼系数ξ;wn为固有频率ωn;A为振动幅度;phi为初相位θclcwn=sqrt(k/m);z=c/2/m/wn;wd=wn*sqrt(1-z^2);fprintf('固有频率为%、3g、rad/s、\n',wd);fprintf('阻尼系数为%、3g、\n',z);fprintf('有阻尼得固有频率为%、3g、rad/s、\n',wd);t=0:tf/1000:tf;if z<1A=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);elseif z==1a1=x0;a2=v0+wn*x0;fprintf('a1=%、3g\n',a1);fprintf('a2=%、3g\n',a2);x=(a1+a2*t)、*exp(-wn*t);elsea1=(-v0+(-z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);a2=(v0+(z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);fprintf('a1=%、3g\n',a1);fprintf('a2=%、3g\n',a2);x=exp(-z*wn*t)、*(a1*exp(-wn*sqrt(z^2-1)*t)+a2*exp(wn*sqrt(z^2-1)*t)); endplot(t,x),gridxlabel('时间(s)')ylabel('位移')title('位移与时间得关系')function jzdd%矩阵迭代法求系统得三阶固有频率与主阵型clear allclose allfid1=fopen('A-1','wt'); %建立主振型文件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; %原始动力矩阵A=ones(3,1); %初始振型for i=1:3 %计算三阶固有频率与主振型pp0=0;iB=D*A;pp=1、0/B(3); %B(3)为B中得最后一个元素A=B/B(3);while abs((pp-pp0)/pp)>1、e-6pp0=pp;B=D*A;pp=1、0/B(3);A=B/B(3);endf=sqrt(pp)/2/pi %固有频率单位转换为Hzfprintf(fid1,'%20、5f',A); %输入主振型数据fprintf(fid2,'%20、5f',f); %输入固有频率数据D=D-A*A'*M/(A'*M*A*pp);endfid1=fopen('A-1','rt'); %打开主振型文件A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵fid2=fopen('B-1','rt'); %打开固有频率文件f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵t=1:3;h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),title('固有频率'),hold on,gridh1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('1阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('2阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('3阶主振型'),hold on,gridpause(0、1)end%chuandijuzhen、m; %传递矩阵方法求固有频率clear all,clear closeJ1=1;J2=1;J3=2;k2=1100000;k3=1200000;k4=100000;fid=fopen('chuandi','wt'); %建立(打开)速度文件M1L=0;for WN=0:0、01:2000shita1R=1;M1R=-WN^2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;shita3R=shita2R+1/k3*M2R;M3R=shita2R*(-WN^2*J3)+(1+(-WN^2*J3)/k3)*M2R; shita4R=shita3R+1/k4*M3R;if abs(shita4R)<0、005WN %搜索到得固有频率(rad/s),并显示shita=[shita1R;shita2R;shita3R;shita4R];%搜索到振型,并显示bar(shita),xlabel('对应得质量块'),ylabel('振型向量')pause(1、0)endfprintf(fid,'%30、15f',shita4R);endfid=fopen('chuandi','rt');x=fscanf(fid,'%f',[1,200001]);t=1:200001;plot(0、01*t,x);grid,xlabel('频率rad/s'),ylabel('第四个质量块得转角(rad/s)'),title('用传递矩阵法求固有频率')function cdjz2%chuandijuzhen、m; %传递矩阵方法求固有频率clear all,clear closeJ1=0、5;J2=1;k2=100000;k3=100000;fid=fopen('chuandi3','wt'); %建立(打开)速度文件M1L=0;for WN=0:0、01:2000shita1R=1;M1R=-WN^2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R; shita3R=shita2R+1/k3*M2R;if abs(shita3R)<0、005WN %搜索到得固有频率(rad/s),并显示shita=[shita1R;shita2R;shita3R;] %搜索到振型,并显示bar(shita),xlabel('对应得质量块'),ylabel('振型向量')pause(1、0)endfprintf(fid,'%30、15f',shita3R);endfid=fopen('chuandi3','rt');x=fscanf(fid,'%f',[1,200001]);t=1:200001;plot(0、01*t,x);grid,xlabel('频率rad/s'),ylabel('第四个质量块得转角(rad/s)'),title('用传递矩阵法求固有频率')function zuoye8%矩阵迭代法求系统得三阶固有频率与主阵型clear allclose allfid1=fopen('A-2','wt'); %建立主振型文件fid2=fopen('B-2','wt'); %建立固有频率文件%输入质量矩阵M(1,1)=1;M(2,2)=1;M(3,3)=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 %计算特征值与特征向量D=inv(K)*M; %原始动力矩阵A=ones(3,1); %初始振型for i=1:3 %计算三阶固有频率与主振型pp0=0;iB=D*A;pp=1、0/B(3); %B(3)为B中得最后一个元素A=B/B(3);while abs((pp-pp0)/pp)>1、e-6pp0=pp;B=D*A;pp=1、0/B(3);A=B/B(3);endf=sqrt(pp)fprintf(fid1,'%20、5f',A); %输入主振型数据fprintf(fid2,'%20、5f',f); %输入固有频率数据D=D-A*A'*M/(A'*M*A*pp);endfid1=fopen('A-2','rt'); %打开主振型文件A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵fid2=fopen('B-2','rt'); %打开固有频率文件f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵t=1:3;h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),title('固有频率'),hold on,gridh1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('1阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('2阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('3阶主振型'),hold on,gridpause(0、1)endfunction zuoye9%矩阵迭代法求系统得三阶固有频率与主阵型clear allclose allfid1=fopen('A-3','wt'); %建立主振型文件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; %原始动力矩阵U=inv(K)A=ones(3,1); %初始振型for i=1:3 %计算三阶固有频率与主振型pp0=0;iB=D*A;pp=1、0/B(1); %B(1)为B中得第一个元素A=B/B(1);while abs((pp-pp0)/pp)>1、e-6pp0=pp;B=D*A;pp=1、0/B(1);A=B/B(1);endf=sqrt(pp)fprintf(fid1,'%20、5f',A); %输入主振型数据fprintf(fid2,'%20、5f',f); %输入固有频率数据D=D-A*A'*M/(A'*M*A*pp);endfid1=fopen('A-3','rt'); %打开主振型文件A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵fid2=fopen('B-3','rt'); %打开固有频率文件f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵t=1:3;h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),title('固有频率'),hold on,gridh1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('1阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('2阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('3阶主振型'),hold on,gridpause(0、1)endfunction cdjz2%chuandijuzhen、m; %传递矩阵方法求固有频率clear all,clear closeJ1=0、5;J2=1;k2=10000;k3=10000;fid=fopen('chuandi3','wt'); %建立(打开)速度文件M1L=0;for WN=0:0、01:2000shita1R=1;M1R=-WN^2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R; shita3R=shita2R+1/k3*M2R;if abs(shita3R)<0、005WN %搜索到得固有频率(rad/s),并显示shita=[shita1R;shita2R;shita3R;] %搜索到振型,并显示bar(shita),xlabel('对应得质量块'),ylabel('振型向量')pause(1、0)endfprintf(fid,'%30、15f',shita3R);endfid=fopen('chuandi3','rt');x=fscanf(fid,'%f',[1,200001]);t=1:200001;plot(0、01*t,x);grid,xlabel('频率rad/s'),ylabel('第三个质量块得转角(rad/s)'),title('用传递矩阵法求固有频率')第3章function vtb3(m,c,k,x0,v0,tf,w,f0,delt)%用欧拉法计算单自由度系统谐迫振动响应wn=sqrt(k/m); %计算固有频率fid1=fopen('disp','wt'); %建立一个位移文件disp、datfor t=0:delt:tf; %delt为时间步长xdd=(f0*sin(w*t)-k*x0-c*v0)/m; %计算加速度xd=v0+xdd*delt; %计算速度x=x0+xd*delt; %计算位移xfprintf(fid1,'%10、4f',x0); %向文件中写数据x0=x;v0=xd;tendfid2=fopen('disp','rt'); %打开disp、dat文件n=tf/delt; %disp、dat文件中位移得个数x=fscanf(fid2,'%f',[1,n]); %将disp、dat文件中文艺写成矩阵t=1:n;plot(t,x),gridxlabel('时间(s)')ylabel('位移(s)')title('位移与时间得关系')function vtb4(m,c,k,x0,v0,tf,w,f0,delt)%用改进得欧拉法计算单自由度系统谐迫振动响应wn=sqrt(k/m); %计算固有频率fid1=fopen('disp','wt'); %建立一个位移文件disp、datfor t=0:delt:tf; %delt为时间步长xdd=(f0*sin(w*t)-k*x0-c*v0)/m; %计算加速度x3d=(f0*w*cos(w*t)-k*v0-c*xdd)/m;xd=v0+xdd*delt+x3d*delt^2/2; %计算速度x=x0+xd*delt+xdd*delt^2/2; %计算位移xfprintf(fid1,'%10、4f',x0); %向文件中写数据x0=x;v0=xd;tendfid2=fopen('disp','rt'); %打开disp、dat文件n=tf/delt; %disp、dat文件中位移得个数x=fscanf(fid2,'%f',[1,n]); %将disp、dat文件中文艺写成矩阵t=1:n;plot(t,x),gridxlabel('时间(s)')ylabel('位移(s)')title('位移与时间得关系')function vtb5(tf,delt) %用线性加速度法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen('disp5','wt'); %建立一个位移文件dip5、datm=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵c=1、5*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵x0=[1 1 1]'; %初始位移v0=[1 1 1]'; %初始速度md=inv(m+delt/2*c+1/6*delt^2*k);m1=inv(m);[E,F]=eig(m1*k);diag(sqrt(F)); %计算固有频率(rad/s)for t=0:delt:tf; %delt为时间步长f=[2、00*sin(3、754*t) -2、00*cos(2、2*t) 1、00*sin(2、8*t)]';if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度elsex=md*(m*(x0+delt*v0+delt^2/3*xdd0)+c*(delt/2*x0+delt^2/3*v0+delt^3/12*x dd0)+delt^2/6*f);%计算位移xdd=6/delt^2*(x-(x0+delt*v0+delt^2/3*xdd0)); %计算加速度xd=v0+delt/2*(xdd0+xdd);%计算速度xdd0=xdd;v0=xd;x0=x;fprintf(fid1,'%10、4f',x0);%向文件中写数据t %显示计算时间步长endendfid2=fopen('disp5','rt'); %打开disp5、dat文件n=tf/delt; %disp5、dat文件中位移得个数x=fscanf(fid2,'%f',[3,n]); %将disp5、dat文件中得位移写成矩阵t=1:n;figure('numbertitle','off','name','自由度-1得位移','pos',[450 180 400 420]);plot(t,x(1,t)),grid,xlabel('时间*0、1秒'),title('自由度-1得位移与时间得关系') figure('numbertitle','off','name','自由度-2得位移','pos',[350 160 400 420]);plot(t,x(2,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系') figure('numbertitle','off','name','自由度-3得位移','pos',[250 140 400 420]);plot(t,x(3,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系') function vtb6(tf,delt) %用线纽马克-β法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen('disp6','wt'); %建立一个位移文件dip6、datm=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵c=1、5*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵x0=[1 1 1]'; %初始位移v0=[1 1 1]'; %初始速度bita=1/6;md=inv(m+delt/2*c+bita*delt^2*k);m1=inv(m);[E,F]=eig(m1*k);diag(sqrt(F)); %计算固有频率(rad/s)for t=0:delt:tf; %delt为时间步长f=[2、00*sin(3、754*t) -2、00*cos(2、2*t) 1、00*sin(2、8*t)]';if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度elsexdd=md*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-bita)*delt^2*xdd0)); %计算加速度xd=v0+delt/2*(xdd0+xdd);%计算速度x=x0+delt*v0+delt^2/2*xdd0+bita*delt^3*(xdd-xdd0)/delt; %计算位移v0=xd;x0=x;fprintf(fid1,'%10、4f',x0);%向文件中写数据t %显示计算时间步长endendfid2=fopen('disp6','rt'); %打开disp6、dat文件n=tf/delt; %disp6、dat文件中位移得个数x=fscanf(fid2,'%f',[3,n]); %将disp6、dat文件中得位移写成矩阵t=1:n;figure('numbertitle','off','name','自由度-1得位移','pos',[450 180 400 420]);plot(t,x(1,t)),grid,xlabel('时间*0、1秒'),title('自由度-1得位移与时间得关系') figure('numbertitle','off','name','自由度-2得位移','pos',[350 160 400 420]);plot(t,x(2,t)),grid,xlabel('时间*0、1秒'),title('自由度-2得位移与时间得关系')420]);plot(t,x(3,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系') function vtb7(tf,delt) %用威尔逊θ法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen('disp7','wt'); %建立一个位移文件dip6、datm=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵c=1、5*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵x0=[1 1 1]'; %初始位移v0=[1 1 1]'; %初始速度theta=1、4;md=inv(k+3*c/theta/delt+6/(theta*delt^2)*m);m1=inv(m);[E,F]=eig(m1*k);diag(sqrt(F)); %计算固有频率(rad/s)for t=0:delt:tf; %delt为时间步长f=[2、00*sin(3、754*t) -2、00*cos(2、2*t) 1、00*sin(2、8*t)]';if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度elsextheta=md*(m*(2*xdd0+6/theta/delt*v0+6/(theta*delt)^2*x0)+c*(theta*delt /2*xdd0+2*v0+3/theta/delt*x0)+f); %计算(t+θdelt)时刻得速度xddtheta=6/(theta*delt)^2*(xtheta-x0)-6/theta/delt*v0-2*xdd0; %计算(t+θdelt)时刻得加速度xdd=(1-1/theta)*xdd0+1/theta*xddtheta; %计算(t+delt)时刻得加速度xd=v0+delt/2*(xdd0+xdd); %计算(t+delt)速度x=x0+delt*v0+delt^2*(2*xdd0+xdd)/6; %计算(t+delt)位移v0=xd;x0=x;xdd0=xdd;fprintf(fid1,'%10、4f',x0);%向文件中写数据t %显示计算时间步长endendfid2=fopen('disp7','rt'); %打开disp6、dat文件n=tf/delt; %disp7、dat文件中位移得个数x=fscanf(fid2,'%f',[3,n]); %将disp7、dat文件中得位移写成矩阵t=1:n;figure('numbertitle','off','name','自由度-1得位移','pos',[450 180 400 420]);plot(t,x(1,t)),grid,xlabel('时间*0、1秒'),title('自由度-1得位移与时间得关系')420]);plot(t,x(2,t)),grid,xlabel('时间*0、1秒'),title('自由度-2得位移与时间得关系') figure('numbertitle','off','name','自由度-3得位移','pos',[250 140 400 420]);plot(t,x(3,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系')。