# Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

### Howard Wilson (view profile)

14 Oct 2002 (Updated )

Companion Software (amamhlib)

rector
```function rector
% Example:  rector
% ~~~~~~~~~~~~~~~~
% This program uses point matching to obtain an
% approximate solution for torsional stresses
% in a Saint Venant beam having a rectangular
% cross section. The complex stress function is
% analytic inside the rectangle and has its
% real part equal to abs(z*z)/2 on the
% boundary. The problem is solved approximately
% using a polynomial stress function which fits
% the boundary condition in the least square
% sense. Surfaces and contour curves describing
% the stress and deformation pattern in the
% beam cross section are drawn.
%
% User m functions required: recstrs

clear; close
fprintf('\n===  TORSIONAL STRESS CALCULATION');
fprintf(' IN A RECTANGULAR  ===');
fprintf('\n===      BEAM USING LEAST SQUARE ');
fprintf('APPROXIMATION      ===\n');
fprintf('\nInput the lengths of the ');
fprintf('horizontal and the vertical sides\n');
fprintf('(make the long side horizontal)\n');
u=input('> ? ','s'); u=eval(['[',u,']']);
a=u(1)/2; b=u(2)/2;

% The boundary conditions are approximated in
% terms of the number of least square points
% used along the sides
nsegb=100; nsega=ceil(a/b*nsegb);
nsega=fix(nsega/2); nsegb=fix(nsegb/2);
fprintf('\nInput the number of terms ');
fprintf('used in the stress function');
fprintf('\n(30 terms is usually enough)\n');
ntrms=input('> ? ');

% Define a grid for evaluation of stresses.
% Include the middle of each side.
nx=41; ny=fix(b/a*nx); ny=ny+1-rem(ny,2);

[c,phi,stres,z] = ...
recstrs(a,nsega,b,nsegb,ntrms,nx,ny);
[smax,k]=max(abs(stres(:))); zmax=z(:);
zmax=zmax(k); xmax=abs(real(zmax));
ymax=abs(imag(zmax));
disp(' '), disp(['The Maximum Shear ',...
'Stress is ',num2str(smax)]);
disp(['at x = ',num2str(xmax),' and y = ',...
num2str(ymax)]);
disp(' '); disp('All Done');

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

function [c,phi,stres,z]=...
recstrs(a,nsega,b,nsegb,ntrms,nxout,nyout)
%
% [c,phi,stres,z]=...
%   recstrs(a,nsega,b,nsegb,ntrms,nxout,nyout)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function uses least square fitting to
% obtain an approximate solution for torsional
% stresses in a Saint Venant beam having a
% rectangular cross section. The complex stress
% function is analytic inside the rectangle
% and has its real part equal to abs(z*z)/2 on
% the boundary. The problem is solved
% approximately using a polynomial stress
% function which fits the boundary condition
% in the least square sense. The beam is 2*a
% wide parallel to the x axis and 2*b deep
% parallel to the y axis. The shear stresses
% in the beam are given by the stress formula:
%
% (tauzx-i*tauzy)/(mu*alpha) = -i*conj(z)+f'(z)
%
% where
%
%   f(z)=i*sum( c(j)*z^(2*j-2), j=1:ntrms )
%
% and c(j) are real.
%
% a,b     - half the side lengths of the
%           horizontal and vertical sides
% nsega,  - numbers of subintervals used to
% nsegb     form the least square equations
% ntrms   - number of terms used in the
%           polynomial stress function
% nxout,  - number of grid points used to
% nyout     evaluate output
% c       - coefficients defining the stress
%           function
% phi     - values of the membrane function
% stres   - array of complex stress values
% z       - complex point array at which
%           stresses are found
%
% User m functions called:  none
%----------------------------------------------

% Generate vector zbdry of boundary points
% for point matching.
zbdry=[a+i*b/nsega*(0:nsega-1)';
i*b+a/nsegb*(nsegb:-1:0)'];

% Determine a scaling parameter used to
% prevent occurrence of large numbers when
% high powers of z are used
s=max(a,b);

% Form the least square equations to impose
% the boundary conditions.
neq=length(zbdry); amat=ones(neq,ntrms);
ztmp=(zbdry/s).^2; bvec=.5*abs(zbdry).^2;
for j=2:ntrms
amat(:,j)=amat(:,j-1).*ztmp;
end

% Solve the least square equations.
amat=real(amat); c=pinv(amat)*bvec;

% Generate grid points to evaluate
% the solution.
xsid=linspace(-a,a,nxout);
ysid=linspace(-b,b,nyout);
[xg,yg]=meshgrid(xsid,ysid);
z=xg+i*yg; zz=(z/s).^2;

% Evaluate the warping function
phi=-imag(polyval(flipud(c),zz));

% Evaluate stresses and plot results
cc=(2*(1:ntrms)-2)'.*c;
stres=-i*conj(z)+i* ...
polyval(flipud(cc),zz)./(z+eps*(z==0));
am=num2str(-a);ap=num2str(a);
bm=num2str(-b);bp=num2str(b);

% Plot results
disp(' '), disp('Press [Enter] to plot')
disp('the warping surface'), pause
[pa,k]=max(abs(phi(:)));
Phi=a/4*sign(phi(k))/phi(k)*phi;
close, colormap('default')
surfc(xg,yg,Phi)
title('Warping of the Cross Section')
xlabel('x axis'), ylabel('y axis')
zlabel('transverse warping'); axis('equal')
shg, disp(' ')
disp('Press [[Enter]] to plot the')
disp('total stress surface'), pause
% print -deps warpsurf

surfc(xg,yg,abs(stres));
title('Total Shear Stress Surface')
xlabel('x axis'); ylabel('y axis')
zlabel('total stress'), axis('equal'), shg
disp(' '), disp('Press [Enter] to plot the')
disp('stress contours'), pause
% print -deps rectorst

contour(xg,yg,abs(stres),20); colorbar
title('Total Stress Contours');
xlabel('x axis'); ylabel('y axis')
shg, disp(' ')
disp('Press [Enter] to plot the maximum')
disp('stress on a rectangle side'), pause
% print -deps torcontu

plot(xsid,abs(stres(1,:)),'k');
grid; ylabel('tangential stress');
xlabel('position on a horizontal side');
title('Stress for y = b/2'); shg
% print -deps torstsid
```