www.pudn.com > 87361026FEA.rar > notch.m, change:2006-02-28,size:7602b


% notch.m
%
% by Jack Chessa
% Northwestern University
%
% A simple 2D explicit nonlinear finite element code (total Lagranian)
% 
% Simulation of an impact on a notched plate assuming hyper-elastic material.
clear
colordef black
state = 0;

%******************************************************************************
%***                             I N P U T                                  ***
%******************************************************************************

% HERE WE DEFINE SOME SYSTEM PARPMETERS

% GMSH FILE PARPMETERS
gmshFileName = 'notch.msh';  % name of gmsh input file
domainID     = 11;           % zone id for domain mesh
impactID     = 12;           % zone id for displacement boundary mesh
fixedID      = 13;            % zone id for traction boundary mesh

% MATERIAL PROPERTIES
rho   = 10;    % density
young = 10000; % youngs modulus
nu    = 0.3;   % poisson ratio

% BOUNDARY CONDITIONS
vImpact=0.25;

% TIME STEPPING PARAMETERS
tstart = 0;       % initial time
tfin = 0.20;      % final time
Cr = 0.75;        % C.F.L. number
hmin=0.1;         % minimal element size factor
deltInit = Cr*hmin/sqrt(young/rho);  % initial time step

%******************************************************************************
%***                    P R E - P R O C E S S I N G                         ***
%******************************************************************************

% READ GMSH FILE
[node,elements,elemType]=msh2mlab(gmshFileName);

% CHECK FOR ONLY TRIANGULAR ELEMENTS

% GET NODE AND CONNECTIVITY MATRICIES
node=node(:,1:2);                % a nodal coordinate matix (reference config)
element  = elements{domainID};   % connectivity matrix for interior
impactBndy = elements{impactID}; % connectivity matrix for impact boundary
fixedBndy = elements{fixedID};   % connectivity matrix for impact boundary

% CHECK FOR BAD JACOBIANS
tricheck(node,element);

% GET NODES ON DISPLACEMENT BOUNDARY
impactBndyNodes=unique(impactBndy);
fixedBndyNodes=unique(fixedBndy);

numnode=size(node,1);    % number of nodes
numelem=size(element,1); % number of elements

% DEFINE SYSTEM DATA STRUCTURES
u=zeros(2*numnode,1);        % nodal displacement vector
v=zeros(2*numnode,1);        % nodal velocity vector
X=node;                      % nodal coordinates in current configureation
f=zeros(2*numnode,1);        % load vector  force external - force internal
fext=zeros(2*numnode,1);     % external load vector
elementStress=cell(numelem); % element stresses
M=zeros(2*numnode,1);        % diagonalized mass matrix (vector)
xs=1:numnode;                % x portion of u and v vectors
ys=(numnode+1):2*numnode;    % y portion of u and v vectors

% DEFINE INITIAL CONDITIONS
v(:)=0;
u(:)=0;

% PLOT MESH
trimesh(element,node(:,1),node(:,2),0*node(:,1))
view(2)
axis equal
disp('paused')
pause

%******************************************************************************
%***                          P R O C E S S I N G                           ***
%******************************************************************************

%%%%%%%%%%%%%%%%%%% CALCULATE MASS MATRIX
Mtmp=sparse(numnode,numnode);
[W,Q]=quadrature( 2, 'TRIANGULAR', 2 );  % define quadrature rule
for e=1:numelem                          % start of element loop
  
  sctr=element(e,:);
  
  for q=1:size(W,1)                      % quadrature loop
    pt=Q(q,:);                           % quadrature point
    wt=W(q);                             % quadrature weight
    [N,dNdxi]=lagrange_basis('T3',pt);   % element shape functions
    J0=X(sctr,:)'*dNdxi;  % element Jacobian matrix w.r.t ref configuration 
    Mtmp(sctr,sctr)=Mtmp(sctr,sctr)+N*N'*det(J0)*wt;
    
  end                                    % of quadrature loop
