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)

vdb
function vdb        
% Example: vdb
% ~~~~~~~~~~~~
%
% This program calculates the shear, moment, 
% slope, and deflection of a variable depth 
% indeterminate beam subjected to complex 
% loading and general end conditions. The 
% input data are defined in the program 
% statements below.
%
% User m functions required:
%   bmvardep, extload, lintrp, oneovrei, 
%   sngf, trapsum

clear all; Problem=1;
if Problem == 1
  Title=['Problem from Arbabi and Li'];
  Printout=10;  % Output frequency
  BeamLength=3; % Beam length
  NoSegs=301;   % # of beam divisions for 
                % integration
  % External concentrated loads and location
  ExtForce= [-1]; ExtForceX=[.5];
  % External ramp loads and range
  %        q1  q2  x1  x2
  ExtRamp=[-1  -1   1   2];
  % Interior supports: initial displacement 
  % and location
  IntSupX=    [1; 2]; IntSupDelta=[0; 0];
  % End (left and right) conditions
  EndCondVal= [0; 0; 0; 0];   % magnitude
  % 1=shear,2=moment,3=slope,4=delta
  EndCondFunc=[3; 4; 3; 4];   
  % 1=left end,2=right end
  EndCondEnd= [1; 1; 2; 2];   
  % EI or beam depth specification
  EIorDepth=1;  % 1=EI values specified 
                % 2=depth values specified
  if EIorDepth == 1
    % Discretize the parabolic haunch for the 
    % three spans
    Width=1; E=1; a=0.5^2; Npts=100;
    h1=0.5; k1=1; x1=linspace(0,1,Npts);
    h2=1.5; k2=1; x2=linspace(1,2,Npts);
    h3=2.5; k3=1; x3=linspace(2,3,Npts);
    y1=(x1-h1).^2/a+k1; y2=(x2-h2).^2/a+k2;
    y3=(x3-h3).^2/a+k3;
    EIx=[x1 x2 x3]'; h=[y1 y2 y3]';
    EIvalue=E*Width/12*h.^3;
    mn=min(EIvalue); EIvalue=EIvalue./mn;
  else
    % Beam width and Young's modulus
    BeamWidth=[]; BeamE=[]; Depth=[]; DepthX=[];
  end
elseif Problem == 2
  Title=['From Timoshenko and Young,', ...
         ' p 434, haunch beam'];
  Printout=12; NoSegs=144*4+1; BeamLength=144; 
  ExtForce=[]; ExtForceX=[]; 
  ExtRamp=[-1 -1 0 108];
  IntSupX=[36; 108]; IntSupDelta=[0; 0];
  EndCondVal=[0; 0; 0; 0]; 
  EndCondFunc=[2; 4; 2; 4];
  EndCondEnd= [1; 1; 2; 2]; EIorDepth=2;  
  if EIorDepth == 1
    EIvalue=[]; EIx=[];
  else
    BeamWidth=[1]; BeamE=[1];
    % Discretize the parabolic sections
    a=36^2/5; k=2.5; h1=0; h2=72; h3=144; 
    N1=36; N2=72; N3=36;
    x1=linspace(  0, 36,N1); y1=(x1-h1).^2/a+k;
    x2=linspace( 36,108,N2); y2=(x2-h2).^2/a+k;
    x3=linspace(108,144,N3); y3=(x3-h3).^2/a+k;
    Depth=[y1 y2 y3]'; DepthX=[x1 x2 x3]';
    % Comparison values
    I=BeamWidth*Depth.^3/12; Imin=min(I); L1=36;
    k1=BeamE*Imin/L1; k2=k1/2; k3=k1;
    t0=10.46/k1; t1=15.33/k1; t2=22.24/k1; 
    t3=27.95/k1;
    fprintf('\n\nValues from reference');
    fprintf('\n  Theta (x=  0): %12.4e',t0);
    fprintf('\n  Theta (x= 36): %12.4e',t1);
    fprintf('\n  Theta (x=108): %12.4e',t2);
    fprintf('\n  Theta (x=144): %12.4e\n',t3);
  end
end

