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


% modal_afm.m
%
% by Jack Chessa
% Northwestern University
%
clear
colordef black
state = 0;

%******************************************************************************
%***                             I N P U T                                  ***
%******************************************************************************
tic;
disp('************************************************')
disp('***          S T A R T I N G    R U N        ***')
disp('************************************************')

disp([num2str(toc),'   START'])

% MATERIAL PROPERTIES
E0  = 160;      % Youngs modulus in GPa
nu0 = 0.27;     % Poisson ratio
rho = 2.330e-9; % density in 10e12 Kg/m^3

% MESH PARAMETERS
quadType='GAUSS';
quadOrder=2;

% GMSH PARAMETERS
fileName='afm.msh';
domainID=50;
fixedID=51;
topID=52;

% EIGENPROBELM SOLUTION PARAMETERS
numberOfModes=8;   % number of modes to compute
consistentMass=0;  % use a consistent mass matrix
fixedBC=1;         % use fixed or free bcs

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ GMSH FILE
disp([num2str(toc),'   READING GMSH FILE'])
[node,elements,elemType]=msh2mlab(fileName); 
[node,elements]=remove_free_nodes(node,elements);
element=elements{domainID};
element=brickcheck(node,element,1);
if ( fixedBC )
  fixedEdge=elements{fixedID};
else
  fixedEdge=[];
end
topSurface=elements{topID};

plot_mesh(node,element,elemType{domainID},'r-')
disp([num2str(toc),'   INITIALIZING DATA STRUCTURES'])
numnode=size(node,1);    % number of nodes
numelem=size(element,1); % number of elements

% GET NODES ON DISPLACEMENT BOUNDARY
fixedNodeX=unique(fixedEdge);
fixedNodeY=fixedNodeX;
fixedNodeZ=fixedNodeX;

uFixed=zeros(size(fixedNodeX));  % displacement for fixed nodes
vFixed=zeros(size(fixedNodeY));
wFixed=zeros(size(fixedNodeZ));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPUTE COMPLIANCE MATRIX
C=zeros(6,6);
C(1:3,1:3)=E0/(1+nu0)/(1-2*nu0)*[ 1-nu0      nu0      nu0; 
                                    nu0    1-nu0      nu0; 
                                    nu0      nu0    1-nu0 ];
C(4:6,4:6)=E0/(1+nu0)*eye(3);                                    
             
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINE SYSTEM DATA STRUCTURES
K=sparse(3*numnode,3*numnode); % stiffness matrix
if ( consistentMass)
  M=sparse(3*numnode,3*numnode); % mass matrix
else
  M=zeros(3*numnode,1); % mass vector
end  
 
%******************************************************************************
%***                          P R O C E S S I N G                           ***
%******************************************************************************

%%%%%%%%%%%%%%%%%%%%% COMPUTE SYSTEM MATRICIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp([num2str(toc),'   COMPUTING STIFFNESS AND MASS MATRIX'])
[W,Q]=quadrature(quadOrder,quadType,3); % define quadrature rule
et=elemType{domainID};
nn=size(element,2);
for e=1:numelem                         % start of element loop
  
  sctr=element(e,:);           % element scatter vector
  sctrB0=[ sctr sctr+numnode sctr+2*numnode ]; % vector that scatters a B matrix
    
  for q=1:size(W,1)                        % quadrature loop
    pt=Q(q,:);                             % quadrature point
    wt=W(q);                               % quadrature weight
    [N,dNdxi]=lagrange_basis(et,pt);     % element shape functions
    
    J0=node(sctr,:)'*dNdxi;                % element Jacobian matrix        
    invJ0=inv(J0);
    dNdx=dNdxi*invJ0;
    detJ0=det(J0);
    
    if (detJ0 <= 0)
      disp(['ERROR: NEGATIVE JACOBIAN IN ELEMENT ',num2str(e)]);
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % COMPUTE B MATRIX
    B0=zeros(6,3*nn);
    B0(1,1:nn)        = dNdx(:,1)';
    B0(2,nn+1:2*nn)   = dNdx(:,2)';
    B0(3,2*nn+1:3*nn) = dNdx(:,3)';
    
    B0(4,2*nn+1:3*nn) = dNdx(:,2)';
    B0(4,nn+1:2*nn)   = dNdx(:,3)';
    
    B0(5,1:nn)        = dNdx(:,3)';
    B0(5,2*nn+1:3*nn) = dNdx(:,1)';
    
    B0(6,nn+1:2*nn)   = dNdx(:,1)';
    B0(6,1:nn)        = dNdx(:,2)';
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % COMPUTE ELEMENT STIFFNESS AT QUADRATURE POINT
    K(sctrB0,sctrB0)=K(sctrB0,sctrB0)+B0'*C*B0*wt*detJ0;
    
    % COMPUTE ELEMENT MASS AT QUADRATURE POINT
    mQPt=N*rho*N'*wt*detJ0;
    if ( ~consistentMass )
      mQPt=sum(mQPt)';
      M(sctr)           = M(sctr)+mQPt;
      M(sctr+numnode)   = M(sctr+numnode)+mQPt;
      M(sctr+2*numnode) = M(sctr+2*numnode)+mQPt;
    else
      M(sctr,sctr)                        = M(sctr,sctr)+mQPt;
      M(sctr+numnode,sctr+numnode)        = M(sctr+numnode,sctr+numnode)+mQPt;
      M(sctr+2*numnode,sctr+2*numnode)    = M(sctr+2*numnode,sctr+2*numnode)+mQPt;
    end 
  end                                      % of quadrature loop
