trs3d.m

6根杆四个节点

%-------    Program TRS3D.M is used to calculate frequencies and eigenvectors
%-------     of space truss. This routine is wrote by Chen Huaihai. All rights
%-------     are reserved. Except permission,no one are allowed to use it.
%-------                             1998.01.02

%------- er=mass density; ea=area; ee=Young'selastic modul
         er=1.;
         ea=1.;
         ee=1.;
%------  ne=elements;nfe=dofs of each element;nf=total dofs of the system
         ne=7;
         nfe=6;
         nf=6;
%------- xyz=3d coordinates of nodes from first to last one
         xyz=[0,0,0;
              1,0,0;
              1,1,0;
              0,1,0;
              1,1,-1;
              1,0,-1];
%-----  ndg=the number of dof of each element,the number of fixed dof must
%-----      great than than nf (nf=the number of dofs of the system)
         ndg=[7,7,7,1,2,3;
              1,2,3,4,5,6;
              4,5,6,7,7,7;
              7,7,7,4,5,6;
              1,2,3,7,7,7;
              4,5,6,7,7,7;
              1,2,3,7,7,7];
%-----  nen=number of the two nodes of each element
         nen=[1 2;
              2 3;
              3 4;
              1 3;
              2 4;
              3 5;
              2 6];
%-----   li,lj,la are for the use of gplot
         li=nen(:,1)
         lj=nen(:,2)
         la=sparse(li,lj,1,nf,nf);
         gplot(la,xyz)
%------  sk=global stiffness matrix;sm=global mass matrix;m=a work matix
         sk=zeros(nf,nf);
         sm=zeros(nf,nf);
         m=zeros(nfe);
%------  loop for every element to form global matrix sk, sm
         for kk=1:ne,
         knd1=nen(kk,1);
         knd2=nen(kk,2);
         xkk1=xyz(knd1,1);
         ykk1=xyz(knd1,2);
         zkk1=xyz(knd1,3);
         xkk2=xyz(knd2,1);
         ykk2=xyz(knd2,2);
         zkk2=xyz(knd2,3);
         lkk=(xkk1-xkk2)^2+(ykk1-ykk2)^2+(zkk2-zkk1)^2;
         lkk=sqrt(lkk)
         cx=(xkk2-xkk1)/lkk;
         cy=(ykk2-ykk1)/lkk;
         cz=(zkk2-zkk1)/lkk;
         cx2=cx*cx;
         cy2=cy*cy;
         cz2=cz*cz;
         cxy=cx*cy;
         cyz=cy*cz;
         cxz=cx*cz;
         ck=ee*ea/lkk;
         cm=er*ea*lkk/2.;
         ek=[cx2,cxy,cxz,-cx2,-cxy,-cxz;cxy,cy2,cyz,-cxy,-cy2,-cyz;
             cxz,cyz,cz2,-cxz,-cyz,-cz2;-cx2,-cxy,-cxz,cx2,cxy,cxz;
             -cxy,-cy2,-cyz,cxy,cy2,cyz;-cxz,-cyz,-cz2,cxz,cyz,cz2];
         ek=ck*ek;
         em=[1,0,0,0,0,0;
             0,1,0,0,0,0;
             0,0,1,0,0,0;
             0,0,0,1,0,0;
             0,0,0,0,1,0;
             0,0,0,0,0,1];
         em=cm*em;
            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

         [v,d]=eig(sk,sm);
         ev=diag(d);
         ev=sqrt(ev);
%------  rearrange eigenvalues and eigenvectors by ascent
         [evs,ser]=sort(ev);
         for j=1:nf,
             for i=1:nf,
             vec(i,j)=v(i,ser(j));
             end
         end
%------  print evs=eigenvalues(rad/s) and eigenvectors
         evs
         vec