% steppedshaft.m Last modified 4/16/09 % This script calculates and plots the load, shear, moment and deflection % diagrams for a stepped, simply supported shaft supported at arbitrary % locations with up to four steps with an arbitrary number of loads of % arbitrary size. % This m-file was written at the University of Wyoming in the Electrical % and Computer Engineering Department and is to be distributed without % cost. clc clear all disp(' This script calculates and plots the load, shear, moment and deflection') disp(' diagrams for a stepped, simply supported shaft supported at arbitrary') disp(' locations with up to four steps with an arbitrary number of loads of') disp(' arbitrary size. The units are arbitrary but must be consistent. For') disp(' instance length in inches and force in pouunds hence bending stiffness') disp(' EI has units of lb-in^2.') disp(' ') L=-1; while L<=0 L=input(' Shaft length = '); end disp(' ') xb=2*L; xa=-1; disp(' Input the locations of the two supports from left to right') while xa<0|xa>L disp(' ') xa=input(' The location of the left support is xa = '); end xaindex=round(xa*200/L)+1; while xbL disp(' ') xb=input(' The location of the right support is xb = '); end xbindex=round(xb*200/L)+1; nstep=-1; disp(' '); disp(' The number of steps in the shaft is less than or equal to four') disp(' If it is a constant cross-section shaft enter 0') disp(' ') while (nstep<0)&((nstep~=1)|(nstep~=2)|(nstep~=3)|(nstep~=4)) nstep=input(' Number of steps between the ends of the shaft = '); end steplocation=zeros(1,4); steplocationindex=zeros(1,5); if nstep>=1 disp(' Input the locations of the steps in ascending order') end if nstep>0 for i=1:nstep disp(' ') while steplocationindex(1,i)==0|steplocationindex(1,i)>=200 steplocation(1,i)=input([' Location of step no. ' num2str(i) ' = ']); steplocationindex(1,i)=round((steplocation(1,i)/L)*200)+1; end end steplocation(1,nstep+1)=L; end if nstep==0 steplocation(1,1)=L; steplocationindex=201; end; nspan=nstep+1; stiffness=zeros(1,nspan); for i=1:nspan disp(' ') stiffness(1,i)=input([' Bending stiffness of span no. ' num2str(i) ' = ']); end EI=zeros(1,201); if nspan==1 EI=stiffness(1,1)*ones(1,201); end if nspan>1 EI(1,1:steplocationindex(1,1)-1)=stiffness(1,1)*ones(1,steplocationindex(1,1)-1); for i=1:nspan-2 EI(1,steplocationindex(1,i)+1:steplocationindex(1,i+1)-1)=stiffness(1,i+1)... *ones(1,steplocationindex(1,i+1)-steplocationindex(1,i)-1); end EI(1,steplocationindex(1,nstep)+1:201)=stiffness(1,nspan)*ones(1,201-steplocationindex(1,nstep)); end for i=1:nstep EI(1,steplocationindex(1,i))=0.5*(EI(1,steplocationindex(1,i)-1)... +EI(1,steplocationindex(1,i)+1)); end EImax=max(stiffness); thk=((EI/EImax))*5; thk(1,1)=0; thk(1,201)=0; disp(' ') nload=-1; while nload<=0 nload=input(' Number of loads = '); end xj=-ones(1,nload); disp(' ') disp(' Loads and deflections downward are positive') disp(' ') for j=1:nload P(1,j)=input([' Load ' num2str(j) ' = ']); xj(1,j)=-1; disp(' ') while xj(1,j)<0|xj(1,j)>L xj(1,j)=input([' Load ' num2str(j) ' location = ']); end disp(' ') end %calculate reactions loadindex=round(200*xj/L)+1; Rb=-sum(P.*((xj-xa*ones(1,nload))))/(xb-xa); Ra=-sum(P)-Rb; disp(' ') x=linspace(0,L,201); forceindex=[]; F=[]; xforce=[]; if xaindex==1 forceindex=[1]; xforce=[0]; F=[Ra]; end if loadindex(1,1)==1 forceindex=[1]; xforce=[0]; F=[P(1,1)]; end % fill F and forceindex vectors for i=2:201 flag=0; for j=1:nload if i==loadindex(1,j)&flag==0 forceindex=[forceindex round(xj(1,j)*200/L)+1]; xforce=[xforce xj(1,j)]; F=[F P(1,j)]; flag=1; elseif i==xaindex&flag==0 forceindex=[forceindex xaindex]; xforce=[xforce xa]; F=[F Ra]; flag=1; elseif i==xbindex&flag==0 forceindex=[forceindex xbindex]; xforce=[xforce xb]; F=[F Rb]; flag=1; end end end %summarize the shaft, load and reaction data clc disp(' SHAFT SUMMARY') disp([' Shaft length L = ' num2str(L)]) if nstep==0 string1=num2str(L); else string1=num2str(steplocation(1,1)); end disp([' Span 1 from 0 to ' string1 ', EI = ' num2str(stiffness(1,1))]); if nstep>0 for i=2:nspan disp([' Span ' num2str(i) ' from ' num2str(steplocation(1,i-1)) ' to ' num2str(steplocation(1,i)) ', EI = ' num2str(stiffness(1,i))]); end end disp(' ') disp(' ') disp(' LOAD SUMMARY') disp(' For loads, downward will be taken as positive') for i=1:nload disp([' Load ' num2str(i) ' = ' num2str(P(1,i)) ' at x = ' num2str(xj(1,i))]) end disp(' ') disp(' ') disp(' REACTION SUMMARY') disp(' For reactions, downward will be taken as positive') disp([' Ra = ' num2str(Ra) ' at x = ' num2str(xa)]) disp([' Rb = ' num2str(Rb) ' at x = ' num2str(xb)]) disp(' ') disp(' Press Enter to Continue') pause figure(1);clf reset axis([ -.1*L 1.1*L -100 100]) hold on plot([0 L],[0 0], 'LineWidth',[2.5]) hold on text(0,90,'Freebody Diagram of Beam---Downward is Positive') xarrow=[0 0 .03*L -.03*L 0]; yarrow=[70 0 20 20 0]; for j=1:nload+2 if F(1,j)>0 plot(xarrow+xforce(1,j)*ones(1,5),yarrow,'LineWidth',[1.5]) hold on text(xforce(1,j)-.05*L, 72,num2str(F(1,j))) elseif F(1,j)<0 plot(xarrow+xforce(1,j)*ones(1,5),-yarrow,'LineWidth',[1.5]) hold on text(xforce(1,j)-.05*L, -72,num2str(F(1,j))) end end if Ra>0 ytextloc=40; else ytextloc=-40; end text(xa+.02*L,ytextloc,'Ra') if Rb>0 ytextloc=40; else ytextloc=-40; end text(xb+.02*L,ytextloc,'Rb') box on plot(x,thk,'k'); hold on plot(x,-thk,'k'); xlabel('Distance from Left End, x') text(0,-90,'Press Enter to Continue') hold off pause % calculate shears shear=zeros(1,201); if forceindex(1,1)==1 % support at left end shear(1,2)=-F(1,1); start=2; else shear(1,2)=0; start=1; end for i=3:200% location index flag=0; for j=start:nload+2% cycle thru the forces if i==forceindex(1,j)&flag==0% if so modify shear shear(1,i)=shear(1,i-1)-F(1,j); flag=1; elseif i~=forceindex(1,j)&flag~=1% if not leave shear the same shear(1,i)=shear(1,i-1); %flag=1; end end end for i=1:nload+2 if forceindex(1,i)>2&forceindex(1,i)<200 shear(1,forceindex(1,i))=(shear(1,forceindex(1,i)-1)+shear(1,forceindex(1,i)+1))/2; end end figure(2);clf reset maxshear=max(shear); minshear=min(shear); axis([-.1*L 1.1*L 1.2*minshear 1.2*maxshear]) hold on plot([0 L],[0 0], 'LineWidth',[2.5]) hold on plot(x,shear) box on xlabel('Distance from Left End, x') ylabel('Shear Force, V(x)') text(0,1.08*maxshear,'Shear Diagram---Up on the Left Edge of an Element is Positive') text(0,1.05*minshear,'Press Enter to Continue') hold off pause moment=zeros(1,201); for i=2:200 m=0; for j=1:nload+2; if iforceindex(1,j) m=m-F(1,j)*(x(1,i)-xforce(1,j)); end moment(1,i)=m; end end maxmoment=max(moment); minmoment=min(moment); figure(3);clf reset Maxmom=max([abs(maxmoment) abs(minmoment)]); if maxmoment>0&minmoment>=0 textloc1=1.1*maxmoment; textloc2=-.1*maxmoment; axis([-.1*L 1.1*L -.2*maxmoment 1.2*maxmoment]) hold on elseif maxmoment<=0&minmoment<0 textloc1=-.1*minmoment; textloc2=1.1*minmoment; axis([-.1*L 1.1*L 1.2*minmoment -.2*minmoment]) hold on elseif maxmoment>0&minmoment<0 range=maxmoment-minmoment; textloc1=maxmoment+.1*range; textloc2=minmoment-.1*range; axis([-.1*L 1.1*L minmoment-.2*range maxmoment+.2*range]) hold on end plot([0 L],[0 0], 'LineWidth',[2.5]) hold on plot(x,moment) xlabel('Distance from Left End, x') ylabel('Bending Moment, M(x)') text(0,textloc1,'Moment Diagram---Compression in Top an Element is Positive') text(0,textloc2,['Press Enter to Continue Max. Absolute Moment = ' num2str(Maxmom)]) box on hold off pause f=moment./EI; deltaL=L/200; slope=zeros(1,201); for i=1:200 slope(1,i+1)=slope(1,i)-f(1,i)*deltaL; end %calculate deflections defl1=zeros(1,101); for i=1:200 defl1(1,i+1)=defl1(1,i)+slope(1,i)*deltaL; end defl2=defl1-defl1(1,xaindex)*ones(1,201); deflection=defl2-defl2(1,xbindex)*((x-xa*ones(1,201))/(xb-xa)); maxdefl=max(deflection); mindefl=min(deflection); hold off figure(4);clf reset Maxdef=max([abs(maxdefl) abs(mindefl)]); if maxdefl>0&mindefl>=0 textloc1=1.1*maxdefl; textloc2=-.1*maxdefl; h1=axes('Ydir','reverse'); hold on axis([-.1*L 1.1*L -.2*maxdefl 1.2*maxdefl]) hold on elseif maxdefl<=0&mindefl<0 textloc1=-.1*mindefl; textloc2=1.1*mindefl; h1=axes('Ydir','reverse'); hold on axis([-.1*L 1.1*L 1.2*mindefl -.2*mindefl]) hold on elseif maxdefl>0&mindefl<0 range1=maxdefl-mindefl; textloc1=maxdefl+.1*range1; textloc2=mindefl-.1*range1; h1=axes('Ydir','reverse'); hold on axis([-.1*L 1.1*L mindefl-.2*range1 maxdefl+.2*range1]) hold on end plot([0 L],[0 0], 'LineWidth',[2.5]) hold on plot(x,deflection) hold on xlabel('Distance from Left End, x') ylabel('Deflection, y(x)') text(0,textloc2,'Deflection Diagram---Downward is Positive') text(0,textloc1,['Press Enter to Continue Maximum Absolute Defl. = ' num2str(Maxdef)]) box on pause