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)

bemimpac
function bemimpac      
% Example: bemimpac
% ~~~~~~~~~~~~~~~~~
% This program analyzes an impact dynamics 
% problem for an elastic Euler beam of 
% constant cross section which is simply 
% supported at each end. The beam is initially 
% at rest when a harmonically varying moment 
% m0*cos(w0*t) is applied to the right end. 
% The resulting transverse displacement and 
% bending moment are computed.  The 
% displacement and moment are plotted as 
% functions of x for the three time values. 
% Animated plots of the entire displacement
% and moment history are also given.
%
% User m functions required: 
%    beamresp, beamanim, sumser, ndbemrsp

fprintf('\nDYNAMICS OF A BEAM WITH AN ');
fprintf('OSCILLATING END MOMENT\n');
ei=1; arho=1; len=1; m0=1; w0=.90*pi^2;
tmin=0; tmax=5; nt=101;
xmin=0; xmax=len; nx=151; ntrms=200;
[t,x,displ,mom]=beamresp(ei,arho,len,m0,w0,...
             tmin,tmax,nt,xmin,xmax,nx,ntrms);
disp(' ')
disp('Press [Enter] to see the deflection') 
disp('for three positions'), pause

np=[3 5 8]; clf; pltsave=1;
dip=displ(np,:); mop=mom(np,:);
plot(x,dip(1,:),'-k',x,dip(2,:),':k',...
     x,dip(3,:),'--k'); 
xlabel('x axis'); ylabel('displacement');
hh=gca; 
r(1:2)=get(hh,'XLim'); r(3:4)=get(hh,'YLim'); 
xp=r(1)+(r(2)-r(1))/10;
dp=r(4)-(r(4)-r(3))/10;
tstr=['Displacement for Nearly Resonant' ...
      ' Moment Acting at Right End'];
title(tstr);
text(xp,dp,['Number of series terms ' ...
            'used = ',int2str(ntrms)]);
legend('t=0.10','t=0.20','t=0.35',3)        
disp(' ')
disp('Press [Enter] to the bending moment') 
disp('for three positions')
shg; pause  
if pltsave, print -deps 3positns, end 

clf; 
plot(x,mop(1,:),'-k',x,mop(2,:),':k',...
     x,mop(3,:),'--k'); 
h=gca; 
r(1:2)=get(h,'XLim'); r(3:4)=get(h,'YLim'); 
mp=r(3)+(r(4)-r(3))/10;
xlabel('x axis'); ylabel('moment');
tstr=['Bending Moment for Nearly Resonant' ...
      ' Moment Acting at Right End'];
title(tstr);
text(xp,mp,['Number of series terms ' ...
            'used = ',int2str(ntrms)]);
legend('t=0.10','t=0.20','t=0.35',2),
disp(' '), disp(...
'Press [Enter] to see the deflections surface') 
shg, pause
if pltsave, print -deps 3moments, end

inct=2; incx=2; 
ht=0.75; it=1:inct:.8*nt; ix=1:incx:nx;
tt=t(it); xx=x(ix); 
dd=displ(it,ix); mm=mom(it,ix);
a=surf(xx,tt,dd);
tstr=['Transverse Deflection as a ' ...
      'Function of Time and Position'];
title(tstr);
xlabel('x axis'); ylabel('time'); 
zlabel('transverse deflection');
disp(' '), disp(['Press [Enter] to ',...
'see the bending moment surface']) 
colormap(gray), brighten(.7)
shg, pause
if pltsave, print -deps dispsrf, end

a=surf(xx,tt,mm);
title(['Bending Moment as a Function ' ...
       'of Time and Position'])
xlabel('x axis'); ylabel('time'); 
zlabel('bending moment'); disp(' ')
disp('Press [Enter] to see animation of');
disp('the beam deflection')
colormap(gray), brighten(.7)
shg, pause 
if pltsave, print -deps momsrf, end
beamanim(x,displ,.1,'Transverse Deflection', ...
        'x axis','deflection'), disp(' ') 
disp('Press [Enter] to see animation');
disp('of the bending moment'); pause
beamanim(x,mom,.1,'Bending Moment History', ...
        'x axis','moment'); 
fprintf('\nAll Done\n'); close;

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

function [t,x,displ,mom]= ...
      beamresp(ei,arho,len,m0,w0,tmin,tmax, ...
               nt,xmin,xmax,nx,ntrms)
