beam3d.m

将梁分成9个单元

        clear all
         E=210e9;
         G=80.8e9; Ro=7.8e3;
         %Iy=2.67e-7;%方截面
         %Iz=4.17e-7;
         %Jx=5.58e-7;
         %A=2e-3;
         %Iy=3.07e-7;% 圆截面
         %Iz=3.07e-7;
         %Jx=6.14e-7;
         %A=1.96e-3;
         Iy=0.17e-8;%2号角钢
         Iz=0.63e-8;
         Jx=360e-12;
         A=1.132e-4;
        

        xyz=[0,0,0,  1,0,0,  0.5, 1,   0;
             1,0,0,  2,0,0,  1.3, 1,   0;
             2,0,0,  3,0,0,  2.3, 1,   0;
             3,0,0,  4,0,0,  3.3, 1,   0;
             4,0,0,  5,0,0,  4.3, 1,   0;
             5,0,0,  6,0,0,  5.3, 1,   0;
             6,0,0,  7,0,0,  6.3, 1,   0;
             7,0,0,  8,0,0,  7.3, 1,   0;
             8,0,0,  9,0,0,  8.3, 1,   0;
             9,0,0,  10,0,0, 9.3, 1,   0]*0.1;
       ndg=[200,200,200,200,200,200,1,2,3,4,5,6;
            1,2,3,4,5,6,7,8,9,10,11,12;
            7,8,9,10,11,12,13,14,15,16,17,18;
            13,14,15,16,17,18,19,20,21,22,23,24;
            19,20,21,22,23,24,25,26,27,28,29,30;
            25,26,27,28,29,30,31,32,33,34,35,36;
            31,32,33,34,35,36,37,38,39,40,41,42;
            37,38,39,40,41,42,43,44,45,46,47,48;
            43,44,45,46,47,48,49,50,51,52,53,54;
            49,50,51,52,53,54,55,56,57,58,59,60];
         
         ne=10;
         nfe=12;%单元自由度数
         nf=60;%总自由度数
         sk=zeros(nf,nf);
         sm=zeros(nf,nf);
         m=zeros(nfe);
         for kk=1:ne,
         ek=zeros(nfe,nfe);
         em=zeros(nfe,nfe);
        
xi=xyz(kk,1);
yi=xyz(kk,2);
zi=xyz(kk,3);
xj=xyz(kk,4);
yj=xyz(kk,5);
zj=xyz(kk,6);
xk=xyz(kk,7);
yk=xyz(kk,8);
zk=xyz(kk,9);
lkk=sqrt((xi-xj)^2+(yi-yj)^2+(zi-zj)^2);

ek(1,1)=E*A/lkk;
ek(7,1)=-ek(1,1);
ek(1,7)=ek(7,1);

ek(2,2)=12*E*Iz/lkk^3;
ek(6,2)=6*E*Iz/lkk^2;
ek(2,6)=ek(6,2);
ek(8,2)=-ek(2,2);
ek(2,8)=ek(8,2);
ek(12,2)=ek(6,2);
ek(2,12)=ek(12,2);
ek(3,3)=12*E*Iy/lkk^3;
ek(5,3)=-6*E*Iy/lkk^2;
ek(3,5)=ek(5,3);
ek(9,3)=-ek(3,3);
ek(3,9)=ek(9,3);
ek(11,3)=ek(5,3);
ek(3,11)=ek(11,3);
ek(4,4)=G*Jx/lkk;
ek(10,4)=-ek(4,4);
ek(4,10)=ek(10,4);

ek(5,5)=4*E*Iy/lkk;
ek(9,5)=-ek(5,3);
ek(5,9)=ek(9,5);
ek(11,5)=2*E*Iy/lkk;
ek(5,11)=ek(11,5);

ek(6,6)=4*E*Iz/lkk;
ek(8,6)=-ek(12,2);
ek(6,8)=ek(8,6);
ek(12,6)=2*E*Iz/lkk;
ek(6,12)=ek(12,6);

ek(7,7)=ek(1,1);

ek(8,8)=ek(2,2);
ek(12,8)=-ek(12,2);
ek(8,12)=ek(12,8);
   
    
ek(9,9)=ek(3,3);
ek(11,9)=-ek(5,3);
ek(9,11)=ek(11,9);
    
ek(10,10)=ek(4,4);

ek(11,11)=ek(5,5);

ek(12,12)=ek(6,6);
   
em(1,1)=1/3;
em(7,1)=1/6;
em(1,7)=em(7,1);

em(2,2)=13/35;
em(6,2)=11/210*lkk;
em(2,6)=em(6,2);
em(8,2)=9/70;
em(2,8)=em(8,2);
em(12,2)=-13/420*lkk;
em(2,12)=em(12,2);

em(3,3)=em(2,2);
em(5,3)=-em(6,2);
em(3,5)=em(5,3);
em(9,3)=em(8,2);
em(3,9)=em(9,3);
em(11,3)=-em(12,2);
em(3,11)=em(11,3);
em(4,4)=Jx/3/A;
em(10,4)=Jx/6/A;
em(4,10)=em(10,4);
em(5,5)=lkk*lkk/105;
em(9,5)=em(12,2);
em(5,9)=em(9,5);
em(11,5)=-lkk*lkk/140;
em(5,11)=em(11,5);
em(6,6)=em(5,5);
em(8,6)=-em(12,2);
em(6,8)=em(8,6);
em(12,6)=em(11,5);
em(6,12)=em(12,6);
em(7,7)=em(1,1);

em(8,8)=em(2,2);
em(12,8)=-em(6,2);
em(8,12)=em(12,8);

em(9,9)=em(2,2);
em(11,9)=em(6,2);
em(9,11)=em(11,9);
em(10,10)=em(4,4);
em(11,11)=em(5,5);
em(12,12)=em(5,5);

em=em*Ro*A*lkk;
           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
%f=zeros(nf,1);
%f(nf-1)=1;
%u=inv(sk)*f
[V,D]=eig(sk,sm);
     fhz=sqrt(diag(D))/6.28;%ANSYS   
     fhz=sort(fhz);
     format short e
     fhz'