% Load input parameters into matrices
Force=[ExtForce,ExtForceX]; 
NoExtForce=length(ExtForce);
[NoExtRamp,ncol]=size(ExtRamp); 
IntSup=[IntSupDelta,IntSupX]; 
NoIntSup=length(IntSupX); 
EndCond=[EndCondVal,EndCondFunc,EndCondEnd];
if EIorDepth == 1
  BeamProp=[]; NoEIorDepths=length(EIx);
  EIdata=[EIvalue EIx];
else
  BeamProp=[BeamWidth BeamE]; 
  NoEIorDepths=length(DepthX);
  EIdata=[Depth DepthX];
end

% more on

% Output input data
label1=['shear     ';'moment    '; ...
        'slope     ';'deflection'];
label2=['left   ';'right  '];
fprintf('\n\nAnalysis of a Variable Depth ');
fprintf('Elastic Beam');
fprintf('\n--------------------------------');
fprintf('---------');
fprintf('\n\n');
disp(['Title: ' Title]);
fprintf...
  ('\nBeam Length:                    %g', ...
  BeamLength);
fprintf...
  ('\nNumber of integration segments: %g', ...
  NoSegs);
fprintf...
  ('\nPrint frequency for results:    %g', ...
  Printout);
fprintf('\n\nInterior Supports: (%g)', ...
  NoIntSup);
if NoIntSup > 0
  fprintf('\n  |   #   X-location   Deflection');
  fprintf('\n  | --- ------------ ------------');
  for i=1:NoIntSup
    fprintf('\n  |%4.0f %12.4e %12.4e', ...
            i,IntSup(i,2),IntSup(i,1)); 
  end
end
fprintf('\n\nConcentrated Forces: (%g)', ...
  NoExtForce);
if NoExtForce > 0 
  fprintf('\n  |   #   X-location        Force');
  fprintf('\n  | --- ------------ ------------');
  for i=1:NoExtForce
    fprintf('\n  |%4.0f %12.4e %12.4e', ...
            i,Force(i,2),Force(i,1)); 
  end
end
fprintf('\n\nRamp loads: (%g)', NoExtRamp);
if NoExtRamp > 0
  fprintf('\n  |   #      X-start         Load');
  fprintf('        X-end         Load');
  fprintf('\n  | --- ------------ ------------');
  fprintf(' ------------ ------------');
  for i=1:NoExtRamp
    fprintf('\n  |%4.0f %12.4e %12.4e ', ...
            i,ExtRamp(i,3),ExtRamp(i,1));
    fprintf('%12.4e %12.4e', ...
            ExtRamp(i,4),ExtRamp(i,2)); 
  end
end
fprintf('\n\nEnd conditions:');
fprintf('\n  | End    Function           Value');
fprintf('\n  ');
fprintf('| ------ ----------  ------------\n');
for i=1:4
  j=EndCond(i,3); k=EndCond(i,2);
  strg=sprintf('  %12.4e',EndCond(i,1));
  disp(['  | ' label2(j,:) label1(k,:) strg]);
end
if EIorDepth == 1
  fprintf('\nEI values are specified');
  fprintf('\n  |   #      X-start     EI-value') 
  fprintf('\n  | --- ------------ ------------');
  for i=1:NoEIorDepths
    fprintf('\n  |%4.0f %12.4e %12.4e', ...
            i,EIdata(i,2),EIdata(i,1));  
  end
else 
  fprintf('\nDepth values are specified for ');
  fprintf('rectangular cross section');
  fprintf('\n  | Beam width:      %12.4e', ...
          BeamProp(1));
  fprintf('\n  | Young''s modulus: %12.4e', ...
          BeamProp(2));
  fprintf('\n  |');
  fprintf('\n  |   #      X-start        Depth') 
  fprintf('\n  | --- ------------ ------------');
  for i=1:NoEIorDepths
    fprintf('\n  |%4.0f %12.4e %12.4e', ...
            i,EIdata(i,2),EIdata(i,1));  
  end
end
disp(' ');
  
% Begin analysis
x=linspace(0,BeamLength,NoSegs)'; t=clock;
[V,M,Theta,Delta,Reactions]= ...
  bmvardep(NoSegs,BeamLength,Force,ExtRamp, ...
           EndCond,IntSup,EIdata,BeamProp);
t=etime(clock,t); 

