# Introducing Mechanical Vibrations with a Book Oscillating on a Coffee Mug

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as pat
from resonance.nonlinear_systems import SingleDoFNonLinearSystem

class BookCupSystem(SingleDoFNonLinearSystem):

def __init__(self):

super(BookCupSystem, self).__init__()

self.constants['d'] = 0.029  # m
self.constants['l'] = 0.238  # m
self.constants['r'] = 0.042  # m
self.constants['m'] = 1.058  # kg
self.constants['g'] = 9.81  # m/s**2
self.coordinates['theta'] = 0.0  # rad
self.speeds['omega'] = 0.0  # rad/s

def rhs(theta, omega, d, l, r, g):
"""Returns the derivatives of the state variables.

Parameters
==========
theta : float
Angle of the book in radians.
omega : float
Angular rate of the book in radians per second.
d : float
Book thickness in meters.
l : float
Book width in meters.
r : float
Cup radius in meters.
g : float
Acceleration due to gravity in meters per squared seconds.

Returns
=======
The angular rate of the book in radians per second.
The angular acceleration of the book in radians per second.

"""
omegadot = (6 * d * g * np.sin(theta) - 12 * g * r * theta * np.cos(theta) -
12 * r**2 * omega**2 * theta) / (4 * d**2 + l**2 + 12 * r**2 * theta**2)

self.diff_eq_func = rhs

def bottom_left_x(r, l, theta):
return r * np.sin(theta) - (r * theta + l / 2) * np.cos(theta)

def bottom_left_y(r, l, theta):
return r + r * np.cos(theta) + (r * theta + l / 2) * np.sin(theta)

def create_plot(time, r, l, d, theta, bottom_left_x, bottom_left_y):
fig, ax = plt.subplots(1, 1)
width = max(l, r * 2)
ax.set_xlim((-width / 2 - width / 10, width / 2 + width / 10))
ax.set_ylim((0.0, 2 * r + 2 * d))
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_aspect('equal')

circ = pat.Circle((0.0, r), radius=r)

rect = pat.Rectangle((bottom_left_x, bottom_left_y),
l, d,
color='black')

title = ax.set_title('Time: {:.2f} [s]'.format(time))

# make sure to return the rectangle, which moves at each time step!
return fig, rect, title

self.config_plot_func = create_plot

def update_frame(time, theta, bottom_left_x, bottom_left_y, rect,
title):
title.set_text('Time: {:.2f} [s]'.format(time))
rect.set_xy((bottom_left_x, bottom_left_y))

self.config_plot_update_func = update_frame
hello world

## Course Learning Objectives

The primary course objectives are (succinctly) that students will be able to:

• characterize vibrations based on measurements from real and simulated dynamic systems
• design vibrating systems with desired characteristics
• create predictive mathematical and computational models of real vibrating systems

In this notebook, you will get a taste of each of these three using a simple real system we can examine in class. Each class session will follow a similar pattern but with new and more complicated systems as we progress. You will be able to complete each step in the entire modeling-analysis-design iterative loop when the course is over.

## Prerequesites

• Familarity using an interpreted scientific programming language (Python, Matlab, R, etc)
• Basic Python skills will be helpful: how to write a function, assign variables, import things
• Intro dynamics: F=ma, ODEs, energy

## Execution Instructions

Execute code cells by clicking on the cell and pressing the "Run" button in the toolbar or by simultaneously pressing the "Shift" and "Enter" keys on your keyboard.

# Introduction

This notebook introduces a single degree of freedom vibratory system: a textbook balancing on a cylindrical coffee mug. The system is implemented as a computational model that you can interact with in order to visualize the free response and compare the computer simulation to a demonstration in the classroom. The video below shows the real system in action:

from IPython.display import YouTubeVideo
YouTubeVideo('B12HbAOKnqI', width=600, height=480)

## What are vibrations?

Here we will study a simple vibratory system. A vibrating mechanical system is typically defined as a collection of rigid and flexible objects that interact in a closed envelope. Thinking about the book and the mug, it will oscillate if initially displaced at a small non-horizontal angle and let go. Note that it oscillates about a horizontal position. This position is called an equilibrium point, equilibrium state, or equilibrium configuration which is a natural position the system comes to when there is no motion. Vibration is formally defined as an oscillation about an equilibrium and occurs if there is a moving inertial object which has restoring forces acting on it.

Vibrations require:

• inertial object that can move
• restoring forces
• equilibrium configuration

During this class, we will examine and explore many different vibratory systems, such as this simple book on a mug system. We will have some live demos, as we are showing now, but in general we will work with computational representations of systems to experiment and learn about the nature of vibration.

## Python imports and setup

