Complex vibration modes arise in experimental research and numerical simulations when non proportional damping is adopted. In such cases a state space formulation of the second order differential dynamic equilibrium equation is the preferred way to adress the problem.
This notebook is inspired in one of Pete Avitabile's Modal Space articles, namely the one discussing the difference between complex modes and real normal modes. Additional information about state space formulation in structural dynamics can be found for example here.
Non proportionally damped system
We will start by setting up the computational environment for this notebook. Furthermore, we will need numpy and scipy for the numerical simulations and matplotlib for the plots:
In [1]:
import sys
import numpy as np
import scipy as sp
import matplotlib as mpl
print('System: {}'.format(sys.version))
print('numpy version: {}'.format(np.__version__))
print('scipy version: {}'.format(sp.__version__))
print('matplotlib version: {}'.format(mpl.__version__))
We will also need a couple of specific modules and a litle "IPython magic" to show the plots:
In [2]:
from numpy import linalg as LA
import matplotlib.pyplot as plt
%matplotlib inline
In structural dynamics the second order differential dynamic equilibrium equation can be written in terms of generalized coordinates (d[isplacement]) and their first (v[elocity]) and second (a[cceleration]) time derivatives:
\begin{equation} \mathbf{M} \times \mathbf{a(t)} + \mathbf{C} \times \mathbf{v(t)} + \mathbf{K} \times \mathbf{d(t)} = \mathbf{F(t)} \end{equation}where:
$\mathbf{M}$ is the mass matrix
$\mathbf{C}$ is the damping matrix
$\mathbf{K}$ is the stiffness matrix
$\mathbf{a(t)}$ is the acceleration vector
$\mathbf{v(t)}$ is the velocity vector
$\mathbf{d(t)}$ is the displacement vector
$\mathbf{F(t)}$ is the force input vector
All these matrices are of size $NDOF \times NDOF$, where $NDOF$ is the number of generalized degrees of freedom of the dynamic system.
In a state space formulation the second order differential dynamic equilibrium equation is turned into a system of first order differential dynamic equilibrium equations:
\begin{equation} \begin{matrix} \mathbf{\dot{x}(t)} = \mathbf{A} \cdot \mathbf{x(t)} + \mathbf{B} \cdot \mathbf{u(t)} \\ \mathbf{y(t)} = \mathbf{C} \cdot \mathbf{x(t)} + \mathbf{D} \cdot \mathbf{u(t)} \end{matrix} \end{equation}where
$\mathbf{A}$ is the system matrix
$\mathbf{B}$ is the imput matrix
$\mathbf{C}$ is the output matrix
$\mathbf{D}$ is the feedthrough matrix
$\mathbf{x(t)}$ is the state vector
$\mathbf{y(t)}$ is the output vector
$\mathbf{u(t)}$ is the input vector
The state vector, of size $2 \times NDOF$ by $1$, has the following form:
\begin{equation} \mathbf{x(t)} = \left[ \begin{matrix} \mathbf{u(t)} \\ \mathbf{\dot{u}(t)} \end{matrix} \right] \end{equation}The system matrix, of size $2 \times NDOF$ by $2 \times NDOF$, is built using the M, C and K matrices:
\begin{equation*} \mathbf{A} = \left[ \begin{matrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{M}^{-1} \cdot \mathbf{K} & -\mathbf{M}^{-1} \cdot \mathbf{C} \end{matrix} \right] \end{equation*}The loading matrix, of size $2 \times NDOF$ by $1$, is composed of 0's and 1's.
In [3]:
MM = np.matrix(np.diag([2,3]))
print(MM)
In [4]:
KK = np.matrix([[2, -1],[-1, 1]])
print(KK)
Let us perform the eigenanalysis of the (undamped) second order differential dynamic equilibrium equation for later comparison of results:
In [5]:
W2, F1 = LA.eig(LA.solve(MM,KK)) # eigenanalysis
ix = np.argsort(np.absolute(W2)) # sort eigenvalues in ascending order
W2 = W2[ix] # sorted eigenvalues
F1 = F1[:,ix] # sorted eigenvectors
print(np.round_(W2, 4))
print(np.round_(F1, 4))
The angular frequencies are computed as the square root of the eigenvalues:
In [6]:
print(np.sqrt(W2))
The modal vectors, the columns of the modal matrix, have unit norm:
In [7]:
print(LA.norm(F1, axis=0))
In [8]:
C0 = np.matrix(np.zeros_like(MM))
print(C0)
The system matrix is the following:
In [9]:
A = np.bmat([[np.zeros_like(MM), np.identity(MM.shape[0])], [LA.solve(-MM,KK), LA.solve(-MM,C0)]])
print(A)
Performing the eigenanalysis on this matrix yields the following complex valued results:
In [10]:
w0, v0 = LA.eig(A)
ix = np.argsort(np.absolute(w0))
w0 = w0[ix]
v0 = v0[:,ix]
print(np.round_(w0, 4))
print(np.round_(v0, 4))
As we can see, the eigenvalues come in complex conjugate pairs. Therefore we can take for instance only the ones in the upper half-plane:
In [11]:
print(np.round_(w0[[0,2]], 4))
In this case, since damping is zero, the real part of the complex eigenvalues is also zero (apart from round-off errors) and the imaginary part is equal to the angular frequency computed previously for the dynamic system:
In [12]:
print(w0[[0,2]].imag)
The columns of the modal matrix, the modal vectors, also come in conjugate pairs. Each vector has unit norm as in the dynamic system:
In [13]:
print(LA.norm(v0[:,[0,2]], axis=0))
Moreover, we can check that the modal matrix is composed of four blocks, each with $NDOF \times NDOF$ dimension. Some column reordering is necessary in order to match both modal matrices:
In [14]:
AA = v0[:2,[0,2]]
AB = AA.conjugate()
BA = np.multiply(AA,w0[[0,2]])
BB = BA.conjugate()
v0_new = np.bmat([[AA, AB], [BA, BB]])
print(np.round_(v0_new[:,[0,2,1,3]], 4))
To help visualize the complex valued modal vectors we will do a polar plot of the corresponding amplitudes and angles:
In [15]:
fig, ax = plt.subplots(1, 2, subplot_kw=dict(polar=True))
for mode in range(2):
ax[mode].set_title('Mode #{}'.format(mode+1))
for dof in range(2):
r = np.array([0, np.absolute(v0[dof,2*mode])])
t = np.array([0, np.angle(v0[dof,2*mode])])
ax[mode].plot(t, r, 'o-', label='DOF #{}'.format(dof+1))
plt.legend(loc='lower left', bbox_to_anchor=(1., 0.))
plt.show()
In a proportionally damped system, the damping matrix is proportional to the mass and stiffness matrices:
\begin{equation*} \mathbf{C} = \alpha \times \mathbf{M} + \beta \times \mathbf{K} \end{equation*}Let us assume $\alpha$ to be 0.1 and $\beta$ to be 0.1. This yields the following damping matrix:
In [16]:
C1 = 0.1*MM+0.1*KK
print(C1)
This damping matrix is orthogonal because the mass and stiffness matrices are also orthogonal:
In [17]:
print(np.round_(F1.T*C1*F1, 4))
The system matrix is the following:
In [18]:
A = np.bmat([[np.zeros_like(MM), np.identity(MM.shape[0])], [LA.solve(-MM,KK), LA.solve(-MM,C1)]])
print(A)
The eigenanalysis yields the eigenvalues and eigenvectors:
In [19]:
w1, v1 = LA.eig(A)
ix = np.argsort(np.absolute(w1))
w1 = w1[ix]
v1 = v1[:,ix]
print(np.round_(w1, 4))
print(np.round_(v1, 4))
As we can see, the eigenvalues come in complex conjugate pairs. Let us take only the ones in the upper half-plane:
In [20]:
print(np.round_(w1[[0,2]], 4))
These complex eigenvalues can be decomposed into angular frequency and damping coefficient:
In [21]:
zw = -w1.real # damping coefficient time angular frequency
wD = w1.imag # damped angular frequency
zn = 1./np.sqrt(1.+(wD/-zw)**2) # the minus sign is formally correct!
wn = zw/zn # undamped angular frequency
print('Angular frequency: {}'.format(wn[[0,2]]))
print('Damping coefficient: {}'.format(zn[[0,2]]))
The columns of the modal matrix, the modal vectors, also come in conjugate pairs, each vector having unit norm:
In [22]:
print(LA.norm(v1[:,[0,2]], axis=0))
Moreover, the modal matrix is composed of four blocks, each with $NDOF \times NDOF$ dimension. Some column reordering is necessary in order to match both modal matrices:
In [23]:
AA = v1[:2,[0,2]]
AB = AA.conjugate()
BA = np.multiply(AA,w1[[0,2]])
BB = BA.conjugate()
v1_new = np.bmat([[AA, AB], [BA, BB]])
print(np.round_(v1_new[:,[0,2,1,3]], 4))
We will visualize again the complex valued modal vectors with a polar plot of the corresponding amplitudes and angles:
In [24]:
fig, ax = plt.subplots(1, 2, subplot_kw=dict(polar=True))
for mode in range(2):
ax[mode].set_title('Mode #{}'.format(mode+1))
for dof in range(2):
r = np.array([0, np.absolute(v1[dof,2*mode])])
t = np.array([0, np.angle(v1[dof,2*mode])])
ax[mode].plot(t, r, 'o-', label='DOF #{}'.format(dof+1))
plt.legend(loc='lower left', bbox_to_anchor=(1., 0.))
plt.show()
In [25]:
C2 = np.matrix([[0.4, -0.1],[-0.1, 0.1]])
print(C2)
Non proportinal damping carries the fact that the damping matrix is not orthogonal anymore:
In [26]:
print(np.round_(F1.T*C2*F1, 4))
The system matrix is the following:
In [27]:
A = np.bmat([[np.zeros_like(MM), np.identity(MM.shape[0])], [LA.solve(-MM,KK), LA.solve(-MM,C2)]])
print(A)
The eigenanalysis yields the eigenvalues and eigenvectors:
In [28]:
w2, v2 = LA.eig(A)
ix = np.argsort(np.absolute(w2))
w2 = w2[ix]
v2 = v2[:,ix]
print(np.round_(w2, 4))
print(np.round_(v2, 4))
As we can see, the eigenvalues come in complex conjugate pairs. Again, let us take only the ones in the upper half-plane:
In [29]:
print(np.round_(w2[[0,2]], 4))
These complex eigenvalues can be decomposed into angular frequency and damping coefficient much like in the propotional damping case:
In [30]:
zw = -w2.real # damping coefficient time angular frequency
wD = w2.imag # damped angular frequency
zn = 1./np.sqrt(1.+(wD/-zw)**2) # the minus sign is formally correct!
wn = zw/zn # undamped angular frequency
print('Angular frequency: {}'.format(wn[[0,2]]))
print('Damping coefficient: {}'.format(zn[[0,2]]))
Again, the columns of the modal matrix, the modal vectors, come in conjugate pairs, and each vector has unit norm:
In [31]:
print(LA.norm(v2[:,[0,2]], axis=0))
Moreover, the modal matrix is composed of four blocks, each with $NDOF \times NDOF$ dimension. Some column reordering is necessary in order to match both modal matrices:
In [32]:
AA = v2[:2,[0,2]]
AB = AA.conjugate()
BA = np.multiply(AA,w2[[0,2]])
BB = BA.conjugate()
v2_new = np.bmat([[AA, AB], [BA, BB]])
print(np.round_(v2_new[:,[0,2,1,3]], 4))
Once more we will visualize the complex valued modal vectors through a polar plot of the corresponding amplitudes and angles:
In [33]:
fig, ax = plt.subplots(1, 2, subplot_kw=dict(polar=True))
for mode in range(2):
ax[mode].set_title('Mode #{}'.format(mode+1))
for dof in range(2):
r = np.array([0, np.absolute(v2[dof,2*mode])])
t = np.array([0, np.angle(v2[dof,2*mode])])
ax[mode].plot(t, r, 'o-', label='DOF #{}'.format(dof+1))
plt.legend(loc='lower left', bbox_to_anchor=(1., 0.))
plt.show()
Several conclusion can be drawn from this very simple example. First of all, damping changes the vibration frequencies and mode shapes. Furthermore, the polar plots show clearly that:
This notebook was created by Paulo Xavier Candeias.