工程振动与控制大作业

《工程振动与控制大作业》由会员分享，可在线阅读，更多相关《工程振动与控制大作业（10页珍藏版）》请在装配图网上搜索。

1、工程振动与控制大作业题目：作业A【作业A试分别用8单元和128单元FEM法求解如图所示两端简支等截面梁（总长度为厶，弯曲刚度EJ,密度和横截面积分别为p和力，m = pAL/16, k = EJ/16I?）的前5阶固有频率的近似值，若边界条件改为两端固定重新进行 计算。（要求用MATLAB编程计算，附程序清单）。编程思想：（1）离散化FEM的皋本思想是将结构离散化，将无限H由度的问题转化多自由的求解问题。以八单元为例，将梁分为8个单元，则共有九个节点，在不考虑边界条件时, 共有18个自rh度，再加上两个质量块的两个位移自由度，则该系统共有20个自 由度。任取某单兀进行分析，在该单兀两端有两个节

2、点，共有四个自由度，有课 本P99知在局部坐标系下梁单元的质量阵和刚度阵分别如下:15622154-13/2224Z213/-3Z25413/156-221.-13 /-312-2214121261-126161412-612”-126112-61.612卩-614Z2.pALMe = 420由于此单元分析具有任意性，因此，该系统的8个单元在其相应的局部坐标系下 的质量阵和刚度阵具有相同的形式。（2总刚度阵和总质量阵的组装系统总的自由度为q = 91，92，不，92oF，系统第一个单元的自由度为血，他，他，94卩，设将单元自由度转换为总体自由度的转换矩阵为S,即S1 =4TQ293Q4.=1

3、* QbQ2/Q3 * Q20卩10000100001000010000000000.004x20则el = SMeSt , Rel = SKeSL o同理可得Me2. Me3.MeQ和心2、*3兀8，则该梁的总质量阵和总刚度阵如下:8 8由于有两个质量块和两个弹赞，需要对总体质量阵和刚度阵进行修正，由于在求 单元刚度的质量阵和刚度阵时是基F系统的动能和势能，因此:E；=扌尬尿9 + |加於0jW _ kqjCjg - jkq9q5 + -kc + y kq3-kq13q2O -+ ykq;0因此,m(19,19)应加m, m(20,20)应加m, k(5,5)应加匕 k(19,19)应加k,

4、 /c(5,19)应减 k, R(19,5)应减k, k(13,13)应加k, k(20,20)应加上,k(13,20)应减上,k(20,13)应减k,则 由以上分析町得系统的总质最阵Af20x20，总刚度阵/20x20(3) 施加边界条件在FEM方法中只考虑位移边界条件，不考虑力边界条件。当两端简支时，山=0, q17 = 0,则将总质量阵和总刚度阵中的第1行、第 1列和第17行、第17列直接划去，故可得系统最终的质量阵M18X18,刚度阵 18x18当两端固定时，Q1 = 0, $2 = 0, 917=0, 918 = 0,则将总质量阵和总刚度 阵中的第1行、第1列、第2行、第2列和第17

5、行、第17列、第18行、第18 列直接划去，故可得系统最终的质量阵M16X16,刚度阵16X160(4) 求特征值由K-bi2M=Q.即B|J MlK-a)2 = 0可求得系统的特征值。(5) 划分128节点根据以上步骤，同理可求出将单元划分为128节点时系统的特征值。未施加 边界条件时，得到系统的总质量阵为M260x260 甩、刚丿叟阵为260x260 实现上述过程MATLAB程序%分别用8单元和128单元FEM对两端简支的等截面梁前五阶固有频率进行计算 clear all;close all;clc;format longn=input(请输入划分单元个数：n=);if (mod(n,4)

6、=0)%判断n能否被4整除errorf3入n值不是4的整数倍，请重新运行程序Jend%各个参数的输入dispf请选择边界条件:)dispC0.自由JdispC1.简支)disp(2.固定JBL=input(请输入左边边界条件：*)；BR=input(请输入右边边界条件：*)；EJ=inputf请输入 E*J(N m.T)值：*)；PA=inpul(irf输入 PA(Kg/m)值：*)；L=inputC请输入梁总长L(m)值:*)；%数据的前处理t=2*(n+l)+2;%总门由度数m=PA*L/16;%附加质量k=EJ/(16*L.A3);%附加刚度L=L/n;%每个单元长度a=eye(4);%

7、4阶单位阵，用于转换矩阵M=zeros(tzt);%总质量初值Ma=156 22壮 54 -13*L；22*L 4*LA2 13*L -3*LA2；.54 13*L 156 -22*L;-13*L -3*L.A2 -22*L4*L.A2; %质量矩阵的中间变换矩阵Me=PA*L/420*Ma;%局部坐标下的质量阵K=zeros( t,t);%总刚度阵初值Ka=|12 6*L -12 6*L;6*L 4*L.A2 -6*L 2*L.A2;.-12 -6*L 12 -6*L;6*L 2*L.A2 -6*L4*L.A2;%刚度矩阵的中间变换矩阵Ke=EJ/(L.A3)*Ka;%局部坐标下的刚度阵%总

