Modes of a Vibrating Building

In this notebook we will find the vibrational modes of a simple model of a building. We will assume that the mass of the floors are much more than the mass of the walls and that the lateral stiffness of the walls can be modeled by a simple linear spring. We will investigate how the building may vibrate under initial conditions that could be caused by a gust of wind and during ground vibration.

In [1]:
from IPython.display import YouTubeVideo
In [2]:
YouTubeVideo('g0cz-oDfUg0', width=600)
Out[2]:
In [3]:
YouTubeVideo('hSwjkG3nv1c', width=600)
Out[3]:
In [4]:
YouTubeVideo('kzVvd4Dk6sw', width=600)
Out[4]:
In [5]:
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 [6]:
np.set_printoptions(precision=5, linewidth=100, suppress=True)
In [7]:
%matplotlib inline

Simulate the four story building

In [8]:
sys = FourStoryBuildingSystem()
In [9]:
sys.constants
Out[9]:
{'m1': 4000,
 'm2': 4000,
 'm3': 4000,
 'm4': 4000,
 'k1': 5000,
 'k2': 5000,
 'k3': 5000,
 'k4': 5000}
In [10]:
sys.coordinates
Out[10]:
{'x1': 0.001, 'x2': 0.01, 'x3': 0.02, 'x4': 0.025}
In [11]:
sys.plot_configuration();
In [12]:
traj = sys.free_response(30, sample_rate=10)
In [13]:
traj[list(sys.coordinates.keys())].plot(subplots=True);
In [14]:
sys.animate_configuration(fps=10)
Out[14]:
<matplotlib.animation.FuncAnimation at 0x7faa6be98820>
In [15]:
M, C, K = sys.canonical_coefficients()
In [16]:
M
Out[16]:
array([[4000,    0,    0,    0],
       [   0, 4000,    0,    0],
       [   0,    0, 4000,    0],
       [   0,    0,    0, 4000]])
In [17]:
C
Out[17]:
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])
In [18]:
K
Out[18]:
array([[10000, -5000,     0,     0],
       [-5000, 10000, -5000,     0],
       [    0, -5000, 10000, -5000],
       [    0,     0, -5000,  5000]])

Exercise

The system can be normalized by the mass matrix and transformed into a symmetric eigenvalue problem by introducing the new coordinate vector:

$$\mathbf{q}=\mathbf{L}^T\mathbf{x}$$

$\mathbf{L}$ is the Cholesky decomposition of the symmetric mass matrix, i.e. $\mathbf{M}=\mathbf{L}\mathbf{L}^T$.

The equation of motion becomes:

$$\ddot{\mathbf{q}} + \tilde{\mathbf{K}} \mathbf{q} = 0$$

Compute $\tilde{\mathbf{K}}$.

In [19]:
L = np.linalg.cholesky(M)
L
Out[19]:
array([[63.24555,  0.     ,  0.     ,  0.     ],
       [ 0.     , 63.24555,  0.     ,  0.     ],
       [ 0.     ,  0.     , 63.24555,  0.     ],
       [ 0.     ,  0.     ,  0.     , 63.24555]])
In [20]:
M**0.5
Out[20]:
array([[63.24555,  0.     ,  0.     ,  0.     ],
       [ 0.     , 63.24555,  0.     ,  0.     ],
       [ 0.     ,  0.     , 63.24555,  0.     ],
       [ 0.     ,  0.     ,  0.     , 63.24555]])
In [21]:
import numpy.linalg as la
In [22]:
from numpy.linalg import inv
In [23]:
K_tilde = inv(L) @ K @ inv(L.T)
In [24]:
K_tilde
Out[24]:
array([[ 2.5 , -1.25,  0.  ,  0.  ],
       [-1.25,  2.5 , -1.25,  0.  ],
       [ 0.  , -1.25,  2.5 , -1.25],
       [ 0.  ,  0.  , -1.25,  1.25]])

Notice that $\tilde{\mathbf{K}}$ is symmetric, so we are guaranteed to get real eigenvalues and orthogonal eigenvectors when solving this system.

Exercise

Find the eigenvalues and eigenvectors. Create the spectral matrix $\mathbf{\Lambda}$ and the matrix $P$ which contains the orthonormal eigenvectors of $\tilde{\mathbf{K}}$.

