Casual Introduction to Numerical Methods – a spring-mass-damper system model – part#3

Here is the third part of a tutorial in both elementary dynamics and numerical methods. It is written at a basic level and it shows you how to set up a dynamic model for numerical solving of simple differential equations.

The dynamic model makes use of an infinite loop, which make the calculations advance in time. Instead of a large table with formulas, this time we have a small number of cells with formulas in the upper part of the table (usually calculations take place in the top most row let’s call it “the current time row” or “the present”) and a large number of cells with historical data (constants) under this active row.

During every loop sequence, a new set of calculations in the “current time row” occur. The input data used in this calculations is taken from the last set of historical data (previous time step). The historical data is organized like a FIFO stack: the oldest time step of historical data is deleted each cycle, the rest of historical data is shifted down one time step per loop cycle to make room for the result of the new calculation which now will represent the latest historical data on top of the stack.

Download the model and hit Reset first and then RunStop.

Enjoy, George

[sociallocker][/sociallocker]

A casual approach to numerical modeling Spring-Mass-Damper-System – part 3.

by George Lungu

A tutorial about the implementation of a dynamic S-M-D model in Excel 2003

Advantages of a dynamic model:

1. Small size spreadsheet since there are just a few formulas (written for only one or a small group of time steps)

2. The model can be run as an infinite loop

3. The model can be paused at any time and restarted from exactly the same exact time step any time later

 

Drawbacks of the dynamic model:

There are no real drawbacks to a dynamic model. Dynamic models
have slightly more VBA code and can sometimes can be a bit slower than
the static ones but there are ways to overcome this by writing mixed models
which are dynamic models nonetheless, but having 10-1000 lines of static
calculation.

These mixed models are in average as fast as their static
counterparts and have all the advantages of the dynamic models.

<excelunusual.com> 2

 

A theoretical review:

Since it is very important to understand the principles of numerical
simulation, let’s again look at the principles behind modeling the spring-mass-
damper system. This was done in the first part of the presentation already.

 

Definitions used in this model:

The mechanics studied in the high school physics class is very simple.
In that class the movement of a body is either uniform or uniformly accelerated.

By definition, the velocity (speed) is the first derivative of position with respect
to time and the acceleration is the first derivative of the speed with respect to
time or the second derivative of position with respect to time.

dx dv d 2x
Velocity formula: v
Acceleration formula: a

dt dt dt
We can also derive the acceleration from Newton’s second law:

<excelunusual.com> 3

In general, elementary formulas are not very useful in a real life situations,
since we almost never encounter a body traveling at a constant speed or with a perfectly
constant acceleration over long periods of time.

The formulas serve two purposes, firstly they have educational value and secondly they
can be applied with good precision in real life situations during very small intervals of time.

The question is, how good a precision we are talking about? And the answer is – any precision,
the shorter the time interval, the better the precision.

The picture below shows the evolution of a general natural phenomenon in time
(blue). A certain cause (the green curve below) produce an effect (the upper blue curve).
Most of the times calculating the effect at moment “t” is either impossible or very difficult
unless we use numerical methods.

Unknown real life curve – this is what we are trying to estimate using numerical simulation.

Unknown effect
Known causes
t = 0
t_current

<excelunusual.com> 4

 

If we divide the modeling time interval in N elementary intervals (A), we can
calculate the effects on each interval using the cause during that time interval (B) .

If the intervals are very short, we can use extremely simple laws to calculate the effects.
elementary effects (on very small intervals) calculated by using linearization and simple
numerical method approximations

t = 0 t_current t = 0

We can approximately calculate the global effect using the algebraic sum of the elementary
effects which is nothing else than numerical integration (C) =>

effects after integration
(algebraic summation of effects)

<excelunusual.com> 5

 

Let’s review our particular system:

L0 = 1m (unstressed)
M = 1Kg
(Mass)
(Spring Constant K = 1N*m)
x = 0 (position from the point of equilibrium)
(Damping Constant = 1N*s/m)

There are a total of 3 forces acting on mass M:

1. Elastic force : Fe = – K * x (the elastic force is proportional and opposite
to the spring deformation)

2. Damping force : Fd = – DC * v (the force of friction is proportional and
opposite to the speed – the proportionality factor is the damping coefficient)

3. Inertia : Fi = – M * a
Because this is a one-dimensional problem we dropped the vector notations. We
assume that positive vectors face to the right while negative ones face to the left.
<excelunusual.com> 6

