NTS-02. On Rayleigh damping coefficients for FE analysis

Note to self. How to compute Rayleigh damping coefficients for given damping ratios $\xi_1$ and $\xi_2$ at frequencies $f_1$ and $f_2$.

This is textbook content, I just need to remind myself too often how this is done and end up re-deriving the equations.

Given the second-order system of differential equations representing the FE model

$$M \ddot{u} + C \dot{u} + K u = F(t)$$

The damping matrix $C$ can be written as a Rayleigh damping matrix:

$$C = a_0 M + a_1 K$$

$a_0$ and $a_1$ are Rayleigh damping coefficients found by solving

$$\left[ \begin{array}{cc} \dfrac{1}{2\pi f_1} & 2 \pi f_1 \ \dfrac{1}{2\pi f_2} & 2 \pi f_2 \end{array} \right] \left[ \begin{array}{c} a_0 \ a_1 \end{array} \right]¶ \left[ \begin{array}{c} \xi_1 \ \xi_2 \end{array} \right]$$

Which I do in the following code:

In [6]:
import sympy  as sym

xi1 = sym.S("xi_1")
xi2 = sym.S("xi_2")
xi = sym.S("xi")

f1 = sym.S("f_1")
f2 = sym.S("f_2")

w1 = 2*sym.pi*f1
w2 = 2*sym.pi*f2

A = sym.Matrix([[ 1/w1, w1],[1/w2, w2]])/2
b = sym.Matrix([[xi1],[xi2] ])

a = A**-1*b

print sym.pretty(sym.simplify(a[0]))
print sym.simplify(a[0])
print sym.pretty(sym.simplify(a[1]))
print sym.simplify(a[1])

4⋅π⋅f₁⋅f₂⋅(f₁⋅ξ₂ - f₂⋅ξ₁)
─────────────────────────
2     2
f₁  - f₂
4*pi*f_1*f_2*(f_1*xi_2 - f_2*xi_1)/(f_1**2 - f_2**2)
f₁⋅ξ₁ - f₂⋅ξ₂
─────────────
⎛  2     2⎞
π⋅⎝f₁  - f₂ ⎠
(f_1*xi_1 - f_2*xi_2)/(pi*(f_1**2 - f_2**2))


$$a_0 = \frac{4 \pi f_{1} f_{2}}{f_{1}^{2} - f_{2}^{2}} \left(f_{1} \xi_{2} - f_{2} \xi_{1}\right)$$ $$a_1 = \frac{f_{1} \xi_{1} - f_{2} \xi_{2}}{\pi \left(f_{1}^{2} - f_{2}^{2}\right)}$$

Damping at other frequencies is given by,

$$\xi(f) = \dfrac{a_0}{2\pi f} + \pi a_1 f$$

For example for $f_1 = 0.3$ Hz and $f_2 = 1.2$ Hz, and $\xi_1 = 0.05$, $\xi_2 = 0.02$.

In [25]:
import scipy as sp
import matplotlib.pyplot as plt

f_1 = 0.3
f_2 = 1.2
xi_1 = 0.05
xi_2 = 0.02

a_0 = 4*pi*f_1*f_2*(f_1*xi_2 - f_2*xi_1)/(f_1**2 - f_2**2)
a_1 = (f_1*xi_1 - f_2*xi_2)/(pi*(f_1**2 - f_2**2))

f = sp.linspace(0.0001, 3, 100)
xi = a_0/(4*pi*f) + pi*a_1*f

plt.figure().set_size_inches([12,6],forward = True)
plt.plot(f,xi*100, linewidth =2)
plt.plot(f_1, xi_1*100, "ro")
plt.plot(f_2, xi_2*100, "ro")
plt.xlabel(r"$f$")
plt.ylabel(r"Damping Ratio, $\xi$, [% of critical]")
plt.ylim([0,10])
plt.grid()