d2_bciz.m

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

function [K,M] = D2_BCIZ (x1, y1, x2, y2, x3, y3, Em, h, miu)

% returns a 2D Kirchhoff Plate bending (KPB) element stiffness matrix
% (Bazeley, Cheung, Irons & Zienkiewicz, 1965) in global coordinates
% INPUT:
%       xi, yi - coordinates of joint 'i' of the element
%       Em     - elastic modulus for isotropic material
%       h      - thickness
%       miu    - Poisson's ratio
% OUTPUT:
%       K, M   - stiffness, mass matrix (M is divided by density !!!)
% ---------------------------------------------------------------------------

K = zeros(9,9); M = zeros(9,9);  % memory allocation
A = (0.5 * det([1 1 1; x1 x2 x3; y1 y2 y3]));    % triangle area
y21=y2-y1; y13=y1-y3; y32=y3-y2; y23=y2-y3; y31=y3-y1; y12=y1-y2;
x13=x1-x3; x21=x2-x1; x32=x3-x2; x23=x2-x3; x31=x3-x1; x12=x1-x2;

% D   - isotropic plate rigidity matrix (force*length):
%       [D11 D12 D13; D12 D22 D23; D13 D23 D33] 
D = Em*h^3/(12*(1-miu^2)).*[ 1  miu  0; miu  1  0; 0  0  (1+miu)/2 ];
% additional coefficients for Ba:
g1=y13-y21+1.5*(-2*y1+y2+y3);
g2=x21-x13+1.5*(2*x1-x2-x3);
g3=y21-y32+1.5*(y1-2*y2+y3);
g4=x32-x21+1.5*(-x1+2*x2-x3);
g5=y32-y13+1.5*(y1+y2-2*y3);
g6=x13-x32+1.5*(-x1-x2+2*x3);
% Gauss integration points 'xi'
Mxi = [ 2/3 1/6 1/6; 1/6 2/3 1/6; 1/6 1/6 2/3 ];

for i=1:3
    xi1=Mxi(i,1); xi2=Mxi(i,2); xi3=Mxi(i,3);
    % B - double differentiated shape functions matrix
    Ba = [6*(xi1+xi2+xi3) 2*(y21*xi2-y13*xi3) 2*(x13*xi3-x21*xi2) 0 0 0 0 0 0;
      0 0 0 6*(xi1+xi2+xi3) 2*(y32*xi3-y21*xi1) 2*(x21*xi1-x32*xi3) 0 0 0;
      0 0 0 0 0 0 6*(xi1+xi2+xi3) 2*(y13*xi1-y32*xi2) 2*(x32*xi2-x13*xi1);
      (6*xi1+2*xi3) (2*y21*xi1+xi3*g1) (-2*x21*xi1+xi3*g2) (6*xi2+2*xi3) ...
         (-2*y21*xi2+xi3*g3) (2*x21*xi2+xi3*g4) (2*xi3) (xi3*g5) (xi3*g6);
      (2*xi1) (xi1*g1) (xi1*g2) (6*xi2+2*xi1) (2*y32*xi2+xi1*g3) (-2*x32*xi2+xi1*g4) ...
          (6*xi3+2*xi1) (-2*y32*xi3+xi1*g5) (2*x32*xi3+xi1*g6);
      (6*xi1+2*xi2) (-2*y13*xi1+xi2*g1) (2*x13*xi1+xi2*g2) (2*xi2) (xi2*g3) (xi2*g4) ...
          (6*xi3+2*xi2) (2*y13*xi3+xi2*g5) (-2*x13*xi3+xi2*g6)];
    Bb = (1/(4*A^2)).*[(y23^2)      (x32^2)      (2*x32*y23);
                       (y31^2)      (x13^2)      (2*x13*y31);
                       (y12^2)      (x21^2)      (2*x21*y12);
                       (2*y23*y31)  (2*x32*x13)  (2*(x32*y31+x13*y23));
                       (2*y31*y12)  (2*x13*x21)  (2*(x13*y12+x21*y31));
                       (2*y12*y23)  (2*x21*x32)  (2*(x21*y23+x32*y12))];
    B = Bb'*Ba;
    K = K + 1/3*(B'*D*B); % Gauss integration

    % shape functions
    N1 = xi1^2*(xi1+3*xi2+3*xi3)+2*xi1*xi2*xi3;
    N2 = xi1^2*(y21*xi2-y13*xi3)+xi1*xi2*xi3*g1;
    N3 = -xi1^2*(x21*xi2-x13*xi3)+xi1*xi2*xi3*g2;
    N4 = xi2^2*(xi2+3*xi1+3*xi3)+2*xi1*xi2*xi3;
    N5 = xi2^2*(y32*xi3-y21*xi1)+xi1*xi2*xi3*g3;
    N6 = -xi2^2*(x32*xi3-x21*xi1)+xi1*xi2*xi3*g4;
    N7 = xi3^2*(xi3+3*xi1+3*xi2)+2*xi1*xi2*xi3;
    N8 = xi3^2*(y13*xi1-y32*xi2)+xi1*xi2*xi3*g5;
    N9 = -xi3^2*(x13*xi1-x32*xi2)+xi1*xi2*xi3*g6;
    % [N] - shape functions matrix for M calculation
    %N = [diag(ones(1,9).*N1)  diag(ones(1,9).*N2) diag(ones(1,9).*N3) diag(ones(1,9).*N4) ...
    %        diag(ones(1,9).*N5) diag(ones(1,9).*N6)  diag(ones(1,9).*N7) ...
    %        diag(ones(1,9).*N8) diag(ones(1,9).*N9)];
   %%%    BN = [h 0 0; 0 h^3/12 0; 0 0 h^3/12];
    %M =  M + h*1/3*(N*N'); % Gauss integration
end;

K = K*A;
M1 = diag([1 1/12 1/12 1 1/12 1/12 1 1/12 1/12]);
M = 1/3*h*A*M1; % equal mass for all nodes 

% displacement vector: [ w1 taux1 tauy1 w2 taux2 tauy2 w3 taux3 tauy3 ] 
%
%          o 3
%         .  .
%        .     .
%       .        .
%      .           .
%     .              .
%  1 o.................o 2
%       ^ Y
%       |
%       o-----> X
%       Z up (towards you)
%