
In this lecture we'll use the techniques discussed in the previous two lectures to calculate damping matrices with specified levels of modal damping. We'll start by demonstrating how to calculate a Rayleigh damping matrix and then see how we can obtain a damping matrix with all modal damping specified. We'll apply these calculations to our 3-storey test structure from previous lectures. We first need to define some previously calculated quantities.

```python
#Define masses and the mass matrix
m1 = 21145.85
m2 = 20773.5
m3 = 20476.5
M = np.matrix([[m1,0,0],[0, m2,0],[0,0,m3]])
```

```python
#Define stiffnesses and the stiffness matrix
k1 = 2494164.77
k2 = 3226337.45
k3 = 4593750.0
K = np.matrix([[k1+k2, -k2, 0],
               [-k2, k2+k3, -k3],
               [0, -k3, k3]])
```

```python
#Define the mass normalised modal matrix
Phi = np.array(
[[ 0.00265937,  0.00577103,  0.00262936],
 [ 0.00423975,  0.00053654, -0.00546579],
 [ 0.00482667, -0.00376177,  0.00337473]])
```

```python
#Define the modal stiffness matrix
modalK = np.array([[27.28, 0, 0],
                   [0, 256.341, 0],
                   [0, 0, 587.694]])
```

### Rayleigh damping example

Let's determine a damping matrix that provides the following levels of modal damping, $\xi_1 = 4\%$, $\xi_2 = 2\%$. Remember, with Rayleigh damping we can only specify damping for two modes.

We start by calculating the modal frequencies.

```python
w1 = math.sqrt(modalK[0,0]) #(rads/s) Angular frequency for mode 1
w2 = math.sqrt(modalK[1,1]) #(rads/s) Angular frequency for mode 2
w3 = math.sqrt(modalK[2,2]) #(rads/s) Angular frequency for mode 3
```

Next we need to solve the simultaneous equations defined previously,

$$
\begin{Bmatrix}
\alpha_0\\
\alpha_1
\end{Bmatrix}=
2
\begin{bmatrix}
\frac{1}{\omega_i} & \omega_i\\
\frac{1}{\omega_j} & \omega_j\\
\end{bmatrix}^{-1}
\begin{Bmatrix}
\xi_i\\
\xi_j
\end{Bmatrix}
\tag{1}
$$

```python
xi1 = 0.04 #Damping ratio specified for mode 1
xi2 = 0.02 #Damping ratio specified for mode 2
dampingVector = np.array([[xi1],[xi2]])
alphaVector = 2*np.mat([[1/w1, w1],[1/w2, w2]]).I*dampingVector

a0 = alphaVector[0].item()
a1 = alphaVector[1].item()
```

Now we can construct the Rayleigh damping matrix.

```python
C_ray = a0*M + a1*K #Rayleigh damping matrix
print('The Rayleigh damping matrix is...')
print('')
print(C_ray)
```

The Rayleigh damping matrix is...

$$
\begin{bmatrix}
13833.82675777 & -3135.10899077  &   0\\
-3135.10899077 & 15728.33075358 & -4463.85635401\\
0 & -4463.85635401 & 12476.99572837
\end{bmatrix}
$$

Note that the resulting level of damping in the third mode can be determined from the expression we derived in the Rayleigh damping lecture,

$$
2\xi_3\omega_3 = \alpha_0 + \alpha_1\omega_3^2
\tag{2}
$$

Therefore,

$$
\xi_3 = \frac{1}{2}\left(\frac{\alpha_0}{\omega_3}+\alpha_1\omega_3\right)
\tag{3}
$$

```python
xi3 = 0.5*((a0/w3)+(a1*w3))
print('The damping ratio in mode 3 is {one}%'.format(one=round(100*xi3,2)))
```

The damping ratio in mode 3 is $1.98\%$.

### General determination example

Next, we will calculate a damping matrix that provides the following levels of modal damping, $\xi_1 = 4\%$, $\xi_2 = 2\%$ and $\xi_3 = 2\%$. The key thing to note here is that we are specifying damping levels in more than 2 modes. As we saw in the previous lecture, this method can be scaled to allow us to specify levels of damping in all $N$ modes of a structure.

We can start by constructing the matrix $\begin{bmatrix}\hat{C}\end{bmatrix}$ that we would ultimately like to end with after transformation.

