% movingload2.m Last modified 7/29/2008 % Animates the response of a simply-supported Bernoulli-Euler beam % with constant moving load (force only) Mg starting at the left end at t=0. % Input is the ratio of transit time to the first natural period of % the beam. The load is a constant force only without mass. % This m-file was written at the University of Wyoming in the Electrical % and Computer Engineering Department and is to be distributed without % cost. clear all string1='The ratio of transit time Tt to first natural period T1, Tt/T1 = '; nmode=5; Tratio=input(string1); sine=[]; omegaoveromega1=0.5/(Tratio); omegajoveromega1=zeros(1,10); xoverL=linspace(0,1,101); xoverL1=linspace(0,4,401); factor=96/(pi^4); for j=1:nmode sine(j,:)=sin(j*pi*xoverL); omegajoveromega1(1,j)=j^2; end xarrow=[0 0 .02 -.02 0]; yarrow=[-.5 -.01 -.1 -.1 -.01]; yoveryst=zeros(1,101); omegat=linspace(0,4*pi,401); for k=1:101%time loop sum=zeros(1,101); for j=1:nmode% series loop if j*omegaoveromega1==omegajoveromega1(1,j) term=0.5*(1/(j*omegaoveromega1))*((1/(j*j))*sin(j*omegat(1,k))-(1/omegaoveromega1)*omegat(1,k).*cos(j*omegat(1,k)))*sine(j,:); else num=(sin(j*omegat(1,k))-(j*omegaoveromega1/omegajoveromega1(1,j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*omegat(1,k))); term=(num/((omegajoveromega1(1,j)^2)-(j*omegaoveromega1)^2))*sine(j,:); end sum=sum+term; end yoveryst(k,:)=factor*sum; end y1=3*xoverL(1,1:51)-(4*(xoverL(1,1:51).^3)); yst=[y1 fliplr(y1(1,1:50))]; % evaluate initial conditions for unforced motion q and qprime are column % vectors omegat1=pi; for j=1:nmode if j*omegaoveromega1==omegajoveromega1(1,j); q(j,1)=0.5*(1/(j*omegaoveromega1))*((1/(j*j))*sin(j*omegat1)-(1/omegaoveromega1)*omegat1*cos(j*omegat1)); qprime(j,1)=0.5*(1/(j*omegaoveromega1))*(((1/j)-(1/omegaoveromega1))*cos(omegat1)+(j/(omegaoveromega1))*sin(omegat1)); else num1=(sin(j*omegat1)-(j*omegaoveromega1/omegajoveromega1(1,j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*omegat1)); q(j,1)=factor*(num1/((omegajoveromega1(1,j)^2)-(j*omegaoveromega1)^2)); num2=j*(cos(j*omegat1)-cos((omegajoveromega1(1,j)/omegaoveromega1)*omegat1)); qprime(j,1)=factor*num2/((omegajoveromega1(1,j)^2)-(j*omegaoveromega1)^2); end end for k=102:401 sum=zeros(1,101); for j=1:nmode sum=sum+(q(j,1)*cos((omegajoveromega1(1,j)/omegaoveromega1)*(omegat(1,k)-pi))+qprime(j,1)*(omegaoveromega1/(j*j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*(omegat(1,k)-pi)))*sine(j,:); end yoveryst=[yoveryst;sum]; end ytip=[]; for k=1:101 ytip=[ytip yoveryst(k,k)]; end ytip=[ytip zeros(1,300)]; figure(1);clf;%Animate the beam deflection as load moves across the beam h1=axes('Ydir','reverse'); axis([-.1 1.1 -2 2]) hold on plot(xoverL, yst); hold on plot([0 1],[0 0],'o'); hold on plot([-.1 0],[0 0]); hold on; plot([1 1.1],[0 0]); hold on text(-.05,1.7,['T_t/T_1 = ' num2str(Tratio)]); hold on A=plot(xarrow,yarrow,'LineWidth',[2],'EraseMode','xor'); Y=plot(xoverL,yoveryst(1,:),'Color','k','LineWidth',[3],'EraseMode','xor'); xlabel('Dimensionless Position, x/L') hold on ylabel('Deflection, 48EIy(x,t)/MgL^3') hold on box on Maximum=max(max(yoveryst)); text(.5,1.7,['y_m_a_x/y_s_t_a_t_i_c_m_a_x = ' num2str(Maximum)]); hold on t1=text(-.05,-1.7,'Press Enter to Continue'); pause set(t1,'String',' '); for k=1:401 set(A,'Xdata',xarrow+xoverL1(1,k)*ones(1,5),'Ydata',yarrow+ytip(1,k)*ones(1,5)); set(Y,'Ydata',yoveryst(k,:)); pause(.06) end hold off set(t1,'String','Press Enter to Continue'); pause MoverMst=[]; for k=1:101%time loop sum1=zeros(1,101); for j=1:5% series loop if j*omegaoveromega1==omegajoveromega1(1,j) term1=0.5*j*j(1/(j*omegaoveromega1))*((1/(j*j))*sin(j*omegat(1,k))-(1/omegaoveromega1)*omegat(1,k).*cos(j*omegat(1,k)))*sine(j,:); else num=(sin(j*omegat(1,k))-(j*omegaoveromega1/omegajoveromega1(1,j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*omegat(1,k))); term1=(j*j*num/((omegajoveromega1(1,j)^2)-(j*omegaoveromega1)^2))*sine(j,:); end sum1=sum1+term1; end MoverMst=[MoverMst;(8/(pi*pi))*sum1]; end for k=102:401 sum1=zeros(1,101); for j=1:5 sum1=sum1+(q(j,1)*cos((omegajoveromega1(1,j)/omegaoveromega1)*(omegat(1,k)-pi))+qprime(j,1)*(omegaoveromega1/(j*j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*(omegat(1,k)-pi)))*sine(j,:); end MoverMst=[MoverMst;sum1*(8/(pi*pi))]; end aoverL=xoverL; Mstatic=[zeros(1,101)]; for k=2:100 m1=4*(1-aoverL(1,k))*xoverL(1,1:k); m2=4*(aoverL(1,k)*ones(1,101-k)-aoverL(1,k)*xoverL(1,k+1:101)); Mstatic=[Mstatic;m1 m2]; end Mstatic=[Mstatic;zeros(1,101)]; MaxMoment=max(max(MoverMst)); figure(4);clf; axis([-.1 1.1 -2 2]) hold on plot([-.1 0],[0 0]); hold on; plot([1 1.1],[0 0]); hold on text(-.05,1.7,['T_t/T_1 = ' num2str(Tratio)]); hold on box on text(.5,1.7,['M_m_a_x/M_s_t_a_t_i_c_m_a_x = ' num2str(MaxMoment)]) B=plot([0 0],[-.5 -.1],'LineWidth',[1],'EraseMode','xor'); Y=plot(xoverL,MoverMst(1,:),'Color','k','LineWidth',[1],'EraseMode','xor'); MM=plot(xoverL,Mstatic(1,:),'EraseMode','xor'); xlabel('Dimensionless Position, x/L') hold on ylabel('Bending Moment, 4M(x,t)/MgL') hold on t1=text(-.05,-1.7,'Press Enter to Continue'); pause set(t1,'String',' '); for k=1:401 set(B,'Xdata',xoverL1(1,k)*ones(1,2)); set(Y,'Ydata',MoverMst(k,:)); if k>101 set(MM,'Ydata',zeros(1,101)) else set(MM,'Ydata',Mstatic(k,:)) end pause(.04) end hold off set(t1,'String','Press Enter to Continue'); pause;