Modes of a Vibrating Building Under Sinusoidal Forcing

The equations of motion for the four story building take this form:

$$ \mathbf{M} \dot{\bar{s}} + \mathbf{K} \bar{c} = \bar{F} $$

where

$$ \bar{c} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{bmatrix}, \bar{s} = \begin{bmatrix} v_1 \\ v_2 \\ v_3 \\ v_4 \end{bmatrix} $$$$ \mathbf{M} = \begin{bmatrix} m_1 & 0 & 0 & 0 \\ 0 & m_2 & 0 & 0 \\ 0 & 0 & m_3 & 0 \\ 0 & 0 & 0 & m_4 \end{bmatrix} $$$$ \mathbf{K} = \begin{bmatrix} k_1+k_2 & -k_2 & 0 & 0 \\ -k_2 & k_2 + k_3 & -k_3 \\ 0 & -k_3 & k_3+k_4 & -k_4 \\ 0 & 0 & -k_4 & k_4 \end{bmatrix} $$$$ \bar{F} = \begin{bmatrix} F_1(t) \\ F_2(t) \\ F_3(t) \\ F_4(t) \end{bmatrix} $$

The forces $F_1, F_2, F_3, F_4$ are the lateral, arbitrary, forces applied to each floor.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from resonance.linear_systems import FourStoryBuildingSystem

This gives a bit nicer printing of large NumPy arrays.

In [2]:
np.set_printoptions(precision=5, linewidth=100, suppress=True)
In [3]:
%matplotlib inline

Determine the model frequiences by calculating the eigenvalues

In [4]:
sys = FourStoryBuildingSystem()
In [5]:
M, C, K = sys.canonical_coefficients()
In [6]:
L = np.linalg.cholesky(M)
In [7]:
K_tilde = np.linalg.inv(L) @ K @ np.linalg.inv(L.T)
In [8]:
evals, evecs = np.linalg.eig(K_tilde)

These are the modal frequencies in radians per second:

In [9]:
ws = np.sqrt(evals)
ws
Out[9]:
array([2.10122, 1.71293, 1.11803, 0.38829])

Forcing the fourth floor at the largest natural frequency

Note that the initial state values are not all zero. Set them all to zero so we see only the effects of forcing.

In [10]:
sys.states
Out[10]:
{'x1': 0.001, 'x2': 0.01, 'x3': 0.02, 'x4': 0.025, 'v1': 0.0, 'v2': 0.0, 'v3': 0.0, 'v4': 0.0}
In [11]:
sys.coordinates['x1'] = 0.0
sys.coordinates['x2'] = 0.0
sys.coordinates['x3'] = 0.0
sys.coordinates['x4'] = 0.0
In [12]:
sys.states
Out[12]:
{'x1': 0.0, 'x2': 0.0, 'x3': 0.0, 'x4': 0.0, 'v1': 0.0, 'v2': 0.0, 'v3': 0.0, 'v4': 0.0}

Create two new constants that describe an amplitude for a sinusoidal forcing.

$$F(t) = A\sin(\omega T)$$
In [13]:
sys.constants['amplitude'] = 100  # N
sys.constants['frequency'] = np.deg2rad(10.0)  # rad/s

Now define a function that takes constants and/or time as inputs and outputs the entries of $\bar{F}$ in the same order as the coordinates and speeds. The following definition applies a sinusoidal force only to the 4th floor.

In [14]:
def push_floor(amplitude, frequency, time):
    F1 = 0.0 if np.isscalar(time) else np.zeros_like(time)
    F2 = 0.0 if np.isscalar(time) else np.zeros_like(time)
    F3 = 0.0 if np.isscalar(time) else np.zeros_like(time)
    F4 = amplitude * np.sin(frequency * time)
    return F1, F2, F3, F4

This function should work with scalar values of time and 1d arrays of time:

In [15]:
push_floor(1.0, 2.0, 3.0)
Out[15]:
(0.0, 0.0, 0.0, -0.27941549819892586)
In [16]:
push_floor(1.0, 2.0, np.ones(5))
Out[16]:
(array([0., 0., 0., 0., 0.]),
 array([0., 0., 0., 0., 0.]),
 array([0., 0., 0., 0., 0.]),
 array([0.9093, 0.9093, 0.9093, 0.9093, 0.9093]))

