function [ek]=shellek(dyhm,jdzb,jdzb1,dybh)
%----------------------------------------------------------------------------------------------
% dyhm: node connectivity for each element
% jdzb: coordinate of the nodes in the top surface of element
% jdzb1: coordinate of the nodes in the bottom surface of element
% dybh: the number of element whose stiffness matrix is calculated in this time of invoking.
temp=0.577350292; % value needed by Gauss integral
gaosi=[ 1 1 1;
1 -1 1;
-1 1 1;
-1 -1 1;
1 1 -1;
1 -1 -1;
-1 1 -1;
-1 -1 -1];
d=[]; % elastic matrix in 7.7-25. five order square matrix.
% Its value is determined by the character of material and isn��t given in this example.
ek=zeros(40);
for i=1:8
for j=1:3
v3i(i,j)=jdzb1(dybh(dyhm,i),j)-jdzb(dybh(dyhm,i),j);
end
end
tempi=[1 0 0];
for i=1:8
temp1=cross(v3i(i,:)',tempi);
xv2i(i,:)=temp1/norm(temp1);
temp1=cross(xv2i(i,:),v3i(i,:));
xv1i(i,:)=temp1/norm(temp1);
end
xv3i=v3i/t;
for gaosii=1:8 % begin of Gauss integral
zb1=temp*gaosi(gaosii,1);
zb2=temp*gaosi(gaosii,2);
zb3=temp*gaosi(gaosii,3);
ni(1)=1/4*(1+zb1)*(1+zb2)*(zb1+zb2-1);
ni(2)=1/4*(1-zb1)*(1+zb2)*(-zb1+zb2-1);
ni(3)=1/4*(1-zb1)*(1-zb2)*(-zb1-zb2-1);
ni(4)=1/4*(1+zb1)*(1-zb2)*(zb1-zb2-1);
ni(5)=1/2*(1-zb1^2)*(1+zb2);
ni(6)=1/2*(1-zb1)*(1-zb2^2);
ni(7)=1/2*(1-zb1^2)*(1-zb2);
ni(8)=1/2*(1+zb1)*(1-zb2^2);
v3=ni(1)*v3i(1,:)'+ni(2)*v3i(2,:)'+ni(3)*v3i(3,:)'+ni(4)*v3i(4,:)'+ni(5)*v3i(5,:)'+ni(6)*v3i(6,:)'+ni(7)*v3i(7,:)'+ni(8)*v3i(8,:)';
temp1=cross(v3,tempi');
tn2=sqrt(temp1(1)^2+temp1(2)^2+temp1(3)^2);
xv2=temp1/tn2;
temp1=cross(xv2,v3);
tn1=sqrt(temp1(1)^2+temp1(2)^2+temp1(3)^2);
xv1=temp1/tn1;
tn3=sqrt(v3(1)^2+v3(2)^2+v3(3)^2);
xv3=v3/tn3;
theta=[xv1,xv2,xv3];
for i=1:8
for j=1:3
zmtemp(i,j)=(jdzb1(dybh(dyhm,i),j)+jdzb(dybh(dyhm,i),j))/2;
end
end
pni(1,1)=1/4*(1+zb2)*(zb1+zb2-1)+(1/4+1/4*zb1)*(1+zb2);
pni(1,2)=(1/4+1/4*zb1)*(zb1+zb2-1)+(1/4+1/4*zb1)*(1+zb2);
pni(1,3)=0;
pni(2,1)=-1/4*(1+zb2)*(-zb1+zb2-1)-(1/4-1/4*zb1)*(1+zb2);
pni(2,2)=(1/4-1/4*zb1)*(-zb1+zb2-1)+(1/4-1/4*zb1)*(1+zb2);
pni(2,3)=0;
pni(3,1)=-1/4*(1-zb2)*(-zb1-zb2-1)-(1/4-1/4*zb1)*(1-zb2);
pni(3,2)=-(1/4-1/4*zb1)*(-zb1-zb2-1)-(1/4-1/4*zb1)*(1-zb2);
pni(3,3)=0;
pni(4,1)=1/4*(1-zb2)*(zb1-zb2-1)+(1/4+1/4*zb1)*(1-zb2);
pni(4,2)=-(1/4+1/4*zb1)*(zb1-zb2-1)-(1/4+1/4*zb1)*(1-zb2);
pni(4,3)=0;
pni(5,1)=-zb1*(1+zb2);
pni(5,2)=1/2-1/2*zb1^2;
pni(5,3)=0;
pni(6,1)=-1/2+1/2*zb2^2;
pni(6,2)=-2*(1/2-1/2*zb1)*zb2;
pni(6,3)=0;
pni(7,1)=-zb1*(1-zb2);
pni(7,2)=-1/2+1/2*zb1^2;
pni(7,3)=0;
pni(8,1)=1/2-1/2*zb2^2;
pni(8,2)=-2*(1/2+1/2*zb1)*zb2;
pni(8,3)=0;
x1=zmtemp(1,1);
y1=zmtemp(1,2);
z1=zmtemp(1,3);
x2=zmtemp(2,1);
y2=zmtemp(2,2);
z2=zmtemp(2,3);
x3=zmtemp(3,1);
y3=zmtemp(3,2);
z3=zmtemp(3,3);
x4=zmtemp(4,1);
y4=zmtemp(4,2);
z4=zmtemp(4,3);
x5=zmtemp(5,1);
y5=zmtemp(5,2);
z5=zmtemp(5,3);
x6=zmtemp(6,1);
y6=zmtemp(6,2);
z6=zmtemp(6,3);
x7=zmtemp(7,1);
y7=zmtemp(7,2);
z7=zmtemp(7,3);
x8=zmtemp(8,1);
y8=zmtemp(8,2);
z8=zmtemp(8,3);
dx1=v3i(1,1);
dy1=v3i(1,2);
dz1=v3i(1,3);
dx2=v3i(2,1);
dy2=v3i(2,2);
dz2=v3i(2,3);
dx3=v3i(3,1);
dy3=v3i(3,2);
dz3=v3i(3,3);
dx4=v3i(4,1);
dy4=v3i(4,2);
dz4=v3i(4,3);
dx5=v3i(5,1);
dy5=v3i(5,2);
dz5=v3i(5,3);
dx6=v3i(6,1);
dy6=v3i(6,2);
dz6=v3i(6,3);
dx7=v3i(7,1);
dy7=v3i(7,2);
dz7=v3i(7,3);
dx8=v3i(8,1);
dy8=v3i(8,2);
dz8=v3i(8,3);
jtemp(1,1)=1/4*(1+zb2)*(zb1+zb2-1)*x1+(1/4+1/4*zb1)*(1+zb2)*x1-1/4*(1+zb2)*(-zb1+zb2-1)*x2-(1/4-1/4*zb1)*(1+zb2)*x2-1/4*(1-zb2)*(-zb1-zb2-1)*x3-(1/4-1/4*zb1)*(1-zb2)*x3+1/4*(1-zb2)*(zb1-zb2-1)*x4+(1/4+1/4*zb1)*(1-zb2)*x4-zb1*(1+zb2)*x5-1/2*(1-zb2^2)*x6-zb1*(1-zb2)*x7+1/2*(1-zb2^2)*x8+1/8*zb3*(1+zb2)*(zb1+zb2-1)*dx1+1/2*zb3*(1/4+1/4*zb1)*(1+zb2)*dx1-1/8*zb3*(1+zb2)*(-zb1+zb2-1)*dx2-1/2*zb3*(1/4-1/4*zb1)*(1+zb2)*dx2-1/8*zb3*(1-zb2)*(-zb1-zb2-1)*dx3-1/2*zb3*(1/4-1/4*zb1)*(1-zb2)*dx3+1/8*zb3*(1-zb2)*(zb1-zb2-1)*dx4+1/2*zb3*(1/4+1/4*zb1)*(1-zb2)*dx4-1/2*zb3*zb1*(1+zb2)*dx5-1/4*zb3*(1-zb2^2)*dx6-1/2*zb3*zb1*(1-zb2)*dx7+1/4*zb3*(1-zb2^2)*dx8;
jtemp(1,2)=1/4*(1+zb2)*(zb1+zb2-1)*y1+(1/4+1/4*zb1)*(1+zb2)*y1-1/4*(1+zb2)*(-zb1+zb2-1)*y2-(1/4-1/4*zb1)*(1+zb2)*y2-1/4*(1-zb2)*(-zb1-zb2-1)*y3-(1/4-1/4*zb1)*(1-zb2)*y3+1/4*(1-zb2)*(zb1-zb2-1)*y4+(1/4+1/4*zb1)*(1-zb2)*y4-zb1*(1+zb2)*y5-1/2*(1-zb2^2)*y6-zb1*(1-zb2)*y7+1/2*(1-zb2^2)*y8+1/8*zb3*(1+zb2)*(zb1+zb2-1)*dy1+1/2*zb3*(1/4+1/4*zb1)*(1+zb2)*dy1-1/8*zb3*(1+zb2)*(-zb1+zb2-1)*dy2-1/2*zb3*(1/4-1/4*zb1)*(1+zb2)*dy2-1/8*zb3*(1-zb2)*(-zb1-zb2-1)*dy3-1/2*zb3*(1/4-1/4*zb1)*(1-zb2)*dy3+1/8*zb3*(1-zb2)*(zb1-zb2-1)*dy4+1/2*zb3*(1/4+1/4*zb1)*(1-zb2)*dy4-1/2*zb3*zb1*(1+zb2)*dy5-1/4*zb3*(1-zb2^2)*dy6-1/2*zb3*zb1*(1-zb2)*dy7+1/4*zb3*(1-zb2^2)*dy8;
jtemp(1,3)=1/4*(1+zb2)*(zb1+zb2-1)*z1+(1/4+1/4*zb1)*(1+zb2)*z1-1/4*(1+zb2)*(-zb1+zb2-1)*z2-(1/4-1/4*zb1)*(1+zb2)*z2-1/4*(1-zb2)*(-zb1-zb2-1)*z3-(1/4-1/4*zb1)*(1-zb2)*z3+1/4*(1-zb2)*(zb1-zb2-1)*z4+(1/4+1/4*zb1)*(1-zb2)*z4-zb1*(1+zb2)*z5-1/2*(1-zb2^2)*z6-zb1*(1-zb2)*z7+1/2*(1-zb2^2)*z8+1/8*zb3*(1+zb2)*(zb1+zb2-1)*dz1+1/2*zb3*(1/4+1/4*zb1)*(1+zb2)*dz1-1/8*zb3*(1+zb2)*(-zb1+zb2-1)*dz2-1/2*zb3*(1/4-1/4*zb1)*(1+zb2)*dz2-1/8*zb3*(1-zb2)*(-zb1-zb2-1)*dz3-1/2*zb3*(1/4-1/4*zb1)*(1-zb2)*dz3+1/8*zb3*(1-zb2)*(zb1-zb2-1)*dz4+1/2*zb3*(1/4+1/4*zb1)*(1-zb2)*dz4-1/2*zb3*zb1*(1+zb2)*dz5-1/4*zb3*(1-zb2^2)*dz6-1/2*zb3*zb1*(1-zb2)*dz7+1/4*zb3*(1-zb2^2)*dz8;
jtemp(2,1)=(1/4+1/4*zb1)*(zb1+zb2-1)*x1+(1/4+1/4*zb1)*(1+zb2)*x1+(1/4-1/4*zb1)*(-zb1+zb2-1)*x2+(1/4-1/4*zb1)*(1+zb2)*x2-(1/4-1/4*zb1)*(-zb1-zb2-1)*x3-(1/4-1/4*zb1)*(1-zb2)*x3-(1/4+1/4*zb1)*(zb1-zb2-1)*x4-(1/4+1/4*zb1)*(1-zb2)*x4+(1/2-1/2*zb1^2)*x5-2*(1/2-1/2*zb1)*zb2*x6-(1/2-1/2*zb1^2)*x7-2*(1/2+1/2*zb1)*zb2*x8+1/2*zb3*(1/4+1/4*zb1)*(zb1+zb2-1)*dx1+1/2*zb3*(1/4+1/4*zb1)*(1+zb2)*dx1+1/2*zb3*(1/4-1/4*zb1)*(-zb1+zb2-1)*dx2+1/2*zb3*(1/4-1/4*zb1)*(1+zb2)*dx2-1/2*zb3*(1/4-1/4*zb1)*(-zb1-zb2-1)*dx3-1/2*zb3*(1/4-1/4*zb1)*(1-zb2)*dx3-1/2*zb3*(1/4+1/4*zb1)*(zb1-zb2-1)*dx4-1/2*zb3*(1/4+1/4*zb1)*(1-zb2)*dx4+1/2*zb3*(1/2-1/2*zb1^2)*dx5-zb3*(1/2-1/2*zb1)*zb2*dx6-1/2*zb3*(1/2-1/2*zb1^2)*dx7-zb3*(1/2+1/2*zb1)*zb2*dx8;
jtemp(2,2)=(1/4+1/4*zb1)*(zb1+zb2-1)*y1+(1/4+1/4*zb1)*(1+zb2)*y1+(1/4-1/4*zb1)*(-zb1+zb2-1)*y2+(1/4-1/4*zb1)*(1+zb2)*y2-(1/4-1/4*zb1)*(-zb1-zb2-1)*y3-(1/4-1/4*zb1)*(1-zb2)*y3-(1/4+1/4*zb1)*(zb1-zb2-1)*y4-(1/4+1/4*zb1)*(1-zb2)*y4+(1/2-1/2*zb1^2)*y5-2*(1/2-1/2*zb1)*zb2*y6-(1/2-1/2*zb1^2)*y7-2*(1/2+1/2*zb1)*zb2*y8+1/2*zb3*(1/4+1/4*zb1)*(zb1+zb2-1)*dy1+1/2*zb3*(1/4+1/4*zb1)*(1+zb2)*dy1+1/2*zb3*(1/4-1/4*zb1)*(-zb1+zb2-1)*dy2+1/2*zb3*(1/4-1/4*zb1)*(1+zb2)*dy2-1/2*zb3*(1/4-1/4*zb1)*(-zb1-zb2-1)*dy3-1/2*zb3*(1/4-1/4*zb1)*(1-zb2)*dy3-1/2*zb3*(1/4+1/4*zb1)*(zb1-zb2-1)*dy4-1/2*zb3*(1/4+1/4*zb1)*(1-zb2)*dy4+1/2*zb3*(1/2-1/2*zb1^2)*dy5-zb3*(1/2-1/2*zb1)*zb2*dy6-1/2*zb3*(1/2-1/2*zb1^2)*dy7-zb3*(1/2+1/2*zb1)*zb2*dy8;
jtemp(2,3)=(1/4+1/4*zb1)*(zb1+zb2-1)*z1+(1/4+1/4*zb1)*(1+zb2)*z1+(1/4-1/4*zb1)*(-zb1+zb2-1)*z2+(1/4-1/4*zb1)*(1+zb2)*z2-(1/4-1/4*zb1)*(-zb1-zb2-1)*z3-(1/4-1/4*zb1)*(1-zb2)*z3-(1/4+1/4*zb1)*(zb1-zb2-1)*z4-(1/4+1/4*zb1)*(1-zb2)*z4+(1/2-1/2*zb1^2)*z5-2*(1/2-1/2*zb1)*zb2*z6-(1/2-1/2*zb1^2)*z7-2*(1/2+1/2*zb1)*zb2*z8+1/2*zb3*(1/4+1/4*zb1)*(zb1+zb2-1)*dz1+1/2*zb3*(1/4+1/4*zb1)*(1+zb2)*dz1+1/2*zb3*(1/4-1/4*zb1)*(-zb1+zb2-1)*dz2+1/2*zb3*(1/4-1/4*zb1)*(1+zb2)*dz2-1/2*zb3*(1/4-1/4*zb1)*(-zb1-zb2-1)*dz3-1/2*zb3*(1/4-1/4*zb1)*(1-zb2)*dz3-1/2*zb3*(1/4+1/4*zb1)*(zb1-zb2-1)*dz4-1/2*zb3*(1/4+1/4*zb1)*(1-zb2)*d
评论0