cst_csq_stress.m

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

function [SIGsys] = CST_CSQ_stress (in_data,resp)

% stress field for CST/CSQ plane FEM model
% ---------------------------------------------------------------------------
% INPUT:
%        in_data -  structure with input data: ND, EL, etc.
%        D       -  computed nodal displacements vector
% ---------------------------------------------------------------------------

fj = size(in_data.EL,1);
SIGsys = zeros(fj(1)*3,1);

for i=1:fj % going thru all elements

   if in_data.EL(i,2)==4  % "4" - CST element (triangle) ----------------------------
      node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
      node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
      node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
      
      if in_data.mater.E~=0   E1=in_data.mater.E;      else E1=in_data.EL(i,6);       end;
      if in_data.mater.h~=0   h_1=in_data.mater.h;     else h_1 = in_data.EL(i,7);    end;
      if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else  miu_1 = in_data.EL(i,8); end;
    
      [Bsys,Esys] = D2_CST (in_data.ND(node1,2),in_data.ND(node1,3),in_data.ND(node2,2),...
          in_data.ND(node2,3),in_data.ND(node3,2),in_data.ND(node3,3),E1,h_1,miu_1);
      
      Dlocal = [resp.static.D(node1*2-1) resp.static.D(node1*2) resp.static.D(node2*2-1) resp.static.D(node2*2) ...
                resp.static.D(node3*2-1) resp.static.D(node3*2)];
      SIGlocal = Esys*Bsys*Dlocal';
      SIGsys((i*3-2):(i*3)) = SIGlocal;    
   end;
   
   if in_data.EL(i,2)==5  % "5" - CSQ element (quadrilateral) -----------------------
      node1 = find(in_data.ND(:,1)==in_data.EL(i,3));
      node2 = find(in_data.ND(:,1)==in_data.EL(i,4));
      node3 = find(in_data.ND(:,1)==in_data.EL(i,5));
      node4 = find(in_data.ND(:,1)==in_data.EL(i,6));
      
      if in_data.mater.E~=0  Em=in_data.mater.E;   else Em=in_data.EL(i,7);   end;
      if in_data.mater.h~=0 
          if length(in_data.mater.h)==1   hT=[1 1 1 1].*in_data.mater.h; end;
          if length(in_data.mater.h)==4   hT=in_data.mater.h;            end;
      else hT = in_data.EL(i,8:11); end; % element thickness at 4 nodes
        
      if in_data.mater.miu~=0 miu_1=in_data.mater.miu; else  miu_1 = in_data.EL(i,12); end;
      
      [B_lc,E_lc] = D2_CSQ (in_data.ND(node1,2),in_data.ND(node1,3), in_data.ND(node2,2),...
          in_data.ND(node2,3),in_data.ND(node3,2),in_data.ND(node3,3),in_data.ND(node4,2),...
          in_data.ND(node4,3), Em, hT, miu_1);
      
      Dlocal = [resp.static.D(node1*2-1) resp.static.D(node1*2) resp.static.D(node2*2-1) resp.static.D(node2*2) ...
              resp.static.D(node3*2-1) resp.static.D(node3*2) resp.static.D(node4*2-1) resp.static.D(node4*2)];
      SIGlocal = E_lc*B_lc*Dlocal';
      SIGsys((i*3-2):(i*3)) = SIGlocal;
   end;

end;