d2_csq.m

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

function [B,E,K,M] = D2_CSQ (x1, y1, x2, y2, x3, y3, x4, y4, Em, hT, miu)

% returns a 2D linear constant stress quadrilateral (CSQ) element stiffness 
% and mass matrix in global coordinates
% INPUT:
%       x1, y1 - coordinates of joint 1 of the quadrilateral element
%       x2, y2 - coordinates of joint 2 of the quadrilateral element
%       x3, y3 - coordinates of joint 3 of the quadrilateral element
%       x4, y4 - coordinates of joint 4 of the quadrilateral element
%       Em     - elastic modulus for isotropic material
%       hT     - element thickness at the nodes [h1 h2 h3 h4]
%       miu    - Poisson's ratio
% OUTPUT:
%       K,M    - stiffness,mass matrix (M is divided by density !!!)
%       B      - differentiated shape functions matrix
%       E      - elastic modulus matrix for plain stress
% ---------------------------------------------------------------------------

% E   - elastic modulus matrix for plain stress:
%       [E11 E12 E13; E12 E22 E23; E13 E23 E33] 
E  = Em/(1-miu^2).*[ 1  miu  0; miu  1  0; 0  0  (1-miu)/2 ];
K  = zeros(8,8); K2 = zeros(8,8);
M  = zeros(8,8); M2 = zeros(8,8);
B  = zeros(3,8); B2 = zeros(3,8);
int_point = [ -sqrt(3/5)  0  sqrt(3/5) ]; % 3-point Gauss integration

% integrating stiffness matrix
for i=1:3
   if i==1 niu=int_point(1); end;
   if i==2 niu=int_point(2); end;
   if i==3 niu=int_point(3); end;
   for k=1:3
      if k==1 xi=int_point(1); end;
      if k==2 xi=int_point(2); end;
      if k==3 xi=int_point(3); end;
      R = 1/4.*[  (-1+niu) (1-niu) (1+niu) (-1-niu) ; (-1+xi) (-1-xi) (1+xi) (1-xi) ];
      % Jacobian of 2D transform connecting (x,y) with (xi,niu) 
      J =  R * [ x1 y1; x2 y2; x3 y3; x4 y4 ];
      dN = inv(J)*R; % dN = [ dNi / dx ; dNi / dy ] 
      % B1 - differentiated shape functions matrix
      B1 = [dN(1,1)    0    dN(1,2)    0    dN(1,3)    0    dN(1,4)    0;
              0    dN(2,1)    0    dN(2,2)    0    dN(2,3)    0    dN(2,4);
           dN(2,1) dN(1,1) dN(2,2) dN(1,2) dN(2,3) dN(1,3) dN(2,4) dN(1,4)];
        
      N1 = 1/4*(1-xi)*(1-niu); % shape function of quadrilateral 
      N2 = 1/4*(1+xi)*(1-niu); % shape function of quadrilateral 
      N3 = 1/4*(1+xi)*(1+niu); % shape function of quadrilateral 
      N4 = 1/4*(1-xi)*(1+niu); % shape function of quadrilateral
      
      N = [diag(ones(1,8).*N1)  diag(ones(1,8).*N2) diag(ones(1,8).*N3) diag(ones(1,8).*N4) ];
        
      hi  = hT(1)*N1 + hT(2)*N2 + hT(3)*N3 + hT(4)*N4; 
      M1 = hi.*N*N'.*det(J) ;
      K1 = hi.*B1'*E*B1.*det(J) ;
      
      if xi==0 & niu==0 B = B1; end;  % B taken at the center of element

      if k==1 M2 = M2 + 5/9.*M1; end; % Gauss integration for M
      if k==2 M2 = M2 + 8/9.*M1; end; % Gauss integration for M
      if k==3 M2 = M2 + 5/9.*M1; end; % Gauss integration for M
      
      if k==1 K2 = K2 + 5/9.*K1; end; % Gauss integration for K
      if k==2 K2 = K2 + 8/9.*K1; end; % Gauss integration for K
      if k==3 K2 = K2 + 5/9.*K1; end; % Gauss integration for K
   end;
   if i==1 M = M + 5/9.*M2; end; % Gauss integration for M
   if i==2 M = M + 8/9.*M2; end; % Gauss integration for M
   if i==3 M = M + 5/9.*M2; end; % Gauss integration for M
   
   if i==1 K = K + 5/9.*K2; end; % Gauss integration for K
   if i==2 K = K + 8/9.*K2; end; % Gauss integration for K
   if i==3 K = K + 5/9.*K2; end; % Gauss integration for K
   
   K2 = K2.*0; M2 = M2.*0;
end;


% ------------------end 
%
%           ^ niu
%           |
%    4      |     3
%      o.........o
%     .     |     .
%    .    -------- . --> xi
%   .       |       . 
%  o.................o
%  1                  2
%
%