m = 2; % mass of the system c = 30; % damping of the system k = 5000; % stiffness of the system P0 = 1e2; % magnitude of the excitation force Omega = 20; % frequency of the excitation force t_final = 3; % final time of calculation n = 5000; % number of data points wn = sqrt(k/m); % natural frequency (rad/sec) zeta = c/(2*wn*m); % viscous damping factor
Input some constants of the system, like mass,
stiffness, etc. Calculate the natural frequency of the system and the damping
coefficient, zeta. P0, Omega are the amplitude and the excitation frequency
of the external harmonic force, repectively.
wd = wn*sqrt(1-zeta^2); disp 'damped natural frequency is' % damped natural frequency wd disp 'damping coefficient is' % viscous damping factor zeta disp 'static displacement is' % static displacement P0/k
Displaying the damped natural frequency, damping
coefficient and the static displacement of the spring.
% % initial conditions (for total response x(t)) % x(t) = x_h(t) + x_p(t) % x_0 = 0; % initial displacement xdot_0 = 0; % initial velocity t = 0:t_final/n:t_final; % make a vector for calculation times % % Particular solution, x_p(t) % psi_s = atan(c*Omega/(k-m*Omega^2)); % phase angle of x_p(t) D = P0/(sqrt((k-m*Omega^2)^2+(c*Omega)^2)); % amplitude of x_p(t) xp = D*sin(Omega*t-psi_s); % particular solution
The above is for the closed form solution for
the particular part. 'psi_s' is the phase angle and 'D' is the amplitude
of the particular solution. If you miss the derivation of those equations
above, click here.
% % Homogeneous solution, x_h(t) % % underdamped case % if zeta < 1 B = x_0+D*sin(psi_s); % initial condition for x_h(0) % note that x_h(0) = x(0)-x_p(0) C = (xdot_0-D*Omega*cos(psi_s)+zeta*wn*B)/wd; A = sqrt(((B*wd)^2+(C+zeta*wn*B)^2)/wd^2); % amplitude of x_h(t) phi_c = atan((C+zeta*wn*B)/(wd*B)); % phase angle of x_h(t) xh = A*cos(wd*t-phi_c).*exp(-zeta*wn*t); % displacement x_h(t) % overdamped case % elseif zeta > 1 D1 = (xdot_0-D*Omega*cos(psi_s)+wn*(zeta+sqrt(zeta^2-1))*... (x_0+D*sin(psi_s)))/(2*sqrt(zeta^2-1)*wn); % initial condition for x_h(0) % note that x_h(0) = x(0)-x_p(0) D2 = x_0+D*sin(psi_s)-D1; xh = (D1.*exp(wn*sqrt(zeta^2-1)*t)+D2.*exp(-wn*sqrt(zeta^2-1)*t))... .*exp(-zeta*wn*t); % displacement x_h(t) % critically damped case % elseif zeta == 1 B = x_0+D*sin(psi_s); % initial condition for x_h(0) C = xdot_0-D*Omega*cos(psi_s)+B*wn; xh = (B+C.*t).*exp(-wn*t); % displacement x_h(t) end
Note that 'x_0' and 'xdot_0' are the initial conditions
of total response. Therefore, when you calculate the homogeneous response,
you should be careful to substitute the correct initial conditions.
The above is the homogeneous solution. The
derivation of the above equations is here
%
% plotting the results
%
figure(1);clf;
subplot(3,1,1),plot(t,xh);
h=legend('homogeneous response');
title('Forced vibration on Mass-Spring-Damper system')
subplot(3,1,2),plot(t,xp,'m');
h=legend('particular response');
axis([0 t_final min(xp)*1.5 max(xp)*2.5]);
subplot(3,1,3),plot(t,xp+xh,'r');
h=legend('total response');
xlabel('time (sec.)');
Plotting the results, the homogeneous solution,
the particular solution and the total solution.
Click here for the entire Matlab source code for the solution including plotting.