home
PDF
Simulink files to download
(a)
project_step.mdl
(b)
project_sin.mdl
(c)
project_impulse.mdl
(d)
project.mdl

Using Simulink to analyze 2 degrees of freedom system

Nasser M. Abbasi

Spring 2009 page compiled on June 29, 2015 at 4:20pm

Abstract

A two degrees of freedom system consisting of two masses connected by springs and subject to 3 different type of input forces is analyzed and simulated using Simulink

Contents

1 Introduction and Theory
2 Analytical solution
 2.1 Finding the homogenous solution
 2.2 Finding particular solutions
  2.2.1 Finding the particular solution for unit impulse input
  2.2.2 Finding the particular solution for unit step input
  2.2.3 Finding the particular solution for sin ωt
3 Simulink simulation and block diagrams
 3.1 Unit step simulink diagram and output
  3.1.1 Verification of result from Simulink by Numerically solving the differential equations
4 Unit impulse simulink diagram and output
  4.0.1 Verification of result from Simulink by Numerically solving the differential equations
 4.1 sin (ωt)  input simulink diagram and output
5 Discussion

1 Introduction and Theory

The system that is being analyzed is show in the following diagram

pict

In the above, F  (t)  is to be taken as each of the following

1.
Unit impulse force.
2.
Unit step force.
3.
sinωt

It is required to find x1 (t)  and x2(t)  analytically and then to use Matlab's Simulink software for the analysis.

The mathematical model of the system is first developed and the equation of motions obtained using Lagrangian formulation then the analytical solution is found by solving the resulting coupled second order differential equations for m
  1   and m
  2   . Next, a simulink model is developed to implement the differential equations and the output x1(t)  and x2 (t)  from Simulink is shown and compared to the output from the analytical solution.

2 Analytical solution

The following is the free body diagram of the above system

pict

Assuming positive is downwards and that x  >  x
  2    1   , force-balance equations for m
  1   results in

m1 ¨x1 = F (t) + k1 (x2 − x1)

And force-balance equations for m2   results in

m2 ¨x2 = − k1(x2 − x1) − k2x2

Hence the EQM for the system become

pict

Or in matrix form

(         ) (   )    (               ) (   )    (     )
  m1    0     ¨x1       k1      − k1      x1       F (t)
                  +                          =
   0   m2     ¨x2      − k1  (k1 + k2)    x2        0

The above can be written in matrix form as

M ¨x + Kx   = F

Where ¨x,x, F  are 2 by 1 vectors and M  and K  are the mass and stiffness matrices. The solution to the above is

x =  xh + xp
(1)

2.1 Finding the homogenous solution

We start by finding xh  from the following

(m1     0 ) ( ¨x1) + (  k1     − k1   ) (x1 ) =  (0)
   0   m2     ¨x2      − k1  (k1 + k2)    x2      0

Now assume x1(t) = A1 cos(ωt + ϕ)  and x2(t) = A2 cos(ωt + ϕ )  , hence x˙1 (t) = − ωA1 sin (ωt + ϕ)  and x˙2 (t) = − ωA2  sin (ωt + ϕ)  and             2
¨x1 (t) = − ω A1 cos(ωt + ϕ)  and            2
˙x2(t) = − ω A2 cos(ωt + ϕ)  . Substituting the above values in the above system results in

(        ) (                    )   (                ) (               )    (  )
  m1   0     − ω2A1 cos (ωt + ϕ)       k1     − k1       A1 cos(ωt + ϕ)       0
                 2                +                                       =
   0  m2     − ω A2 cos (ωt + ϕ)      − k1  (k1 + k2)    A2 cos(ωt + ϕ)       0

Divide by cos(ωt + ϕ)  since not zero (else no solution exist) we obtain

(        ) (        )   (                ) (   )    (  )
  m1   0     − ω2A1        k1     − k1       A1       0
  0   m      − ω2A    +   − k   (k +  k )    A    =   0
        2          2         1    1    2      2

Rewrite the above as

pict

From the last equation above, we see that to obtain a solution we must have

