二维梁的固有频率和振型.docVIP

一、综合实验题目和要求 题目：求一二维梁的固有振型和频率。 要求：用有限元理论，求一二维梁的固有振型和频率： 用二维梁有限元对梁进行分析数值计算求出其主振型向量和频率； 求出其理论精确解，精确主振型向量和频率； 将理论结果和计算结果进行比较。 二、程序流程图 读取梁的相关参数,输入梁单元数量 读取梁的相关参数,输入梁单元数量 读取梁单元的相关参数，刚度、质量矩阵等 组集总刚度阵、总质量阵 约束处理 求特征值、特征向量 计算圆频率，特征向量归一化 读取理论特征值，计算理论圆频率和振型函数 画出理论和计算主振型图 主程序流程图 （1） （2） （3） 子程序流程图 读取梁单元的刚度和质量矩阵 与下一梁单元的刚度和质量矩阵相加 是否叠加完毕 形成总刚阵，总质量阵 否 是 （1） 子程序流程图 子程序流程图 特征向量归一化 提取与纵向位移有关的向量元素 是否归一化完毕 形成新的特征向量矩阵 否 是 （2） 初始计算的特征向量 子程序流程图 理论主振型的绘制 绘制计算节点处的位移 是否画图完毕 画出理论和计算主振型进行比较 否 是 （3） 理论主振型的归一化 三、实验结果 1.前六阶振型 同一有限元数不同阶数比较（以有限元20为例）如下图所示： 一阶 二阶 三阶 四阶 五阶 六阶 四、实验分析 对于二维梁有限元的划分（以下只对二维梁而言），要根据需求精度进行合理划分，既兼顾精度，同时也兼顾计算量（随着计算精度的提高，单元数量增加，相应计算量也会增加，计算时间也会增加），经过试验随着单元数量增加，其计算精度也不段提高，当将梁分到七单元时，通过计算得到的主振型和频率和理论值吻合的非常好。当梁取一单元时（elementno=1），由于梁总体只有两自由度，故只能得出前两阶主振型；当梁取二单元时（elementno=2），由于梁总体有四自由度，故只能得出前四阶主振型；对于梁取三单元（elementno=3）以及三单元以上（elementno3）时，梁总体有六自由度以及更高自由度，这里只画出前六阶主振型图。下六图是在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*