裄架结构有限元计算程序

function U=FEA_Plane_Truss_GeneralFlexuralMember(Nodes,Elements,Force)
    % 平面裄架结构有限元求解函数
    % Nodes存储各节点的坐标，每行为一个节点，第一列为x坐标，第二列为y坐标
    %         第3-5列存储未知量的编号
    % Elements存储各单元信息，每行为一个单元，前两列为该单元的两个节点编号
    %         第三列为单元长度，第四列为单元抗拉刚度，
    %         第五列为截面面积，第六列为截面惯性矩。
    % Force为载荷列阵，每三行为一组，对应一个节点的载荷，分别为fx，fy，M弯矩


    % 初始化内变量，预分配内存
    Dofs=3*size(Nodes,1);      % 结构总自由度，一个节点对应三个自由度
    EleCount=size(Elements,1); % 单元个数
    Conarr=reshape(Elements(:,7:12),6*size(Elements,1),1);
    Conarr=unique(Conarr);
    ConCount=size(Conarr,1)-1;
    K=zeros(ConCount,ConCount);        % 总刚度矩阵
    % 计算每个单元的刚度矩阵并汇总于总刚度矩阵
    for k=1:EleCount
        % 提取单元信息
        N1=Elements(k,1);
        N2=Elements(k,2);
        L =Elements(k,3);
        E =Elements(k,4);
        A =Elements(k,5);
        I =Elements(k,6);m
        % 提取节点信息
        x1=Nodes(N1,1);
        x2=Nodes(N2,1);
        y1=Nodes(N1,2);
        y2=Nodes(N2,2);
        u1=Elements(k,7);
        u2=Elements(k,10);
        v1=Elements(k,8);
        v2=Elements(k,11);
        r1=Elements(k,9);
        r2=Elements(k,12);
        % 计算单元坐标系到全局坐标系的旋转矩阵
        c=(x2-x1)/L;
        s=(y2-y1)/L;
        R=[ c s 0  0 0 0;...
           -s c 0  0 0 0;...
            0 0 1  0 0 0;...
            0 0 0  c s 0;...
            0 0 0 -s c 0;...
            0 0 0  0 0 1];
        
        % 计算单元刚度矩阵
        k1=E*A/L;
        k2=E*I/L^3;
        k3=E*I/L^2;
        k4=E*I/L;
        ke=[k1   0      0     -k1 0      0    ;...
             0   12*k2  6*k3  0   -12*k2 6*k3 ;...
             0   6*k3   4*k4  0   -6*k3  2*k4 ;...
           -k1   0      0     k1  0      0    ;...
             0   -12*k2 -6*k3 0   12*k2  -6*k3;...
             0   6*k3   2*k4  0   -6*k3  4*k4];
        % 将单元坐标系下的单元刚度矩阵旋转至全局坐标系下
        Globalke=R'*ke*R;
        % 将全局坐标系的单元刚度矩阵组装到总刚度矩阵的对应位置
        index1=[u1 v1 r1];index1(:,find(index1==0))=[];
        index2=[u2 v2 r2];index2(:,find(index2==0))=[];
        for i1=1:size(index1,2)
            for i2=1:size(index1,2)
                K(index1(i1),index1(i2))=K(index1(i1),index1(i2))+Globalke(i1,i2);
            end
        end
        for i1=1:size(index2,2)
            for i2=1:size(index2,2)
                K(index2(i1),index2(i2))=K(index2(i1),index2(i2))+Globalke(i1+3,i2+3);
            end
        end
    end
    
    
    
    % solve
    U=K\Force;
end