小项目论文

1 机械振动小项目论文 班 级： 94040102 组 员 一： 2009040401039 韩姣 组 员 二： 2009040401040 马驰 指导教师： 王 克 明 2012年 05 月 22 日 2 目 录 一．项目原题 ...................................................................................... 1 二．两种方法的模块简单分析 ........................................................... 2 1.预备程序段说明 ............................................................................................................................... 2 2.传递矩阵法 ....................................................................................................................................... 2 3.矩阵迭代法 ....................................................................................................................................... 3 三．完整的程序 .................................................................................. 3 1.传递矩阵法 ....................................................................................................................................... 3 2 矩阵迭代法 ....................................................................................................................................... 4 四．程序运行结果及简要分析 ........................................................... 5 1.传递矩阵法 ....................................................................................................................................... 5 2.矩阵迭代法 ....................................................................................................................................... 7 五．体会 ............................................................................................. 8 1 一．项目原题 Score:(039) (040) Project Calculation of Natural Frequencies and Normal Modes of Torsional Vibration of a Multi-Disk Rotor System (Using matrix iteration and transfer matrix methods) System description: A simplified rotor model is as shown in the figure. The parameters are given bellow: Disk（盘） Geometry dimensions: diameter（直径）×thickness = 0.40×0.04 (m); Material properties: Young’s modulus（杨氏模量） E=2×1011 (N/m2), shear modulus（剪切模量） G=7.69×1010 (N/m2), mass densityρ =7800 (kg/m3); All the disks are identical and they are assumed rigid. Shaft(杆) Geometrical dimensions: diameter（直径） = 0.04 (m), length is as shown in the figure, a = 0.12 (m); Material properties: Young’s modulus E=2×1011 (N/m2), shear modulus G=7.69× 1010 (N/m2), mass densityρ =7800 (kg/m3); Neglect shaft mass and suppose the shaft stiffness is not affected by the disks. Task: Write computer programs to calculate all the 8 natural frequencies and the corresponding normal modes of torsional vibration of the rotor system; using matrix iteration method and transfer matrix method respectively. Requirements: Your report should include method description and calculation results (frequencies in Hz and normal modes in figures), and the two computer programs. a a a a a a 6a a k1 k2 k3 k4 k5 k6 k7 k8 2 二．两种方法的模块简单分析 1.预备程序段说明 （1）转动惯量：利用项目原题中提供的关于系统材料的原始数据，可以利用转动惯量 公式 J=m*^2/2 得出各段轴的转动惯量，从而得到对应的刚度矩阵。 （2）转动惯量矩阵（对应质量矩阵）：从示意图不难看出各段质量的分布情况，因其质 量是离散的，所以可以得到对角阵形式的转动惯量矩阵。 （3）扭转刚度矩阵（对应刚度矩阵）：扭转刚度公式为 k=π *d^4*G/(32*l)，由此获得 各阶的扭转刚度，从而得到扭转刚度矩阵。 （4）参考数据（reference data）：通过列动力矩阵可以用特征值特征向量的办法得到这 八自由度系统的各阶频率和主振型，作为检验这两种方法的参照。 2.传递矩阵法 传递矩阵即使按照顺序把前一阶段的圆盘和轴的转角和扭矩组成一个状态矢量（点矩阵 与场矩阵的乘积）传递到后面待计算的轴段上。因此，其方法的核心在于传递矩阵的构建。 单独研究轴，有：θ i L =θ i-1 R +Ti L / ki=Ri-1 R+ Ti-1 R / ki （a）； Ti L=Ri-1 R （b）；二者组合有： [θ ,T]L i=[1,1/k;0,1] [θ ,T] R i-1;类似，由盘的左右两侧转角相同也 能列出等式，整理得， [θ ,T]R i=[1,0;- ω 2J,1] [θ ,T]L i 把两矩阵相乘，得到程序中的 T=[1 1/k(i);-J(i)*omtr^2 1-J(i)*omtr^2/k(i)]； 再找解的过程中，才用频率尝试的办法，循环带入看各个被尝试的频率是否满足传递矩 阵使其为零，分三种情况： （i）若刚好为零，则该频率为其中一个固有频率 （ii）若相近的两个尝试频率的乘积小于零，则说明某固有频率在它们极小的数值差之 间，此时，采用 相似三角形的办法，近似地找到该阶固有频率，因为各尝试频率的差异很 小，所以，算得的估计值可以满足精确度。 （iii）固有频率数与自由度数相等，所以当找到的固有频率数大于自由度数时，循环自 动跳出。 算主振型时，需要把边界条件带入，本题（左端）为[1 0]'，并结合算出的各阶频率， 带入传递矩阵，便可得出各阶频率对应的主振型。 对于相似三角形的办法，用下图给予说明： 3 3.矩阵迭代法 本题中，边界非约束，所以柔度矩阵不存在，为了采用矩阵迭代法，需要把刚度矩阵做 变化让其对应的柔度矩阵存在，来满足矩阵迭代法的格式。 首先假设一个经过归一化的主振型，用动力矩阵前乘它，对通过乘法运算得到的振型矢 量进行归一化，将初始假设振型示为系统主振型的线性组合。每次迭代都把前一阶阵型除 去，从而消除了前阶振型对所求当前阶的影响，如果在规定的精确度内迭代前后的阵型相等， 则可认为得到了该阶主振型。最后把最初已经变化的刚度矩阵还原，固有频率还原。即可。 各动力矩阵可按递推公式求得： [D]s =[D]s-1 -｛A (s)｝T｛A(s)｝[M] /(ω 2 s M s) 三．完整的程序 1.传递矩阵法 %transfer matrix method for natural frequencies and normal modes clc clear %initial n=8; p=7800;%(kg/m3) th=0.04;%m d=0.04;%m r=0.2; g=7.69*10^10; %N/m^2 a=0.12;%单位是米 %moment of inertia转动惯量 for i=1:n J(i)=(1/2)*p*pi*th*r^4; end Ip=pi*d^4/32; k0=g*Ip/a; k00=g*Ip/(6*a); k=[k0 k0 k0 k0 k0 k0 k00 k0]; %reference data M=[J(1) 0 0 0 0 0 0 0 0 J(2) 0 0 0 0 0 0 0 0 J(3) 0 0 0 0 0 0 0 0 J(4) 0 0 0 0 0 0 0 0 J(5) 0 0 0 0 0 0 0 0 J(6) 0 0 0 0 0 0 0 0 J(7) 0 0 0 0 0 0 0 0 J(8)];%只要是离散系统直 接是对角阵 K=[k(2) -k(2) 0 0 0 0 0 0 -k(2) k(2)+k(3) -k(3) 0 0 0 0 0 0 -k(3) k(3)+k(4) -k(4) 0 0 0 0 0 0 -k(4) k(4)+k(5) -k(5) 0 0 0 4 0 0 0 -k(5) k(5)+k(6) -k(6) 0 0 0 0 0 0 -k(6) k(6)+k(7) -k(7) 0 0 0 0 0 0 -k(7) k(7)+k(8) -k(8) 0 0 0 0 0 0 -k(8) k(8)]; D=inv(M)*K; [A lam]=eig(D); for i=1:n omt(i)=sqrt(abs(lam(i,i))); A(:,i)=A(:,i)/A(1,i); end omt A %transfer matrix method for natural frequencies J=[J(1) J(2) J(3) J(4) J(5) J(6) J(7) J(8)]; k=[k0 k0 k0 k0 k0 k0 k00 k0]; step=1;%频率搜索 步长尽量小 几个自由 度几个根 j=1;%模态数 T21p=0; for ii=1:100000 omtr=step*(ii-1); Tp=eye(2);%单位阵 for i=1:n T=[1 1/k(i);-J(i)*omtr^2 1-J(i)*omtr^2/k(i)]*Tp;%组合单元连续矩阵 Tp=T; end if T(2,1)==0 omega(j)=omtr; j=j+1; end if T(2,1)*T21p<0 t=step/(1+abs(T21p/T(2,1)));%用相 似三角形完成 omega(j)=omtr-t;%-step/2,omega(j)=omtr 图 形有特点 o(∩_∩)o 哈哈 j=j+1; end if j>n break end T21p=T(2,1); end omega %normal modes for j=1:n SV(:,1)=[1 0]';%状态矢量 转角是假设 的 1 for i=1:n SV(:,i+1)=[1 1/k(i); -J(i)*omega(j)^2 1-J(i)*omega(j)^2/k(i)]*SV(:,i); end AM(:,j)=SV(1,:)'; end AM plot(AM) f=omega/(2*pi) 2 矩阵迭代法 % Matrix iteration clc clear %initial n=8; p=7800;%(kg/m3) th=0.04;%m d=0.04;%m r=0.2; g=7.69*10^10; %N/m^2 a=0.12;%单位是米 %moment of inertia转动惯量 for i=1:n J(i)=(1/2)*p*pi*th*r^4; end Ip=pi*d^4/32; k0=g*Ip/a; k00=g*Ip/(6*a); k=[k0 k0 k0 k0 k0 k0 k00 k0]; %reference data M=[J(1) 0 0 0 0 0 0 0 0 J(2) 0 0 0 0 0 0 0 0 J(3) 0 0 0 0 0 0 0 0 J(4) 0 0 0 0 5 0 0 0 0 J(5) 0 0 0 0 0 0 0 0 J(6) 0 0 0 0 0 0 0 0 J(7) 0 0 0 0 0 0 0 0 J(8)];%只要是离散系统直 接是对角阵 K=[k(2) -k(2) 0 0 0 0 0 0 -k(2) k(2)+k(3) -k(3) 0 0 0 0 0 0 -k(3) k(3)+k(4) -k(4) 0 0 0 0 0 0 -k(4) k(4)+k(5) -k(5) 0 0 0 0 0 0 -k(5) k(5)+k(6) -k(6) 0 0 0 0 0 0 -k(6) k(6)+k(7) -k(7) 0 0 0 0 0 0 -k(7) k(7)+k(8) -k(8) 0 0 0 0 0 0 -k(8) k(8)]; D=inv(M)*K; [A lam]=eig(D); for i=1:n omt(i)=sqrt(abs(lam(i,i))); A(:,i)=A(:,i)/A(1,i); end omt A alf=1;%(设定值) K=K+alf*M; D=inv(K)*M; for j=1:n A=[1 0 0 0 0 0 0 0]'; if j>1 Mg=B'*M*B; S=eye(n)-B*B'*M/Mg; D=D*S; end for i=1:100 B=D*A; a=B(1); B=B/a; if max((abs(B-A)))<0.0000000001 break end A=B; end AM(:,j)=A; oms(j)=1/sqrt(a); om(j)=sqrt(abs((oms(j))^2-alf)); end om AM plot(AM) f=om/(2*pi) 四．程序运行结果及简要分析 1.传递矩阵法 omt = 0.0000 128.3627 261.3835 463.0231 640.9295 660.3980 787.1890 876.0350 A = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9198 0.6674 -0.0438 -1.0000 -1.1233 -2.0169 -2.7364 1.0000 0.7658 0.1127 -1.0419 -1.0000 -0.8614 1.0511 3.7514 1.0000 0.5503 -0.4794 -0.9525 1.0000 1.2296 0.9480 -3.7776 1.0000 0.2908 -0.9120 0.1311 1.0000 0.7098 -2.0152 2.8079 1.0000 0.0078 -1.0413 1.0778 -1.0000 -1.3172 1.1013 -1.0981 1.0000 -1.6934 0.2612 0.0078 -1.0000 3.3019 -0.1355 0.0830 1.0000 -1.8411 0.3914 -0.1787 1.0000 -2.9394 0.0672 -0.0303 6 omega = 0 128.3601 261.3808 463.0230 640.9322 660.3844 787.1863 876.0339 AM = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9198 0.6674 -0.0438 -1.0000 -1.1233 -2.0169 -2.7364 1.0000 0.7658 0.1128 -1.0419 -1.0000 -0.8615 1.0511 3.7514 1.0000 0.5504 -0.4794 -0.9525 1.0000 1.2295 0.9481 -3.7775 1.0000 0.2908 -0.9120 0.1311 1.0000 0.7100 -2.0152 2.8078 1.0000 0.0079 -1.0413 1.0778 -1.0001 -1.3170 1.1012 -1.0979 1.0000 -1.6933 0.2611 0.0078 -0.9994 3.2988 -0.1341 0.0807 1.0000 -1.8410 0.3913 -0.1787 0.9996 -2.9361 0.0646 -0.0243 f = 0 20.4291 41.6000 73.6924 102.0075 105.1034 125.2846 139.4251 振型图 7 简要分析：运用传递矩阵法得到的固有频率和主振型与特征值特征向量方法得到的值大体 一致，但是在传递矩阵法运算的过程中采用了估计值的办法，比如用相似三角形取得固有频 率，导致固有频率值的取得略有偏差，进而用这些值算出来的各阶主振型也有一定偏差。虽 然有些许瑕疵，但是这种方法在自由度很多的但只想得到前几阶固有频率和主振型时就显示 了很好的灵活性与使用性。 2.矩阵迭代法 omt = 0.0000 128.3627 261.3835 463.0231 640.9295 660.3980 787.1890 876.0350 A = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9198 0.6674 -0.0438 -1.0000 -1.1233 -2.0169 -2.7364 1.0000 0.7658 0.1127 -1.0419 -1.0000 -0.8614 1.0511 3.7514 1.0000 0.5503 -0.4794 -0.9525 1.0000 1.2296 0.9480 -3.7776 1.0000 0.2908 -0.9120 0.1311 1.0000 0.7098 -2.0152 2.8079 1.0000 0.0078 -1.0413 1.0778 -1.0000 -1.3172 1.1013 -1.0981 1.0000 -1.6934 0.2612 0.0078 -1.0000 3.3019 -0.1355 0.0830 1.0000 -1.8411 0.3914 -0.1787 1.0000 -2.9394 0.0672 -0.0303 om = 0.0000 128.3627 261.3835 463.0231 640.9446 660.3980 787.1890 876.0350 AM = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9198 0.6674 -0.0438 -1.0001 -1.1237 -2.0169 -2.7364 1.0000 0.7658 0.1127 -1.0419 -0.9999 -0.8611 1.0511 3.7514 1.0000 0.5503 -0.4794 -0.9525 1.0002 1.2302 0.9480 -3.7776 1.0000 0.2908 -0.9120 0.1311 0.9998 0.7090 -2.0152 2.8079 1.0000 0.0078 -1.0413 1.0778 -1.0002 -1.3180 1.1013 -1.0981 1.0000 -1.6934 0.2612 0.0078 -0.9967 3.3135 -0.1355 0.0830 1.0000 -1.8411 0.3914 -0.1787 0.9970 -2.9499 0.0672 -0.0303 f = 0.0000 20.4296 41.6005 73.6924 102.0095 105.1056 125.2850 139.4253 8 振型图 简要分析：运用矩阵迭代法算出的固有频率和主振型与特征值特征向量法算出来的值在精 度允许的范围内完全一致，只要在判断迭代前后的两个主振型的值的差值允许范围足够小， 完全可以达到精度要求。 五．体会 这是我们第一次用 Matlab 软件解决学科问题，对于机械振动、线性代数以及计算机编 程都有了更深刻的认识。之前学软件只是单纯的学习，对于它的工程实际应用没有过多的了 解，而起初对于机械振动的计算方法也很模糊。这次这个题目让这些知识得到了实际应用， 让我们也对自己解决实际问题能力的信心有所提升，对于得出的结论也养成了想办法检验的 好习惯。每条语句的来龙去脉都弄得有理有据，对于理论方法也更加明确和清楚。这种训练 的不在于次数的多少，而在于深入研究了多少，对于自己整个学习态度与习惯的改善提升了 多少，我觉得这是我感受最深的地方。 另外，对于《机械振动》课，给我留下更深印象的是“双语”的氛围与思维很年轻的老 师，老师对于课堂上程式化的问题的观点是先进的，人性化的，在这样的氛围里面更适宜我 们专注于思考问题而不怕犯错误，不怕问浅显的问题，老师给我们的“学生，不要顾及太多， 学知识是根本”的态度让我们更敢于尝试新的上课状态和与老师交流的方式。 还有一些不想名状的东西，已经深深扎根在我们心里。虽然目前还不能知晓结果，但是 我已深深的受到了教育，得到了洗礼，那将会是我最重要的人生财富，指引我正确的人生的 一盏明灯。