% Output results
disp(' ');
disp(['Solution time was ',num2str(t),' secs.']);
if NoIntSup > 0
  fprintf('\nReactions at Internal Supports:');
  fprintf('\n  |   X-location     Reaction');
  fprintf('\n  | ------------ ------------');
  for i=1:NoIntSup
    fprintf('\n  | %12.8g %12.4e', ...
            IntSup(i,2),Reactions(i));
  end
end
fprintf('\n\nTable of Results:');
fprintf('\n  |  X-location        Shear');
fprintf('       Moment');
fprintf('        Theta        Delta');
fprintf('\n  | ----------- ------------ ');
fprintf('------------');
fprintf(' ------------ ------------');
if Printout > 0
  for i=1:Printout:NoSegs
    fprintf('\n  |%12.4g %12.4e %12.4e', ...
            x(i),V(i),M(i));
    fprintf(' %12.4e %12.4e',Theta(i),Delta(i));
  end
  disp(' ');
else
  i=1; j=NoSegs;
  fprintf('\n  |%12.4g %12.4e %12.4e', ...
          x(i),V(i),M(i));
  fprintf(' %12.4e %12.4e',Theta(i),Delta(i));
  fprintf('\n  |%12.8g %12.4e %12.4e', ...
          x(j),V(j),M(j));
  fprintf(' %12.4e %12.4e',Theta(j),Delta(j));
end
fprintf('\n\n');
subplot(2,2,1); 
  plot(x,V,'k-'); grid; xlabel('x axis');
  ylabel('Shear'); title('Shear Diagram');
subplot(2,2,2); 
  plot(x,M,'k-'); grid; xlabel('x axis');
  ylabel('Moment'); title('Moment Diagram')
subplot(2,2,3); 
  plot(x,Theta,'k-'); grid; xlabel('x axis');
  ylabel('Slope'); title('Slope Curve');
subplot(2,2,4); 
  plot(x,Delta,'k-'); grid; xlabel('y axis');
  ylabel('Deflection'); 
  title('Deflection Curve'); subplot
drawnow; subplot ; figure(gcf)
%print -deps vdb

% more off

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

function [V,M,Theta,Delta,Reactions]= ...
  bmvardep(NoSegs,BeamLength,Force,ExtRamp, ...
  EndCond,IntSup,EIdata,BeamProp)
% [V,M,Theta,Delta,Reactions]=bmvardep ...
% (NoSegs,BeamLength,Force,ExtRamp,EndCond, ...
% IntSup,EIdata,BeamProp)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function computes the shear, moment, 
% slope, and deflection in a variable depth 
% elastic beam having specified end conditions, 
% intermediate supports with given 
% displacements, and general applied loading, 
% allowing concentrated loads and linearly 
% varying ramp loads. 
%
% NoSegs     - number of beam divisions for
%              integration
% BeamLength - beam length
% Force      - matrix containing the magnitudes
%              and locations for concentrated
%              loads
% ExtRamp    - matrix containing the end
%              magnitudes and end locations
%              for ramp loads
% EndCond    - matrix containing the type of
%              end conditions, the magnitudes,
%              and whether values are for the
%              left or right ends
% IntSup     - matrix containing the location
%              and delta for interior supports
% EIdata     - either EI or depth values
% BeamProp   - either null or beam widths
%
% V          - vector of shear values
% M          - vector of moment values
% Theta      - vector of slope values
% Delta      - vector of deflection values
% Reactions  - reactions at interior supports
%
% User m functions required: 
%    oneovrei, extload, sngf, trapsum
%----------------------------------------------

if nargin < 8, BeamProp=[]; end
% Evaluate function value coordinates and 1/EI
x=linspace(0,BeamLength,NoSegs)'; 
kk=oneovrei(x,EIdata,BeamProp);

% External load contributions to shear and 
% moment interior to span and at right end
[ve,me]=extload(x,Force,ExtRamp);
[vv,mm]=extload(BeamLength,Force,ExtRamp);

% Deflections and position of interior supports
ns=size(IntSup,1);
if ns > 0 
  ysprt=IntSup(:,1); r=IntSup(:,2); 
  snf=sngf(x,r,1);
else 
  ysprt=[]; r=[]; snf=zeros(NoSegs,0); 
end

% Form matrix governing y''(x)
smat=kk(:,ones(1,ns+3)).* ...
     [x,ones(NoSegs,1),snf,me];