|                                |
|− ω2m   + k          − k        |
||      1    1            1       || = 0
|    − k1      − ω2m2 +  (k1 + k2)|

since if we had (A1A2 ) =  0  then no solution will exist. Therefore, taking the determinant and setting it to zero results in

pict

Let  4    2
ω  = λ   , hence the above becomes

 2                                     ( 2           2)
λ m1m2  + λ (− k1m1 − k2m1  − k1m2 ) +  k1 + k1k2 − k1 =  0

Solving for λ  gives

           (                                                                                   )
λ  =  --1--- k m  + k  m  + k m   + ∘k2-m2--+-k2m2-+--k2m2-+-2k--k-m2-+-2k2-m-m---−-2k-k-m--m--
 1    2m1m2 (  1  1    1 2    2  1   ∘ -1--1----1--2----2--1-----1-2--1-----1--1-2-----1-2--1--2)
λ2 =  --1--- k1m1 + k1m2  + k2m1  −   k2m2  + k2m2 +  k2m2 + 2k1k2m2  + 2k2 m1m2  − 2k1k2m1m2
      2m1m2                             1  1    1  2    2  1          1     1
(3)

For each of the above solutions, we obtain a different (    )
  A1

  A2 from equation (2) as follows

For λ1   , (2) becomes

pict

From the first equation above, we have

pict

Similarly for λ2   ,

− λ m  + k     (A  ) (2)
---2--1---1-=   --2
    k1          A1

Let

pict

Hence now xh  can be written as

(   )     (  (1)                  (2)             )
  x1       A 1  cos(ω1t + ϕ1) + A1  cos(ω2t + ϕ2)
  x    =   A (1) cos(ω t + ϕ ) + A(2)cos(ω t + ϕ )
   2  h      2       1     1     2       2     2

But   (1)      (1)
A 2  = r1A 1   and  (2)      (2)
A2  = r2A 1   , hence the above becomes

(   )     (    (1)                  (2)               )
  x1         A 1  cos(ω1t + ϕ1) + A1  cos(ω2t + ϕ2)
  x    =   r A (1) cos(ω t + ϕ ) + rA (2)cos(ω t + ϕ )
   2  h     1  1       1     1     2 1       2     2
(5)

Now, given numerical values for k1,k2,m1, m2   we can find ω1,ω2   from (3) above, and next find r1,r2   from (4).  Hence (5) contains 4 unknowns,   (1)  (2)
A 1 ,A 1 ,ϕ1,ϕ2   which now can be found from initial conditions (after we find the particular solution) which we will now proceed to do.

2.2 Finding particular solutions

There are 3 different F (t)  which we are asked to consider

1.
Unit impulse force.
2.
Unit step force.
3.
sinωt

For each of the above, we find xp  and then add it to xh  found above in (5) to obtain (1).

2.2.1 Finding the particular solution for unit impulse input

Using the standard response for a unit impulse which for a single degree of freedom system is          1
x (t) = mωn-sin ωnt  , then we write xp  as

      (   )    (  1             1       )
       x1        m-ω1 sinω1t + m-ω2 sinω2t
xp =   x     =              0
         2 p

Hence, the general solution becomes

(   )    (     (1)                  (2)                )   (                        )
  x1         A 1 cos (ω1t + ϕ1 ) + A 1 cos (ω2t + ϕ2 )      m1ω1 sin ω1t + m1ω2 sinω2t
  x    =   r A (1)cos (ω  t + ϕ ) + r A (2)cos (ω t + ϕ ) +              0
    2       1  1       1     1    2  1       2     2
(6)

2.2.2 Finding the particular solution for unit step input

Since unit step is 1  for t > 0  , then, using convolution we write

pict

Then, since now we have 2 natural frequencies, we can write xp  as

     (   )     (                                       )
       x1        m1ω2[1 − cos(ω1t)] + m-1ω2[1 − cos(ω2t)]
xp =         =      1                   2
       x2  p                       0

Hence, the general solution becomes

