shrstr2.m

% Example: shrstr2
% ~~~~~~~~~~~~~~~~
% This example determines the shear stress
% distribution in 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.
% V   - shear on the cross section
% No_planes - number of horizontal planes to
%             use for clipping
%
% 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:
%    prop, shftprop, genprint, clip, inout,
%    intrsect, genarea, findt
%----------------------------------------------

clear;
%...Input definitions
Problem=3;
if Problem == 1
  %...Rectangle, bottom left corner at (0,0)
  x=[ 0  1  1  0]; y=[ 0  0  5  5];
  V=1; No_planes=10; 
elseif Problem == 2
  %...T-shaped cross section
  x=[0 .5 .5 2.25 2.25 -1.75 -1.75 0];
  y=[0 0 2 2 2.5 2.5 2 2];
  V=1.5; No_planes=20; 
elseif Problem == 3
  %...I beam
  x=[3 3 0.5 0.5 3 3 -3 -3 -.5 -.5 -3 -3];
  y=[-5 -4.5 -4.5 4.5 4.5 5 5 ...
     4.5 4.5 -4.5 -4.5 -5];
  V=2000; No_planes=80; 
end

%...Initialize
No_pts=length(x); n1=No_planes+1;
Q=zeros(1,n1); q=zeros(1,n1);
t=zeros(1,n1); tau=zeros(1,n1);

%...Determine geometrical properties
[A,Qx,Qy,Ix,Iy,Ixy,xbar,ybar]=prop(x,y);
%...Shift to centroidal axis
[Ixbar,Iybar,Ixybar]= ...
  shftprop(A,xbar,ybar,Ix,Iy,Ixy);
xs=x-xbar; ys=y-ybar;
%...Extreme coordinates
ymax=max(ys); ymin=min(ys);

%...Horizontal planes to find shear stresses at
y_shear=linspace(ymin,ymax,n1);

%...Calculate shear at each plane
No_pts_clip=No_pts; x_clip=xs; y_clip=ys;
for i=2:No_planes
  %...Generate the clipped area  
  [No_pts_clip,x_clip,y_clip]=genarea( ...
    No_pts_clip,x_clip,y_clip,y_shear(i));
  %...Find the first moment of area
  [D1,Q(i),D2,D3,D4,D5,D6,D7]= prop( ...
    x_clip,y_clip);
  %...Determine the t at the plane of interest
  [t(i)]=findt(x_clip,y_clip,y_shear(i));
end
%...Compute shear flow and stress
q=(V/Ixbar)*Q; 
n=n1-1; tau(2:n)=q(2:n)./t(2:n);

%...Shift coordinates to original system
y_shear=y_shear+ybar;

fprintf('\n\nShear Stresses in Cross Section');
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  Qy:    %g',Qy);
fprintf(  '\n  Ix:    %g',Ix);
fprintf(  '\n  Iy:    %g',Iy);
fprintf(  '\n  Ixy:   %g',Ixy);
fprintf(  '\n  x-bar: %g',xbar);
fprintf(  '\n  y-bar: %g',ybar);
fprintf('\n\nGeometrical Properties About ');
fprintf(    'Centroidal Axis');
fprintf(  '\n  Ix-bar:  %g',Ixbar);
fprintf(  '\n  Iy-bar:  %g',Iybar);
fprintf(  '\n  Ixy-bar: %g',Ixybar);
fprintf('\n\nShear (V): %g',V);
for i=1:No_planes+1
  fprintf('\n\nClipping plane number %g',i);
  fprintf('\n  y:                %g', ...
          y_shear(i));
  fprintf('\n  First moment (Q): %g',Q(i));
  fprintf('\n  Thickness:        %g',t(i));
  fprintf('\n  Shear flow (q):   %g',q(i));
  fprintf('\n  Shear stress:     %g',tau(i));
end
fprintf('\n\n');

xp=[x(:);x(1)]; yp=[y(:);y(1)]; 
clf;
ax=subplot(1,2,1);
  fill(xp,yp,'y'); hold on
  plot(xp,yp,'-',xp,yp,'o');
  title('Cross Section Analyzed');
  xlabel('x'); ylabel('y');
  axis('equal'); hold off; drawnow; 
  yminmax=ylim;

subplot(1,2,2);
  plot(tau,y_shear,'-')
  title( ...
  'Shear Stress Distribution');
  xlabel('Shear Stress'); ylabel('y');
  ylim([yminmax(1) yminmax(2)]);
  drawnow;
% genprint('shrstr2');