Now add the forcing function to the system:

In [17]:
sys.forcing_func = push_floor

The forced_response() function works like the free_response() function but it will apply the forcing in the simulation.

In [18]:
traj = sys.forced_response(100)
In [19]:
traj[sys.coordinates.keys()].plot(subplots=True);
In [20]:
sys.animate_configuration(fps=10, repeat=False)
Out[20]:
<matplotlib.animation.FuncAnimation at 0x7f356f1e1eb0>

Exercise: Forcing at the modal frequencies

Update the frequency value to simulate the fourth floor being forced at each of the four natural frequencies and note your observations. Compare the obeserved motion to the mode shapes associated with that modal frequency. Use the animation to help visualize what is happening.

In [21]:
sys.constants['frequency'] = ws[0]  # rad/s
traj = sys.forced_response(100)
traj[sys.coordinates.keys()].plot(subplots=True)
sys.animate_configuration(fps=10, repeat=False)
Out[21]:
<matplotlib.animation.FuncAnimation at 0x7f356ef02f40>
In [22]:
sys.constants['frequency'] = ws[1]  # rad/s
traj = sys.forced_response(100)
traj[sys.coordinates.keys()].plot(subplots=True)
sys.animate_configuration(fps=10, repeat=False)
Out[22]:
<matplotlib.animation.FuncAnimation at 0x7f356ebab880>
In [23]:
sys.constants['frequency'] = ws[2]  # rad/s
traj = sys.forced_response(100)
traj[sys.coordinates.keys()].plot(subplots=True)
sys.animate_configuration(fps=10, repeat=False)
Out[23]:
<matplotlib.animation.FuncAnimation at 0x7f356e973730>
In [24]:
sys.constants['frequency'] = ws[3]  # rad/s
traj = sys.forced_response(100)
traj[sys.coordinates.keys()].plot(subplots=True)
sys.animate_configuration(fps=10, repeat=False)
Out[24]:
<matplotlib.animation.FuncAnimation at 0x7f356e9ad190>

Exercise: Forcing at a node

Recall that the 3rd mode shape has an eigenvector component that is zero:

In [25]:
evecs
Out[25]:
array([[ 0.42853,  0.65654, -0.57735,  0.22801],
       [-0.65654, -0.22801, -0.57735,  0.42853],
       [ 0.57735, -0.57735, -0.     ,  0.57735],
       [-0.22801,  0.42853,  0.57735,  0.65654]])

This is called a "node". This node is associated with $x_3$, the third floor and it tells us that there is no motion at floor 3 if this mode is excited.

Adjust the forcing function to apply sinusoidal forcing at the third floor. Use the third modal frequency to apply forcing to the third floor, then use one of the other modal frequencies to force the third floor. Compare the results and discuss your observations.

In [26]:
sys.constants['frequency'] = ws[2]  # rad/s
def push_floor(amplitude, frequency, time):
    F1 = 0.0 if np.isscalar(time) else np.zeros_like(time)
    F2 = 0.0 if np.isscalar(time) else np.zeros_like(time)
    F3 = amplitude * np.sin(frequency * time)
    F4 = 0.0 if np.isscalar(time) else np.zeros_like(time)
    return F1, F2, F3, F4
sys.forcing_func = push_floor
traj = sys.forced_response(100)
traj[sys.coordinates.keys()].plot(subplots=True)
sys.animate_configuration(fps=10, repeat=False)
Out[26]:
<matplotlib.animation.FuncAnimation at 0x7f356e5a2730>
In [27]:
sys.constants['frequency'] = ws[0]  # rad/s
def push_floor(amplitude, frequency, time):
    F1 = 0.0 if np.isscalar(time) else np.zeros_like(time)
    F2 = 0.0 if np.isscalar(time) else np.zeros_like(time)
    F3 = amplitude * np.sin(frequency * time)
    F4 = 0.0 if np.isscalar(time) else np.zeros_like(time)
    return F1, F2, F3, F4
sys.forcing_func = push_floor
traj = sys.forced_response(100)
traj[sys.coordinates.keys()].plot(subplots=True)
sys.animate_configuration(fps=10, repeat=False)
Out[27]:
<matplotlib.animation.FuncAnimation at 0x7f356e64f310>