This is boiler plate import code that makes all of the necessary libraries available for use and sets up some display setting for the notebook. In Python, all special functionality comes from importing different Python modules and packages.

import numpy as np
import matplotlib.pyplot as plt

%matplotlib widget
hello world

# Model Description

For the purposes of the introduction we will assume you have come up with a simplified version of the conceptual mechanics that govern the system's motion. In this case there are several assumptions made.

### Exercise

If we want to learn specifically about how the book oscillates on the cup, what might be some good assumptions to make to simplify reality so that we can create a mathematical model?

Here are the main assumptions I started with:

• The book and the cup surfaces are perfectly shaped geometric entities.
• The book rolls without slip on the cup's surface.
• The book is a uniformly dense rectangle that only has planar motion.

The essence of modeling is coming up with appropriate assumptions that let you describe the system in the simplest way such that the model's behavior is sufficient to answer your questions about the system. It is a skill that takes a lot of practice and you will get to practice this a fair amount in this class. A model may be good or it may be bad, depending on what we are hoping to predict. We will check its "goodness" in the analysis section.

Below you will find a pre-determined free body diagram that can likely capture the essential vibratory motion of this system. # Objective 1: Numerical Analysis of the Motion

In this section we will make use of a premade system to handle the numerical evaluation of the equations of motion. The following cells load this system and name it sys so we can make use of it.

sys = BookCupSystem()
hello world

sys is now a new system object that you can interact with. This system has many variables and functions associated with it. You can see what the system has and can do by using the Python dot notation. Type sys. and press the tab key to see a list of all the variables and functions that are associated with this system.

## Analysis Step 1: Visualize the System's Configuration

It is often very helpful to visualize a system's configuration. In this case we need a two dimensional drawing similar to the diagram above. The plot_configuration() function let's us see a simple visualization.

sys.plot_configuration();
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-2d1ddf8a9d4f> in <module> ----> 1 sys.plot_configuration(); NameError: name 'sys' is not defined

### System constants

One thing that systems have are different constants, for example this system has geometry, such as the book's thickness and length and the cup's radius. The book also has a mass and, in this case, an underlying assumption is that the book is uniformly dense. Note that all of these do not change with time, i.e. they are constant. You can view all of the constants, which are stored in a Python dictionary by typing:

sys.constants
hello world

A Python dictionary maps keys, in this case the constant's names, to values, the numerical values you'd like to assign to the constant. For example the key 'l', for "length", is associated with a value 0.029. An individual constant's value can be accessed by using square brackets:

sys.constants['l'] = 0.1  # meters
sys.plot_configuration();
hello world

Set the length constant back to it's original value before moving on.

sys.constants['l'] = 0.238  # meters
hello world

### Coordinates¶

There are other system values of interest too. Another very important type are those that vary with time.

There are are an infinite number of time varying parameters, but it is often preferable to choose a uniquely simple set of time varying parameters, often called generalized coordinates. These coordinates define the configuration of the system. In our case, the vertical and horizontal location of the book's mass center could uniquely describe the configuration of the system (if the book can't slip on the cup). But a better choice would be to use the single time varying angle of the books surface relative to horizontal to define the configuration.

The angle of the book is thus a generalized coordinate because no fewer number of time varying parameters can possibly be used to describe the configuration. For simple systems, the number of generalized coordinates corresponds to the number of degrees of freedom of a system. The degrees of freedom are the number of independent parameters that define the configuration. The non-slipping book on a cup has 1 degree of freedom which is described by the single generalized coordinate, the book's angle. The system's generalized coordinates can be accessed as such:

sys.coordinates['theta'] = np.deg2rad(10)  # radians
sys.plot_configuration();

hello world

## Free Response¶

Now that we have a system with a defined equation of motion. There are two ways to create this motion: apply perturbing forces to the system or set the coordinate to an initial angle other than the equilibrium angle. We will do the later here. The resulting motion is called the free response of the system or the solution to the initial value problem, meaning that no external forces are causing the motion. To simulate the free response of the system, some values of time are needed. In this case a final time value, effectively the duration, is passed into the free_response() function. First, set the initial angle of the book and then call free_response(), storing the returned result in a variable named trajectories :

trajectories = sys.free_response(5.0)
trajectories
hello world

The data frames have useful plotting functions associated with them, so it is rather easy to plot the various coordinates and measurements versus time:

trajectories.plot(subplots=True);
hello world

## Analysis Step 3: Animating The Motion

Now that we we have a time varying response, we can animate the configuration figure to visualize how the system moves. Use the animate_configuration() function to do so:

sys = BookCupSystem()
sys.animate_configuration(fps=10, repeat=False)