d3_lb.m

结构力学中的有限元例子,包含了7个分类文件夹

function [K,B,E,V] = D3_LB (x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4,...
   x5, y5, z5, x6, y6, z6, x7, y7, z7, x8, y8, z8, Em, miu) 

% returns a 3D linear constant stress brick (LB) element stiffness matrix
% in global coordinates
% INPUT:
%       x1, y1 - coordinates of joint 1 ... of the brick element
%       Em     - elastic modulus for isotropic material
%       miu    - Poisson's ratio
% OUTPUT:
%       K      - stiffness matrix
%       B      - differentiated shape functions matrix
%       E      - elastic modulus matrix for 3D stress
%       V      - element volume
% ---------------------------------------------------------------------------

% shape functions for 8-node 3D brick:
%         N1 = 1/8*(1-xi1)*(1-xi2)*(1-xi3);
%         N2 = 1/8*(1+xi1)*(1-xi2)*(1-xi3);
%         N3 = 1/8*(1+xi1)*(1+xi2)*(1-xi3);
%         N4 = 1/8*(1-xi1)*(1+xi2)*(1-xi3);
%         N5 = 1/8*(1-xi1)*(1-xi2)*(1+xi3);
%         N6 = 1/8*(1+xi1)*(1-xi2)*(1+xi3);
%         N7 = 1/8*(1+xi1)*(1+xi2)*(1+xi3);
%         N8 = 1/8*(1-xi1)*(1+xi2)*(1+xi3);

% E   - elastic modulus matrix for plain stress:
%       [E11 : E16; E21 : E26; ...; E61 : E66] 
s=(1-miu);
E  = Em*(1-miu)/((1+miu)*(1-2*miu)).*...
         [   1   miu/s  miu/s       0                0                0;
           miu/s   1    miu/s       0                0                0;
           miu/s miu/s    1         0                0                0;
             0     0      0   (1-2*miu)/(2*s)        0                0;
             0     0      0         0         (1-2*miu)/(2*s)         0;
             0     0      0         0                0        (1-2*miu)/(2*s)];
             
K  = zeros(24,24); K2 = zeros(24,24); K3 = zeros(24,24);
B  = zeros(6,24);
V2 = 0; V3 = 0; V = 0;

int_point = [ -sqrt(3/5)  0  sqrt(3/5) ]; % 3-point Gauss integration

