# Be Careful with Modal Damping

Modal damping is kind of the it-spell in the dark art that is modeling viscous damping in structures. Although modal damping is pretty straightforward, you should be aware of an important aspect of its implementation in OpenSees.

The issue, which is described in section 9 of this paper, is that OpenSees assembles the dynamic tangent into the matrix storage scheme you choose via the ‘system’ command. This is a problem for modal damping because the damping matrix is full even when the matrix storage, based on the nodal connectivity, is banded, profile, or sparse. As a result, some terms of the dynamic tangent will be zero even though they are non-zero in the modal damping matrix. Convergence problems will ensue–even for linear problems–because the dynamic tangent will be inconsistent with the residual.

Let’s describe this issue with a canonical shear frame model.

Using floor mass $m=1.0352$ and story stiffness $k=610$, the model has the following mass and stiffness matrices ${\bf M}=\left[ \begin{array}{cccc} 1.0352 & 0 & 0 & 0 \\ 0 & 1.0352 & 0 & 0 \\ 0 & 0 & 1.0352 & 0 \\ 0 & 0 & 0 & 0.5176 \end{array}\right] \hspace{0.1in} {\bf K}=\left[ \begin{array}{cccc} 1220 & 0 & -610 & 0 \\ 0 & 1220 & -610 & -610 \\ -610 & -610 & 1220 & 0 \\ 0 & -610 & 0 & 610 \end{array}\right]$

Note that the equations are numbered 1-3-2-4 (use the ‘Plain’ equation numberer) so that the stiffness matrix is not tridiagonal–its bandwidth is 3, the profile dips inside the band, and there’s sparsity within the profile. Accordingly, the ‘BandSPD’, ‘ProfileSPD’, and ‘UmfPack’ systems will store different numbers of non-zeros.

Now suppose we want to use 5% damping ( $\zeta=0.05$) in the first mode. Based on the mass and stiffness, the first mode frequency is $\omega_1=9.4715$ rad/sec and the corresponding eigenvector and the resulting modal damping matrix are $\phi_1=\left[ \begin{array}{c} 0.26596 \\ 0.64208 \\ 0.49143 \\ 0.69498 \end{array} \right] \hspace{0.1in} {\bf C}=2\zeta\omega_1({\bf M}\phi_1)({\bf M}\phi_1)^T=\left[ \begin{array}{cccc} 0.07179465 & 0.17332762 & 0.13265922 & 0.09380423 \\ 0.17332762 & 0.4184499 & 0.32026769 & 0.22646345 \\ 0.13265922 & 0.32026769 & 0.24512228 & 0.17332762 \\ 0.09380423 & 0.22646345 & 0.17332762 & 0.12256114 \end{array} \right]$

(Update April 11, 2022: corrected modal damping matrix to include mass matrix)

Assuming Newmark time integration with constant average acceleration ( $\gamma=0.5$ and $\beta=0.25$) and time step $\Delta t=0.01$ sec, the consistent dynamic tangent is ${\bf K}_T = {\bf K} + \frac{\gamma}{\beta \Delta t}{\bf C} + \frac{1}{\beta(\Delta t)^2}{\bf M} = \left[ \begin{array}{cccc} 4.26423589e+04 & 3.46655246e+01 & -5.83468156e+02 & 1.87608467e+01 \\ 3.46655246e+01 & 4.27116900e+04 & -5.45946463e+02 & -5.64707309e+02 \\ -5.83468156e+02 & -5.45946463e+02 & 4.26770245e+04 & 3.46655246e+01 \\ 1.87608467e+01 & -5.64707309e+02 & 3.46655246e+01 & 2.13385122e+04 \end{array} \right]$

