Mass-Spring-Damper system : Forced Vibration

Problem Statement

The following problem is for the system consists of mass, spring and damper. This is exactly same as a previous example, Damped Free Vibration except we now have an external force. The procedure will be same as before: FBD first, then the equation of motion using FBD and force balance, Matlab coding and finally some results.

The below is the drawing of the mass-spring-damper system.

Free Body Diagram and Equation of Motion

We can write down the equation of motion for this 1 degree of freedom system with FBD above using force balance.

From the free body digram we can write the force balance equations as follow. Consider the right-going direction as the positive x-direction to get

\begin{displaymath} -c\dot{x}-kx+F(t)=m \ddot{x}\end{displaymath} (1)

where m is the mass of the system, k is the stiffness of the spring, and c is the viscous damping coefficient F(t) is the excitation force acting on the mass. In the current problem, the external force is a harmonic function $F(t)=P_0\mbox{sin}(\Omega t)$.Thus the above equation becomes

\begin{displaymath} m \ddot{x} +c\dot{x} +kx = P_0 \mbox{sin}(\Omega t) \end{displaymath} (2)

Closed Form Solution

We will solve for the homogeneous solution and the particular solution of the system. The particular solution is necessary here since the external force is not zero as in the free vibrations examples. Both the homogeneous and particular solutions together will be the closed form solution.

The total response of the system with the external force is the sum of the homogeneous and particular solutions.

x(t)=xh(t)+xp(t)

(1)

We already have the closed form solution to the homogeneuos problem from the previous MATLAB tutorial, therefore we will focus on the particular solution here.

Any function that satisfies the equation of motion is the particular soluton. Since the external force is harmonic, it is reasonable to guess that the particular solution is also harmonic.

\begin{displaymath} x_p(t)=D\mbox{sin}(\Omega t - \psi_p)\end{displaymath} (2)

Here D and $\psi_p$ will be determined to satisfy the equation of motion. substituting the assumed solution form into the governing equation yields

\begin{displaymath} \mbox{tan}\psi_p = \frac{c\Omega}{k-m\Omega^2}\end{displaymath} (3)
\begin{displaymath} D = \frac{P_0}{\sqrt{\left(k-m\Omega^2\right)^2+\left(c\Omega\right)^2}}\end{displaymath} (4)

Since the amplitude of the particular solution is constant, the particular solution is also called the steady-state solution. After the transient part of the solution (homogeneous) is damped out, the system will oscillate according to xp(t) so long as the driving force $P_0\mbox{sin}(\Omega t)$ is applied. Thus the steady state solution can be written as

\begin{displaymath} x_p(t) = \frac{P_0}{\sqrt{\left(k-m\Omega^2\right)^2+\left(c... ...ox{tan}^{-1} \left( \frac{c\Omega}{k-m\Omega^2} \right)\right]\end{displaymath} (5)

The closed form solution can be coded into Matlab. The first solution consists of both the homogeneous (transient) solution plus the particular (steady-state) solution. The second is just the steady-state solution

Integrating EOM with Matlab

Once again we will also use Matlab to intgrate the equation of motion. The only difference from the free vibration problems is that the state space equation for position also contains the force term on the right hand side. The following is the MATLAB code for ODE45 solver solution of mass-spring system. If you need to see the basic operation of ODE45 solver, see Matlab Assistance. 

m = 1;                      %  mass of the system in meter
k = 50;                     %  stiffness of the spring (N/m)
c = 20;                     %  damping coefficient
wn = sqrt(k/m);             %  natural frequency (rad/sec)
zeta = c/(2*wn*m)           %  viscous damping factor

t_final = 2;                %  calculation time
The initial conditions are shown. tol is the tolerance for ODE45 solver. The higher it is, the more accurate the answer will be, but the calculation takes more time. 
x_0 = 1;                    %  initial displacement
xdot_0 = 0;                %  initial velocity
X_0 = [x_0,xdot_0];        %  form a vector(array) of initial conditions
tol = 1e-6;                 %  tolerance of error for ODE45 solver
Using ODE45 solver, we can solve for the displacement x(t) and velocity xdot(t). 
t_final = 3;                %  calculation time

[t,y] = ode45('forcedsub',0,t_final,X_0,tol);

x = y(:,1);                 %  displacement x(t)
x_dot  = y(:,2);            %  velocity x_dot(t)
Here's the subprogram forcedsub.m. 
function ydot = forcedsub(t,y)

m = 2;			%  mass of the system in meter
k = 5000;		%  stiffness of the spring (N/m)
c = 30;			%  damping coefficient
P = 100;		%  magnitude of force
Omega = 20;		%  frequency of force
	
yd1 = y(2);
yd2 = (-k*y(1)-c*y(2)+P*sin(Omega*t))/m;

ydot = [yd1;yd2];

NOTE that we did not discretize the critically damped case from underdamped case and overdamped case since this ODE45 solver takes care of that automatically, i.e. we did not have singular matrix as we had with closed form solution.

Here are the complete list of the main program and subprogram.

Plotting Results

Results shows the homogeneous response (transient response), the particular response (steady-state response) and the sum of those two ( total response).
Also, the results for the frequency-amplitude and frequency-phase angle solution.

Return to Example Table