% integrating stiffness matrix
for i=1:3
   if i==1 xi1=int_point(1); end;
   if i==2 xi1=int_point(2); end;
   if i==3 xi1=int_point(3); end;
   for j=1:3
      if j==1 xi2=int_point(1); end;
      if j==2 xi2=int_point(2); end;
      if j==3 xi2=int_point(3); end;
      for k=1:3
         if k==1 xi3=int_point(1); end;
         if k==2 xi3=int_point(2); end;
         if k==3 xi3=int_point(3); end;
         R = 1/8.*[  (-1+xi3+xi2-xi2*xi3) (1-xi3-xi2+xi2*xi3) (1-xi3+xi2-xi2*xi3) ...
               (-1+xi3-xi2+xi2*xi3) (-1-xi3+xi2+xi2*xi3) (1+xi3-xi2-xi2*xi3) ...
               (1+xi3+xi2+xi2*xi3) (-1-xi3-xi2-xi2*xi3);
            (-1+xi3+xi1-xi1*xi3) (-1+xi3-xi1+xi1*xi3) (1-xi3+xi1-xi1*xi3) ...
               (1-xi3-xi1+xi1*xi3) (-1-xi3+xi1+xi1*xi3) (-1-xi3-xi1-xi1*xi3) ...
               (1+xi3+xi1+xi1*xi3) (1+xi3-xi1-xi1*xi3);
            (-1+xi2+xi1-xi1*xi2) (-1+xi2-xi1+xi1*xi2) (-1-xi2-xi1-xi1*xi2) ...
               (-1-xi2+xi1+xi1*xi2) (1-xi2-xi1+xi1*xi2) (1-xi2+xi1-xi1*xi2) ...
               (1+xi2+xi1+xi1*xi2) (1+xi2-xi1-xi1*xi2)];
         % Jacobian of 3D transform connecting (x,y,z) with (xi1,xi2,xi3) 
         J =  R * [ x1 y1 z1; x2 y2 z2; x3 y3 z3; x4 y4 z4; x5 y5 z5;
                    x6 y6 z6; x7 y7 z7; x8 y8 z8];
         dN = inv(J)*R; % dN = [ dNi / dx ; dNi / dy ; dNi / dz ] 
         % B1 - differentiated shape functions matrix
         B1 = [dN(1,1) 0  0  dN(1,2) 0  0  dN(1,3) 0  0 dN(1,4) 0  0 dN(1,5) 0  0 ...
                 dN(1,6) 0  0 dN(1,7) 0  0 dN(1,8) 0  0;
               0 dN(2,1) 0  0 dN(2,2) 0  0 dN(2,3) 0 0 dN(2,4) 0 0 dN(2,5) 0 0 ...
                 dN(2,6) 0  0 dN(2,7) 0  0 dN(2,8) 0;
               0 0 dN(3,1) 0 0 dN(3,2) 0  0 dN(3,3) 0 0 dN(3,4) 0 0 dN(3,5) 0 0 ...
                 dN(3,6) 0  0 dN(3,7) 0  0 dN(3,8);
               dN(2,1) dN(1,1) 0 dN(2,2) dN(1,2) 0 dN(2,3) dN(1,3) ...
               0 dN(2,4) dN(1,4) 0 dN(2,5) dN(1,5) 0 dN(2,6) dN(1,6) 0 ...
                 dN(2,7) dN(1,7) 0 dN(2,8) dN(1,8) 0;
               0 dN(3,1) dN(2,1) 0 dN(3,2) dN(2,2) 0 dN(3,3) dN(2,3) ...
               0 dN(3,4) dN(2,4) 0 dN(3,5) dN(2,5) 0 dN(3,6) dN(2,6) 0 ...
                 dN(3,7) dN(2,7) 0 dN(3,8) dN(2,8);
               dN(3,1) 0 dN(1,1) dN(3,2) 0 dN(1,2) dN(3,3) 0 dN(1,3) ...
               dN(3,4) 0 dN(1,4) dN(3,5) 0 dN(1,5) dN(3,6) 0 dN(1,6) ...
                 dN(3,7) 0 dN(1,7) dN(3,8) 0 dN(1,8)];
         
         Vo = det(J);
         K1 = B1'*E*B1.*Vo ;

         if k==1 V2 = V2 + 5/9.*Vo; end; % volume integration
         if k==2 V2 = V2 + 8/9.*Vo; end;
         if k==3 V2 = V2 + 5/9.*Vo; end;

         if xi1==0 & xi2==0 & xi3==0 B = B1; end;  % B taken at the center of element
         
         if k==1 K2 = K2 + 5/9.*K1; end; % Gauss integration
         if k==2 K2 = K2 + 8/9.*K1; end; % Gauss integration
         if k==3 K2 = K2 + 5/9.*K1; end; % Gauss integration
      end;

      if j==1 V3 = V3 + 5/9.*V2; end; % volume integration
      if j==2 V3 = V3 + 8/9.*V2; end;
      if j==3 V3 = V3 + 5/9.*V2; end;
      
      if j==1 K3 = K3 + 5/9.*K2; end; % Gauss integration
      if j==2 K3 = K3 + 8/9.*K2; end; % Gauss integration
      if j==3 K3 = K3 + 5/9.*K2; end; % Gauss integration
      K2 = K2.*0; V2 = 0;
   end;
   
   if i==1 V = V + 5/9.*V3; end; % volume integration
   if i==2 V = V + 8/9.*V3; end;
   if i==3 V = V + 5/9.*V3; end;
   
   if i==1 K = K + 5/9.*K3; end; % Gauss integration
   if i==2 K = K + 8/9.*K3; end; % Gauss integration
   if i==3 K = K + 5/9.*K3; end; % Gauss integration
   K3 = K3.*0; V3 = 0;
end;


% ------------------end