% stringspring.m Last modified 8/4/2008 % A taut string with tension tau and mass per unit length mu is % terminated at x=L with a transverse spring of stiffness k. The solution % is accomplished via separation of variables and generalized Fourier % series. Note that the eigenvalues are not the same as those for % the string fixed ends, they are smaller.Input is the ratio % K = k*L/Tension. % 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 disp('A taut string with tension tau and mass per unit length mu is') disp('terminated at x=L with a transverse spring of stiffness k.') disp('The solution is accomplished via separation of variables') disp('and generalized Fourier series. Note that the eigenvalues') disp('are not the same as those for the string fixed ends, they are') disp('smaller. Input is the ratio K = k*L/Tension.') disp(' ') clear all set(0,'DefaultAxesFontSize',12); set(0,'DefaultTextFontSize',12); K=-1; while (K<.2)|(K>5) disp('Choose a value between 0.2 and 5') K=input('Input the stiffness parameter k*L/Tension = '); end omega1t=linspace(0,20,201); Nmode=20; lambdaL=zeros(1,Nmode); z0=[]; for i=1:2:2*Nmode+1 z0=[z0 i*(pi/2)+(.2/i)]; end for j=1:Nmode e=1; z=z0(1,j); while (abs(e)>1e-5) z1=z-((K*tan(z)+z)*((cos(z))^2))/(K+(cos(z)^2)); e=z-z1; z=z1; end lambdaL(1,j)=z1; end phi=zeros(Nmode,101); xoverL=linspace(0,1,101); alpha=zeros(1,Nmode); for i=1:Nmode phi(i,:)=sin(lambdaL(1,i)*xoverL); %Normalization Const. alpha(1,i)=((2*lambdaL(1,i))-(sin(2*lambdaL(1,i))))/(4*lambdaL(1,i)); end % expand Initial condition b=zeros(1,Nmode); for i=1:Nmode b(1,i)=(xoverL*phi(i,:)'-0.5*xoverL(1,101)*phi(i,101))*.01/alpha(1,i); end %Calculate the response array omega1t=linspace(0,20,201); omegaioveromega1=lambdaL/lambdaL(1,1); framedata=zeros(202,101); framedata(1,:)=zeros(1,101); framedata(2,:)=xoverL; for t=3:202 sum=zeros(1,101); for i=1:Nmode sum=sum+b(1,i)*phi(i,:)*(cos(omegaioveromega1(1,i)*omega1t(1,t-1))); end framedata(t,:)=sum; end bpatchx=[0 0 -.15 -.1 -.16 0]; bpatchy=[1.4 -1.4 -.3 .2 .7 1.4]; rpatchy=[1.2 1.2 1.6 1.8 1.9 1.2]; rpatchx=[1.01 1.04 1.03 1.04 1.01 1.01]; spatchx=[.9 1.1 1.05 .97 .9]; spatchy=[6 6 6.2 6.3 6]; lpatchy=rpatchy; lpatchx=[.99 .96 .97 .96 .99 .99]; xspring=[ 1 1 1.1 .9 1.1 .9 1 1]; y1=framedata(:,101); y2=framedata(:,101)*.875; y3=framedata(:,101)*.625; y4=framedata(:,101)*.375; y5=framedata(:,101)*.125; ylabeltext='Displacement, y(x,t)/y_0'; figure(1);clf h1=axes('Ydir','reverse'); axis([-.2 1.2 -1.5 6.5]); hold on box on xlabel('Distance, x/L'); ylabel(ylabeltext); patch(bpatchx,bpatchy,'r'); hold on patch(rpatchx,rpatchy,'r'); hold on patch(lpatchx,lpatchy,'r'); hold on patch(spatchx,spatchy,'r') hold on t=1; P=plot(xoverL,framedata(1,:),'k','EraseMode','xor','LineWidth',2); yspring=[y1(1,1) 3+y1(1,1) 3.25+y2(1,1) 3.75+y3(1,1) 4.25+y4(1,1) 4.75+y5(1,1) 5 6]; S=plot(xspring, yspring,'k','EraseMode','xor'); text(0,3,['\omega_1 = (' num2str(lambdaL(1,1)) '/L)*sqrt(\tau/\mu)']); text(0,3.4,'\tau = string tension'); text(0,3.8,'\mu = mass per unit length'); text(0,4.2,['kL/\tau = ' num2str(K)]); text(.81,4.5,'k'); message='Press Enter to Set Initial Condition'; T=text(0,5,message); pause set(P,'Ydata',framedata(2,:)); yspring=[y1(2,1) 3+y1(2,1) 3.25+y2(2,1) 3.75+y3(2,1) 4.25+y4(2,1) 4.75+y5(2,1) 5 6]; set(S, 'Ydata',yspring); message='Press Enter to Animate '; set(T,'String',message); pause message=' '; set(T,'String', message); for t=3:202 set(P,'Ydata',framedata(t,:)); yspring=[y1(t,1) 3+y1(t,1) 3.25+y2(t,1) 3.75+y3(t,1) 4.25+y4(t,1) 4.75+y5(t,1) 5 6]; set(S, 'Ydata',yspring); pause(.04) end message='Press Enter to Continue'; set(T,'String',message); hold off pause figure(2);clf h1=axes('Ydir','reverse'); axis([0 20 -1.5 1.5]); hold on box on k=0; for j=21:20:101 k=k+1; plot(omega1t,framedata(2:202,j)); hold on plot([.1 1.4],[framedata(3,j) framedata(3,j)+.3]); text(1.45,framedata(3,j)+.35,['x/L = ' num2str(.2*k)]); end grid on text(12,1.2,['kL/\tau = ' num2str(K)]); xlabel('Dimensionless Time, \omega_1t'); ylabel('Displacemant, y(x,t)/y_0'); hold off text(3,-1.3,'Press Enter to Continue'); pause figure(3);clf h1=axes('Ydir','reverse'); axis([0 10 -1.5 1.5]); hold on box on k=0; for j=21:20:101 k=k+1; plot(omega1t(1,1:101),framedata(2:102,j)); hold on plot([.1 1.4],[framedata(3,j) framedata(3,j)+.3]); text(1.45,framedata(3,j)+.35,['x/L = ' num2str(.2*k)]); end grid on text(6.1,1.2,['kL/\tau = ' num2str(K)]); xlabel('Dimensionless Time, \omega_1t'); ylabel('Displacemant, y(x,t)/y_0'); hold off text(2,-1.3,'Press Enter to Continue'); pause figure(4);clf h1=axes('Ydir','reverse'); axis([0 1 -1.5 1.5]) hold on box on grid on for j=2:5:42 plot(xoverL, framedata(j,:)); plot([.77 .84],[framedata(j,78) 1.4-(j-2)*.04]); text(.85, 1.4-(j-2)*.04,['\omega_1t = ' num2str(.01*(j-2))]); xlabel('Distance Down String, x/L'); ylabel('Displacement, y(x,t)/y_0'); end text(.1,1.2,'Press Enter to Continue'); hold off pause figure(5) [X,Y]=meshgrid(xoverL,omega1t(1,1:101)); mesh(X,Y,framedata(2:102,:)) xlabel('Distance, x/L','Rotation',-34) ylabel('Time, \omega_1t','Rotation',12) zlabel('Displacement, y(x,t)/y_0') axis([0 1 0 10 -1.2 1.2]) view(60,30) text(0,2,.9,'Press Enter to Continue','Rotation',12); pause