8、质量阵和总刚度阵的组装for i=l: nT=zeros(4,t);T(:,(2*i-l):(2*i+2)=a;Mei=T*Me*T;M=M+Mei;%求总质量阵Kei=T*Ke*T;K=K+Kei;%求总刚度阵end%附加质星阵和附加刚度阵的组装M(t-l/t-l)=M(t-l/t-l)+m;%第一质量块的质量阵修正M(t/t)=M(t/t)+m;%第二质最块的质量阵修正K(n/2+l,n/2+l)=K(n/2+l,n/2+l)+k;K(t-l/t-l)=K(t-l,t-l)+k;K(n/2+l,t-l)=K(n/2+l/t-l)-k;K(t-l,n/2+l)=K(t-l,n/2+l)-k;

9、 %第一质量块的刚度阵修正K(3*n/2+lz3*n/2+l)=K(3*n/2+l/3*n/2+l)+k;K(tzt)=K(tzt)+k;K(3*n/2+l,t)=K(3*n/2+l/t)-k;K(tz3*n/2+l)=K(t3*n/2+l)-k;%第二质量块的刚度阵修正%根据边界条件消减自ftl度if(BL=l)M(1,:)=;M(:/1)=；K(1Z:)=;K(:Z1)=;else if (BL=2)K(1:2/:)=;K(:/1：2)=；else if(BL=O)else%根据左边界条件划行划列%简支时消去一行一列%固定时消去两行两列%1；1由时不作处理dispr左边界条件有误，请重新运

10、行程仔，）endendendif(BR=l)%根据右边界条件划行划列M(2*n+l-BL)/:)=;M(:/(2*n+l-BL)=;K(2*n+l-BL)/:)=;K(:/(2*n+l-BL)=;else if (BR=2)M(2*n+l-BL):(2*n+2-BL)/:)=;M(:/(2*n+l-BL):(2*n+2-BL)=;K(2*n+l-BL):(2*n+2-BL)/:)=;K(:/(2*n+l-BL):(2*n+2-BL)=； else if (BR=O)elsedispr右边界条件有误，请重新运行程序jendendend%数据的后处理%求特征值和特征向量%对特征值进行排序ww=ei

11、g( K*inv( M);w=sort(sqrt(ww);w(6:end/)=；%取前五阶固有频率dispC前五阶固有频率为：Jdisp(w=)forj=l:5%特征值太小时置0fprintf(%.6fn：w(j)；if(abs(w(j)le-2)w(j)=O;endend程序说明及运行结果程序设计中设置的参数有：划分单元格数n一一必须为4的倍数左右边界条件BL, BR0为|由端，1为简支端，2为固定端梁的属性参数EJ (E是弹性模最，J为抗弯刚度)一一单位是Nm2梁属性参数PA (P为密度，A为截面积)单位是Kg/m梁的总长度L单位是m程序中总白由度的排序方式为：吗，&i，a)2f &2,

12、33, %，3卄1，% + 1，wn+2/ Wn+3)T在程序运行时，为方便起见，取EJTN-m?, PA=lKg/m, L=lm运行结果如下：以8单元的两端简支情况为例，在matlab窗口输入和输出数据如下，其它 情况类似，汇总在后面表格中：请输入划分单元个数：n=8请选择边界条件：0.自由1. 简支2. 固定请输入左边边界条件：1请输入右边边界条件：1请输入E*J（N m八2）值:1请输入PA（Kg/m）值：1请输入梁总长L（m）值：1前五阶固有频率为：w =0.9993430.9999199.87616339.49183988.941428汇总表格:1两端简支，不同单元数前五阶固有频率近

13、似解对比（rad/s）模态阶数i8单元32单元128单元10.9993430.9993430.99934320.9999190.9999190.99991939.8761639.8760019.8760001439.49183939.48162739.481586588.94142888.82760688.8271452.两端固定，不冋单元数前五阶固有频率近似解对比（rad/s）模态阶数i8单元32单元128单元10.9998980.9998980.99989820.9999640.9999640.999964322.37726022.37537822375371461.71390861.67

14、509461.6749405121.190871120.905530120.904368结果分析（1）由于有限元方法是将无限自由度的问题转化为多门由度问题，因此相当于 给系统附加了刚度，单元数越多，越接近于真实解，附加刚度也就越小, 因此FEM求的固有频率均大丁该系统真实的固有频率，并且随着单元数 目的增加，固有频率逐渐减小，并逐渐收敛到真实的固有频率。（2）单元数目增加后，用于描述位移场的节点数增加，因此位移场的描述就更 准确，从而根据位移场求岀的质量阵和刚度阵也就更好的反应真实情况和 高阶的振动形态，因而结果就更加精确。（3）边界条件有简支变为两端固定后，限制了梁两端的转角，相当于给系统增 加了附加的刚度，因此系统的各阶固有频率都有相应的提髙。