bilinearquadelementcaculate.m

四节点矩形单元matlab有限元程序

function [x,y]=BilinearQuadElementCaculate(P,U,K,N)
%     LinearTriangleElementCaculate的作用是:由刚度矩阵、力的边界条件和位移边界条件求解线性方程组
%     此程序有三种方法：刚度矩阵分块求解法、结构刚度矩阵的主元素置大数法、结构刚度矩阵主元素置1法。

%    刚度矩阵分块求解法.最终把刚度矩阵、结点力矩阵、位移矩阵重新组合成如下形式，其中UA为位移边界条件给出的已知位移初值，PB为已知的已知结点力
%    初值，同时相应的调整刚度矩阵。
% [PA]  [KA   KB] [UA]
% [  ] =[       ]*[  ]
% [PB]  [KC   KD] [UB]
i=1;
KB=K;
PB=P;
KA=[];
for I=1:N
    if (abs(U(I)-888))>1.0e-5
        KA=[KA,K(1:N,I)];
        KB(:,I-i+1)=[];   %前面的已经消除了，所以要用 I－i＋1
        PB(I-i+1)=[];
        UA(i)=U(I);
        i=i+1;
    end
end
j=1;
KC=KA;
KD=KB;
KA=[];
KB=[];
for I=1:N
    if (abs(U(I)-888))>1.0e-5
        KA=[KA;KC(I-j+1,:)];
        KC(I-j+1,:)=[];
        KB=[KB;KD(I-j+1,:)];
        KD(I-j+1,:)=[];   %前面的已经消除了，所以要用 I－i＋1
        j=j+1;
    end
end
for I=1:(N-i+1)
    if (abs(PB(I)-888))<1.0e-5
        PB(I)=0;
    end
end
UB=inv(KD)*(PB'-KC*UA');

s=1;
for I=1:N
    if (abs(U(I)-888))<1.0e-5
        U(I)=UB(s);
        s=s+1;
    end
end
P=K*U';

%   刚度矩阵分块求解法


%   结构刚度矩阵的主元素置大数法。是一种近似算法
% for I=1:N
%     if (abs(U(I)-888))>1.0e-5
%         K(I,I)=1.0e30;
%         P(I)=U(I)*1.0e30;
%         if (abs(U(I)-0))<1.0e-5
%             U(I)=1.0e-5;
%         end
%     end
% end
% U=inv(K)*P';
% P=K*U;
%   结构刚度矩阵的主元素置大数法


%   结构刚度矩阵主元素置1法。是一种准确的算法
% KK=K;
% UU=U;
% for I=1:N
%     if (abs(U(I)-888))>1.0e-5
%         if (abs(P(I)-888))<1.0e-5
%             P(I)=0;
%         end
%         P=P-K(I,:)*U(I);
%         K(:,I)=0;
%         K(I,:)=0;
%         K(I,I)=1;
%     end 
% end
% U=inv(K)*P';
% P=KK*U;
%   结构刚度矩阵主元素置1法      


x=P;
y=U;