bilinearquadelementstiffness.m

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

function y = BilinearQuadElementStiffness(E,NU,t,x1,y1,x2,y2,x3,y3,x4,y4,p)
%此子程序计算4节点矩形单元刚度矩阵%

a=x3-x1;
b=y4-y2;

C1=b/(3*a)+(1-NU)*a/(6*b);
C2=NU/4+(1-NU)/8;
C3=a/(3*b)+(1-NU)*b/(6*a);
C4=(1-NU)*a/(12*b)-b/(3*a);
C5=NU/4-(1-NU)/8;
C6=a/(6*b)-(1-NU)*b/(6*a);
C7=b/(6*a)-(1-NU)*a/(6*b);
C8=(1-NU)*b/(12*a)-a/(3*b);
%k=zeros[8,8];
k(1,1)=C1;
k(2,1)=C2;
k(1,2)=C2;
k(2,2)=C3;
k(3,1)=C4;
k(1,3)=C4;
k(3,2)=-1*C5;
k(2,3)=-1*C5;
k(3,3)=C1;
k(4,1)=C5;
k(1,4)=C5;
k(4,2)=C6;
k(2,4)=C6;
k(4,3)=-1*C2;
k(3,4)=-1*C2;
k(4,4)=C3;
k(5,1)=C1/(-2);
k(1,5)=C1/(-2);
k(5,2)=-1*C2;
k(2,5)=-1*C2;
k(5,3)=C7;
k(3,5)=C7;
k(5,4)=-1*C5;
k(4,5)=-1*C5;
k(5,5)=C1;
k(6,1)=-1*C2;
k(1,6)=-1*C2;
k(6,2)=C3/(-2);
k(2,6)=C3/(-2);
k(6,3)=C5;
k(3,6)=C5;
k(6,4)=C8;
k(4,6)=C8;
k(6,5)=C2;
k(5,6)=C2;
k(6,6)=C3;
k(7,1)=C7;
k(1,7)=C7;
k(7,2)=C5;
k(2,7)=C5;
k(7,3)=C1/(-2);
k(3,7)=C1/(-2);
k(7,4)=C2;
k(4,7)=C2;
k(7,5)=C4;
k(5,7)=C4;
k(7,6)=-1*C5;
k(6,7)=-1*C5;
k(7,7)=C1;
k(8,1)=-1*C5;
k(1,8)=-1*C5;
k(8,2)=C8;
k(2,8)=C8;
k(8,3)=C2;
k(3,8)=C2;
k(8,4)=C3/(-2);
k(4,8)=C3/(-2);
k(8,5)=C5;
k(5,8)=C5;
k(8,6)=C6;
k(6,8)=C6;
k(8,7)=-1*C2;
k(7,8)=-1*C2;
k(8,8)=C3;
if p==1
    D=(E*t)/(1-NU*NU);
elseif p==2
    D=(E*t)/((1+NU)*(1-2*NU));
end
k=D*k;
y=k;