The spring and mass

You will get experimental data of a mass-spring system, and write a Python program to solve the equation of motion. You will then use this simulation data to plot a visualization of position vs. time and the phase plot.  You will also calculate the energy and plot it.  You will have to discuss the output, and eventually improve the code.

This lab contains explicit questions for you to answer labelled with Q. However, you still need to maintain a notebook for the lab, make plots and write Python code.

Background knowledge for exercise 4

Python lists, arrays, numpy, pyplot

Physics Harmonic oscillator, Hooke’s law. Review these from your first year text.

1 First Session

1.1 Introduction

Physics of springs is based on Hooke’s Law

Fspring  = .ky ,

where y is the vertical displacement, and k is the spring constant mea- sured in N/m.  Hooke’s law can be used with Newton’s second law to derive the equation of motion:

d2 y k

dt2        m

Q1 Show how equation (2) was derived. Point out the main approxima- tions involved (if any).

The most common method used to find numerical solutions to equa- tions of motion is by setting up a pair of coupled  ordinary  differential equations.

Given a mass, m, coordinate -q, momentum p-, and force ,

=

d-q p-

=

dt      m .

Figure 1:  A vertical spring-mass system

1.2 Numerical Methods

The solution to equation (2) is well-known, but in general, differential equations can be difficult to solve. In practice, we can calculate approximate solutions by approximating the derivatives as

dt  s ∆t ┌y(t + ∆t) . y(t)┐ .

(4)

Therefore, equations (3) can be approximated with

p-(t + ∆t) = p-(t) + ∆t F (-q(t))

-q(t + ∆t) = -q(t) + p t) .

(5a)

(5b)

Equations (5) are update formulae.  Given initial position, momentum, and starting time t = t0 , we can determine a numerical approximation of position and momentum at times t0 + ∆t, 2∆t, . . .. The numerical solution will approach the actual solution as ∆t → 0.

This time-stepping scheme is called the “Forward Euler” method, because it uses the derivative at the start of an interval to extrapolate forward. We will see other time-stepping schemes in the next session.

The scheme can also be written in a recursive form,

p-i+1  = p-i + ∆t F (q-i )

∆t

m

where p-i  = p-(ti ), q-i  = -q(ti ), and ti  = t0 + i∆t.

(6a) (6b)

1.3 Numerical methods for the mass-spring system

The equation of motion for the vertical spring-mass system can be written in the coupled form as:

= .Ω0(2)y

dt

dy

(7a)

(7b)

where Ω0  = ′k/m is the angular frequency of oscillation.

The initial conditions will be y0  = A (the amplitude measured in experiment), and v0  = 0 m/s. Using (7), the numerical approximation can be written as

yi+1  = yi + ∆t vi

vi+1  = vi . ∆t Ω0(2)yi

(8a) (8b)

for i = 0, 1, 2, . . . . This is the update scheme you will use in your first program.

Q2 Show the derivation of equations (8) using the equation of motion (2) and the definition of the Forward Euler method (5).

1.4 The experiment

You should have a mass hanging from a spring.  Underneath the mass is a motion sensor, which is plugged into  a  data  acquisition  device  (DAQ).  The  output  of the  DAQ  is  analyzed  by  a  LabView  application (~\Desktop\2nd  Year  lab  Files\MotionSensor.vi)

The goal of the lab exercise is to determing the period of oscillation of the bob, T (and therefore the frequency Ω0 ). This will be used in our later analysis.

Do the following:

1.  Measure the mass of the bob.

2.  Open MotionSensor.vi.

3.  Switch to the “Setup/Callibrate” panel, and change the samples/second to 100.

4.  Switch to the “Collect” panel. Position the bob at about 20 cm above the motion sensor.

5.  Start the bob oscillation by either raising or lowering it from equilibrium  and letting go;  a few centimeters should be enough.

6.  Click “start detector” the motion sensor should start buzzing now. Click “Collect Data” to begin the data collection.

7.  Let the bob oscillate for about 10 seconds, then stop the detector.

8.  If the  motion sensor is  not  aligned well,  the  data may  be poor.  As  long  as  the graph shows  a fairly regular motion, it should be possible to measure the period.

9.  Switch to the “Analyze” panel. Adjust the yellow cursors to the begining and end of 5 to 10 oscilla- tions.

10.  You may need to  change  the  bounds  of the plot to get a  closer view.   You  do  this  by highlighting the last value on either scale and replacing it by your value.

11.  Read the time between cursors, and divide by the number of oscillations observed to get the period.

12.  Record the period.

13.  Open File/Position Data Set. Save the data file on your memory stick.

14.  Plot Distance vs.  Time with labelled axes.  Ignore the  ”Dist.  Error” for now..  Calculate the period from your plot.

Q3 What is the value of the spring constant, k?

1.5 Python programming

Scripts for simulations typically follow the same pattern

1.  define constants/parameters

2.  set the initial conditions

3.  use the numerical approximation (8) to step forward in time

4.  loop until the final time

5. plot the graph and interpret the result

Additional steps that may be useful include

● wrapping it in a function for repeated simulations,

●  saving the data for later analysis.    The outline for our specific program is,

1.  import required modules (eg. pylab, or numpy and matplotlib.pyplot)

2.  define constants ∆t and Ω0 .

3.  calculate the time values, ti  for the plot using ∆t = 0.01 s, t0  = 0 s, T = 10.0 s.

4. preallocate the array (with the zeros() function)

5.  set initial conditions

6. write the time-stepping loop

7. plot the data (position vs. time) with labelled axes