end   % of element loop
% lump mass matrix
Mtmp=sum(Mtmp')';
M=[Mtmp;Mtmp];
clear Mtmp;

%%%%%%%%%%%%%%%%%%% CALCULATE COMPLIENCE MATRIX
mu=young/2/(1+nu);
lambda=nu*young/(1+nu)/(1-2*nu);
C=[lambda+2*mu      lambda     0;
        lambda lambda+2*mu     0;
             0           0     mu];
ind2voit=[1;4;2];  % the are used to switch between indicie and voigt form
voit2ind=[1 3;3 2];

% ***************************** TEMPORAL LOOP *********************************
tn=tstart;
delt=deltInit;
step=0;

while ( tn <= tfin ) 
  tn=tn+delt;
  step=step+1;
  
  % ADD EXTERNAL FORCE 
  f=fext;
  
  %%%%%%%%%%%%%%%%%%%%% COMPUTE INTERNAL FORCE FOR TIME n-1
  [W,Q]=quadrature( 1, 'TRIANGULAR', 2 );  % define quadrature rule
  for e=1:numelem                          % start of element loop
    
    sctr=element(e,:);
    sctrx=sctr;
    sctry=sctr+numnode;
    sctrB0=[sctr(1) sctr(1)+numnode sctr(2) ...
                  sctr(2)+numnode sctr(3) sctr(3)+numnode ];
    
    for q=1:size(W,1)                      % quadrature loop
      pt=Q(q,:);                           % quadrature point
      wt=W(q);                             % quadrature weight
      [N,dNdxi]=lagrange_basis('T3',pt);   % element shape functions
      
      J0=X(sctr,:)'*dNdxi;  % element Jacobian matrix w.r.t ref configuration       
      invJ0=inv(J0);
      dNdX=dNdxi*invJ0;
      B0=[dNdX(1,1)         0 dNdX(2,1)         0 dNdX(3,1)         0;
                  0 dNdX(1,2)         0 dNdX(2,2)         0 dNdX(3,2);
          dNdX(1,2) dNdX(1,1) dNdX(2,2) dNdX(2,1) dNdX(3,2) dNdX(3,1)];
      
      H=[u(sctrx)';u(sctry)']*dNdX;       % definitions from TB's book
      F=eye(2)+H;                         % Box 4.6
      E=1/2*(H+H'+H'*H);
      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % COMPUTE ELEMENT STRESS (add your constitiutive law here)
      S = C*E(ind2voit);
      elementStress{e}=S;
      P=S(voit2ind)*F';  
      
      % subtact internal force from force vector
      finte = -B0'*P(ind2voit)*det(J0)*wt;
      f(sctrB0) = f(sctrB0) + finte;
      
    end                                    % of quadrature loop
  end   % of element loop
  
  % UPDATE NODAL VELOCITIES
  v = v + delt*(f./M);
  
  % APPLY ESSENTIAL BCs TO NODAL VELOCITIES
  v(fixedBndyNodes+numnode)=0.0;     % set y-velocity to zero
  v(fixedBndyNodes)=0.0;             % set x-velocity to zero
  v(impactBndyNodes)=vImpact;        % set x-velocity to impact velocity
  v(impactBndyNodes+numnode)=0.0;    % set y-velocity to impact velocity (0.0)
  
  % UPDATE NODAL DISPLACEMENTS
  u = u + delt*v;

  % POST-PROCESS STRESS VALUES ( get nodal average values )
  count=zeros(numnode,1);
  Sxx=zeros(numnode,1);
  Syy=zeros(numnode,1);
  Sxy=zeros(numnode,1);
  for e=1:numelem
    sctr=element(e,:);
    Sxx(sctr)=Sxx(sctr)+elementStress{e}(1);
    Syy(sctr)=Syy(sctr)+elementStress{e}(2);
    Sxy(sctr)=Sxy(sctr)+elementStress{e}(3);
    count(sctr)=count(sctr)+1;
  end
  Sxx=Sxx./count;
  Syy=Syy./count;
  Sxy=Sxy./count;
  S1=(Sxx+Syy)/2+sqrt( ((Sxx-Syy)/2).^2 - Sxy.^2 );
  S2=(Sxx+Syy)/2-sqrt( ((Sxx-Syy)/2).^2 - Sxy.^2 );
  Svm=sqrt(S1.^2+S2.^2);

  % PLOT RESULTS
  hold off
  trisurf(element,X(:,1),X(:,2),0*X(:,1),Svm)
  axis equal
  axis([-12 10 -7 7])
  caxis([0 80])
  shading interp
  view(2)
  pause(0.1)
  film(step)=getframe;
  
  if ( step==1 )
    colorbar
  end
  
  % SAVE RESULTS

  
end % of time loop
% *****************************************************************************
state=1;

% play back the movie
movie(film)

% break
% %%%%%%%%%%%%%%%%%%%%%%%%
% % patch test
% clear
% e=1;
% X=[0 0;
%    1 0;
%    1 1;
%    0 1];
%  element=[1 2 4;3 4 2];
%  numnode=4; numelem=2;
%  u=[0 0 0 0 0 0 .01 .01]';
%  f=[0 0 0 0 0 0 0 0]';
%  C=[2 1 0;
%     1 2 0;
%     0 0 0.5];
% ind2voit=[1;4;2];  
% voit2ind=[1 3;3 2];