cbeamex.m

% Example:cbeamex 
% ~~~~~~~~~~~~~~~
% This program calculates the stress 
% distribution in a curved beam having a 
% polygonal cross section.
%
% Data is defined in the declaration statements
% below, where:
%
% x,y - vectors containing the coordinates for
%       the polygon nodes.  CCW fashion for
%       positive contributions, CW fashion
%       for negative contributions.
% M   - moment on the cross section
% No_planes - number of planes to calculate
%       stresses at 
%
% NOTE: To make the plotting perform 
%       as desired an undocumented MATLAB
%       capability has been used (Renderlimits).
%       MathWorks does not promise that it
%       will be supported in the future.
%       Renderlimits allows the two side-by-
%       side plots to maintain the same 
%       y-axes. (no longer applicable!)
%
% User m functions required:
%    cbprop, genprint
%----------------------------------------------

clear;
%...Input definitions
Problem=2;
if Problem == 1
  %...Rectangular
  x=[0    0.04  0.04 0];
  y=[0.04 0.04  0.1  0.1];
  M=-1162; No_planes=10;
elseif Problem == 2
  %...Trapezoidal
  x=[-1.57/2 1.57/2 0.4/2 -0.4/2];
  y=[3.0     3.0    1.5    1.5];
  M=-17000; No_planes=10;
elseif Problem == 3
  %...T-shape cross section
  x=[-40 40 40 10 10 -10 -10 -40];
  y=[ 30 30 50 50 90  90  50  50];
  x=x/1000; y=y/1000; % metric conversion
  M=940.5; No_planes=10;
end

%...Initialize
No_pts=length(x); n1=No_planes+1;
sigma=zeros(1,n1);

%...Planes to find stresses at
ymin=min(y); ymax=max(y);
yi=linspace(ymin,ymax,n1);

%...Geometrical properties of cross section
[A,Qx,K,ybar]=cbprop(x,y);

%...Radius of curvature and the distance 
%...between neutral axis and centroidal axis
rho=A/K; e=ybar-rho; 

%...Define stress equation constants and
%...solve for bending stresses 
alpha=M*A/(A^2-Qx*K); beta=-M*K/(A^2-Qx*K);
sigma=alpha./yi+beta;

fprintf( ...
'\n\nStresses in Cross Section of Curved Beam');
fprintf( ...
'\n----------------------------------------');
fprintf( ...
  '\n\nTable of Nodes Forming Cross Section');
fprintf( ...
    '\n\n  Node #         x             y');
for i=1:No_pts
  fprintf('\n  %3.0f %12.5g  %12.5g', ...
          i,x(i),y(i));
end
fprintf('\n\nGeometrical Properties');
fprintf(  '\n  A:     %g',A);
fprintf(  '\n  Qx:    %g',Qx);
fprintf(  '\n  K:     %g',K);
fprintf(  '\n  y-bar: %g',ybar);
fprintf(  '\n  e:     %g',e);
fprintf(  '\n  rho:   %g',rho);
fprintf('\n');
fprintf('\nBending moment M:         %g',M);
fprintf('\n\nStress Equation Coefficients');
fprintf(  '\n  alpha: %g',alpha);
fprintf(  '\n  beta:  %g',beta);
fprintf('\n');

for i=1:No_planes+1
  fprintf('\nPlane number %g',i);
  fprintf('\n  y:                      %g',yi(i));
  fprintf('\n  Bending stress:         %g', ...
    sigma(i));
end
fprintf('\n\n');

xp=x; yp=y; 
xp(No_pts+1)=xp(1); yp(No_pts+1)=yp(1);
xr1=[min(xp) max(xp)]; 
xr2=[min(sigma) max(sigma)]; 
yybar=[ybar ybar]; yrho=[rho rho];
clf;
ax=subplot(1,2,1);
  fill(xp,yp,'y'); hold on;
  plot(xp,yp,'-',xp,yp,'o');
  plot(xr1,yybar,'r:');
  plot(xr1,yrho,'g-');
  title('Cross Section Analyzed');
  xlabel('x'); ylabel('y');
  set(gca,'YDir','reverse');
  axis('equal'); hold off; drawnow; 
  yminmax=ylim;
subplot(1,2,2);
  plot(sigma,yi,'-y',xr2,yrho,'g-');
  title('Bending Stress Distribution');
  xlabel('Stress'); ylabel('y');
  set(gca,'YDir','reverse');
  ylim([yminmax(1) yminmax(2)]);
  drawnow;
  %genprint('cbtotal');