$$ \mathbf{P} = \left[ \mathbf{v}_1, \ldots, \mathbf{v}_4 \right] $$
In [25]:
evals, evecs = np.linalg.eig(K_tilde)
evals
Out[25]:
array([4.41511, 2.93412, 1.25   , 0.15077])
In [26]:
evecs
Out[26]:
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]])
In [27]:
Lambda = np.diag(evals)
Lambda
Out[27]:
array([[4.41511, 0.     , 0.     , 0.     ],
       [0.     , 2.93412, 0.     , 0.     ],
       [0.     , 0.     , 1.25   , 0.     ],
       [0.     , 0.     , 0.     , 0.15077]])
In [28]:
P = evecs

Exercise

Prove that the eigenvectors in $\mathbf{P}$ are orthonormal.

In [29]:
np.dot(P[:, 0], P[:, 1])
Out[29]:
-2.7755575615628914e-16
In [30]:
np.linalg.norm(P[:, 0])
Out[30]:
1.0
In [31]:
P[:, 0].T @ P[:, 1]
Out[31]:
-2.7755575615628914e-16
In [32]:
P[:, 0].T @ P[:, 0]
Out[32]:
1.0

An orthonormal matrix has the property that its transpose multiplied by itself is the identity matrix.

In [33]:
P.T @ P
Out[33]:
array([[ 1., -0.,  0., -0.],
       [-0.,  1.,  0.,  0.],
       [ 0.,  0.,  1., -0.],
       [-0.,  0., -0.,  1.]])

Exercise

Find the natural freqencies of the system in both radians per second and Hertz, store them in an array in the order of the eigenvalues with names ws and fs.

In [34]:
ws = np.sqrt(evals)
ws
Out[34]:
array([2.10122, 1.71293, 1.11803, 0.38829])
In [35]:
fs = ws / 2 / np.pi
fs
Out[35]:
array([0.33442, 0.27262, 0.17794, 0.0618 ])

Exercise

Transform the eigenvectors back into the coordinate system associated with $\mathbf{x}$.

$$ \mathbf{S} = \left[ \mathbf{u}_1, \ldots, \mathbf{u}_4 \right] $$
In [36]:
S = np.linalg.inv(L.T) @ P
S
Out[36]:
array([[ 0.00678,  0.01038, -0.00913,  0.00361],
       [-0.01038, -0.00361, -0.00913,  0.00678],
       [ 0.00913, -0.00913, -0.     ,  0.00913],
       [-0.00361,  0.00678,  0.00913,  0.01038]])
In [37]:
sys.coordinates
Out[37]:
{'x1': 0.001, 'x2': 0.01, 'x3': 0.02, 'x4': 0.025}

Exercise: visualize the modeshapes

The eigenmodes (mode shapes) are contained in each column of $\mathbf{S}$. Create a plot for each mode shape with these specifications:

  • The title of each plot should be the frequency of the corresponding modeshape in Hz.
  • The y axis should be made up of the values [0, 3, 6, 9, 12] meters.
  • The x axis should plot the five values. The first should be zero and the remaining values should be the components of the mode shape in order of the component associated with the lowest floor to the highest.
  • Plot lines with small circles at each data point.
In [38]:
S[:, 0]
Out[38]:
array([ 0.00678, -0.01038,  0.00913, -0.00361])
In [39]:
np.hstack((0, S[:, 0]))
Out[39]:
array([ 0.     ,  0.00678, -0.01038,  0.00913, -0.00361])
In [40]:
u1 = S[:, 0]
u1
Out[40]:
array([ 0.00678, -0.01038,  0.00913, -0.00361])
In [41]:
u1[::-1]
Out[41]:
array([-0.00361,  0.00913, -0.01038,  0.00678])
In [42]:
S[:, 2]
Out[42]:
array([-0.00913, -0.00913, -0.     ,  0.00913])
In [43]:
fig, axes = plt.subplots(1, 4)

for i in range(4):
    axes[i].plot(np.hstack((0, S[:, i])), [0, 3, 6, 9, 12], marker='o')
    axes[i].set_title('{:1.2f} Hz'.format(fs[i]))
    
plt.tight_layout()
In [44]:
fs[0]
Out[44]:
0.3344190049004654
In [45]:
S[:, 0]
Out[45]:
array([ 0.00678, -0.01038,  0.00913, -0.00361])
In [46]:
sys.coordinates['x1'] = S[0, 2]
sys.coordinates['x2'] = S[1, 2]
sys.coordinates['x3'] = S[2, 2]
sys.coordinates['x4'] = S[3, 2]
In [47]:
traj = sys.free_response(30, sample_rate=10)
In [48]:
traj[list(sys.coordinates.keys())].plot(subplots=True)
Out[48]:
array([<matplotlib.axes._subplots.AxesSubplot object at 0x7faa6be0ec70>,
       <matplotlib.axes._subplots.AxesSubplot object at 0x7faa6bfc3fd0>,
       <matplotlib.axes._subplots.AxesSubplot object at 0x7faa6be080d0>,
       <matplotlib.axes._subplots.AxesSubplot object at 0x7faa6bdcfd90>], dtype=object)
