空间Timoshenko梁单元.rar_Timoshenko梁_matlab_shouthx7_梁_空间梁

E=2.1e11; A=4e-4;%横截面积 rou=7860; l=1;%梁长 v=0.3; a=0.02;%梁宽（横截面） h=0.02;%梁厚 Iz=a*h^3/12;%横截面惯性矩 Iy=h*a^3/12; J=Iz+Iy; k=5/6; element=40; node=41; Nonel=2; Ndof=6;%节点自由度 nc=zeros(41,2);%节点坐标 nc(1,1)=0; nc(1,2)=0; for ig=1:40 nc(ig+1,1)=nc(ig,1)+l/element; nc(ig+1,2)=nc(1,2); end % connode=[1,1,1 1 1 1 1; % 41,1,1 1 1 1 1]; % conval=[1,0,0 0 0 0 0; % 41,0,0 0 0 0 0]; k0=zeros(Ndof*Nonel,Ndof*Nonel); m0=zeros(Ndof*Nonel,Ndof*Nonel); kk=zeros(node*Ndof,node*Ndof); mm=zeros(node*Ndof,node*Ndof); index=zeros(12,1); % bcdof=zeros(node*Ndof,1); % bcval=zeros(node*Ndof,1); % [n1,n2]=size(connode); % for ni=1:n1 % ki=connode(ni,1); % bcdof((ki-1)*Ndof+1:ki*Ndof)=connode(ni,2:Ndof+1); % bcval((ki-1)*Ndof+1:ki*Ndof)=conval(ni,2:Ndof+1); % end for iel=1:element nd(1)=iel; nd(2)=iel+1; x(1)=nc(nd(1),1);y(1)=nc(nd(1),2); x(2)=nc(nd(2),1);y(2)=nc(nd(2),2); L=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); [k0,m0]=FrameElement(E,v,rou,A,Iz,Iy,J,L,k); index=femEldof(nd,Nonel,Ndof); [kk]=femAssemblyel(kk,k0,index); [mm]=femAssemblyelMass(mm,m0,index); end % [kk,mm]=femApplybc(kk,mm,bcdof); fix=[1,2,3,4,5,6,241,242,243,244,245,246]; dofall=1:246; doffree=setdiff(dofall,fix); V=zeros(246,246); % [v(doffree,:),d]=eigs(kk(doffree,doffree),mm(doffree,doffree),10,'sm');(doffree,:) [V,D]=eig(kk(doffree,doffree),mm(doffree,doffree)); [lambda,ki]=sort(diag(D)); omega=sqrt(lambda); omega1=sqrt(lambda)/(2*pi); V=V(:,ki);