rod01.m

对杆进行有限元分析

         clear all %清内存
         close all %关闭所有图形窗口
         ea=1;  %截面积
         ee=1;  %弹性模量
         er=1;%mass density
         xy=[0,0,1,0;%单元两个节点的坐标
             1,0,2,0;
             2,0,3,0;
             3,0,4,0;
             4,0,5,0;
             5,0,6,0;
             6,0,7,0];%单元坐标信息
         ndg=[8,8,1,8;
             1,8,2,8;
             2,8,3,8;
             3,8,4,8;
             4,8,5,8;
             5,8,6,8;
             6,8,7,8];%单元自由度信息
         ne=7;%单元数
         nfe=4;%单元自由度数
         nf=7;%系统自由度数
         sk=zeros(nf,nf);%定义总刚度矩阵
                  sm=zeros(nf,nf);%定义总矩阵

         m=zeros(nfe);%定义工作数组
         for kk=1:ne,%---------------------形成总刚度矩阵--------------
         ek=zeros(4,4);
         em=zeros(4,4);
         xkk1=xy(kk,1);
         ykk1=xy(kk,2);
         xkk2=xy(kk,3);
         ykk2=xy(kk,4);
         lkk=(xkk1-xkk2)^2+(ykk1-ykk2)^2;
         lkk=sqrt(lkk);
         ca=(xkk2-xkk1)/lkk;
         sa=(ykk2-ykk1)/lkk;
         ck=ee*ea/lkk;
         T1=[ca,sa,0,0;
             0,0,ca,sa];
         T2=T1.';
         ek1=ck*[1,-1;-1,1];
         em1=er*ea*lkk*[1/3,1/6;1/6,1/3];
         ek=T2*ek1*T1;
         em=T2*em1*T1;
         
         %ek=[ca2,csa,-ca2,-csa;csa,sa2,-csa,-sa2;-ca2,-csa,ca2,csa;-csa,-sa2,csa,sa2];
         %ek=ck*ek;
            for i=1:nfe,%组装总刚度矩阵
                 m(i)=ndg(kk,i);
             end
             for i=1:nfe,
                 mi=m(i);
                 if(mi<=nf),
                 for j=1:nfe,
                     mj=m(j);
                      if(mj<=nf),
                        sk(mi,mj)=sk(mi,mj)+ek(i,j);
                        sm(mi,mj)=sm(mi,mj)+em(i,j);

                      end
                   end
                  end
             end
         end %for kk-----------------------------------------总刚度、总量矩阵形成完毕
         [v,s]=eig(sk,sm);%求解广义特征值问题
         s=sqrt(s)%rad/s
         sn=diag(s);%取对角元
         sn=sn'%转置运算
         s1=1/2*pi/7;%杆纵向振动理论解
         sa=[s1,3*s1,5*s1,7*s1,9*s1,11*s1,13*s1]