1121.m

用matlab编写的

% 空间梁单元有限元程序
    %本程序为第 1 子程序 单元的刚度矩阵
   
%刚度矩阵子程序开始
    
    %定义一些函数
    syms felem_i felem_e felem_g felem_as felem_l;
    stiffness_fai=12*felem_e*felem_i/(felem_g*felem_as*felem_l^2);
    %t=simple(stiffness_fai)   
    syms felem_i felem_e felem_g felem_as felem_l stiffness_fai;
    stiffness_a=12*felem_e*felem_i/(felem_l^3*(1+stiffness_fai));
    stiffness_c=6*felem_e*felem_i/(felem_l^2*(1+stiffness_fai));
    stiffness_e=(4+stiffness_fai)*felem_e*felem_i/(felem_l*(1+stiffness_fai));
    stiffness_f=(2-stiffness_fai)*felem_e*felem_i/(felem_l*(1+stiffness_fai));
    %t=simple(stiffness_22)
    
    %子程序开始
    Beam4.mat.ex=
    Beam4.mat.g=
    Beam4.sect.area=
    Beam4.sect.iz=
    Beam4.sect.iy=
    Beam4.sect.ix=
    Beam4.sect.asz=
    Beam4.sect.asy=    
    elem.longth=
    
    
    %函数中边梁与梁单元的数据接口
    felem_e=Beam4.mat.ex
    felem_g=Beam4.mat.g
    felem_l=elem.longth
        %------------结束------------
    %辅助函数az，cz，ez，fz
    felem_i=Beam4.sect.iz;
    felem_as=Beam4.sect.asy;
    stiffness_fai=eval(stiffness_fai);
    temp.faz=eval(stiffness_a);
    temp.fcz=eval(stiffness_c);
    temp.fez=eval(stiffness_e);    
    temp.ffz=eval(stiffness_f); 
        %------------结束------------
    %辅助函数ay，cy，ey，fy
    felem_i=Beam4.sect.iy;
    felem_as=Beam4.sect.asz;
    stiffness_fai=eval(stiffness_fai);
    temp.fay=eval(stiffness_a);
    temp.fcy=eval(stiffness_c);
    temp.fey=eval(stiffness_e);
    temp.ffy=eval(stiffness_f); 
        %------------结束------------
        
        
    %--------------刚度集成部分---------------    
    %刚度矩阵
    Beam4.stiffness(1,1)=0.5*elem.mat.ex*elem.sect.area/elem.longth ;
    Beam4.stiffness(2,2)=0.5*temp.faz;
    Beam4.stiffness(3,3)=0.5*temp.fay;
    Beam4.stiffness(4,4)=0.5*elem.mat.g*Beam4.sect.ix/elem.longth ;
    Beam4.stiffness(5,5)=0.5*temp.fey;
    Beam4.stiffness(6,6)=0.5*temp.fez;  
    
    Beam4.stiffness(7,1)=-elem.mat.ex*elem.sect.area/elem.longth;  
    Beam4.stiffness(6,2)=temp.fcz;
    Beam4.stiffness(8,2)=-temp.faz;
    Beam4.stiffness(12,2)=temp.fcz;    
    Beam4.stiffness(5,3)=-temp.fcy;
    Beam4.stiffness(9,3)=-temp.fay;
    Beam4.stiffness(11,3)=-temp.fcy;      
    Beam4.stiffness(10,4)=-elem.mat.g*Beam4.sect.ix/elem.longth ;   
    Beam4.stiffness(9,5)=temp.fcy;
    Beam4.stiffness(11,5)=temp.ffy;
    Beam4.stiffness(8,6)=-temp.fcz;
    Beam4.stiffness(12,6)=temp.ffz;  
    for i=1:6;
        for j=1:6;
            Beam4.stiffness(6+i,6+j)=Beam4.stiffness(i,j);
        end;
    end;
    Beam4.stiffness=Beam4.stiffness+Beam4.stiffness';