Write and run the program. Save the program and its output on your memory stick.

1.6 Programming - Analysis

Q4 What do you see happening in the plot? Is this what you expect? Can you interpret the graph?

There are other graphs that are useful for analysis:

1.  ‘velocity vs. time’ (time on horizontal axis, velocity on vertical axis)

2.  ‘velocity vs. position’ (position on horizontal axis, velocity on vertical axis)

The latter is called a phase-plot.

Create these plots with labelled axes and save them.

1.6.1 Energy

Energy (and other conserved quantities) are very useful for analyzing dynamical systems. For the spring-mass system, the mechanical energy is conserved,

Etot  = K(y˙) + U(y) ,                                                                  (9)

where K is kinetic energy (K = mv2 ) and U is potential energy (both elastic and gravitational).  The energy expression we will use is

Etot  = mv2 + ky2 .                                                                (10)

The term for gavitational potential energy can be omitted by wisely choosing the reference for potential energy.  You should be able to prove this statement.

●  Modify your program to calculate the energy at each step. Note that the energy is not zero at t = 0. The mass (m) and spring constant (k) will need to be added to your code.

●  Plot ‘Energy vs. time’

Q5 What does the energy plot suggest?  Does this explain the strange appearance of the velocity vs.  time and velocity vs. position plots?

Q6 For a spring-mass system, the phase plot should be an ellipse. Explain why, using energy conservation. Give an explanation for the appearance of your phase plot.

Q7 (bonus) Determine the leading error in our numerical method. Hint: use a Taylor expansion of y(t + ∆t) and find the terms we ignored in equation (5).

1.7 More on numerical integration schemes

If you properly implemented the numerical approximation in Equation (8), you should have found that the simulation was bad: the amplitude increased, and energy was not conserved. The reason for the failure was the integration method.  We will discuss other numerical integrators, and implement a symplectic scheme, which is suitable for conservative systems.

In the first program, we used the most primitive numerical integration method called the Forward Euler or Explicit Euler method,

y[i + 1] = y[i] + ∆t v[i]

v[i + 1] = v[i] . ∆t Ω0(2)y[i] .

Forward Euler extrapolates the position through an interval ∆t using only the  (approximate) derivative at the begining of the interval.  It  can  be proved [1] that Forward Euler is  unstable,  which  means  that  the amplitude  and  total  energy  increases  in  time.   In phase space, the increase in energy looks like a spiral outward. The accuracy of the Forward Euler method can be improved by decreasing the time step, but the energy will always increase.

A similar method to Forward Euler is Backward Euler or Implicit Euler. This method uses the derivative at the end of the interval (making an implicit scheme) to calculate the next point,

y[i + 1] = y[i] + ∆t v[i + 1]

v[i + 1] = v[i] . ∆t Ω2 y[i + 1] .

(11a) (11b)

Backward Euler is absolutely stable, which means that longer timesteps can be used (especially for “stiff” or rapidly changing ODEs). However, it artificially dissipates energy.

A simple remedy is to combine the Forward and Backward Euler methods to create the Euler-Cromer or

Symplectic Euler Method,

y[i + 1] = y[i] + ∆t v[i]

v[i + 1] = v[i] . ∆t    y[i + 1] .

The result is an explicit method, which conserves approximate energy. Do the following:

(12a) (12b)

2 Second session

In the first part of this exercise, we analyzed the equation of motion of a spring-mass system and wrote Python code to numerically integrate the equation.  In this session, we will add a physical damping term, and qualitatively compare the numerical simulation to the experiment.

2.1 Adding a damping term

In the real world, mechanical energy cannot be perfectly conserved. So far, the equations of motion you have used model the system without damping  (physical dissipation).  In general, the damping force exerted on a body depends on the velocity -v of the body.  Drag is one source of damping caused by velocity of a body relative to the surrounding medium (air or water) commonly modelled as quadratic drag,

d  = . CρAl-vl-v ,                                                                   (13)

where

C  is the drag coefficient (dimensionless),

ρ  is the density of the medium (kg/m3 ),

A  is the cross-sectional area perpendicular to the flow (m2 ), and

-v is the velocity of the body relative to the medium (m/s).

The direction of the damping force is always opposite to the direction of velocity. The drag coefficient is not a universal constant: it depends on viscosity of the medium, shape of the body, and roughness of its surface. It also depends on the velocity through the Reynolds number.

The Reynolds  number  Re,  is a very useful dimensionless quantity used to characterize flows in fluid mechanics. In our case, it will characterize the dependence of the drag coefficient on velocity. The Reynolds number is the ratio of the inertial force of the medium to the viscous force:

ρvl

Re =

η

where

(14)

ρ  is the density of the medium (kg/m3 ),

v  is the velocity of the body relative to the medium (m/s),

l  is the characteristic length in the direction of flow (m), and

η  is the dynamic viscosity (N ·s/m2 ).

If the flow is laminar, the Reynolds number takes small values (Re < 2300) and the drag coefficient is inversely proportional to the velocity. This makes the drag force directly proportional to velocity:

d  = .γ-v ,                                                                         (15)

where γ is the damping coefficient.

When the flow is turbulent, the Reynolds number is large (Re > 4000) the drag coefficient is approximately constant and the drag force is dependent on the square of the velocity.

Q9 What is the Reynolds number for the spring-mass system?  Hint:  this may depend on the  amplitude  of the oscillation.

We expect that the motion and speed of a typical mass-spring system in air correspond to small Reynolds numbers. Therefore, the equation of motion of our system would be written as:

+ γ      + Ω0(2)y = 0 .

(16)