工程振动与控制大作业.docxVIP

工程振动与控制大作业题目：作业A【作业A】试分别用8单元和128单元FEM法求解如图所示两端简支等截面梁（总长度为，弯曲刚度EJ，密度和横截面积分别为和，，）的前5阶固有频率的近似值，若边界条件改为两端固定重新进行计算。（要求用MATLAB编程计算，附程序清单）。编程思想：（1）离散化FEM的基本思想是将结构离散化，将无限自由度的问题转化多自由的求解问题。以八单元为例，将梁分为8个单元，则共有九个节点，在不考虑边界条件时，共有18个自由度，再加上两个质量块的两个位移自由度，则该系统共有20个自由度。任取某单元进行分析，在该单元两端有两个节点，共有四个自由度，有课本P99知在局部坐标系下梁单元的质量阵和刚度阵分别如下：由于此单元分析具有任意性，因此，该系统的8个单元在其相应的局部坐标系下的质量阵和刚度阵具有相同的形式。（2）总刚度阵和总质量阵的组装系统总的自由度为，系统第一个单元的自由度为，设将单元自由度转换为总体自由度的转换矩阵为，即得则，。同理可得，则该梁的总质量阵和总刚度阵如下：由于有两个质量块和两个弹簧，需要对总体质量阵和刚度阵进行修正，由于在求单元刚度的质量阵和刚度阵时是基于系统的动能和势能，因此：因此，，，，，应减，应减，，，应减，应减，则由以上分析可得系统的总质量阵，总刚度阵。（3）施加边界条件在FEM方法中只考虑位移边界条件，不考虑力边界条件。当两端简支时，，，则将总质量阵和总刚度阵中的第1行、第1列和第17行、第17列直接划去，故可得系统最终的质量阵，刚度阵。当两端固定时，，，，，则将总质量阵和总刚度阵中的第1行、第1列、第2行、第2列和第17行、第17列、第18行、第18列直接划去，故可得系统最终的质量阵，刚度阵。（4）求特征值由=0，即即可求得系统的特征值。（5）划分128节点根据以上步骤，同理可求出将单元划分为128节点时系统的特征值。未施加边界条件时，得到系统的总质量阵为，总刚度阵为。实现上述过程MATLAB程序%分别用8单元和128单元FEM对两端简支的等截面梁前五阶固有频率进行计算clear all;close all;clc;format longn=input(请输入划分单元个数：n=);if (mod(n,4)~=0) %判断n能否被4整除 error(输入n值不是4的整数倍，请重新运行程序)end%各个参数的输入disp(请选择边界条件：)disp( 0. 自由)disp( 1. 简支)disp( 2. 固定)BL=input(请输入左边边界条件： );BR=input(请输入右边边界条件： );EJ=input(请输入E*J(N·m.^2)值 : );PA=input(请输入PA(Kg/m)值： );L=input(请输入梁总长L(m)值： );%数据的前处理t=2*(n+1)+2; %总自由度数m=PA*L/16; %附加质量k=EJ/(16*L.^3); %附加刚度L=L/n; %每个单元长度a=eye(4); %4阶单位阵，用于转换矩阵M=zeros(t,t); %总质量初值Ma=[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]; %质量矩阵的中间变换矩阵Me=PA*L/420*Ma; %局部坐标下的质量阵K=zeros(t,t); %总刚度阵初值Ka=[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]; %刚度矩阵的中间变换矩阵Ke=EJ/(L.^3)*Ka; %局部坐标下的刚度阵%总质量阵和总刚度阵的组装for i=1:n T=zeros(4,t); T(:,(2*i-1):(2*i+2))=a; Mei=T*Me*T; M=M+Mei; %求总质量阵 Kei=T*Ke*T; K=K+Kei; %求总刚度阵end%附加质量阵和附加刚度阵的组装M(t-1,t-1)=M(t-1,t-1)+m; %第一质量块的质量阵修正M(t,t)=M(t,t)+m; %第二质量块的质量阵修正K(n/2+1,n/2+1)=K(n/