% Integrate twice to get slope and deflection 
% matrices
smat=trapsum(0,BeamLength,smat); 
ymat=trapsum(0,BeamLength,smat);

% External load contributions to 
% slope/deflection at the right end
ss=smat(NoSegs,ns+3); yy=ymat(NoSegs,ns+3);

% Equations to solve for left end conditions 
% and internal reactions
ns4=ns+4; j=1:4; a=zeros(ns4,ns4); 
b=zeros(ns4,1); js=1:ns; js4=js+4;

% Account for four independent boundary 
% conditions.  Usually two conditions will be 
% imposed at each end.
for k=1:4
  val=EndCond(k,1); typ=EndCond(k,2); 
  wchend=EndCond(k,3);
  if wchend==1
    b(k)=val; row=zeros(1,4); row(typ)=1; 
    a(k,j)=row;
  else
    if typ==1     % Shear
      a(k,j)=[1,0,0,0]; b(k)=val-vv;
      if ns>0 
        a(k,js4)=sngf(BeamLength,r,0); 
      end
    elseif typ==2 % Moment
      a(k,j)=[BeamLength,1,0,0]; b(k)=val-mm;
      if ns>0 
        a(k,js4)=sngf(BeamLength,r,1); 
      end
    elseif typ==3 % Slope
      a(k,j)=[smat(NoSegs,1:2),1,0]; 
      b(k)=val-ss;
      if ns>0 
        a(k,js4)=smat(NoSegs,3:ns+2); 
      end
    else          % Deflection
      a(k,j)=[ymat(NoSegs,1:2),BeamLength,1]; 
      b(k)=val-yy;
      if ns>0 
        a(k,js4)=ymat(NoSegs,3:ns+2); 
      end
    end
  end
end

% Interpolate to assess how support deflections 
% are affected by end conditions, external 
% loads, and support reactions.
if ns>0
  a(js4,1)=interp1(x,ymat(:,1),r);
  a(js4,2)=interp1(x,ymat(:,2),r);
  a(js4,3)=r; a(js4,4)=ones(ns,1);
  for j=1:ns-1 
    a(j+5:ns+4,j+4)= ...
      interp1(x,ymat(:,j+2),r(j+1:ns));
  end
end
b(js4)=ysprt-interp1(x,ymat(:,ns+3),r);

% Solve for unknown reactions and end conditions
c=a\b; v0=c(1); m0=c(2); s0=c(3); y0=c(4); 
Reactions=c(5:ns+4);

% Compute the shear, moment, slope, deflection
% for all x
if ns > 0
  V=v0+ve+sngf(x,r,0)*Reactions;
  M=m0+v0*x+me+sngf(x,r,1)*Reactions;
  Theta=s0+smat(:,ns+3)+smat(:,1:ns+2)* ...
        [v0;m0;Reactions];
  Delta=y0+s0*x+ymat(:,ns+3)+ ...
        ymat(:,1:ns+2)*[v0;m0;Reactions];
else
  Reactions=[]; V=v0+ve; M=m0+v0*x+me;
  Theta=s0+smat(:,ns+3)+smat(:,1:2)*[v0;m0];
  Delta=y0+s0*x+ymat(:,ns+3)+ ...
        ymat(:,1:2)*[v0;m0];
end

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

function [V,M,EITheta,EIDelta]=extload ...
         (x,Force,ExtRamp)
% [V,M,EITheta,EIDelta]=extload ...
%                       (x,Force,ExtRamp)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function computes the shear, moment, 
% slope, and deflection in a uniform depth  
% Euler beam which is loaded by a series of 
% concentrated loads and ramp loads. The values 
% of shear, moment, slope and deflection all 
% equal zero when x=0.
%
% x       - location along beam
% Force   - concentrated force matrix
% ExtRamp - distributed load matrix
%
% V       - shear
% M       - moment
% EITheta - slope
% EIDelta - deflection
%
% User m functions required: sngf
%----------------------------------------------

nf=size(Force,1); nr=size(ExtRamp,1); 
nx=length(x); V=zeros(nx,1); M=V; 
EITheta=V; EIDelta=V;
% Concentrated load contributions
if nf > 0
  F=Force(:,1); f=Force(:,2); V=V+sngf(x,f,0)*F; 
  M=M+sngf(x,f,1)*F;
  if nargout > 2
    EITheta=EITheta+sngf(x,f,2)*(F/2);
    EIDelta=EIDelta+sngf(x,f,3)*(F/6);
  end