With banded storage, the dynamic tangent that OpenSees will use in the analysis is ${\bf K}_T = \left[ \begin{array}{cccc} 4.26423589e+04 & 3.46655246e+01 & -5.83468156e+02 & 0 \\ 3.46655246e+01 & 4.27116900e+04 & -5.45946463e+02 & -5.64707309e+02 \\ -5.83468156e+02 & -5.45946463e+02 & 4.26770245e+04 & 3.46655246e+01 \\ 0 & -5.64707309e+02 & 3.46655246e+01 & 2.13385122e+04 \end{array} \right]$

Notice the zeros at the upper right and lower left corners of the matrix, outside the bandwidth of 3. Because of these zeros, the dynamic tangent used in the analysis will be inconsistent with the residual and your analyses won’t converge in one iteration for linear problems.

Consider the banded matrix solver ‘BandSPD’ in a free vibration analysis of the shear frame. The analysis takes four iterations to find equilibrium at each time step.

Although not shown here, the analysis will take three or four iterations per time step if you switch the solver to ‘ProfileSPD’ or ‘UmfPack’, or any solver that does not store the full matrix.

Because it stores the entire matrix, you can use the ‘FullGeneral’ solver to get the convergence in one iteration that you expect for a linear problem.

Don’t worry though, the analyses will converge to the same solution regardless of which solver you use. However, without getting into spectral radii, it’s possible that the analysis will not converge at all due to inconsistent zeros in the dynamic tangent. But this is unlikely for realistic values of mass, stiffness, and damping.

For large models, using a solver with full matrix storage will be dreadfully slow. You’ll be better off taking the extra iterations with a fast sparse solver like ‘UmfPack’ than using the ‘FullGeneral’ solver to converge in one iteration. When using a fast sparse solver, don’t use the ‘NewtonRaphson’ algorithm–you don’t want to factorize the same inconsistent tangent over and over. Consider ‘KrylovNewton’ or ‘BFGS’ instead. And don’t use the ‘Linear’ algorithm because it won’t iterate to equilibrium.

Recognizing that one size does not fit all, OpenSees was designed to give you total control over the analysis. Therefore, it’s important you are aware of issues like this one with modal damping so that you can choose analysis options wisely.

Here is an OpenSees Tcl script of modal damping, without all the equation numbering and solver shenanigans.

#
# Free vibration analysis of shear frame with modal damping
# Using zero length elements for the model
#

wipe
model basic -ndm 1 -ndf 1

set k 610.0 ;# Story stiffness
set m 1.0352 ;# Floor mass

node 0 0.0; fix 0 1
node 1 0.0
node 2 0.0
node 3 0.0
node 4 0.0

mass 1 $m mass 2$m
mass 3 $m mass 4 [expr 0.5*$m]

uniaxialMaterial Elastic 1 $k element zeroLength 1 0 1 -mat 1 -dir 1 element zeroLength 2 1 2 -mat 1 -dir 1 element zeroLength 3 2 3 -mat 1 -dir 1 element zeroLength 4 3 4 -mat 1 -dir 1 # Impose displacement at roof timeSeries Constant 1 pattern Plain 1 1 { sp 4 1 1.0 } # Do a static analysis for initial condition integrator LoadControl 0.0 constraints Transformation analysis Static analyze 1 wipeAnalysis # Remove displacement so it doesn't affect eigen calculations remove loadPattern 1 set N 2 ;# Number of modes for modal damping eigen$N

modalDamping 0.05 0.02 ;# 5% in mode 1, 2% in mode 2

# Switch to dynamic analysis
analysis Transient

# Record roof displacement
recorder Node -file roofDisp.out -time -node 4 -dof 1 disp

# Analysis time, step, and num steps
set Tf 5.0
set dt 0.01
set Nsteps [expr int($Tf/$dt)]

analyze $Nsteps$dt

## 31 thoughts on “Be Careful with Modal Damping”

1. MSB says:

There are articles on the effect of viscous damping modeling methods on seismic performance using OpenSees.
As an example, https://www.sciencedirect.com/science/article/abs/pii/S2352012418300730 .

Like

1. Positive Definite says:

Hello MSB, thank you for the article!

