ex10102.m

《The finite element method using Matlab》源程序及电子书 用Matlab进行有限元分析

v2(1)=v3(2)*v1(3)-v3(3)*v1(2);   % v2 is cross product of v3 and v1
v2(2)=v3(3)*v1(1)-v3(1)*v1(3);
v2(3)=v3(1)*v1(2)-v3(2)*v1(1);
sum=sqrt(v2(1)^2+v2(2)^2+v2(3)^2); 
v2(1)=v2(1)/sum;	% make vector v2 a unit vector
v2(2)=v2(2)/sum;
v2(3)=v2(3)/sum;

%---------------------------------------------------------
%   construct nodal variable transformation matrix 
%---------------------------------------------------------

a(1,1)=v1(1);
a(2,1)=v1(2);
a(3,1)=v1(3);
a(1,2)=v2(1);
a(2,2)=v2(2);
a(3,2)=v2(3);
a(1,3)=v3(1);
a(2,3)=v3(2);
a(3,3)=v3(3);

ainv=inv(a);

for i=1:nnel
i1=(i-1)*ndof+1;
i2=i1+1;
i3=i2+1;
i4=i3+1;
i5=i4+1;
i6=i5+1;

trsh(i1,i1)=1.0;
trsh(i2,i2)=1.0;
trsh(i3,i3)=1.0;
trsh(i4,i4)=ainv(1,1);
trsh(i4,i5)=ainv(1,2);
trsh(i4,i6)=ainv(1,3);
trsh(i5,i4)=ainv(2,1);
trsh(i5,i5)=ainv(2,2);
trsh(i5,i6)=ainv(2,3);
trsh(i6,i4)=ainv(3,1);
trsh(i6,i5)=ainv(3,2);
trsh(i6,i6)=ainv(3,3);
end
  
%------------------------------------------------------
%  strain transformation matrix
%------------------------------------------------------

rot6(1,1)=v1(1)^2;
rot6(1,2)=v1(2)^2;
rot6(1,3)=v1(3)^2;
rot6(1,4)=v1(1)*v1(2);
rot6(1,5)=v1(2)*v1(3);
rot6(1,6)=v1(1)*v1(3);
rot6(2,1)=v2(1)^2;
rot6(2,2)=v2(2)^2;
rot6(2,3)=v2(3)^2;
rot6(2,4)=v2(1)*v2(2);
rot6(2,5)=v2(2)*v2(3);
rot6(2,6)=v2(1)*v2(3);
rot6(3,1)=v3(1)^2;
rot6(3,2)=v3(2)^2;
rot6(3,3)=v3(3)^2;
rot6(3,4)=v3(1)*v3(2);
rot6(3,5)=v3(2)*v3(3);
rot6(3,6)=v3(1)*v3(3);
rot6(4,1)=2.0*v1(1)*v2(1);
rot6(4,2)=2.0*v1(2)*v2(2);
rot6(4,3)=2.0*v1(3)*v2(3);
rot6(4,4)=v1(1)*v2(2)+v2(1)*v1(2);
rot6(4,5)=v1(2)*v2(3)+v2(2)*v1(3);
rot6(4,6)=v1(3)*v2(1)+v2(3)*v1(1);
rot6(5,1)=2.0*v2(1)*v3(1);
rot6(5,2)=2.0*v2(2)*v3(2);
rot6(5,3)=2.0*v2(3)*v3(3);
rot6(5,4)=v2(1)*v3(2)+v3(1)*v2(2);
rot6(5,5)=v2(2)*v3(3)+v3(2)*v2(3);
rot6(5,6)=v2(3)*v3(1)+v3(3)*v2(1);
rot6(6,1)=2.0*v3(1)*v1(1);
rot6(6,2)=2.0*v3(2)*v1(2);
rot6(6,3)=2.0*v3(3)*v1(3);
rot6(6,4)=v3(1)*v1(2)+v1(1)*v3(2);
rot6(6,5)=v3(2)*v1(3)+v1(2)*v3(3);
rot6(6,6)=v3(3)*v1(1)+v1(3)*v3(1);

%------------------------------------------------------
%  material property matrix transformation
%------------------------------------------------------

dmtxt=rot6'*dmtx*rot6;

%------------------------------------------------------
%  numerical integration loop
%------------------------------------------------------

for intx=1:nglx
r=point1(intx);                    % sampling point in x-axis
wtx=weight1(intx,1);               % weight in x-axis
for inty=1:ngly
s=point1(inty,2);                  % sampling point in y-axis
wty=weight1(inty,2) ;              % weight in y-axis
for intz=1:nglz
t=pointz(intz);                  % sampling point in z-axis
wtz=weightz(intz) ;              % weight in z-axis

[shape2,dhdr,dhds]=feisoq4(r,s);     % compute 2-D shape functions and
                                    % derivatives at sampling point

%[shape1,dhdz]=feisol2(t);     % compute 1-D shape functions and
                                    % derivatives at sampling point
shape1(1)=0.5*(1-t);
shape1(2)=0.5*(1+t);
dhdt(1)=-0.5;
dhdt(2)=0.5;

%--------------------------------------------------
%   compute jacobian matrix and its inverse
%--------------------------------------------------

hz=thick*0.5*t;
hzdt=thick*0.5;

aj(1,1)=dhdr(1)*xcoord(1)+dhdr(2)*xcoord(2)+dhdr(3)*xcoord(3)+ ...
   dhdr(4)*xcoord(4)+dhdr(1)*hz*v3(1)+dhdr(2)*hz*v3(1)+ ...
   dhdr(3)*hz*v3(1)+dhdr(4)*hz*v3(1);
