二维梁的固有频率和振型

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

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

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

4、009;%截面惯性矩elementno=input('输入有限元的数量:');%有限元的数量rodno=elementno+1;%节点数alldimension=rodno*2;l=L/elementno;%单元刚度矩阵ke=E*I/lA3*12-6*l-12-6*l;- 6*l4平26*l2平2;- 126*l126*l;- 6*l2*lA26*l4*lA2;%单元质量矩阵me=rou*A*l/420*156-22*l5413*l;-22*l4*lA2-13*l-3*lA2;54-13*l15622*l;13*l-3*lA222*l4*lA2;K=zeros(alldimen

5、sion,alldimension);M=zeros(alldimension,alldimension);fori=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);K(2*i,2*i)=ke(2,2)+K(2*i,2*i);K(2*i,

6、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)=ke(4,2)+K(2*i+2,2*i);K(2*i+2,2*i+1)=ke(

7、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(2,2)+M(2*i,2*i);M(2*i,2*i+1)=me(2,3)+M(2*i,2*i

8、+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*i+2,2*i);M(2*i+2,2*i+1)=me(4,3)+M(2*i+2,2*i+1);x=0:

9、0.01:1;M(2*i+2,2*i+2)=me(4,4)+M(2*i+2,2*i+2);endKcantalever=K(3:alldimension,3:alldimension);Mcantalever=M(3:alldimension,3:alldimension);V,D=eig(Kcantalever,Mcantalever);lanbuda=diag(D);%特征值omega=sqrt(diag(D)vstandard=zeros(elementno,elementno);%数值计算主振型的生成fori=1:(alldimension-2)Vstandard(:,i)=V(:,i

10、)/V(alldimension-3),i);forj=1:elementnovstandard(j,i)=Vstandard(2*j-1,i);%提取有纵向位移有关的向量元素形成新向量endendlanbudar=1.8754.6947.85510.99614.13717.279'omegaexact=lanbudar.A2*sqrt(E*I/rou/A)xrod=linspace(1/elementno,1,elementno);y1=zeros(101,6);a=(sinh(lanbudar)-sin(lanbudar)./(cosh(lanbudar)+cos(lanbudar

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