end     % of element loop
%%%%%%%%%%%%%%%%%%% END OF SYSTEM MATRIX COMPUTATION %%%%%%%%%%%%%%%%%%%%%%

% ELIMINATE FIXED DOFS FROM EIGENVALUE COMUTATION
disp([num2str(toc),'   FINDING ACTIVE DOFS'])
activeDof=setdiff([1:numnode]',[fixedNodeX;fixedNodeY;fixedNodeZ]);
activeDof=[activeDof;activeDof+numnode;activeDof+2*numnode];

% SOLVE SYSTEM
disp([num2str(toc),'   SOLVING EIGEN PROBLEM'])
if ( consistentMass )
  [modeShape,freq]=eigs(K(activeDof,activeDof),M(activeDof,activeDof),...
                                numberOfModes,0);
else
  Minv=spdiags(1./M,0,3*numnode,3*numnode);
  K=Minv*K;
  [modeShape,freq]=eigs(K(activeDof,activeDof),numberOfModes,0);
end                  
freq=diag(freq)/(2*pi);  % frequency in kHz

%******************************************************************************
%***                     P O S T  -  P R O C E S S I N G                    ***
%******************************************************************************
disp([num2str(toc),'   POST-PROCESSING'])
disp(['THE MODE FREQUENCIES ARE:'])
for m=1:length(freq)
  disp(['   MODE: ',num2str(m),' ',num2str(freq(m))])
  
  % PLOT MODE SHAPE
  figure(m); clf;
  U=zeros(numnode,1);
  U(activeDof)=modeShape(:,m);
  scaleFactor=20/max(abs(U));
  plot_field(node+[U(1:numnode) U(numnode+1:2*numnode) U(2*numnode+1:3*numnode)]*scaleFactor,topSurface,elemType{topID},ones(3*numnode,1));  
  hold on
  plot_mesh(node+[U(1:numnode) U(numnode+1:2*numnode) U(2*numnode+1:3*numnode)]*scaleFactor,topSurface,elemType{topID},'k-'); 
  plot_mesh(node,topSurface,elemType{topID},'r-');
  title(['MODE ',num2str(m),', FREQUENCY = ',num2str(freq(m)),' [kHz]'])
  view(37,36) 
  axis off
  print(m, '-djpeg90', ['afm_mode_',num2str(m),'.jpg']);
  
end

% ANIMATE MODE
nCycles=5; % number of cycles to animate
fpc=10;    % frames per cycle
fact=sin(linspace(0,2*pi,fpc));
m=input('What mode would you like to animate (type 0 to exit) ');
while ( m~=0 )
  
  U=zeros(numnode,1);
  U(activeDof)=modeShape(:,m);
  wt=20/max(abs(U));
  for i=1:fpc
    scaleFactor=fact(i)*wt;
    figure(length(freq+1));
    clf;
    plot_field(node+[U(1:numnode) U(numnode+1:2*numnode) U(2*numnode+1:3*numnode)]*scaleFactor,topSurface,elemType{topID},ones(3*numnode,1));  
    hold on
    plot_mesh(node+[U(1:numnode) U(numnode+1:2*numnode) U(2*numnode+1:3*numnode)]*scaleFactor,topSurface,elemType{topID},'k-'); 
    plot_mesh(node,topSurface,elemType{topID},'w-');
    hold on 
    view(37,36) 
    axis([70 240 30 160 -10 10])
    title(['MODE ',num2str(m),', FREQUENCY = ',num2str(freq(m)),' [kHz]'])
    axis off
    film(i)=getframe;
  end
  
  movie(film,nCycles);
  
  m=input('What mode would you like to animate (type 0 to exit) ');
  if ( m > length(freq) )
    disp(['mode must be less than ',num2str(length(freq))])
  end
end

disp([num2str(toc),'   RUN FINISHED'])
% ***************************************************************************
% ***                     E N D   O F    P R O G R A M                    ***
% ***************************************************************************
disp('************************************************')
disp('***            E N D    O F    R U N         ***')
disp('************************************************')

% compute uexact