《有限单元法》--王勖成，习题2.12 MATLAB 程序

节点、单元编号

只计算到节点位移，应力没计算

 

%% 《有限单元法》--王勖成，习题2.12 MATLAB 程序

% 仅供参考

function fem()

    P=[-10   0    10   0    10    0      -10    0]';%结构节点载荷列向量
    [K,P]=xQiYi(P,make_K());
    solu(K,P);
    
    %% 存储各单元的节点 x、y 值
        function [a,b]=get_xy(m)
            x=[0    12     0;
                0       12     12;];
            y=[4     0       0;
                4      4       0];
            a=x(m,:);
            b=y(m,:);
        end
    
    %% 求解并显示
        function solu(K,P)
            n=1000;
            JieDianWeiYi=inv(K)*P
            %         YingLi=WeiYi(3,1)/12*2e5  %x向拉、压时可用
            plot([0 12 12 0],[4 4 0 0],'b');hold on;xlim([-1 14]);ylim([-4 8]);
            plot(JieDianWeiYi([1 3 5 7],1)'*n+[0 12 12 0],JieDianWeiYi([2 4 6 8],1)'*n+[4 4 0 0],'r');
        end
    
    %%  用“对角元素改 1 法” ，引入位移约束并使 总的K 不奇异
        function [b,a]=xQiYi(P,K)
           
            %---------------------------------------------------------------------------------%
            d=[0  0  1  1  1   1   0  0]';%左侧两结点铰支
            %         共 4 个节点，顺序为 ：左上角开始1 ，顺时针排序
            %         共 8个位移分量，[dx1 dy1  dx2  dy2  ...  ]
            %         位移为 0的分量，对应项填 0，其他填非 0 值
            %         例：右侧两点铰支，则为 d=[ 1   1   0   0    0   0   1   1]'
           
            for i=1:1:8
                if d(i,1)==0
                    for  j=1:1:8
                        K(i,j)=0;
                        K(j,i)=0;
                    end
                    P(i,1)=0;
                    K(i,i)=1;
                end
            end
            a=P;
            b=K;
        end
    
    %%  组装总的结构刚度矩阵
        function K=make_K()
            K=[su([1 1],[1 1])     su(0,[1 2])      su([1 2],[1 3])         su([1 3],0);
                su(0,[2 1])             su(0,[2 2])      su(0,[2 3])              su(0,0);
                su([2 1],[3 1])       su(0,[3 2])      su([2 2],[3 3])         su([2 3],0);
                su([3 1],0)             su(0,0)           su([3 2],0)                su([3 3],0)];
        end
    
    %% 分块计算矩阵的和
        function b=su(a,b)
            K1=get_K(1);
            K2=get_K(2);
            if a==0
                K1=zeros(6,6);
                a=[1,1];
            end
            if b==0
                K2=zeros(6,6);
                b=[1,1];
            end
            h11=a(1)*2-1; h12=a(1)*2; l11=a(2)*2-1; l12=a(2)*2;
            h21=b(1)*2-1; h22=b(1)*2; l21=b(2)*2-1; l22=b(2)*2;
            b=[K1(h11,l11)+K2(h21,l21)   K1(h11,l12)+K2(h21,l22);
                K1(h12,l11)+K2(h22,l21)   K1(h12,l12)+K2(h22,l22)];
        end
    
    %%  计算 m 单元的 刚度矩阵
        function K=get_K(m)
            %-------------------------------------------------------%
            v=0.3;
            D=2e5/(1-v^2)*[ 1        ,v        ,0         ;
                v        ,1        , 0        ;
                0        ,0        ,(1-v)/2];
            %弹性矩阵
           
            [x,y]=get_xy(m);
            B=get_B(x,y);
            A=get_A(x,y);
            K=B'*D*B*A;
        end
    
    %%  计算单元应变矩阵  B
        function B=get_B(x,y)
            [~,b,c]=get_abc(x,y);
            B=zeros(3,6);
            [x,y]=get_xy(1);
            A=get_A(x,y);
            for i=1:1:3
                B(:,i*2-1:i*2)=[b(i)   0;
                    0    c(i);
                    c(i)  b(i)];
            end
            B=B/(2*A);
        end
    
    %% 计算 ai   bi   ci
        function [a,b,c]=get_abc(x , y)
            a=zeros(2,3);
            b=a;
            c=a;
            for l=0:1:2
                i=l+1;
                j=rem(l+1,3)+1;
                m=rem(l+2,3)+1;
               
                a(i)=det([x(j)    y(j);
                    x(m)    y(m)]);
                b(i)= y(j) - y(m);
                c(i)= -x(j) + x(m);
            end
        end
    
    %% 计算单元面积
        function A=get_A(x,y)
            A=-det([1   x(1)      y(1);  1   x(2)      y(2) ; 1     x(3)      y(3)])/2;
        end
    end