Matlab tutorial

Mass-Spring-Damper system

The following problem is for the system consists of mass, spring and damper. This is exactly same as previous example, Matlab tutorial part 5 except we now have a damper. 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.

Problem Statement

Free Body Diagram and Equation of Motion

Closed Form Solution

Integrating EOM with Matlab

Plotting Results

Problem Statement

The below is the drawing of the mass-spring-damper system. We will solve for the homogeneous solution of the system, which means the external excitation force is zero. We will get the closed form of the solution.

In this example, we will use both closed form solution and ODE45 solver to compare the results. (If you have no idea what ODE45 solver is, then click here.)

Free Body Diagram and Equation of Motion

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, c is the viscous damping coefficient, and F(t) is the external force acting on the mass.

However, in the current problem we are focusing on the homogeneous solution, thus the external force is zero, i.e. F(t)=0. Then the above equation becomes

$\begin{displaymath}m \ddot{x} +c\dot{x} +kx = 0 \end{displaymath}$ (2)

and we have an equation of motion for the mass-spring system for the homogeneous solution. We can define the natural frequency of the system, $\omega_n$ , and the viscous damping factor, $\zeta$ , as

$\begin{displaymath}\omega_n^2 = \frac{k}{m}\end{displaymath}$ (3)

$\begin{displaymath}\zeta = \frac{c}{2m\omega_n}\end{displaymath}$ (4)

Then the equation of motion becomes

$\begin{displaymath}\ddot{x} + 2\zeta\omega_n \dot{x} + \omega_n^2 x = 0 \end{displaymath}$ (5)

Closed Form Solution

As in the mass-spring system we will solve for the closed form solution.

Now we need to find the closed form solution for the equation above. The solution will be assumed to be of the form

x(t)=X₀ e^st (1)

where X₀ is the amplitude. Note that

e ^ix = cos(x)+ i sin(x) (2)

Therefore, e^x is another way of representing sinusoidal functions.

Inserting equation (1) into the equation of motion gives

$\begin{displaymath}s^2 +2\zeta\omega_n s +\omega_n^2=0\end{displaymath}$ (3)

which has the roots

$\begin{displaymath}s_{1,2} = \left(-\zeta \pm \sqrt{\zeta^2-1}\right) \omega_n\end{displaymath}$ (4)

Clearly, the roots depend on the value of $\zeta$ .

If $\zeta \neq 1$ the solution form is

x(t) = X₁ e^s₁^t + X₂ e^s₂^t (5)

$\begin{displaymath}= \left[ X_1 e^{\omega_n t\sqrt{\zeta^2-1}}+ X_2 e^{-\omega_n t\sqrt{\zeta^2-1}} \right] e^{-\zeta \omega_n t}\end{displaymath}$ (6)

When $0<\zeta<1$ , the underdamped case, x(t) will become complex since there is nothing preventing the roots from being complex. The imaginary part of the solution is as important as the real part. The solution will the be of the form

$\begin{displaymath}\mbox{solution}=Re\{x(t)\}+i\;Im\{x(t)\}\end{displaymath}$ (7)

To calculate the amplitudes from the initial conditions we can use

x₀ = X₁ + X₂ (8)

$\begin{displaymath}\dot{x}_0 = s_1 X_1 + s_2 X_2\end{displaymath}$ (9)

Since the initial conditions are known, these two equations can be solved for the amplitudes.

When $\zeta=1$ (critical damping) the solution form is

X₁ = x₀ (10)

$\begin{displaymath}X_2 = X_1\omega_n =\dot{x}_0 + x_0 \omega_n\end{displaymath}$ (11)

Thus the above equations yield the complete solution for any value of damping.

The following is the MATLAB code for closed form solution of mass-spring-damper system. It is quite simple and straightforward.

First, as we did in mass-spring system, determine the constant input values and calculate some relevant variables.

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 = 3;                %  calculation time
n = 500;                    %  number of data points

s1 = (-zeta+sqrt(zeta^2-1))*wn;
s2 = (-zeta-sqrt(zeta^2-1))*wn;

Initial conditions

x_0 = 1;                    %  initial displacement
x_dot_0 = 0;                %  initial velocity
x0 = [x_0,x_dot_0];

We need to check whether viscous damping coefficient zeta becomes 1 (critically damped case) or not since the above equations will give singular matrix of A when zeta = 1. We first determine the amplitudes of x(t), X_1, X_2 with given initial conditions, x_0, x_dot_0. The reason we separate x(t) in real and imaginary number is explained above deriving EOM.

t = 0:t_final/n:t_final;

if zeta ~= 1 

  A = [[1,1];[s1,s2]];
  X = A^(-1)*x0';

  X_1 = X(1);
  X_2 = X(2);

  x = real(X_1*exp(s1*t)+X_2*(exp(s2*t)))+...
      imag(X_1*exp(s1*t)+X_2*(exp(s2*t)));               %  displacement x(t)
  x_dot = real(s1*X_1*exp(s1*t)+s2*X_2*(exp(s2*t)))+...
          imag(s1*X_1*exp(s1*t)+s2*X_2*(exp(s2*t)));     %  velocity x_dot(t)

If zeta = 1, then we use the above equations to get displacement and velocity.

elseif zeta == 1

  X_1 = x_0;
  X_2 = x_dot_0+X_1*wn;
  x = (X_1+X_2.*t).*exp(-wn*t);
  x_dot = (-wn*(X_1+X_2.*t)+X_2).*exp(-wn*t);

end

Click here for the entire MATLAB source code for closed form solution.

Integrating EOM with Matlab

Here we will use the ODE45 solver in Matlab to integrate the equations of motion. 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.

Some constant inputs need to be fixed before we proceed.

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

And the initial conditions are shown. tol is the tolerance for ODE45 solver. Higher it is, more accurate the answer will be, more time the calculation takes.

x_0 = 1;                    %  initial displacement
x_dot_0 = 0;                %  initial velocity
X_0 = [x_0,x_dot_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 x_dot(t). y is used as a dummy variable.

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

x = y(:,1);                 %  displacement x(t)
x_dot  = y(:,2);            %  velocity x_dot(t)

Here's the subprogram pmsd_sol.m.

function ydot = pmsd_sol(t,y)

m = 1;                      %  mass of the system in meter
k = 50;                     %  stiffness of the spring (N/m)
c = 20;                     %  damping coefficient

yd1 = y(2);
yd2 = -(k*y(1)+c*y(2))/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 (PC,Mac) , main program (Unix) and subprogram.

Plotting Results

Results are given in 2 parts. The first gives the plotting commands and results using the closed form solution, the second gives the plotting commands and results using the ODE solver. Just as you may have guessed, both methods should give the same answers. Check them out.

Return to Example Table