v1_003.m

% v1_003.m % % This is a script file to solve a sdof system, % given the mass, damping and stiffness terms % in dimensionless units, % % The method uses a Runge-Kutta approach without the use % of any matlab functions (ODE23 or ODE45). % % Reference: Theory of Vibration with Applications % 4th Edition, 1993 % William T. Thomson % % Consider the second order differential equation for a single % degree of freedom vibrating system % % mass * x'' + damp * x' + stiff * x = f % % where the prime ' denotes the derivative operator. % We can rewrite this as a system of first order differential % equations: % % A = x'' % V = x' % X = x % % V' = x'' = 1/mass( f - damp * V - stiff * X ) %********************************************************************** % Author: Randall J. Allemang % Date: 25-Aug-9€7 % Structural Dynamics Research Lab % University of Cincinnati % Cincinnati, Ohio 45221-0072 % TEL: 513-556-2725 % FAX: 513-556-3390 % E-MAIL: randy.allemang@uc.edu %********************************************************************* % clg, clear plt=input('Store plots to file (Yes=1): (0)'); if isempty(plt),plt=0;end; mass=input('Mass value: (20)'); if isempty(mass), mass=20;end; stiff=input('Stiffness value: (9000)'); if isempty(stiff),stiff=9000;end; damp=input('Damping value: (130)'); if isempty(damp),damp=130;end; % Define time limits and increment tinit = 0; tfinal = 2; t=linspace(tinit,tfinal,501); dt=t(2)-t(1) % Define initial conditions xinit = 0; vinit = 0; % Clear force function for ii=1:501 f(ii)=0; end % Define transient excitation (sine pulse) time_dur=0.1; numb_int=round(time_dur/dt) for ii=1:numb_int f(ii)=15*sin(2.*pi./(2.*time_dur).*t(ii)); end % Define transient excitation (constant-ramp pulse) % for ii=1:21 % f(ii)=5; % end % for ii=22:46 % f(ii)=5-(ii-21)*dt*50; % end fig1=figure(1); plot(t,f), title('SDOF time history (force)'),grid,pause if plt==1, print -ddeskjet -f1 v1_003a.eps, end; % Start Runge-Kutta approximation x(1)=xinit; v(1)=vinit; a(1)=(f(1)-stiff.*x(1)-damp.*v(1))./mass; for ii=2:501 T1=t(ii-1); X1=x(ii-1); V1=v(ii-1); A1=(f(ii-1)-stiff.*X1-damp.*V1)./mass; T2=t(ii-1)+dt/2; X2=x(ii-1)+V1*dt/2; V2=v(ii-1)+A1*dt/2; A2=(((f(ii-1)+f(ii))./2)-stiff.*X2-damp.*V2)./mass; T3=t(ii-1)+dt/2; X3=x(ii-1)+V2*dt/2; V3=v(ii-1)+A2*dt/2; A3=(((f(ii-1)+f(ii))./2)-stiff.*X3-damp.*V3)./mass; T4=t(ii); X4=x(ii-1)+V3*dt; V4=v(ii-1)+A3*dt; A4=(f(ii)-stiff.*X3-damp.*V3)./mass; x(ii)=x(ii-1) + (V1+2*V2+2*V3+V4)*dt/6; v(ii)=v(ii-1) + (A1+2*A2+2*A3+A4)*dt/6; a(ii)=(f(ii)-stiff.*x(ii)-damp.*v(ii))./mass; end % Plot results fig2=figure(2); plot(t,x), title('SDOF time history (displacement)'),grid,pause if plt==1, print -ddeskjet -f2 v1_003b.eps, end; fig3=figure(3); plot(t,v), title('SDOF time history (velocity)'),grid,pause if plt==1, print -ddeskjet -f3 v1_003c.eps, end; fig4=figure(4); plot(t,a), title('SDOF time history (acceleration)'),grid,pause if plt==1, print -ddeskjet -f4 v1_003d.eps, end; fig5=figure(5); plot(x,v), title('SDOF equation - phase plane plot'),grid,pause if plt==1, print -ddeskjet -f5 v1_003e.eps, end;