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)

[c,phi,stres,z]=...
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 along the x axis and 2*b deep along 
% the y axis. The shear stresses in the beam 
% are given by the complex 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 in 
% nsegb       formation of least square 
%             equations
% ntrms   - number of polynomial terms used in 
%             polynomial stress function
% nxout,  - number of grid points used to
% nyout       evaluate output
% c       - coefficeints 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 

Contact us