end
% Ramp load contributions
if nr > 0
  P=ExtRamp(:,1); Q=ExtRamp(:,2); 
  p=ExtRamp(:,3); q=ExtRamp(:,4);
  S=(Q-P)./(q-p); sp2=sngf(x,p,2); 
  sq2=sngf(x,q,2); sp3=sngf(x,p,3); 
  sq3=sngf(x,q,3); sp4=sngf(x,p,4); 
  sq4=sngf(x,q,4);
  V=V+sngf(x,p,1)*P-sngf(x,q,1)* ... % Shear
    Q+(sp2-sq2)*(S/2); 
  M=M+sp2*(P/2)-sq2*(Q/2)+ ...       % Moment
    (sp3-sq3)*(S/6);
  if nargout > 2
    EITheta=EITheta+sp3*(P/6)- ...  % EI*Theta
            sq3*(Q/6)+(sp4-sq4)*(S/24);
    EIDelta=EIDelta+sp4*(P/24)- ... % EI*Delta
            sq4*(Q/24)+(sngf(x,p,5)- ...
            sngf(x,q,5))*(S/120); 
  end
end

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

function val=oneovrei(x,EIdata,BeamProp)
% [val]=oneovrei(x,EIdata,BeamProp) 
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% This function computes 1/EI by piecewise 
% linear interpolation through a set of data 
% values.
%
% x        - location along beam
% EIdata   - EI or depth values
% BeamProp - null or width values
%
% val      - computed value for 1/EI
%
% User m functions required: none
%----------------------------------------------

if size(EIdata,1) < 2  % uniform depth case
  v=EIdata(1,1); 
  EIdata=[v,min(x);v,max(x)]; 
end
if ( nargin > 2 ) & ( sum(size(BeamProp)) > 0)
  % Compute properties assuming the cross 
  % section is rectangular and EIdata(:,1) 
  % contains depth values
  width=BeamProp(1); E=BeamProp(2);
  EIdata(:,1)=E*width/12*EIdata(:,1).^3; 
end
val=1./lintrp(EIdata(:,2),EIdata(:,1),x);

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

function y=sngf(x,x0,n)
% y=sngf(x,x0,n) 
% ~~~~~~~~~~~~~~
%
% This function computes the singularity
% function defined by 
%    y=<x-x0>^n for n=0,1,2,...
%
% User m functions required: none
%----------------------------------------------

if nargin < 3, n=0; end
x=x(:); nx=length(x); x0=x0(:)'; n0=length(x0); 
x=x(:,ones(1,n0)); x0=x0(ones(nx,1),:); d=x-x0; 
s=(d>=zeros(size(d))); v=d.*s;
if n==0 
  y=s;
else 
  y=v; 
  for j=1:n-1; y=y.*v; end
end

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

function v=trapsum(a,b,y,n)
%
% v=trapsum(a,b,y,n)
% ~~~~~~~~~~~~~~~~~~
%
% This function evaluates:
%
%   integral(a=>x, y(x)*dx) for a<=x<=b
%
% by the trapezoidal rule (which assumes linear
% function variation between succesive function
% values).
%
% a,b - limits of integration
% y   - integrand that can be a vector-valued
%       function returning a matrix such that
%       function values vary from row to row.
%       It can also be input as a matrix with
%       the row size being the number of
%       function values and the column size
%       being the number of components in the
%       vector function.
% n   - the number of function values used to
%       perform the integration.  When y is a
%       matrix then n is computed as the number
%       of rows in matrix y.
%
% v   - integral value
%
% User m functions called:  none
%----------------------------------------------

if isstr(y)
  % y is an externally defined function
  x=linspace(a,b,n)'; h=x(2)-x(1);
  Y=feval(y,x); % Function values must vary in 
                % row order rather than column 
                % order or computed results 
                % will be wrong.
  m=size(Y,2);
else
  % y is column vector or a matrix
  Y=y; [n,m]=size(Y); h=(b-a)/(n-1);
end
v=[zeros(1,m); ...
  h/2*cumsum(Y(1:n-1,:)+Y(2:n,:))];

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

% function y=lintrp(xd,yd,x)
% See Appendix B

Contact us