(   )    (    (1)                  (2)               )    ( -1--                -1--             )
  x1   =     A1  cos(ω1t + ϕ1) + A 1 cos (ω2t + ϕ2 )    +   mω21 [1 − cos (ω1t)] + m ω22 [1 − cos(ω2t)]
  x2       r1A(1)cos(ω1t + ϕ1) + r2A (2)cos(ω2t + ϕ2)                         0
              1                      1

2.2.3 Finding the particular solution for sin ωt

In this case, we guess that x1p = c1cosωt + c2 sin ωt  , and since there is no forcing function being applied directly on m2   then x2p = 0  hence

      (      )    (                   )
       x1 (t)       c1cos ωt + c2sin ωt
xp =            =
       x2 (t) p              0

Then ˙x1p(t) = − ωc1 sinωt + ωc2 cosωt  and ¨x1p (t) = − ω2c1 cosωt + ω2c2 sin ωt  and now we substitute these into the original ODE for x1   which is

pict

We obtain the following

pict

Hence by comparing coefficients, we obtain

pict

or

pict

c
 1   must be zero since k  − ω2m   = 0
  1       1  only when ω  = ω
      n  and we assume that this is not the case here. Hence

pict

Therefore xp  becomes

      (     )     (               )
       x1 (t)       (m1ω12+k1) sinωt
xp =            =
       x2 (t) p            0

And the general solution becomes

(   )    (    (1)                 (2)               )    ( ---1----sin ωt)
  x1  =     A 1  cos(ω1t + ϕ1) + A1  cos(ω2t + ϕ2)    +    (m1ω2+k1)
  x2      r1A (11)cos(ω1t + ϕ1) + r2A(12)cos(ω2t + ϕ2)             0
(8)

3 Simulink simulation and block diagrams

In simulink, we will directly solve the system from the original formulation

pict

or

pict

Hence

pict

3.1 Unit step simulink diagram and output

The simulink block diagram will be as follows for the unit step input

PIC

For an initial run with parameters m1  = m2  = k1 = k2 = 1  I get this warning below

  EDU>> simulink
  Warning: Using a default value of 0.2 for maximum step size.  The simulation
  step size will be equal to or less than this value.  You can disable this
  diagnostic by setting 'Automatic solver parameter selection' diagnostic to
  'none' in the Diagnostics page of the configuration parameters dialog.

And this is the output for x1(t)  and x2 (t)  for the unit step response

PIC

3.1.1 Verification of result from Simulink by Numerically solving the differential equations

To verify the above output from Simulink, I solved the same coupled differential equations for zero initial conditions numerically (using a numerical differential equation solver) and plotted the solution for x1(t)  and x2 (t)  and the result matches that shown above by simulink. Here is the code the plot as a result of this verification

pict

pict

4 Unit impulse simulink diagram and output

pict

And the output for x1 (t)  and x2 (t)  is as follows

pict

4.0.1 Verification of result from Simulink by Numerically solving the differential equations

To verify the above output from Simulink, The same coupled differential equations were solved numerically for zero initial conditions numerically and the solution plotted for x1 (t)  and x2 (t)  and the result was found to match that shown above by simulink. Here is the code used to do the verification

verification

pict

4.1 sin (ωt )  input simulink diagram and output

The simulink block diagram will be as follows for the sin ωt  input

PIC

For an initial run with parameters m1 =  m2 =  k1 = k2 = 1  this is the output for x1 (t)  and x2(t)  and showing the input signal at the same time

pict

5 Discussion

A coupled system of two masses and springs was analyzed using Simulink. The simulation was done for one set of parameters (masses and stiffness). Simulink made the simulation of this system under different loading conditions easy to do. The 2 masses response were recorded using simulink scope and the signals captured on the same plot to make it easy to compare the response of the first mass to the second mass.

The analytical analysis was more time consuming than actually making the simulation in simulink. The ability to easily change different sources to the system was useful as well as the ability to change the frequency of the input and immediately see the effect on the response.

This was my first project using Simulink, and I can see that this tool will be useful to learn more as it allows one to easily analyze engineering problems.