%
% [t,x,displ,mom]=beamresp(ei,arho,len,m0, ...
%           w0,tmin,tmax,nt,xmin,xmax,nx,ntrms)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function evaluates the time dependent 
% displacement and moment in a constant 
% cross section, simply supported beam which 
% is initially at rest when a harmonically 
% varying moment is suddenly applied at the 
% right end.  The resulting time histories of 
% displacement and moment are computed.
%
% ei        - modulus of elasticity times 
%             moment of inertia
% arho      - mass per unit length of the 
%             beam
% len       - beam length
% m0,w0     - amplitude and frequency of the 
%             harmonically varying right end 
%             moment
% tmin,tmax - minimum and maximum times for
%             the solution
% nt        - number of evenly spaced 
%             solution times
% xmin,xmax - minimum and maximum position 
%             coordinates for the solution. 
%             These values should lie between 
%             zero and len (x=0 and x=len at 
%             the left and right ends).
% nx        - number of evenly spaced solution 
%             positions
% ntrms     - number of terms used in the 
%             Fourier sine series
% t         - vector of nt equally spaced time 
%             values varying from tmin to tmax
% x         - vector of nx equally spaced 
%             position values varying from 
%             xmin to xmax
% displ     - matrix of transverse 
%             displacements with time varying 
%             from row to row, and position 
%             varying from column to column
% mom       - matrix of bending moments with 
%             time varying from row to row, 
%             and position varying from column 
%             to column
%
% User m functions called:  ndbemrsp
%----------------------------------------------

tcof=sqrt(arho/ei)*len^2; dcof=m0*len^2/ei;
tmin=tmin/tcof; tmax=tmax/tcof; w=w0*tcof;
xmin=xmin/len; xmax=xmax/len;
[t,x,displ,mom]=...
ndbemrsp(w,tmin,tmax,nt,xmin,xmax,nx,ntrms);
t=t*tcof; x=x*len; 
displ=displ*dcof; mom=mom*m0;

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

function beamanim(x,u,tpause,titl,xlabl,ylabl)
%
% beamanim(x,u,tpause,titl,xlabl,ylabl,save)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function draws an animated plot of data 
% values stored in array u.  The different 
% columns of u correspond to position values 
% in vector x.  The successive rows of u 
% correspond to different times. Parameter 
% tpause controls the speed of animation.
%
%  u      - matrix of values to animate plots 
%           of u versus x
%  x      - spatial positions for different 
%           columns of u
%  tpause - clock seconds between output of 
%           frames. The default is .1 secs 
%           when tpause is left out. When 
%           tpause=0, a new frame appears 
%           when the user presses any key.
%  titl   - graph title
%  xlabl  - label for horizontal axis
%  ylabl  - label for vertical axis
%
% User m functions called:  none
%----------------------------------------------
 
if nargin<6, ylabl=''; end; 
if nargin<5, xlabl=''; end
if nargin<4, titl=''; end; 
if nargin<3, tpause=.1; end;

[ntime,nxpts]=size(u); 
umin=min(u(:)); umax=max(u(:));
udif=umax-umin; uavg=.5*(umin+umax); 
xmin=min(x); xmax=max(x); 
xdif=xmax-xmin; xavg=.5*(xmin+xmax);
xwmin=xavg-.55*xdif; xwmax=xavg+.55*xdif;
uwmin=uavg-.55*udif; uwmax=uavg+.55*udif; clf;
axis([xwmin,xwmax,uwmin,uwmax]); title(titl);
xlabel(xlabl); ylabel(ylabl); hold on;

for j=1:ntime
  ut=u(j,:); 
  plot(x,ut,'-k'); axis('off'); figure(gcf);
  if tpause==0 
    pause; 
  else 
    pause(tpause); 
  end
  if j==ntime, break, else, cla; end
end
% print -deps cntltrac
hold off; clf;

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

function [u,t,x] = sumser(a,b,c,funt,funx, ...
                   tmin,tmax,nt,xmin,xmax,nx)
