mxnorex.m

% Example: mxnorex
% ~~~~~~~~~~~~~~~~
% This example determines the angles for
% moment loading required to produce a
% maximum or minimum stress in a polygonal
% cross section.
%
% Data is defined in the declaration statements
% below.
%
% x,y - vectors containing the coordinates for
%       the polygon nodes.  CCW fashion for
%       positive contributions, CW fashion
%       for negative contributions.
% P   - Applied axial force
% M   - Applied moment
%
% User m functions required:
%    prop, shftprop, genprint
%----------------------------------------------

clear;
%...Input definitions
Problem=1;
if Problem == 1
  %...Complex cross section
  x=[0 .4 .4 1.66 1.66 2.7 2.7 3.1 3.1 1 1 0];
  y=[0 0 .7 2 3 3 1 1 3.4 3.4 2 1];
  P=1000; M=3000;
elseif Problem == 2
  %...Symmetric angle
  x=[0 10 10 2 2 0]; y=[0 0 2 2 10 10];
  P=-2000; M=1000;
elseif Problem == 3
  %...Z section
  x=[-1 7 7 5 5 1 1 -7 -7 -5 -5 -1];
  y=[-10 -10 -6 -6 -8 -8 10 10 6 6 8 8];
  P=-2000; M=-2000;
else
  %...Channel 
  x=[0 10 10 9 9 1 1 0];
  y=[0 0 10 10 1 1 10 10];
  P=0; M=5000;
end

%...Initialize
No_pts=length(x); raddeg=180/pi; Small=-1e20;

%...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;

%...Constants for stress function
Denom=Ixbar*Iybar-Ixybar^2;
C1=Ixbar/Denom; C2=Iybar/Denom;
C3=Ixybar/Denom;

%..................
%...Find worst case
%..................
Sigma_mxmn=Small;
K1=C2*ys-C3*xs; K2=C3*ys-C1*xs;
Beta=atan2(K2,K1);
cosine=cos(Beta); sine=sin(Beta);
Sigma=(C2*M*cosine+C3*M*sine).*ys- ...
      (C1*M*sine+C3*M*cosine).*xs;

%...Total stress at each corner
Sigma_min=P/A-abs(Sigma); 
Sigma_max=P/A+abs(Sigma);
%...Save the worst
[Sigma_mxmn,Index]=max(abs(Sigma));

for i=1:No_pts
  %...Get the angles right
  if Sigma(i) < 0
    %...Minimum stress
    Angle_min(i)=Beta(i);
    Angle_max(i)=rem(Beta(i)+pi,2*pi);
  else
    %...Maximum stress
    Angle_max(i)=Beta(i);
    Angle_min(i)=Beta(i)-pi;
  end
end
Angle_max=Angle_max*raddeg;
Angle_min=Angle_min*raddeg;

fprintf('\n\n     Angle for Moment Loading');
fprintf('\nRequired to Produce Max/Min Stress');
fprintf('\n   in a Polygonal 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\nMax/Min Stresses');
fprintf('\n  Axial force:    %g',P);
fprintf('\n  Applied moment: %g',M);
fprintf('\n\n  Node %g controls with:',Index);
fprintf('\n    Maximum stress: %g', ...
        Sigma_max(Index));
fprintf( ...
  '\n      at an angle of: %g degrees', ...
  Angle_max(Index));
fprintf('\n    Minimum stress: %g', ...
        Sigma_min(Index));
fprintf( ...
  '\n      at an angle of: %g degrees', ...
  Angle_min(Index));
fprintf('\n\nNote: 1) Angles measured from');
fprintf(' positive x-axis,');
fprintf(  '\n         + for counter-clockwise');
fprintf('\n      2) + for tensile stress,');
fprintf( ...
  '\n         - for compressive stress\n\n');

clf; x=[x(:);x(1)]; y=[y(:);y(1)];
fill(x,y,'y'); hold on;
  plot(x,y,'-',x,y,'o');
  title('Cross Section Analyzed');
  xlabel('x'); ylabel('y');
  axis('equal'); hold off; drawnow;
% genprint('mxnmsect');
  disp('Press key to continue'); pause

clf; x=1:No_pts; 
subplot(2,1,1);
  plot(x,Sigma_max,'-',x,Sigma_max,'o', ...
       x,Sigma_min,':',x,Sigma_min,'*');
  title('Max/Min Stress at Each Node')
  xlabel('Node Number');
  ylabel('Stress');
subplot(2,1,2);
  plot(x,Angle_max,'-',x,Angle_max,'o', ...
       x,Angle_min,':',x,Angle_min,'*');
  title('Moment Angle versus Node')
  xlabel('Node Number');
  ylabel('Moment Angle (degrees)'); drawnow;
% genprint('mxstress');