二维梁的固有频率和振型0001

1、一、综合实验题目和要求题目：求一二维梁的固有振型和频率。要求：用有限元理论，求一二维梁的固有振型和频率：(1)用二维梁有限元对梁进行分析数值计算求出其主振型向量和频率;(2)求出其理论精确解，精确主振型向量和频率；(3)将理论结果和计算结果进行比较。二、程序流程图读取梁的相关参数，输入梁单元数量读取梁单元的相关参数，刚度、质量矩阵等(1)r1组集总刚度阵、总质量阵k约束处理F-1求特征值、特征向量计算圆频率，特征向量归一化主程序流程图(2)读取理论特征值，计算理论圆频率和振型函数(3)画出理论和计算主振型图读取梁单元的刚度和质量矩阵形成总刚阵，总质量阵子程序流程图初始计算的特征向量形成新的特征

2、向量矩阵子程序流程图理论主振型的归一化理论主振型的绘制画出理论和计算主振型进行比较绘制计算节点处的位移子程序流程图三、实验结果1.前六阶振型同一有限元数不同阶数比较（以有限元 20为例）如下图所示:1三阶四阶五阶六阶四、实验分析对于二维梁有限元的划分(以下只对二维梁而言)，要根据需求精度进行合 理划分，既兼顾精度，同时也兼顾计算量(随着计算精度的提高，单元数量增加， 相应计算量也会增加，计算时间也会增加)，经过试验随着单元数量增加，其计 算精度也不段提高，当将梁分到七单元时，通过计算得到的主振型和频率和理论 值吻合的非常好。当梁取一单元时(elementno=1),由于梁总体只有两自由度， 故

3、只能得出前两阶主振型；当梁取二单元时(elementno=2),由于梁总体有四 自由度，故只能得出前四阶主振型；对于梁取三单元( elementno=3)以及三单 元以上(elementno>3)时，梁总体有六自由度以及更高自由度，这里只画出前 六阶主振型图。下六图是在 elementno=20的情况下，通过计算，画出前六阶的 主振型图(其中红线部分为理论主振型图，绿色五角星是计算在梁各单元节点处的振型，数量取决于梁单元划分的数目)。五、源程序清单clear all close all %各参数的设置 rou=2.7e3; %?度 A=1e-3;%横截面积E=72e9;%单性模量L=1;

4、藏长I=8.3333e-009; %1面惯性矩elementno=input('输入有限元的数量：);%有限元的数量 rodno=elementno+1; %点数alldimension=rodno*2;l=L/elementno;%单元刚度矩阵ke=E*I/lA3*12 -6*l -12 -6*l;- 6*l 4*lA2 6*l 2*lA2;- 12 6*l 12 6*l;- 6*l 2*lA2 6*l 4*lA2;%单元质量矩阵me=rou*A*l/420*156 -22*l 54 13*l;-22*l 4*lA2 -13*l-3*lA2;54 -13*l 156 22*l;13*

5、l-3平2 22*l 4平2;K=zeros(alldimension,alldimension); M=zeros(alldimension,alldimension);for i=1:elementno %总刚度矩阵和总 质量矩阵K(2*i-1,2*i-1)=ke(1,1)+K(2*i-1,2*i-1);K(2*i-1,2*i)=ke(1,2)+K(2*i-1,2*i);K(2*i-1,2*i+1)=ke(1,3)+K(2*i-1,2*i+1) ;K(2*i-1,2*i+2)=ke(1,4)+K(2*i-1,2*i+2) ;K(2*i,2*i-1)=ke(2,1)+K(2*i,2*i-1)

6、;K(2*i,2*i)=ke(2,2)+K(2*i,2*i);K(2*i,2*i+1)=ke(2,3)+K(2*i,2*i+1);K(2*i,2*i+2)=ke(2,4)+K(2*i,2*i+2); K(2*i+1,2*i-1)=ke(3,1)+K(2*i+1,2*i-1);K(2*i+1,2*i)=ke(3,2)+K(2*i+1,2*i);K(2*i+1,2*i+1)=ke(3,3)+K(2*i+1,2*i+1 );K(2*i+1,2*i+2)=ke(3,4)+K(2*i+1,2*i+2 );K(2*i+2,2*i-1)=ke(4,1)+K(2*i+2,2*i-1) ;K(2*i+2,2*i

7、)=ke(4,2)+K(2*i+2,2*i);K(2*i+2,2*i+1)=ke(4,3)+K(2*i+2,2*i+1 );K(2*i+2,2*i+2)=ke(4,4)+K(2*i+2,2*i+2 );M(2*i-1,2*i-1)=me(1,1)+M(2*i-1,2*i-1) ;M(2*i-1,2*i)=me(1,2)+M(2*i-1,2*i);M(2*i-1,2*i+1)=me(1,3)+M(2*i-1,2*i+ 1);M(2*i-1,2*i+2)=me(1,4)+M(2*i-1,2*i+ 2);M(2*i,2*i-1)=me(2,1)+M(2*i,2*i-1);M(2*i,2*i)=me(

8、2,2)+M(2*i,2*i);M(2*i,2*i+1)=me(2,3)+M(2*i,2*i+1);M(2*i,2*i+2)=me(2,4)+M(2*i,2*i+2); M(2*i+1,2*i-1)=me(3,1)+M(2*i+1,2*i- 1);M(2*i+1,2*i)=me(3,2)+M(2*i+1,2*i);M(2*i+1,2*i+1)=me(3,3)+M(2*i+1,2*i+ 1);M(2*i+1,2*i+2)=me(3,4)+M(2*i+1,2*i+ 2);M(2*i+2,2*i-1)=me(4,1)+M(2*i+2,2*i- 1);M(2*i+2,2*i)=me(4,2)+M(2*

9、i+2,2*i);M(2*i+2,2*i+1)=me(4,3)+M(2*i+2,2*i+ 1);M(2*i+2,2*i+2)=me(4,4)+M(2*i+2,2*i+2)；endKcantalever=K(3:alldimension,3:alldime nsion);Mcantalever=M(3:alldimension,3:alldim ension);V,D=eig(Kcantalever,Mcantalever); lanbuda=diag(D);%!征信omega=sqrt(diag(D)vstandard=zeros(elementno,elementno);%数值计算主振型的生

10、成for i=1:(alldimension-2)Vstandard(:,i)=V(:,i)/V(alldimension-3) ,i);for j=1:elementnovstandard(j,i)=Vstandard(2*j-1,i); 湖取有纵向位移有关的向量元素形 成新向量endendlanbudar=1.875 4.694 7.85510.996 14.137 17.279'omegaexact=lanbudar.A2*sqrt(E*I/rou/A)xrod=linspace(1/elementno,1,elementno );y1=zeros(101,6);a=(sinh(

11、lanbudar)-sin(lanbudar)./(cosh( lanbudar)+cos(lanbudar);if (elementno<=3) p=2*elementno;else p=6; endfor i=1:py1(:,i)=cosh(lanbudar(i)*x')-cos(lanbuda r(i)*x')-a(i).*(sinh(lanbudar(i)*x')-sin( l a nbudar(i)*x'); %理论主振型函数 ystandardexact= y1(:,i)/ y1(101,i); %理 论主振型归一化figure(i)plot(x,ystandardexact,'r') %理论主 振型函数的绘制hold onplot(xrod,vstandard(:,i),'g*') grid onend六、对本实验课的几点建议首先，通过几周的MATLAB课程的学习，让我掌握了其基本的操作和命令， 以及其Simulink仿真模块，对以