As a parenthesis, I like to express the damping (friction) force function of the damping
ratio rather than damping coefficient. This way the damping time is only related to the
resonance period of the system.

The formula for the damping ratio is:
The damping force becomes:
We can apply Newton’s second law of dynamics to calculate the acceleration:
elastic damping

This way the acceleration becomes:

There will be only three formulas in our spreadsheet for “a”, “v” and “x”. These
quantities will be calculated for each time step using data from the previous time step.

<excelunusual.com> 7

Let’s review the calculations used :

Since the force is the cause of change in movement, our algorithm will calculate the
current acceleration first, using the velocity and position from the previous time step
since the velocity and position of the current step are not yet available:

The velocity is calculated second, using the acceleration from the current time step and
the velocity from the previous time step:
v  v a dt
current previous current
The position is calculated third, using the velocity from the current time step and the position
from the previous time step:

current previous

<excelunusual.com> 8

 

The outline of a dynamic model:

The dynamic model makes use of an infinite loop, which make the
calculations advance in time.

Instead of a large table with formulas, this time we have a small number of
cells with formulas in the upper part of the table (usually calculations take place
in the top most row let’s call it “the current time row” or “the present”) and a large
number of cells with historical data (constants) under this active row.

During every loop sequence, a new set of calculations in the “current time
row” occur. The input data used in this calculations is taken from the last set of
historical data (previous time step).

The historical data is organized like a FIFO stack: the oldest time step of historical
data is deleted each cycle, the rest of historical data is shifted down one time step
per loop cycle to make room for the result of the new calculation which now will
represent the latest historical data on top of the stack.

<excelunusual.com> 9

 

A snapshot of the dynamic model worksheet layout:

Initial conditions (constants)
Current row or present – active calculations
Latest History – constants

The initial conditions idle in the blue area most of the time.

During the Reset only, the “Reset” macro will copy them and paste below the active row” green.

During the next cycle, the active (red) formulas will use them once to calculate (a,v,x) for
the second time step (0+dt).

After that operation, the initial condition data becomes “history” and is shifted down the stack
until it reaches the end and is discarded.

 

Moving history Stack (past) – constants

Copy the workbook “Osc_2” and rename the copy “Osc_3”

In this worksheet we will add the following features:

1. We will move away from the static large table of formula simulation into a
dynamic style of simulation (loop). The table will look almost the same but now most of
it is filled with historical data (constant numbers) except one row of active formulas

2. Upgrade the macros and add a reset button

3. Generate a more detailed animation of the system

<excelunusual.com> 10

 

The worksheet structure:

Column B contains a time series, column C will contain a series of accelerations, column D
a series of speeds and column E the position of the oscillating mass with respect with the neutral
point.

There will be a total of only three formulas in the worksheet placed in row “9”:

C9: “ =-(E10*B$2+D10*B$3*2*SQRT(B$1*B$2))/B$1”
D9: “ =D10+C9*B$4”

The acceleration is equal to the sum of elastic and friction forces divided by the proof mass.
The velocity is calculated as the integral of acceleration.

E9: “ =E10+D9*B$4”

The position is the integral of the velocity (a first order approximation was used).

We used the following recursive, first order approximation of the integral of a function f(t):

f(t)dt  f(t)dt  f (t)t

<excelunusual.com> 11

Macros used in this worksheet:

Besides straightforward spinner button macros (which are not shown in the here)
the following two macros were used: Reset and StartStop

Sub reset()
DoEvents
Range(“B9”) = 0
Range(“B10:E1010”).Clear
Range(“B10:E10”) = Range(“B8:E8”).Value
End Sub

The reset macro achieves three different functions:

1. Resets the time in the current time row to zero
2. Clears the history (deletes everything up the current time row)
3. After the history is cleared, the macro pastes the initial conditions one row
below the active row so that the active formulas can use these initial conditions to
calculate the values of the first time step (acceleration, speed, position)

<excelunusual.com> 12

Macros used in this worksheet (continuation):

Let’s see the Reset macro in action:
Initial conditions (constants)
Current row or present – active calculations
Latest
History – constants
Moving history
Stack (past) – constants
Before “Reset”
After “Reset”

<excelunusual.com> 13

 

The StartStop macro:

Macro declaration statement

It logically “flips” the Boolean control variable RunSim so that we can use the same button to start or stop the macro

Conditional Do loop statement declaration. The loop will run until either RunSim is False or the time (in cell B9) becomes larger than 31 seconds

The macro hogs the processor, this statement allows other processes to happen during each loop cycle such as updating the chart.