```python
xi1 = 0.04 #Damping ratio specified for mode 1
xi2 = 0.02 #Damping ratio specified for mode 2
xi3 = 0.02 #Damping ratio specified for mode 3

c1_hat = 2*xi1*w1 #(Ns/m) Damping coefficient for damper 1
c2_hat = 2*xi2*w2 #(Ns/m) Damping coefficient for damper 2
c3_hat = 2*xi3*w3 #(Ns/m) Damping coefficient for damper 3

C_hat = np.array([[c1_hat, 0, 0],
                  [0, c2_hat, 0],
                  [0, 0, c3_hat]])
```

Next we can construct the required (pre-transformation) damping matrix using the equation derived in the previous lecture,

$$
\begin{bmatrix}
C
\end{bmatrix} =
\begin{bmatrix}
M
\end{bmatrix}
\begin{bmatrix}
\Phi
\end{bmatrix}
\begin{bmatrix}
\hat{C}
\end{bmatrix}
\begin{bmatrix}
\Phi
\end{bmatrix}^T
\begin{bmatrix}
M
\end{bmatrix}
\tag{4}
$$

```python
C_gen = M*Phi*C_hat*Phi.T*M

print('The damping matrix is...')
print('')
print(C_gen)
```

`The damping matrix is...`

$$
\begin{bmatrix}
13856.35354784 & -3181.13438688  &  28.00019939\\
-3181.13438688 & 15822.2691131  & -4521.03912885\\
28.00019939 & -4521.03912885 & 12511.79848761
\end{bmatrix}
$$

We can see the similarity between both damping matrices which is to be expected. However the second method obviously gives us more flexibility.

### Confirm damping matrix orthogonality

Finally we will confirm that this damping matrix is indeed diagonalised by a modal coordinate transformation (as per the mass and stiffness matrices).

```python
modalC = Phi.T*C_gen*Phi
print('The damping matrix after modal transformation is, ')
print('')
print(np.round(modalC,3))
```

`The damping matrix after modal transformation is, `

$$
\begin{bmatrix}
0.418 & 0  &  0\\
0 & 0.64 & 0\\
0 & 0 & 0.97
\end{bmatrix}
$$

And of course we can recover our originally specified damping ration directly from this modal damping matrix.

```python
xi1m = modalC[0,0]/(2*w1) xi2m = modalC[1,1]/(2*w2) xi3m = modalC[2,2]/(2*w3)

print("The 'recovered' modal damping ratios, obtained from the transformed modal matrix are:")
print(f"Mode 1: {round(xi1m,3)}")
print(f"Mode 2: {round(xi2m,3)}")
print(f"Mode 3: {round(xi3m,3)}")

```

`The 'recovered' modal damping ratios, obtained from the transformed modal matrix are:`

$$
\begin{align}
\text{Mode 1:}\: & 0.04\\
\text{Mode 2:}\: & 0.02\\
\text{Mode 3:}\: & 0.02\\
\end{align}
$$

Just to reinforce the point, let's say we naively tried to diagonalise a damping matrix of the form,

$$
C =
\begin{bmatrix}
c_1+c_2 & -c_2 & 0\\
-c_2 & c_2 +c_3 & -c_3\\
0 &-c_3 & c_3
\end{bmatrix}
\tag{5}
$$

where, $c_1=c_2=c_3=1000\:Ns/m$, without taking precautions to convert it into a diagonalisable format;

```python
#Define a damping matrix
C = np.matrix([[2000, -1000, 0],
 [-1000, 2000,  -1000],
 [0, -1000, 1000]])

#Try to diagonalise with modal coordinate transformation
Cm = np.round(Phi.T*C*Phi,3)

print('The naively transformed damping matrix is:')
print("")
print(Cm)
```

`The naively transformed damping matrix is:`

$$
\begin{bmatrix}
0.01 &  0.005 & -0.001\\
0.005 & 0.079 & 0.02\\
-0.001 & 0.02  & 0.151
\end{bmatrix}
$$

Clearly we can see that only after taking precautions to construct a _diagonalisable_ matrix, do we end up with an uncoupled damping matrix after transformation to modal coordinates.


Please log in or enroll to continue

Please log in or enroll to continue