beam3d_shear.m

对三维空间梁进行有限元建模和分析

        clear all
         E=210e9;
         G=E/(2*(1+0.3)); 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;
         As=0.8*A; %横截面有效剪切面积

        

        xyz=[0,0,0,  1,0,0,  0.5, 0.5, 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);

    rz=sqrt(Iz/A);
    ry=sqrt(Iy/A);
    phi_y=0;
    phi_z=0; %phi取0，表示不考虑剪切

    az=(12*E*Iz)/lkk^3/(1+phi_y);
    ay=(12*E*Iy)/lkk^3/(1+phi_z);
    cz=(6*E*Iz)/lkk^2/(1+phi_y);
    cy=(6*E*Iy)/lkk^2/(1+phi_z);
    ez=(4+phi_y)*E*Iz/lkk/(1+phi_y);
    ey=(4+phi_z)*E*Iy/lkk/(1+phi_z);
    fz=(2-phi_y)*E*Iz/lkk/(1+phi_y);
    fy=(2-phi_z)*E*Iy/lkk/(1+phi_z);


    Az=(13/35+7/10*phi_y+1/3*phi_y*phi_y+6/5*(rz/lkk)^2)/(1+phi_y)^2;
    Ay=(13/35+7/10*phi_z+1/3*phi_z*phi_z+6/5*(ry/lkk)^2)/(1+phi_z)^2;
    Bz=(9/70+3/10*phi_y+1/6*phi_y*phi_y-6/5*(rz/lkk)^2)/(1+phi_y)^2;
    By=(9/70+3/10*phi_z+1/6*phi_z*phi_z-6/5*(ry/lkk)^2)/(1+phi_z)^2;
    Cz=(11/210+11/120*phi_y+1/24*phi_y^2+(1/10-1/2*phi_y)*(rz/lkk)^2)*lkk/(1+phi_y)^2;
    Cy=(11/210+11/120*phi_z+1/24*phi_z^2+(1/10-1/2*phi_z)*(ry/lkk)^2)*lkk/(1+phi_z)^2;
    Dz=(13/420+3/40*phi_y+1/24*phi_y^2-(1/10-1/2*phi_y)*(rz/lkk)^2)*lkk/(1+phi_y)^2;
    Dy=(13/420+3/40*phi_z+1/24*phi_z^2-(1/10-1/2*phi_z)*(ry/lkk)^2)*lkk/(1+phi_z)^2;
    Ez=(1/105+1/60*phi_y+1/120*phi_y^2+(2/15+1/6*phi_y+1/3*phi_y^2)*(rz/lkk)^2)*lkk*lkk/(1+phi_y)^2;
    Ey=(1/105+1/60*phi_z+1/120*phi_z^2+(2/15+1/6*phi_z+1/3*phi_z^2)*(ry/lkk)^2)*lkk*lkk/(1+phi_z)^2;
    Fz=(1/140+1/60*phi_y+1/120*phi_y^2+(1/30+1/6*phi_y-1/6*phi_y^2)*(rz/lkk)^2)*lkk*lkk/(1+phi_y)^2;
    Fy=(1/140+1/60*phi_z+1/120*phi_z^2+(1/30+1/6*phi_z-1/6*phi_z^2)*(ry/lkk)^2)*lkk*lkk/(1+phi_z)^2;


    ek(1,1)=E*A/lkk;
    ek(7,1)=-ek(1,1);
    ek(1,7)=ek(7,1);
    ek(2,2)=az;
    ek(6,2)=cz;
    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)=ay;
    ek(5,3)=-cy;
    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)=ey;
    ek(9,5)=-ek(5,3);
    ek(5,9)=ek(9,5);
    ek(11,5)=fy;
    ek(5,11)=ek(11,5);
    ek(6,6)=ez;
    ek(8,6)=-cz;
    ek(6,8)=ek(8,6);
    ek(12,6)=fz;
    ek(6,12)=ek(12,6);
    ek(7,7)=ek(1,1);
    ek(8,8)=az;
    ek(12,8)=-cz;
    ek(8,12)=ek(12,8);
    ek(9,9)=ay;
    ek(11,9)=cy;
    ek(9,11)=ek(11,9);
    ek(10,10)=ek(4,4);
    ek(11,11)=ey;
    ek(12,12)=ez;



    em(1,1)=1/3;
    em(7,1)=1/6;
    em(1,7)=em(7,1);
    em(2,2)=Az;
    em(6,2)=Cz;
    em(2,6)=em(6,2);
    em(8,2)=Bz;
    em(2,8)=em(8,2);
    em(12,2)=-Dz;
    em(2,12)=em(12,2);
    em(3,3)=Ay;
    em(5,3)=-Cy;
    em(3,5)=em(5,3);
    em(9,3)=By;
    em(3,9)=em(9,3);
    em(11,3)=Dy;
    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)=Ey;
    em(9,5)=-Dy;
    em(5,9)=em(9,5);
    em(11,5)=-Fy;
    em(5,11)=em(11,5);
    em(6,6)=Ez;
    em(8,6)=Dz;
    em(6,8)=em(8,6);
    em(12,6)=-Fz;
    em(6,12)=em(12,6);
    em(7,7)=1/3;
    em(8,8)=Az;
    em(12,8)=-Cz;
    em(8,12)=em(12,8);
    em(9,9)=Ay;
    em(11,9)=Cy;
    em(9,11)=em(11,9);
    em(10,10)=Jx/3/A;
    em(11,11)=Ey;
    em(12,12)=Ez;
    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'