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



14 Oct 2002 (Updated )

Companion Software (amamhlib)

function [c,phi,stres,z]=...
% [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.
% Determine a scaling parameter used to 
% prevent occurrence of large numbers when 
% high powers of z are used

% 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 

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

% Generate grid points to evaluate 
% the solution.
z=xg+i*yg; zz=(z/s).^2; 

% Evaluate the warping function

% Evaluate stresses and plot results
stres=-i*conj(z)+i* ...

% Plot results
disp(' '), disp('Press [Enter] to plot')
disp('the warping surface'), pause 
close, colormap('default') 
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

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 
grid; ylabel('tangential stress');
xlabel('position on a horizontal side');
title('Stress for y = b/2'); shg
% print -deps torstsid 

Contact us