k_cst.m

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

function [Kgl, M] = K_CST (ND, EL, rhoX, rhoY)

% assemble global stiffness & mass matrix for CST FEM plane model -----------
% INPUT:
%        ND   -  nodes matrix
%        dof  -  # of system dof's 
%        EL   -  elements matrix
%        rhoX -  density of the element: X direction
%        rhoY -  density of the element: Y direction
% ---------------------------------------------------------------------------
dof = size(ND)*2;
Kgl = zeros(dof(1),dof(1)); % global stiffness matrix: dof = x y xy
M   = zeros(dof(1),dof(1));

for i=1:size(EL)
   node1 = find(ND(:,1)==EL(i,3));
   node2 = find(ND(:,1)==EL(i,4));
   node3 = find(ND(:,1)==EL(i,5));
   E1    = EL(i,6);
   h_1   = EL(i,7);
   miu_1 = EL(i,8);
   if EL(i,10)==4
    [Klc,Bsys,Esys] = D2_CST (ND(node1,2),ND(node1,3),ND(node2,2),ND(node2,3),...
         ND(node3,2),ND(node3,3),E1,h_1,miu_1);
   end;
   t1 = node1*2;
   t2 = node2*2;
   t3 = node3*2;
   x1 = ND(node1,2);
   y1 = ND(node1,3);
   x2 = ND(node2,2);
   y2 = ND(node2,3);
   x3 = ND(node3,2);
   y3 = ND(node3,3);
   A = abs(0.5 * det([1 1 1; x1 x2 x3; y1 y2 y3]) ); % triangle area
   fM = A*EL(i,6)/3.*[rhoX rhoY rhoX rhoY rhoX rhoY]';
   Kgl((t1-1):t1,(t1-1):t1) = Kgl( (t1-1):t1,(t1-1):t1) + Klc(1:2,1:2);
   Kgl((t1-1):t1,(t2-1):t2) = Kgl( (t1-1):t1,(t2-1):t2) + Klc(1:2,3:4);
   Kgl((t1-1):t1,(t3-1):t3) = Kgl( (t1-1):t1,(t3-1):t3) + Klc(1:2,5:6);
   Kgl((t2-1):t2,(t1-1):t1) = Kgl( (t2-1):t2,(t1-1):t1) + Klc(3:4,1:2);
   Kgl((t2-1):t2,(t2-1):t2) = Kgl( (t2-1):t2,(t2-1):t2) + Klc(3:4,3:4);
   Kgl((t2-1):t2,(t3-1):t3) = Kgl( (t2-1):t2,(t3-1):t3) + Klc(3:4,5:6);
   Kgl((t3-1):t3,(t1-1):t1) = Kgl( (t3-1):t3,(t1-1):t1) + Klc(5:6,1:2);
   Kgl((t3-1):t3,(t2-1):t2) = Kgl( (t3-1):t3,(t2-1):t2) + Klc(5:6,3:4);
   Kgl((t3-1):t3,(t3-1):t3) = Kgl( (t3-1):t3,(t3-1):t3) + Klc(5:6,5:6);
   M(t1-1,t1-1) = M(t1-1,t1-1) + fM(1);  
   M(t1,t1)     = M(t1,t1) + fM(2); 
   M(t2-1,t2-1) = M(t2-1,t2-1) + fM(3); 
   M(t2,t2)     = M(t2,t2) + fM(4);
   M(t3-1,t3-1) = M(t3-1,t3-1) + fM(5);
   M(t3,t3)     = M(t3,t3) + fM(6);
end;