It appears to address only Rayleigh damping and which stiffness matrix to use. There are additional papers on this topic:
https://doi.org/10.1002/eqe.541
https://doi.org/10.1061/(ASCE)0733-9445(2008)134:4(581)
https://doi.org/10.1002/eqe.2357

The following paper also addresses modal damping:
https://doi.org/10.1002/eqe.2622

Like

1. MSB says:

Thanks.

Like

2. Kris says:

Hi, can you show me an example with tcl script to explain how can I use modal Damping in OpenSees?
Thank you.

Like

1. Positive Definite says:

Hello Kris,
2. Perform an eigenvalue analysis
eigen(N) ;# N is the number of eigenmodes to use in modal damping
3. Issue the modal damping command
modalDamping(0.02) ;# 2% damping in all N modes
# or
modalDamping(0.02,0.05,…,0.04) ;# 2% damping in mode 1, 5% in mode 2, …, 4% in mode N
For the second option, you need to give N values

PD

Like

1. Kris says:

Thanks.

Like

2. ghxtju says:

Hi,when I use the Modaldamping in transient analysis with the KRAlphaExplicit integrator,
But it warning me with the code “IncrementalIntegrator::getVel() – not implemeneted for this integrator”. How can I fix this problem?
Thank you!

Like

3. Positive Definite says:

Please post errors and possible bugs on the OpenSees GitHub page.

Thanks!

PD

Like

3. Hamed Hasani says:

Hi,
When I use modal damping I can see the following error:
WARNING PenaltySp_FE::getM_Force() – not yet implemented
How I can fix this? The same model does not have any issue with rayleigh damping.
Thank you,

Like

1. Positive Definite says:

Hello Hamed,
Thanks for pointing out the issue. Please post this on the OpenSees GitHub “Issues” page: https://github.com/OpenSees/OpenSees/issues
In the mean time, try using a different constraint handler.
PD

Like

4. Jinghui Huang says:

Hi,

As the codes listed below, we need to remove loadPattern before the eigenanalysis, but will the removal of the loadPattern affect the transient analysis in the next step, as the loadPattern removed is often the gravity loads, which should be the static loads applied to the model through the analysis. Thanks.

# Remove displacement so it doesn’t affect eigen calculations

set N 2 ;# Number of modes for modal damping
eigen \$N

Like

1. Positive Definite says:

Hello Jinghui,
The load pattern being removed is for the initial displacement (sp command) used for free vibration. You are correct, if it was a gravity load, you would not want to remove the pattern.
PD

Like

1. Jinghui Huang says:

As you said, if the loadPattern defined before the eigenanalysis was for the gravity loads, it should not be removed before doing the eigenanalysis, as it will be used for the transient analysis, but if it was not removed, will the displacement caused by the gravity loads affect the eigenanalysis and then the modalDamping assigned? Thanks.

Like

2. Positive Definite says:

The gravity loads can cause a slight change in stiffness (especially for RC structures) that will in turn affect the eigenvalues. So, best to apply gravity loads first, then do eigenvalues and modal damping.

Like

5. Xavier says:

Hi,

Thank you for the post. Why modal damping matrix equation does not include multiplication of mass matrix terms? The manually calculated damping matrix does not match with it. Does OpenSees also not include multiplication of mass matrix terms while calculating damping forces ?

Liked by 1 person

1. Michael H. Scott says:

I should have been more careful with modal damping! I wrote this post before implementing the GimmeMCK integrator which would have made it easier to check what I was calculating against what OpenSees computes.
https://portwooddigital.com/2022/04/03/gimme-all-your-modal-damping/

Like

1. Xavier says:

Thank you for clarifying my doubt. I checked with GimmeMCK integrator, it is matching with manual calculation. I also watched the youtube video explaining the modal damping code. However I did not understand where the does this multiplication was occurring.

Like

2. Michael H. Scott says:

It’s done in IncrementalIntegrator::setupModal and IncrementalIntegator::doMv

Like

This site uses Akismet to reduce spam. Learn how your comment data is processed.