Departments of Civil Engineering and Mathematics
CE 109: Computing for Engineering --- Mathematica Session 3

Damped and Forced Vibrations

In Session 2, you generated some animations for undamped, free vibrating systems. In this session, we'll look at the characteristic behaviours of systems with damping (in section 1) and forcing (in section 2). 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. We will develop on this model by adding, first, a "damping term", and then a "forcing term" to the model equation.

This session contains a small piece of ASSESSED WORK, in Exercises and Explorations 2 and 3. Your answers to this should be printed out (preferably as a Mathematica "notebook" document) and submitted to the "black box" in the General Office by Friday 5 February.

1. The Damped Oscillator

Recall the equation of motion for a simple harmonic oscillator (SHO); if v represents the displacement of such an oscillator from its equilibrium position, then by Newton's second law we have:

m(((d^( 2)) v)/(d (t^( 2))))=-k v ,    or, rearranging:    m(((d^( 2)) v)/(d (t^( 2)))) + k v =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)) v)/(d (t^( 2)))) +c(( d v)/(d t)) + k v =0

This will be our mathematical model for damped oscillations. The physical system which it models can be treated as a light beam, simply supported at each end and carrying a concentrated mass m at the mid-span. We will ignore beam deflections due to gravity, so the equilibrium position of the beam is horizontal:

Now suppose that the beam is deflected from equilibrium by a displacement v, and is moving with a velocity v' . The mass will be subject to two forces: a "spring force", kv, and a damping force, cv'. The following schematic diagram shows (with the vertical scale rather exaggerated) the forces acting when the beam is negatively deflected (v < 0), and moving back towards the horizontal (v' > 0):

The structural damping constant, c, encapsulates several different damping effects: hysteretic damping (cycles of tension and compression inside the material of the beam), air resistance, and friction at the support contacts; to a good approximation, they can all be modelled as a single force that is proportional to velocity. The "spring stiffness", k, can be related to the flexural stiffness of the beam; for this simple end-supported configuration it is given by k= 48EI/(L^3), so typically k is of order (10^6) N (m^(-1)). However, for simplicity we are going to work with non-physical quantities in this task, with m, k and c all taking values of order 1.

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 stiffness, k, damping constant, c, and mass, m, are varied.

Type in the following procedure:

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

    solution = 
      NDSolve[{m\ \(\(v'\)'\)[t] + c\ \(v'\)[t] + k\ v[t] == 0, v[0] == 1.0, 
          \(v'\)[0] == 0.0}, v[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  v and v'. Second, it expects the thing that has to be solved for (v 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).

NDSolve returns a special Mathematica object called an "interpolating  function", which is essentially a list of coordinates corresponding to the solution, with some additional helpful data thrown in. By default, Mathematica keeps hidden the details of the interpolating function, which is fine, but we do need a way of getting actual values out of the "solution" object -- the following command defines a new function "vsol" which does that:

vsol[t_] = First[v[t] /. solution]

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

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

or determine its value at some point in time:

vsol[3.1416]

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\ \(\(v'\)'\)[t] + c\ \(v'\)[t] + k\ v[t] == 0, v[0] == 1.0, 
          \(v'\)[0] == 0.0}, v[t], {t, 10}]; 

    vsol[t_] = First[v[t] /. solution]; 
Plot[vsol[t], {t, 0, 10}]\)

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 "vsol" for the case of m = 1, c = 0 and k = 1. Plot a graph of vsol, and check the following particular values:

vsol[0]vsol[π/2]vsol[π]vsol[2  π]vsol[3  π]

Allowing for the possibility of numerical errors, is our numerical procedure generating the correct motion for this undamped case? What happens if you change the initial conditions to v(0) = 0.0, v'(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 vsol, remember that you can use the PlotRange option to bring out details by specifying the v-range (or both t- and v-ranges).

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 pendulums or masses on springs moving with air resistance: the motion is an oscillation, but the amplitude gradually decreases to zero. This is called under-damped behaviour.

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. Applying more damping beyond this critically damped behaviour turns out to be less effective at bringing the motion to a stop! Hence, the behaviour above the critical value of c is called over-damped.  (Note that, in general for the model we are looking at, (Subscript[c, critical]) = 2 (Sqrt[ k m])).

Critical damping is a very important aspect of mechanical engineering design (e.g. for optimal control of unwanted vibrations in machinery), but very rarely has significance for civil engineering structures -- typically, structural damping constants take values between just 1% and 10% of the critical case.

2. The Forced, Damped Oscillator: Resonance

Now we are going to add another term to the model differential equation, which will describe the effects of  forcing -- applying an external "driving force" to a structure. Important practical examples of this phenomenon are the pulsating forces due to gusts of wind, or vehicle movement across a bridge.

Our model for a (idealised) 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)) v)/(d (t^( 2)))) +c(( d v)/(d t)) +k v =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>).  (This equation is also soluble using symbolic techniques that you'll encounter in your maths course.) Here is a revised Mathematica procedure incorporating the new forcing term:

\(m = 1; 
c = 0.15; 
k = 1; 
A = 2; 
Ω = 1.4; 

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

    Plot[vsol[t], {t, 0, 120}, PlotPoints -> 50]\)

(The plots may show spurious jumps in the solutions as Ω gets large; the PlotPoints option here helps to smooth these out by increasing the resolution of the plot.)

Exercises and explorations 2

1.	Plot the "vsol" given by the procedure above for t from 0 to 120. 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.15, 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. N.B. you will need these results for the ASSESSED WORK below.

({{Ω,     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.

Resonance

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.

Exercises and explorations 3

1.	Use the command ListPlot to make a graph of your table of steady state amplitudes. First, enter your amplitude data into a list with the following "list of lists" structure:

ampData={{0.2,...},{0.4,...},...,{2.0,...}}

Then,

dataPlot = ListPlot[ampData, PlotStyle -> PointSize[0.02]]

2.	Theory gives the following result for the steady state amplitude as a function of Ω and the other parameters:

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

Make a plot of "amplitude" for the parameter values used in Exercises and Explorations 2. Use the Show command (check the Help system) to show a combined plot of dataPlot and the theoretical curve.

ASSESSED WORK: Please submit (preferably in the form of a Mathematica notebook)
      (i)  your list (ampData) of steady state amplitudes (use TableForm[ampData] to get tabular output), and
     (ii)  your data/theory comparison plot.

3. A Structural Simulation Movie [OPTIONAL]

In Session 2, you built up animations with standard Mathematica plots (Plot and Plot3D). Another possible approach is to use a set of functions that are called "graphical primitives" to construct pictures using basic objects such as points, lines, circles and rectangles.

Suppose that we want to create movies to illustrate the damped and forced vibrations of the loaded beam introduced in section 1. There will be four objects in our movie: two fixed supports, a moving mass and a moving beam; we'll take the beam length to be 10 m. The supports are easy to define:

supportA = 
      Line[{{0, 0}, {0.4, -0.4}, {-0.4, -0.4}, {0, 0}}];

supportB = 
      Line[{{10, 0}, {10.4, -0.4}, {9.6, -0.4}, {10, 0}}];

The mass object will be a "blob" whose y-coordinate is equal to the displacement v:

mass[v_] := {PointSize[0.02], Point[{5.0, v}]}

The special plotting construction Show[Graphics[{...}]] is used to display these graphical objects:

Show[Graphics[{supportA, supportB, mass[0.5]}], 
      AspectRatio -> Automatic, PlotRange -> {-0.5, 1.0}]

We'll define the beam as a plot of a quadratic curve in x going through the points (0, 0), (10, 0) and (5, v) :

beam[v_] := 
      Plot[((-v)/25) (x^2) + ((2  v)/5) x, {x, 0, 10}, 
        Axes -> False, DisplayFunction -> Identity]

The DisplayFunction option is required here to turn off graphical output, since we only want to see a single combined picture for each frame of the movie, such as:

Show[beam[0.5], Graphics[{supportA, supportB, mass[0.5]}], 
      DisplayFunction -> $DisplayFunction, AspectRatio -> Automatic, 
      PlotRange -> {-0.5, 1.0}]

(Note that $DisplayFunction is the opposite of Identity -- it turns on graphical output to the screen.)

Finally, generate a solution for the unforced, damped vibration corresponding to the following configuration:
                      m = 1,  c = 0.15,  k =0.5, v(0) =1.0,  v'(0)  = 0.0,  time duration = 35.
Feed the resulting displacement function "vsol" into the following procedure to create the simulation movie:

\(Table[
	
      Show[beam[vsol[t]], Graphics[{supportA, supportB, mass[vsol[t]]}], 
        DisplayFunction -> $DisplayFunction, AspectRatio -> Automatic, 
        PlotRange -> {\(-1.2\), 1.2}], {t, 0, 35, 1}]\)

Try out this procedure for some different damped and forced vibrations. Of course, you will need to adjust the time duration for whatever it is for the vsol you input, and perhaps also the PlotRange. And you may need to specify quite a small timestep to get a smooth animation of the vibrations.