In [49]:
sys.animate_configuration(fps=10)
Out[49]:
<matplotlib.animation.FuncAnimation at 0x7faa6bcdc280>

Simulating the trajectory

The trajectory of building's coordinates can be found with:

$$ \mathbf{x}(t) = \sum_{i=1}^n c_i \sin(\omega_i t + \phi_i) \mathbf{u}_i $$

where

$$ \phi_i = \arctan \frac{\omega_i \mathbf{v}_i^T \mathbf{q}_0}{\mathbf{v}_i^T \dot{\mathbf{q}}_0} $$

and

$$ c_i = \frac{\mathbf{v}^T_i \mathbf{q}_0}{\sin\phi_i} $$

$c_i$ are the modal participation factors and reflect what proportion of each mode is excited given specific initial conditions. If the initial conditions are the eigenmode, $\mathbf{u}_i$, the all but the $i$th $c_i$ will be zero.

Exercise

Show that if $\mathbf{q}_0 = \mathbf{v}_i$ then $c_i = 1$ all other modal participation factors are 0. Also, report all of the phase angles, $\phi_i$, in degrees.

In [50]:
for i in range(4):
    x0 = S[:, i]
    xd0 = np.zeros(4)
    print(x0)

    q0 = L.T @ x0
    qd0 = L.T @ xd0

    phis = np.arctan2(ws * P.T @ q0, P.T @ xd0)
    print(np.rad2deg(phis))

    cs = P.T @ q0 / np.sin(phis)
    print(cs)
    print('=' * 40)

[ 0.00678 -0.01038  0.00913 -0.00361]
[90. 90. 90. 90.]
[ 1. -0.  0. -0.]
========================================
[ 0.01038 -0.00361 -0.00913  0.00678]
[ 90.  90. -90. -90.]
[-0.  1. -0. -0.]
========================================
[-0.00913 -0.00913 -0.       0.00913]
[ 90. -90.  90. -90.]
[ 0. -0.  1.  0.]
========================================
[0.00361 0.00678 0.00913 0.01038]
[ 90. -90. -90.  90.]
[-0. -0.  0.  1.]
========================================

Exercise

Create a function called simulate() that returns the trajectories of the coordinates given an array of monotonically increasing time values and the initial conditions of the system.

It should look like:

def simulate(t, x0, xd0):
    """Returns the state trajectory.

    Parameters
    ==========
    t : ndarray, shape(m,)
        Monotonic values of time.
    x0 : ndarray, shape(n,)
        The initial conditions of each coordinate.
    xd0 : ndarray, shape(n,)
        The initial conditions of each speed.

    Returns
    =======
    x : ndarray, shape(m, n)
        The trajectories of each state.

    """

    # your code here

    return x
In [51]:
def simulate(t, x0, xd0):
    
    q0 = L.T @ x0
    qd0 = L.T @ xd0
    phis = np.arctan2(ws * P.T @ q0, P.T @ xd0)
    cs = P.T @ q0 / np.sin(phis)
    
    x = np.zeros((len(x0), len(t)))
    for ci, wi, phii, ui in zip(cs, ws, phis, S.T):
        x += ci * np.sin(wi * t + phii) * np.tile(ui, (len(t), 1)).T
    
    return x

Exercise

Using the plotting function below, show that the results found here are the same as the simulations from the FourStoryBuildingSystem given the same initial conditions.

In [52]:
def plot_trajectories(t, x):
    
    fig, axes = plt.subplots(4, 1)
    
    for i, ax in enumerate(axes.flatten()):
        ax.plot(t, x[i])
        ax.set_ylabel(r'$x_{}$ [m]'.format(i + 1))
    ax.set_xlabel('Time [s]')
    
    plt.tight_layout()
In [53]:
t = np.linspace(0, 50, num=50 * 60)
x0 = np.array([0.001, 0.010, 0.020, 0.025])
xd0 = np.zeros(4)
x = simulate(t, x0, xd0)
plot_trajectories(t, x)

This shows the plot of a single mode:

In [54]:
x = simulate(t, S[:, 0], np.zeros(4))
plot_trajectories(t, x)