aj(2,1)=dhds(1)*xcoord(1)+dhds(2)*xcoord(2)+dhds(3)*xcoord(3)+ ...
   dhds(4)*xcoord(4)+dhds(1)*hz*v3(1)+dhds(2)*hz*v3(1)+ ...
   dhds(3)*hz*v3(1)+dhds(4)*hz*v3(1);
aj(3,1)=hzdt*v3(1);
aj(1,2)=dhdr(1)*ycoord(1)+dhdr(2)*ycoord(2)+dhdr(3)*ycoord(3)+ ...
   dhdr(4)*ycoord(4)+dhdr(1)*hz*v3(2)+dhdr(2)*hz*v3(2)+ ...
   dhdr(3)*hz*v3(2)+dhdr(4)*hz*v3(2);
aj(2,2)=dhds(1)*ycoord(1)+dhds(2)*ycoord(2)+dhds(3)*ycoord(3)+ ...
   dhds(4)*ycoord(4)+dhds(1)*hz*v3(2)+dhds(2)*hz*v3(2)+ ...
   dhds(3)*hz*v3(2)+dhds(4)*hz*v3(2);
aj(3,2)=hzdt*v3(2);
aj(1,3)=dhdr(1)*zcoord(1)+dhdr(2)*zcoord(2)+dhdr(3)*zcoord(3)+ ...
   dhdr(4)*zcoord(4)+dhdr(1)*hz*v3(3)+dhdr(2)*hz*v3(3)+ ...
   dhdr(3)*hz*v3(3)+dhdr(4)*hz*v3(3);
aj(2,3)=dhds(1)*zcoord(1)+dhds(2)*zcoord(2)+dhds(3)*zcoord(3)+ ...
   dhds(4)*zcoord(4)+dhds(1)*hz*v3(3)+dhds(2)*hz*v3(3)+ ...
   dhds(3)*hz*v3(3)+dhds(4)*hz*v3(3);
aj(3,3)=hzdt*v3(3);

ajinv=inv(aj);
det3=det(aj);

%-----------------------------------------------------------
%   compute global derivatives
%-----------------------------------------------------------

for i=1:nnel
derivg(1,i)=ajinv(1,1)*dhdr(i)+ajinv(1,2)*dhds(i);
derivg(2,i)=ajinv(2,1)*dhdr(i)+ajinv(2,2)*dhds(i);
derivg(3,i)=ajinv(3,1)*dhdr(i)+ajinv(3,2)*dhds(i);
dhdz(1,i)=ajinv(1,3)*hzdt;
dhdz(2,i)=ajinv(2,3)*hzdt;
dhdz(3,i)=ajinv(3,3)*hzdt;
end

%-----------------------------------------------------------
%   compute strain-displacement matrix called bmtx
%-----------------------------------------------------------

bmtx=zeros(ndof,edof);

for i=1:nnel
i1=(i-1)*ndof+1;
i2=i1+1;
i3=i2+1;
i4=i3+1;
i5=i4+1;
i6=i5+1;

gk1=derivg(1,i)*hz+shape2(i)*dhdz(1,i);
gk2=derivg(2,i)*hz+shape2(i)*dhdz(2,i);
gk3=derivg(3,i)*hz+shape2(i)*dhdz(3,i);

bmtx(1,i1)=derivg(1,i);	% elements associated with epsilon_x
bmtx(1,i4)=gk1*(-v2(1)); 
bmtx(1,i5)=gk1*v1(1);
bmtx(2,i2)=derivg(2,i); % elements related to epsilon_y
bmtx(2,i4)=gk2*(-v2(2));
bmtx(2,i5)=gk2*v1(2);
bmtx(3,i3)=derivg(3,i); % elements related to epsilon_z
bmtx(3,i4)=gk3*(-v2(3));
bmtx(3,i5)=gk3*v1(3);
bmtx(4,i1)=derivg(2,i); % elements related to gamma_xy
bmtx(4,i2)=derivg(1,i);
bmtx(4,i4)=gk2*(-v2(1))+gk1*(-v2(2)); 
bmtx(4,i5)=gk2*v1(1)+gk1*v1(2);
bmtx(5,i2)=derivg(3,i); % elements related to gamma_yz
bmtx(5,i3)=derivg(2,i);
bmtx(5,i4)=gk3*(-v2(2))+gk2*(-v2(3)); 
bmtx(5,i5)=gk3*v1(2)+gk2*v1(3);
bmtx(6,i1)=derivg(3,i); % elements related to gamma_xz
bmtx(6,i3)=derivg(1,i);
bmtx(6,i4)=gk3*(-v2(1))+gk1*(-v2(3)); 
bmtx(6,i5)=gk3*v1(1)+gk1*v1(3);
end

detwt=det3*wtx*wty*wtz;

k=k+bmtx'*dmtxt*bmtx*detwt;  % element matrix in local axes

end
end
end	% end of numerical integration loop
   
kt=trsh'*k*trsh;   % element matrix in global axis


index=feeldof(nd,nnel,ndof);% extract system dofs associated with element

kk=feasmbl1(kk,kt,index);  % assemble element matrices 

end

%-----------------------------
%   apply boundary conditions
%-----------------------------

[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);

%----------------------------
%  solve the matrix equation
%----------------------------

disp=kk\ff;   

num=1:1:sdof;
displace=[num' disp]                % print nodal displacements
     			     
%------------------------------------------------------------------------