No BSD License  

Highlights from
Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

image thumbnail

Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

by

 

14 Oct 2002 (Updated )

Companion Software (amamhlib)

trusvibs
function trusvibs      
% Example: trusvibs
% ~~~~~~~~~~~~~~~~~
%
% This program analyzes natural vibration modes 
% for a general plane pin-connected truss. The 
% direct stiffness method is employed in 
% conjunction with eigenvalue calculation to 
% evaluate the natural frequencies and mode 
% shapes. The truss is defined in terms of a 
% set of nodal coordinates and truss members 
% connected to different nodal points. Global 
% stiffness and mass matrices are formed. Then 
% the frequencies and mode shapes are computed 
% with provision for imposing zero deflection 
% at selected nodes. The user is then allowed 
% to observe animated motion of the various 
% vibration modes.
%
% User m functions called: 
%        eigsym, crossdat, drawtrus, eigc,
%        assemble, elmstf, cubrange

global x y inode jnode elast area rho idux iduy
kf=1; idux=[]; iduy=[]; close; disp(' ')
disp(['Modal Vibrations for a Pin ', ...
      'Connected Truss']); disp(' ');

% A sample data file defining a problem is 
% given in crossdat.m
disp(['Give the name of a function which ', ...
      'creates your input data']);
disp(['Do not include .m in the name ', ...
      '(use crossdat as an example)']);
filename=input('>? ','s');
eval(filename); disp(' ');

% Assemble the global stiffness and 
% mass matrices
[stiff,masmat]= ...
  assemble(x,y,inode,jnode,area,elast,rho);

% Compute natural frequencies and modal vectors 
% accounting for the fixed nodes
ifixed=[2*idux(:)-1; 2*iduy(:)];
[modvcs,eigval]=eigc(stiff,masmat,ifixed);
natfreqs=sqrt(eigval); 

% Set parameters used in modal animation
nsteps=31; s=sin(linspace(0,6.5*pi,nsteps));
x=x(:); y=y(:); np=2*length(x);
bigxy=max(abs([x;y])); scafac=.05*bigxy;
highmod=size(modvcs,2); hm=num2str(highmod); 

% Show animated plots of the vibration modes
while 1
  disp('Give the mode numbers to be animated?');
  disp(['Do not exceed a total of ',hm, ...
        ' modes.']); disp('Input 0 to stop');
  if kf==1, disp(['Try 1:',hm]); kf=kf+1; end  
  str=input('>? ','s'); 
  nmode=eval(['[',str,']']); 
  nmode=nmode(find(nmode<=highmod));
  if sum(nmode)==0; break; end
  % Animate the various vibration modes 
  hold off; clf; ovrsiz=1.1; 
  w=cubrange([x(:),y(:)],ovrsiz);
  axis(w); axis('square'); axis('off'); hold on;
  for kk=1:length(nmode)  % Loop over each mode
    kkn=nmode(kk); 
    titl=['Truss Vibration Mode Number ', ...
          num2str(kkn)];
    dd=modvcs(:,kkn); mdd=max(abs(dd));
    dx=dd(1:2:np); dy=dd(2:2:np); 
    clf; pause(1);
    % Loop through several cycles of motion 
    for jj=1:nsteps           
      sf=scafac*s(jj)/mdd; 
      xd=x+sf*dx; yd=y+sf*dy; clf;
      axis(w); axis('square'); axis('off');
      drawtrus(xd,yd,inode,jnode); title(titl);
      drawnow; figure(gcf);
    end
  end
end
disp(' '); hold off

%=============================================

function crossdat
% [inode,jnode,elast,area,rho]=crossdat
% This function creates data for the truss
% vibration program. It can serve as a model
% for other configurations by changing the
% function name and data quantities       
% Data set: crossdat
% ~~~~~~~~~~~~~~~~~~
%
% Data specifying a cross-shaped truss.
%
%----------------------------------------------

global x y inode jnode elast area rho idux iduy

% Nodal point data are defined by:
%   x - a vector of x coordinates
%   y - a vector of y coordinates
x=10*[.5 2.5 1 2 0 1 2 3 0 1 2 3 1 2];
y=10*[ 0   0 1 1 2 2 2 2 3 3 3 3 4 4];

% Element data are defined by:
%   inode - index vector defining the I-nodes
%   jnode - index vector defining the J-nodes
%   elast - vector of elastic modulus values
%   area  - vector of cross section area values
%   rho   - vector of mass per unit volume 
%           values
inode=[1 1 2 2 3 3 4 3 4 5 6 7 5 6 6 6 7 7 7 ...
       8 9 10 11 10 11 10 11 13];
jnode=[3 4 3 4 4 6 6 7 7 6 7 8 9 9 10 11 10 ...
       11 12 12 10 11 12 13 13 14 14 14];
elast=3e7*ones(1,28);
area=ones(1,28); rho=ones(1,28);

