二维梁的固有频率和振型

一、综合实验题目和要求
题目：求一二维梁的固有振型和频率。

要求：用有限元理论，求一二维梁的固有振型和频率：
（1） 用二维梁有限元对梁进行分析数值计算求出其主振型向量和频率； （2） 求出其理论精确解，精确主振型向量和频率； （3） 将理论结果和计算结果进行比较。

二、程序流程图
三、实验结果
1.前六阶振型
同一有限元数不同阶数比较（以有限元20为例）如下图所示：
00.10.20.30.40.50.60.70.80.9
一阶
-0.8
-0.6-0.4-0.200.20.40.60.81
二阶
-0.8
-0.6-0.4-0.200.20.40.60.81
三阶
-0.8
-0.6-0.4-0.200.20.40.60.8
四阶
-0.8
-0.6-0.4-0.200.20.40.60.81
五阶
-0.8
-0.6-0.4-0.200.20.40.60.81
六阶 四、实验分析
对于二维梁有限元的划分（以下只对二维梁而言），要根据需求精度进行合理划分，既兼顾精度，同时也兼顾计算量（随着计算精度的提高，单元数量增加，相应计算量也会增加，计算时间也会增加），经过试验随着单元数量增加，其计算精度也不段提高，当将梁分到七单元时，通过计算得到的主振型和频率和理论值吻合的非常好。

当梁取一单元时（elementno=1），由于梁总体只有两自由度，故只能得出前两阶主振型；当梁取二单元时（elementno=2），由于梁总体有四自由度，故只能得出前四阶主振型；对于梁取三单元（elementno=3）以及三单元以上（elementno>3）时，梁总体有六自由度以及更高自由度，这里只画出前六阶主振型图。

下六图是在elementno=20的情况下，通过计算，画出前六阶的主振型图（其中红线部分为理论主振型图，绿色五角星是计算在梁各单元节点处的振型，数量取决于梁单元划分的数目）。

五、源程序清单
clear all
close 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/l^3*[12 -6*l -12 -6*l;
-6*l 4*l^2 6*l 2*l^2;
-12 6*l 12 6*l;
-6*l 2*l^2 6*l 4*l^2];
%单元质量矩阵
me=rou*A*l/420*[156 -22*l 54 13*l; -22*l 4*l^2 -13*l -3*l^2;
54 -13*l 156 22*l;
13*l -3*l^2 22*l 4*l^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);
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)=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(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*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);
end
Kcantalever=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); %数值计算主振型的生成
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=lanbudar.^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(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 on
plot(xrod,vstandard(:,i),'g*')
grid on
end
六、对本实验课的几点建议
首先，通过几周的MATLAB课程的学习，让我掌握了其基本的操作和命令，以及其Simulink仿真模块，对以后的学习增加了一个实用的数学工具。

其次，在完成作业的过程中进一步让我熟悉了该软件的应用，最关键的是应用到专业相关的题目求解中，对于车辆动态仿真的领域将于相关软件（如ADMAS）结合进行联合仿真，让我受益匪浅。

最后，在上课的过程中，由于机器较少，希望能增加一些机器，希望老师能提供一些除幻灯片外的其他关于MATLAB的学习资料。

总之，感谢游霞老师的悉心指导。