二维梁的固有频率和振型

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

2、，总质量阵否是（1）子程序流程图特征向量归一化提取与纵向位移有关的向量元素是否归一化完毕形成新的特征向量矩阵否是（2）初始计算的特征向量子程序流程图理论主振型的绘制绘制计算节点处的位移是否画图完毕画出理论和计算主振型进行比较否是（3）理论主振型的归一化三、实验结果1.前六阶振型 同一有限元数不同阶数比较（以有限元20为例）如下图所示： 一阶二阶 三阶 四阶 五阶 六阶四、实验分析对于二维梁有限元的划分（以下只对二维梁而言），要根据需求精度进行合理划分，既兼顾精度，同时也兼顾计算量（随着计算精度的提高，单元数量增加，相应计算量也会增加，计算时间也会增加），经过试验随着单元数量增加，其计算精度也不

3、段提高，当将梁分到七单元时，通过计算得到的主振型和频率和理论值吻合的非常好。当梁取一单元时（elementno=1），由于梁总体只有两自由度，故只能得出前两阶主振型；当梁取二单元时（elementno=2），由于梁总体有四自由度，故只能得出前四阶主振型；对于梁取三单元（elementno=3）以及三单元以上（elementno>3）时，梁总体有六自由度以及更高自由度，这里只画出前六阶主振型图。下六图是在elementno=20的情况下，通过计算，画出前六阶的主振型图（其中红线部分为理论主振型图，绿色五角星是计算在梁各单元节点处的振型，数量取决于梁单元划分的数目）。五、 源程序清单clea

4、r allclose all%各参数的设置rou=2.7e3; %密度A=1e-3; %横截面积E=72e9; %弹性模量L=1; %梁长I=8.3333e-009; %截面惯性矩 elementno=input('输入有限元的数量:'); %有限元的数量 rodno=elementno+1; %节点数 alldimension=rodno*2; l=L/elementno; %单元刚度矩阵 ke=E*I/l3*12 -6*l -12 -6*l; -6*l 4*l2 6*l 2*l2; -12 6*l 12 6*l; -6*l 2*l2 6*l 4*l2; %单元质量矩阵me=

5、rou*A*l/420*156 -22*l 54 13*l; -22*l 4*l2 -13*l -3*l2; 54 -13*l 156 22*l;13*l -3*l2 22*l 4*l2; 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*

6、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,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)=

7、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(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)=

8、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+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*

9、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); 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(

10、D); %特征值 omega=sqrt(diag(D) vstandard=zeros(elementno,elementno); %数值计算主振型的生成 for i=1:(alldimension-2) Vstandard(:,i)=V(:,i)/V(alldimension-3),i); for j=1:elementno vstandard(j,i)=Vstandard(2*j-1,i); %提取有纵向位移有关的向量元素形成新向量 end end lanbudar=1.875 4.694 7.855 10.996 14.137 17.279' omegaexact=lanbuda

11、r.2*sqrt(E*I/rou/A) x=0:0.01:1; xrod=linspace(1/elementno,1,elementno); y1=zeros(101,6); a=(sinh(lanbudar)-sin(lanbudar)./(cosh(lanbudar)+cos(lanbudar); if (elementno<=3) p=2*elementno; else p=6; end for i=1:p y1(:,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') %理论主振型函数的绘制hold on plot(xrod,vstandard(:,i),'g*') grid onend 六、 对本实验课的几点建议首先，通过几周的MATLAB课程的学习，让我掌握了其基本的操作和命令，以及其Simuli