Departments of Mechanical Engineering and Mathematics
Mathematica Session 3

Damped and Forced Vibrations

In Session 2, you generated animations for undamped, free vibrations in 1 and 2-dimensions. In reality, damping, in the form of friction and air resistance, is rarely absent from a mechanical system, and in fact, the deliberate use of damping to control and, if possible, eliminate vibration in a system is a vital technique in engineering design.

You should be already familiar with the simple harmonic oscillator, which models small vibrations of springs (where the restoring force is proportional to displacement) and pendulums, in the absence of frictional effects. In this session we're going to do some "experiments" on oscillators where friction is modelled as a damping force that is proportional to velocity -- which is actually a very good model for mechanical devices like car shock absorbers. Then, in part 2, we'll look briefly at another aspect of real vibrating systems,  forced oscillations.

1. The Damped Oscillator

Recall that the equation of motion for a simple harmonic oscillator (SHO) is:

m(((d^( 2)) x)/(d (t^( 2))))=-k x ,    or, rearranging:    m(((d^( 2)) x)/(d (t^( 2)))) + k x =0 .

We are going to add to this equation a damping force that is proportional to velocity and opposes the motion, that is:

m(((d^( 2)) x)/(d (t^( 2)))) +c(( d x)/(d t)) + k x =0

This will be our mathematical model for damped oscillations. The physical system which it models can be treated as a  mass (m) sliding on a smooth floor and connected to a fixed support by a spring (with spring constant k) and a dashpot (idealised shock absorber), with "damping constant" c:

The above differential equation is soluble using analytic techniques that you'll see later on in your mathematics course. Here, though, we want to use Mathematica's numerical solver, NDSolve, to experiment with the behaviour of the damped oscillator as the relative values of the spring constant, k, damping constant, c, and mass, m, are varied.

Type in the following procedure:

\(m = 1; 
c = 0; 
k = 1; 

    solution = 
      NDSolve[{m\ \(\(x'\)'\)[t] + c\ \(x'\)[t] + k\ x[t] == 0, x[0] == 1.0, 
          \(x'\)[0] == 0.0}, x[t], {t, 10}]\)

NDSolve expects as its first argument a list containing the differential equation and as many extra conditions as are needed to specify a particular solution. (A list is a sequence of expressions separated by commas and enclosed within curly brackets). Here, since the differential equation is second-order, two conditions are needed, in this case the initial values of  x and x'. Second, it expects the thing that has to be solved for (x as a function of t). Finally, it expects the independent variable (time, t) and the duration of the numerical integration (10 time units in this case).

We need a way of getting actual values out of the "hidden" function that comes out of the call to NDSolve -- the following command defines a new function "xsol" which does that:

xsol[t_] = First[x[t] /. solution]

You can use xsol like any other function, make a plot of it, say:

Plot[xsol[t], {t, 0, 10}]

Since you are going to have to run the solution procedure many times, edit your previous inputs so that you have the whole procedure set up in one piece, as follows:

\(m = 1; 
c = 0; 
k = 1; 

    solution = 
      NDSolve[{m\ \(\(x'\)'\)[t] + c\ \(x'\)[t] + k\ x[t] == 0, x[0] == 1.0, 
          \(x'\)[0] == 0.0}, x[t], {t, 10}]; 

    xsol[t_] = First[x[t] /. solution]\)

To use this procedure efficiently, we suggest that, if you have not already done so,  you get the hang of the "Copy" (Ctrl-C) and "Paste"  (Ctrl-V)  features (in the Edit menu), and "Copy Input from Above" (Ctrl-L) in the Input menu.

Exercises and explorations 1

1.	Run the solution procedure to set up an "xsol" for the case of m = 1, c = 0 and k = 1. Plot a graph of xsol, and check the following particular values:

xsol[0]xsol[π/2]xsol[π]xsol[2  π]xsol[3  π]

Allowing for the possibility of numerical errors, is our numerical procedure generating the correct motion for this undamped (free vibration) case? What about when k = 2? Or when m = 2? What happens if you change the initial conditions to x(0) = 0.0, x'(0) = 1.0 ?

2.	Run the solution procedure for the case of fixed mass and spring constant, m = 1 and k = 1, and c taking the values 0.25, 0.5, 1.0, 1.5, 2.0, 2.5 and 3.0. Characterise the motion in each case: is it oscillatory? What's happening to the amplitude of the motion, and how fast is it happening -- i.e. how quickly does the motion tend to a stop? You will need to adjust the time duration in some cases. In your graphs of xsol, remember that you can use the PlotRange option to bring out details by specifying the y-range (or both x- and y-ranges).

More damping isn't necessarily better

There are several surprising things that happen as you increase c up to 2 and beyond. For c < 2, the motions behave as you'd expect from having seen pendulumsor masses on springs moving with air resistance: the motion is an oscillation, but the amplitude gradually decreases to zero.

The first surprising result is that for c = 2 and above, there are no oscillations at all. The second surprising result -- which you may not have spotted -- is that the most rapid decay of these non-oscillatory motions occurs just when c = 2 -- more damping is less effective at bringing the motion to a stop!

The motions can be characterised,then, in terms of three regimes: under-damped (oscillating, c < 2), critically damped (non-oscillating, most rapid decay, (Subscript[c, critical])= 2) and over-damped (non-oscillating, but slower decay, c > 2).

Exercises and explorations 2

Show several solution graphs for different values of c on the same axes. Your figure should illustrate what happens as the value of c goes through the "critical" value 2. (One approach would be to create several different solution functions-- xsol1,  xsol2, etc -- and display them using Plot with an appropriate PlotStyle setting).

2. The Forced, Damped Oscillator: Resonance

Another feature of real vibrating systems that we would like to model mathematically is forcing. This may occur intentionally, in a device such as a driven pendulum in a clock mechanism, where you want the oscillation to continue indefinitely in spite of the damping in the system, or unintentionally, such as in a car where the suspension system is made to vibrate by the regular vibrations of the engine.

Our model for a driving force will assume two things: it is an external force that does not depend on the position and velocity of the oscillating body, and it varies periodically in time, let's say (for simplicity) in terms of a sine function. Our revised equation of motion will look like this:

m(((d^( 2)) x)/(d (t^( 2)))) +c(( d x)/(d t)) +k x =A sin(Ω t)

where A is the amplitude of the driving force and Ω is its frequency (that's a capital omega, for which there is a symbol on the basic input palette, or use <esc>O<esc>).  Here is a revised Mathematica procedure incorporating the new force term:

\(m = 1; 
c = 0.25; 
k = 1; 
A = 2; 
Ω = 2; 

    solution = 
      NDSolve[{m\ \(\(x'\)'\)[t] + c\ \(x'\)[t] + k\ x[t] == 
            A\ Sin[Ω\ t], x[0] == 1.0, \(x'\)[0] == 0.0}, 
        x[t], {t, 50}]; 
xsol[t_] = First[x[t] /. solution]\)

Exercises and explorations 3

1.	Plot the "xsol" given by the procedure above for t from 0 to 50. How would you describe what you see? Can you identify transient and steady-state behaviour in this picture? What is the amplitude of the steady-state vibration?

2.	For the same system configuration (i.e. m =1, c = 0.25, k = 1, A = 2), generate solutions for the values of Ω given in the table below. Estimate (to about 1 decimal place)* the amplitude of the steady state vibration in each case, and write it in the box provided. Note that you may need to run for a longer time duration to see the steady state appear.

({{Ω,     0.2      ,   0.4   ,       0.6    ,       0.7      ,      0.8      ,       0.9      ,      1.0      }, {\(Amplitude\ of
	steady\ state\)	,  ,  ,  ,  ,  ,  ,  }, {Ω, 1.1,       1.2    ,    1.3   ,   1.4      ,  1.6     , 1.8  , 2.0 }, {\(Amplitude\ of
	steady\ state\)	,  ,  ,  ,  ,  ,  ,  }})

* There are several ways that Mathematica can assist in the estimation process: (i) enlarge the graphic output: click on the graphic, move the mouse to the "resize" blob at the bottom right-hand corner, click (and hold) on the blob and drag the mouse to the right; (ii) use the option GridLines->Automatic in Plot; (iii) click on the plot graphic and hold down the <Ctrl> key -- the cursor's shape will change to a "+" -- as you move the cursor around, its coordinate position can be read off in the bottom left corner of the Mathematica window.

The phenonenon you have just observed, where the steady-state amplitude goes through a maximum as Ω is varied, is called resonance. It turns out that the resonant forcing frequency is (for small c) very close to the natural frequency of unforced vibrations -- in this case, (Sqrt[k/m]) (actually, it's a little less than this -- the difference increases as c increases). Whilst the mathematics behind this is fairly sophisticated, it is interesting that even young children "derive" this result when they try to get a swing to ride as high as possible.

3.	[OPTIONAL] Use the command ListPlot to make a graph of your table of data from question 2 (check the syntax in the Help system). Theory gives the following result for the steady state amplitude as a function of Ω and the other parameters -- first, we'll use Clear to remove any values given to these variables previously:

Clear[Ω, c, k, m, A]

amplitude[Ω_] := 
      (A/(Sqrt[(((k - m (Ω^2)))^2) + 
            (((c Ω))^2)]))

Make a plot of "amplitude" for the parameter values used in question 1. Use the Show command (again, check the Help system) to show a combined plot of your data points and the theoretical curve. Is the agreement a good one?

3. Mechanical Simulation Movies

This section rounds off the discussion from Session 2 about making movies of mechanical simulations. We'll combine the ideas from there about using graphical "objects" to construct pictures, with using numerical solutions to describe the motions of the objects.

We'll set up a situation like the one shown in the diagram on page 1 -- a mass sliding on a smooth floor and connected to a wall by a spring (represented by a sine wave) and a dashpot (which, for simplicity, we'll represent just as a thick line). Our movie will contain four objects:

wallAndFloor = Line[{{0, 1.5}, {0, 0}, {5, 0}}]

mass[x_] := {GrayLevel[0.8], Rectangle[{x - 0.5, 0}, {x + 0.5, 1}]}

dashpot[x_] := {Thickness[0.022], Line[{{0, 0.33}, {x - 0.5, 0.33}}]}

spring[x_] := 
      Line[Table[{a, 0.66 + 0.1  Sin[32  a/(x - 0.5)]}, {a, 0, x - 0.5, 
            (x - 0.5)/60}]]

And the following procedure will build the movie -- xeqm is the equilibrium position and xsol will be the displacement as a function of time:

\(xeqm = 2.5; 

    Table[
		
      Show[Graphics[{
		\ wallAndFloor, 
		\ mass[xeqm + xsol[t]], 

            		\ spring[xeqm + xsol[t]], 
		\ dashpot[xeqm + xsol[t]]}], 
        AspectRatio -> Automatic, PlotRange -> {{\(-0.1\), 5}, {0, 1}}], 

      	{t, 0, 50, 1}]\)

Of course, you will need to adjust the time duration (here set to 50) for whatever it is for the xsol you input. You may need to specify quite a small timestep to get a smooth animation of the vibrations.

Exercises and explorations 4

1.	Make a movie of an (unforced) underdamped vibration. For example, this configuration produces a good result:
             m = 2, c = 0.25, k = 1, x(0) = 2.0, x'(0) = 0.0, time duration = 50.

2.	Now change the value of c in your configuration to make a movie of an overdamped vibration.

3.	Make a movie of a forced, damped vibration.