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

Howard Wilson

 

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