Complex vibration modes arise in experimental research and numerical simulations with non proportional damping. In such cases the eigenproblem is nonlinear, thus requiring additional work in order to obtain the vibration frequencies and mode shapes of the dynamic system.
The solution to the free vibration response of the dynamic equilibrium equation involves some form of linearization of the quadratic eigenvalue problem it yields. For more information on this subject see for example this report.
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
Considering a dynamic system with $NDOF$ is the number of generalized degrees of freedom, the vectors will have dimensions of $NDOF \times 1$ and the matrices $NDOF \times NDOF$.
When the system is undamped, the damping matrix will be null. In this case the eigenproblem is linear:
\begin{equation} \left[ -\mathbf{M} \times \mathbf{\omega^2} + \mathbf{K} \right] \times \mathbf{v} = \mathbf{0} \end{equation}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}where $\alpha$ and $\beta$ are mass and stiffness proportionality coefficients. These are typically very small positive numbers. The resulting eigenproblem is still linear because the damping matrix can be decomposed by the modal vectors.
In a system with non proportional damping, the damping matrix will not be proportional neither to the mass nor to the stiffness matrices.
In [3]:
MM = np.matrix(np.diag([1., 2.]))
print(MM)
In [4]:
KK = np.matrix([[20., -10.], [-10., 10.]])
print(KK)
We will also setup two damping matrices, one proportional to the mass and stiffness matrices (C1) and the other non proportional (C2):
In [5]:
C1 = 0.1*MM+0.04*KK
print(C1)
In [6]:
C2 = np.matrix([[0.1, 0.2], [0.2, 0.2]])
print(C2)
In [7]:
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 [8]:
print(np.sqrt(W2))
The modal vectors, the columns of the modal matrix, have unit norm:
In [9]:
print(LA.norm(F1, axis=0))
Contrary to what is normaly done, we will visualize the modal vectors in a polar plot of the corresponding amplitudes and angles of equivalent complex values:
In [10]:
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(F1[dof,mode])])
t = np.array([0, np.angle(F1[dof,mode])])
ax[mode].plot(t, r, 'o-', label='DOF #{}'.format(dof+1))
plt.legend(loc='lower left', bbox_to_anchor=(1., 0.))
plt.show()
This damping matrix is orthogonal because the mass and stiffness matrices are also orthogonal:
In [11]:
print(np.round_(F1.T*C1*F1, 4))
The system and input matrices are the following:
In [12]:
A = np.bmat([[np.zeros_like(MM), MM], [MM, C1]])
print(A)
In [13]:
B = np.bmat([[MM, np.zeros_like(MM)], [np.zeros_like(MM), -KK]])
print(B)
The eigenanalysis yields the eigenvalues and eigenvectors:
In [14]:
w1, v1 = LA.eig(LA.solve(A,B))
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 [15]:
print(np.round_(w1[::2], 4))
These complex eigenvalues can be decomposed into angular frequency and damping coefficient:
In [16]:
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 [17]:
print(LA.norm(v1[:,::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 [18]:
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 [19]:
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()
Non proportinal damping carries the fact that the damping matrix is not orthogonal anymore:
In [20]:
print(np.round_(F1.T*C2*F1, 4))
The system and input matrices are the following:
In [21]:
A = np.bmat([[np.zeros_like(MM), MM], [MM, C2]])
print(A)
In [22]:
B = np.bmat([[MM, np.zeros_like(MM)], [np.zeros_like(MM), -KK]])
print(B)
The eigenanalysis yields the eigenvalues and eigenvectors of the system matrix:
In [23]:
w2, v2 = LA.eig(LA.solve(A,B))
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 [24]:
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 [25]:
zw = -w2.real # damping coefficient times 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 [26]:
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 [27]:
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 [28]:
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.