% Any points constrained against displacement 
% are defined by:
%   idux - indices of nodes having zero 
%          x-displacement
%   iduy - indices of nodes having zero 
%          y-displacement
idux=[1 2]; iduy=[1 2];

%=============================================

function drawtrus(x,y,i,j)
%
% drawtrus(x,y,i,j)
% ~~~~~~~~~~~~~~~~~
%
% This function draws a truss defined by nodal
% coordinates defined in x,y and member indices 
% defined in i,j.
%
% User m functions called: none
%----------------------------------------------

hold on;
for k=1:length(i)
  plot([x(i(k)),x(j(k))],[y(i(k)),y(j(k))]);
end

%=============================================

function [vecs,eigvals]=eigc(k,m,idzero)
%
% [vecs,eigvals]=eigc(k,m,idzero)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function computes eigenvalues and 
% eigenvectors for the problem 
%            k*x=eigval*m*x 
% with some components of x constrained to 
% equal zero. The imposed constraint is
%            x(idzero(j))=0 
% for each component identified by the index 
% matrix idzero.
%
% k       - a real symmetric stiffness matrix 
% m       - a positive definite symmetric mass 
%           matrix
% idzero  - the vector of indices identifying 
%           components to be made zero
%
% vecs    - eigenvectors for the constrained 
%           problem. If matrix k has dimension 
%           n by n and the length of idzero is 
%           m (with m<n), then vecs will be a 
%           set on n-m vectors in n space
% eigvals - eigenvalues for the constrained 
%           problem. These are all real.
%
% User m functions called:  eigsym
%----------------------------------------------

n=size(k,1); j=1:n; j(idzero)=[]; 
c=eye(n,n); c(j,:)=[];
[vecs,eigvals]=eigsym((k+k')/2, (m+m')/2, c);

%=============================================

function [evecs,eigvals]=eigsym(k,m,c)
%
% [evecs,eigvals]=eigsym(k,m,c)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function solves the constrained
% eigenvalue problem
%    k*x=(lambda)*m*x, with c*x=0.
% Matrix k must be real symmetric and matrix
% m must be symmetric and positive definite;
% otherwise, computed results will be wrong.
%
% k       - a real symmetric matrix
% m       - a real symmetric positive 
%           definite matrix
% c       - a matrix defining the constraint 
%           condition c*x=0. This matrix is
%           omitted if no constraint exists.
%
% evecs   - matrix of eigenvectors orthogonal
%           with respect to k and m. The
%           following relations apply:
%           evecs'*m*evecs=identity_matrix
%           evecs'*k*evecs=diag(eigvals).
% eigvals - a vector of the eigenvalues
%           sorted in increasing order 
%
% User m functions called: none
%----------------------------------------------

if nargin==3
  q=null(c); m=q'*m*q; k=q'*k*q;
end
u=chol(m); k=u'\k/u; k=(k+k')/2;
[evecs,eigvals]=eig(k);
[eigvals,j]=sort(diag(eigvals));
evecs=evecs(:,j); evecs=u\evecs;
if nargin==3, evecs=q*evecs; end

%=============================================

function [stif,masmat]= ...
  assemble(x,y,id,jd,a,e,rho)
%
% [stif,masmat]=assemble(x,y,id,jd,a,e,rho)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function assembles the global 
% stiffness matrix and mass matrix for a 
% plane truss structure. The mass density of 
% each element equals unity.
%
% x,y    - nodal coordinate vectors
% id,jd  - nodal indices of members
% a,e    - areas and elastic moduli of members
% rho    - mass per unit volume of members
%
% stif   - global stiffness matrix
% masmat - global mass matrix
%
% User m functions called: elmstf
%----------------------------------------------

numnod=length(x); numelm=length(a); 
id=id(:); jd=jd(:);
stif=zeros(2*numnod); masmat=stif;
ij=[2*id-1,2*id,2*jd-1,2*jd];
for k=1:numelm, kk=ij(k,:);
  [stfk,volmk]= ...
    elmstf(x,y,a(k),e(k),id(k),jd(k));
  stif(kk,kk)=stif(kk,kk)+stfk;  
  masmat(kk,kk)=masmat(kk,kk)+ ...
                rho(k)*volmk/2*eye(4,4); 
end

%=============================================

function [k,vol]=elmstf(x,y,a,e,i,j)
%
% [k,vol]=elmstf(x,y,a,e,i,j)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function forms the stiffness matrix for 
% a truss element. The member volume is also 
% obtained.
%
% User m functions called: none
%----------------------------------------------

xx=x(j)-x(i); yy=y(j)-y(i); 
L=norm([xx,yy]); vol=a*L;
c=xx/L; s=yy/L; k=a*e/L*[-c;-s;c;s]*[-c,-s,c,s];

%=============================================

% function range=cubrange(xyz,ovrsiz)
% See Appendix B

Contact us