%
% [u,t,x] = sumser(a,b,c,funt,funx,tmin, ...
%                  tmax,nt,xmin,xmax,nx)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function evaluates a function U(t,x) 
% which is defined by a finite series. The 
% series is evaluated for t and x values taken 
% on a rectangular grid network. The matrix u 
% has elements specified by the following 
% series summation:
%
% u(i,j)  =   sum(    a(k)*funt(t(i)*b(k))*...
%           k=1:nsum
%                          funx(c(k)*x(j))
%
% where nsum is the length of each of the 
% vectors a, b, and c.
%
% a,b,c        - vectors of coefficients in 
%                the series
% funt,funx    - handles of functions accepting 
%                matrix argument.  funt is 
%                evaluated for an argument of 
%                the form funt(t*b) where t is 
%                a column and b is a row. funx
%                is evaluated for an argument 
%                of the form funx(c*x) where 
%                c is a column and x is a row.
% tmin,tmax,nt - produces vector t with nt 
%                evenly spaced values between 
%                tmin and tmax
% xmin,xmax,nx - produces vector x with nx 
%                evenly spaced values between 
%                xmin and xmax
% u            - the nt by nx matrix 
%                containing values of the 
%                series evaluated at t(i),x(j),
%                for i=1:nt and j=1:nx
% t,x          - column vectors containing t 
%                and x values. These output 
%                values are optional.
%
% User m functions called:  none.
%----------------------------------------------

tt=(tmin:(tmax-tmin)/(nt-1):tmax)';
xx=(xmin:(xmax-xmin)/(nx-1):xmax); a=a(:).';
u=a(ones(nt,1),:).*feval(funt,tt*b(:).')*...
  feval(funx,c(:)*xx);
if nargout>1, t=tt; x=xx'; end

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

function [t,x,displ,mom]= ...
    ndbemrsp(w,tmin,tmax,nt,xmin,xmax,nx,ntrms)
%
% [t,x,displ,mom]=ndbemrsp(w,tmin,tmax,nt,...
%                         xmin,xmax,nx,ntrms)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function evaluates the nondimensional 
% displacement and moment in a constant 
% cross section, simply supported beam which 
% is initially at rest when a harmonically 
% varying moment of frequency w is suddenly 
% applied at the right end. The resulting 
% time history is computed.
%
% w          - frequency of the harmonically 
%              varying end moment
% tmin,tmax  - minimum and maximum 
%              dimensionless times
% nt         - number of evenly spaced 
%              solution times
% xmin,xmax  - minimum and maximum 
%              dimensionless position 
%              coordinates. These values 
%              should lie between zero and 
%              one (x=0 and x=1 give the 
%              left and right ends).
% nx         - number of evenly spaced 
%              solution positions
% ntrms      - number of terms used in the 
%              Fourier sine series
% t          - vector of nt equally spaced 
%              time values varying from 
%              tmin to tmax
% x          - vector of nx equally spaced 
%              position values varying 
%              from xmin to xmax
% displ      - matrix of dimensionless 
%              displacements with time 
%              varying from row to row, 
%              and position varying from 
%              column to column
% mom        - matrix of dimensionless 
%              bending moments with time 
%              varying from row to row, and
%              position varying from column 
%              to column
%
% User m functions called:  sumser
%----------------------------------------------

if nargin < 8, w=0; end; nft=512; nh=nft/2;
xft=1/nh*(0:nh)'; 
x=xmin+(xmax-xmin)/(nx-1)*(0:nx-1)';
t=tmin+(tmax-tmin)/(nt-1)*(0:nt-1)'; 
cwt=cos(w*t);

% Get particular solution for nonhomogeneous 
% end condition
if w ==0 % Case for a constant end moment
  cp=[1 0 0 0; 0 0 2 0; 1 1 1 1; 0 0 2 6]\ ...
     [0;0;0;1];
  yp=[ones(size(x)), x, x.^2, x.^3]*cp; yp=yp';
  mp=[zeros(nx,2), 2*ones(nx,1), 6*x]*cp; 
  mp=mp';
  ypft=[ones(size(xft)), xft, xft.^2, xft.^3]*cp;

% Case where end moment oscillates 
% with frequency w
else  
  s=sqrt(w)*[1, i, -1, -i]; es=exp(s);
  cp=[ones(1,4); s.^2; es; es.*s.^2]\ ...
     [0; 0; 0; 1];
  yp=real(exp(x*s)*cp); yp=yp';
  mp=real(exp(x*s)*(cp.*s(:).^2)); mp=mp';
  ypft=real(exp(xft*s)*cp);
end

% Fourier coefficients for 
% particular solution
yft=-fft([ypft;-ypft(nh:-1:2)])/nft;

% Sine series coefficients for 
% homogeneous solution
acof=-2*imag(yft(2:ntrms+1));
ccof=pi*(1:ntrms)'; bcof=ccof.^2;

% Sum series to evaluate Fourier 
% series part of solution. Then combine 
% with the particular solution.
displ=sumser(acof,bcof,ccof,@cos,@sin,...
                  tmin,tmax,nt,xmin,xmax,nx);
displ=displ+cwt*yp; acof=acof.*bcof;
mom=sumser(acof,bcof,ccof,'cos','sin',...
                 tmin,tmax,nt,xmin,xmax,nx);
mom=-mom+cwt*mp;

Contact us