Increments time in cell B9 by delta t (which is the time step from cell B4).

End of “Do” loop statement.

End of macro statement.

Sub StartStop()
RunSim = Not (RunSim)
Do While RunSim = True And [B9] =< 31
DoEvents
Range(“B10:E1010”) = Range(“B9:E1009”).Value
Range(“B9”) = Range(“B9”) + Range(“B4”)
Loop
End Sub

The StartStop macro achieves three different functions:

1. If the calculation loop is stopped it will start it and vice versa

2. If the calculation loop is running it will increment the time

3. If the calculation loop is running it will delete the earliest (bottom) entry in
the history, shift the whole history stack down and write the active data values on the
top of the history stack (see the bold blue line within the macro).
This way the history stack will act like a FIFO (First-In-First-Out type of shift register).

<excelunusual.com> 14

 

Let’s observe the StartStop macro and the conditional “Do” loop in action:

t a v x t a v x t a v x t a v x t a v x
0 -0.4 0 -0.4 0 -0.4 0 -0.4 0 -0.4
0 1 0.1 -0.39 0.1 0.956026 0.195603 -0.37044 0.2 0.888986 0.284501 -0.34199 0.3 0.800994 0.364601 -0.30553 0.4 0.694646 0.434065 -0.262123
0 -0.4 0 1 0.1 -0.39 0.1 0.956026 0.195603 -0.37044 0.2 0.888986 0.284501 -0.34199 0.3 0.800994 0.364601 -0.30553
0 -0.4 0 1 0.1 -0.39 0.1 0.956026 0.195603 -0.37044 0.2 0.888986 0.284501 -0.34199
0 -0.4 0 1 0.1 -0.39 0.1 0.956026 0.195603 -0.37044
0 -0.4 0 1 0.1 -0.39
0 -0.4

After first loop
After second loop
After “Reset” After third loop After fourth loop iteration

loop iteration

t a v x t a v x t a v x t a v x t a v x
0 -0.4 0 -0.4 0 -0.4 0 -0.4 0 -0.4
0.5 0.572949 0.49136 -0.212987 0.6 0.439238 0.535284 -0.159459 0.7 0.297084 0.564992 -0.102959 0.8 0.150199 0.580012 -0.044958 0.9 0.002346 0.580247 0.013067
0.4 0.694646 0.434065 -0.262123 0.5 0.572949 0.49136 -0.212987 0.6 0.439238 0.535284 -0.159459 0.7 0.297084 0.564992 -0.102959 0.8 0.150199 0.580012 -0.044958
0.3 0.800994 0.364601 -0.30553 0.4 0.694646 0.434065 -0.262123 0.5 0.572949 0.49136 -0.212987 0.6 0.439238 0.535284 -0.159459 0.7 0.297084 0.564992 -0.102959
0.2 0.888986 0.284501 -0.34199 0.3 0.800994 0.364601 -0.30553 0.4 0.694646 0.434065 -0.262123 0.5 0.572949 0.49136 -0.212987 0.6 0.439238 0.535284 -0.159459
0.1 0.956026 0.195603 -0.37044 0.2 0.888986 0.284501 -0.34199 0.3 0.800994 0.364601 -0.30553 0.4 0.694646 0.434065 -0.262123 0.5 0.572949 0.49136 -0.212987
0 1 0.1 -0.39 0.1 0.956026 0.195603 -0.37044 0.2 0.888986 0.284501 -0.34199 0.3 0.800994 0.364601 -0.30553 0.4 0.694646 0.434065 -0.262123
0 -0.4 0 1 0.1 -0.39 0.1 0.956026 0.195603 -0.37044 0.2 0.888986 0.284501 -0.34199 0.3 0.800994 0.364601 -0.30553
0 -0.4 0 1 0.1 -0.39 0.1 0.956026 0.195603 -0.37044 0.2 0.888986 0.284501 -0.34199
0 -0.4 0 1 0.1 -0.39 0.1 0.956026 0.195603 -0.37044

After fifth loop iteration
After sixth loop iteration
After seventh iteration
After eighth loop iteration
After ninth loop iteration

Observe the data shifting and renewal that takes place.
The macro can be started at any time and restarted from the place it was turned off.

<excelunusual.com> 15

 

The next tutorial will explain the animation used in the simple mass-spring-damper system model

L1 L2, L0 / L-equilibrium, L3, L4, L5, L6, L7, L8

by George Lungu <excelunusual.com>