哈工大结构动力学大作业2012春

结构动力学大作业对于如下结构，是研究质量块的质量变化和在简支梁上位置的变化对整个系统模态的影响。1以上为一个简支梁结构。集中质量块放于梁上，质量块距简支梁的左端点距离为L.将该简支梁简化为欧拉伯努利梁，并离散为N个单元。每个单元有两个节点，四个自由度。单元的节点位移可表示为：1122,,,evv则单元内一点的挠度可计作：带入边界条件：1332210)(xaxaxaaxv01)0(avxv3322102)(LaLaLaavLxv110ddaxvx2321232ddLaLaaxvLx10va1234NNNNN建立了单元位移模式后，其动能势能均可用节点位移表示。单元的动能为：00111()222llTTTkeeeeeyEdxqNNdxqqmqt其中m为单元质量阵，并有：0lTmNNdx带入公式后积分可得：2222156225413224133541315622420133224lllllllmllllll单元势能可表示为2220011()()222TllTTepeeeeqyEEIdxEINNdxqqKqx其中K为单元刚度矩阵，并有0()lTKEINNdx2232212612664621261266264llllllEIklllllll以上为单元类型矩阵，通过定义全局位移矩阵，可以得到系统刚度矩阵和系统质量矩11a)2(1)(3211222LvvLa)(1)(22122133LvvLa1232133222231)(LxLxxvLxLxxv22232332223LxLxvLxLx24231211)()()()()(xNvxNxNvxNxv阵。当集中质量块加到简支梁上时，可以认为系统的刚度矩阵不变，而质量矩阵在有基础上有所增加。当集中质量块位于第i个单元内，距单元左节点的位移为0x，则其动能可表示为：000()()2TTeemVNxNx则质量块附加到总体刚度矩阵的子矩阵尾为：000()()TmNxNx例如当质量为0m的质量块位于单元的右节点上时，其质量阵可表示为00000000000100000mm根据以上理论，编写matlab程序。（相关参数：密度7860KG/M3，长度L=1m，截面尺寸0.02m*0.02m）以下为matlab源代码：主程序clear;clc;Beam_InputData541;%输入相关参数，划分了40个单%元k=zeros(No_nel*No_dof,No_nel*No_dof);%初始化单元刚度阵和质量阵m=zeros(No_nel*No_dof,No_nel*No_dof);kk=zeros(Sys_dof,Sys_dof);mm=zeros(Sys_dof,Sys_dof);%初始化总体质量阵和刚度阵index=zeros(No_nel*No_dof,1);foriel=1:No_elnd(1)=iel;nd(2)=iel+1;leng=0.025;[k,m]=BeamElement11(prop,leng);%调用子函数形成单元刚度矩阵%质量矩阵index=femEldof(nd,No_nel,No_dof);kk=femAssemble1(kk,k,index);mm=femAssemble1(mm,m,index);%调用子函数形成总体刚度矩阵end%和质量矩阵Mzhiliangkuai=1/5*prop(3)*prop(6)*1;%将集中质量块加进总体质量矩Ms=zeros(82,82);%之中Ms(41,41)=Mzhiliangkuai;%此为质量块为于梁中央时的处mm=Ms+mm;理方法kk(1,:)=[];%引入边界条件，划去相应的kk(80,:)=[];%行和列kk(:,1)=[];kk(:,80)=[];mm(1,:)=[];mm(80,:)=[];mm(:,1)=[];mm(:,80)=[];[V,D]=eig(kk,mm);%求解特征值和特征向量%或使用雅克比迭代求解特征值与特征向量，以下为代码%function[aa,v]=jac(a)%此程序用jacobi方法求实对称矩阵的全部特征值和特征向量%输入x：nxn矩阵%----------------------------------------------------%%n=length(a);%aa=zeros(1,n);%u=zeros(n,n);%l=0;%v是n阶单位矩阵%v=eye(n,n);%while(1)%l=l+1;%fm=0.0;%fori=1:n%forj=1:n%ifi==j%continue;%else%d=abs(a(i,j));%ifdfm%fm=d;%p=i;%q=j;%end%end%end%end%iffm0.0001%break;%end%x=2*a(p,q)*(sign(a(p,p)-a(q,q)));%y=abs(a(p,p)-a(q,q));%ify~=0%cn=sqrt((1+y/sqrt(x^2+y^2))/2.0);%else%cn=1/sqrt(2);%end%ify~=0%sn=x/sqrt(x^2+y^2)/(2*cn);%elseify==0&a(p,q)0%sn=1/sqrt(2);%elseify==0&a(p,q)0%sn=-1/sqrt(2);%end%%更新a%a_t=a;%fori=1:n%ifi~=p&i~=q%a(i,p)=a_t(i,p)*cn+a_t(i,q)*sn;%a(i,q)=-a_t(i,p)*sn+a_t(i,q)*cn;%a(p,i)=a(i,p);%a(q,i)=a(i,q);%end%end%a(p,p)=a_t(p,p)*cn*cn+2*a_t(p,q)*cn*sn+a_t(q,q)*sn*sn;%a(q,q)=a_t(p,p)*sn*sn-2*a_t(p,q)*sn*cn+a_t(q,q)*cn*cn;%a(p,q)=0.0;%a(q,p)=a(p,q);%%更新v%v_t=v;%fori=1:n%v(i,p)=v_t(i,p)*cn+v_t(i,q)*sn;%v(i,q)=-v_t(i,p)*sn+v_t(i,q)*cn;%end%end%fori=1:n%aa(1,i)=a(i,i);%end%end%fori=1:10;%w(i,1)=sqrt(D(i,i))/(2*pi);%end%E=prop(1);Iz=prop(8);rho=prop(3)*prop(6);L=1;%i=(1:10)';%omega2=i.*i*pi^2*sqrt(E*Iz/(rho*L^4));%omega3=omega2/(2*pi);%fori=1:39;%x1(i)=i*0.025;%y1(i)=V(2*i,3);%end%plot(x1,y1);各子函数：Beam_InputData541No_el=40;numberofelementsNo_nel=2;No_dof=2;No_node=(No_nel-1)*No_el+1;Sys_dof=No_node*No_dof;prop(1)=2.1e11;%elasticmodulusprop(3)=7860;prop(4)=0.02;prop(5)=0.02;prop(6)=0.02*0.02;prop(7)=0.02*0.02^3/12;prop(8)=0.02*0.02^3/12;BeamElement11function[k,m]=BeamElement11(prop,leng);%,Opt_mass)E=prop(1);u=prop(2);rho=prop(3);A=prop(6);Iz=prop(8);G=E/(2*(1+u));c=E*Iz/(leng^3);k0=[126*leng-126*leng;6*leng4*leng^2-6*leng2*leng^2;-12-6*leng12-6*leng;6*leng2*leng^2-6*leng4*leng^2];k=c*k0;mass=rho*A*leng;m0=[15622*leng54-13*leng;22*leng4*leng^213*leng-3*leng^2;5413*leng156-22*leng;-13*leng-3*leng^2-22*leng4*leng^2];m=mass/420*m0;femAssemble1function[kk]=femAssemble1(kk,k,index)eldof=length(index);fori=1:eldofii=index(i);forj=1:eldofjj=index(j);kk(ii,jj)=kk(ii,jj)+k(i,j);endendfunction[index]=femEldof(nd,No_nel,No_dof)k=0;fori=1:No_nelstart=(nd(i)-1)*No_dof;forj=1:No_dofk=k+1;index(k)=start+j;endend首先我分析了无集中质量块时梁的模态特性，提取了前十阶固有频率和前四阶振型。然后分析了当集中质量块分别为简支梁质量1/40,1/20.1/10,1/5,1/2,1倍和2倍时这八种工况。对每种工况又进行了质量块距简支梁左端1/8L,1/4L,3/8L和1/2L时每种情形的模态特性。一，简支梁（无质量块的固有频率和振型）的模态分析前十阶固有频率：阶数固有频率（Hz）146.87862187.50713421.89184750.033351171.938761687.620972297.105983000.434093797.6669104688.8936前四阶振型：一弯：二弯：三弯：四弯：经过验证，计算结果与连续体的解析解吻和良好，这充分说明了简化模型的科00.10.20.30.40.50.60.70.80.91-4-3-2-10123400.10.20.30.40.50.60.70.80.91-4-3-2-10123400.10.20.30.40.50.60.70.80.91-4-3-2-10123400.10.20.30.40.50.60.70.80.91-4-3-2-101234学性和算法的精确性。二，集中质量块在简支梁上的位置对模态的影响：集中质量块的质量为0.5M时（M代表简支梁的质量）的前四阶固有频率：1/8L1/4L3/8L1/2LN=143.65511137.95589334.26243533.085684N=2152.28186144.69197166.19866187.50714N=3337.04577387.0983407.1014341.09605N=4640.46341750.03328630.31814750.03328前四阶固有频率随位置变化的曲线：可以看到，第一阶固有频率随着从梁一端而接近梁中央而逐渐下降，二至四阶固有频率变化规律不明显。质量块在简支梁上不同位置处的振型：1/8L处：11.522.533.54010020030040050060070080000.10.20.30.40.50.60.70.80.91-4-3-2-10123400.10.20.30.40.50.60.70.80.91-4-3-2-1012341/4L处：1/2L处：集中质量为2M时的前四阶固有频率：1/8L1/4L3/8L1/2MN=136.33822926.21302822.01278620.865716N=2116.70027123.76571155.61304187.50714N=3303.81482