Modal dynamic analysis

The modal dynamic procedure provides time history analysis of linear systems.

The following topics are discussed:

Solution to the uncoupled system
Solution to the coupled system
Integral coefficients
Structural damping in modal dynamic analysis
Response of nodal and element variables
Initial conditions
Base motion definition

Solution to the uncoupled system

The solution to the uncoupled equations is obtained readily as a particular integral for the loading and a solution to the homogeneous equation (with no right-hand side). These solutions can be combined and written in the general form

\begin{matrix} (3) & {\begin{matrix} q_{t + Δ t} \\ {\dot{q}}_{t + Δ t} \end{matrix}} = [\begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix}] {\begin{matrix} q_{t} \\ {\dot{q}}_{t} \end{matrix}} + [\begin{matrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{matrix}] {\begin{matrix} f_{t} \\ f_{t + Δ t} \end{matrix}}, \end{matrix}

where $a_{i j}$ and $b_{i j}$ , $i, j = 1, 2$ are constants, since we have assumed that the loading only varies linearly over the time increment (that is, $Δ f / Δ t$ is constant).

There are three cases of this solution for nonrigid body motion ( $ω \neq 0$ ), depending on whether the damping in the modal equilibrium equation is greater than, equal to, or less than critical damping (that is, depending on whether $(ξ^{2} - 1)$ is positive, zero, or negative).

Solution to the coupled system

Equation 3 can be generalized to address full coupling in the projected damping matrix. Let the matrix $C$ be split into its diagonal and off-diagonal parts so that

C = C_{d i a g} + C_{o f f} .

Then, with the additional assumption that the off-diagonal damping force varies linearly over a time increment, the equation for the uncoupled system can be rewritten as

\begin{matrix} (4) & \begin{matrix} {\begin{matrix} q_{t + Δ t} \\ \bar{C} {\dot{q}}_{t + Δ t} \end{matrix}} = [\begin{matrix} a_{11} & a_{12} - b_{11} C_{o f f} \\ a_{21} & a_{22} - b_{21} C_{o f f} \end{matrix}] {\begin{matrix} q_{t} \\ {\dot{q}}_{t} \end{matrix}} + [\begin{matrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{matrix}] {\begin{matrix} f_{t} \\ f_{t + Δ t} \end{matrix}} \\ + [\begin{matrix} 0 & - b_{12} C_{o f f} \\ 0 & 0 \end{matrix}] {\begin{matrix} q_{t + Δ t} \\ {\dot{q}}_{t + Δ t} \end{matrix}}, \end{matrix} \end{matrix}

where $\bar{C}$ is given by

\bar{C} = I + b_{22} C_{o f f} .

$\bar{C}$ is fully populated but is a function only of $Δ t$ , so it needs to be factored only a single time for a given analysis.

Integral coefficients

Damping less than critical

This case is the most common. With

\bar{ω} = ω \sqrt{1 - ξ^{2}},

we have

\begin{matrix} a_{11} = & \exp (- ξ ω Δ t) (ξ \frac{ω}{\bar{ω}} \sin \bar{ω} Δ t + \cos \bar{ω} Δ t), \\ a_{12} = & \exp (- ξ ω Δ t) \frac{1}{\bar{ω}} \sin \bar{ω} Δ t, \\ a_{21} = & - \exp (- ξ ω Δ t) \frac{ω}{\sqrt{1 - ξ^{2}}} \sin \bar{ω} Δ t, \\ a_{22} = & \exp (- ξ ω Δ t) (\cos \bar{ω} Δ t - \frac{ξ ω}{\bar{ω}} \sin \bar{ω} Δ t), \end{matrix}

\begin{matrix} b_{11} = & - \exp (- ξ ω Δ t) {(\frac{ξ}{ω \bar{ω}} + \frac{2 ξ^{2} - 1}{ω^{2} \bar{ω} Δ t}) \sin \bar{ω} Δ t + (\frac{1}{ω^{2}} + \frac{2 ξ}{ω^{3} Δ t}) \cos \bar{ω} Δ t} + \frac{2 ξ}{ω^{3} Δ t}, \\ b_{12} = & \exp (- ξ ω Δ t) {\frac{2 ξ^{2} - 1}{ω^{2} \bar{ω} Δ t} \sin \bar{ω} Δ t + \frac{2 ξ}{ω^{3} Δ t} \cos \bar{ω} Δ t} + \frac{1}{ω^{2}} - \frac{2 ξ}{ω^{3} Δ t}, \\ b_{21} = & - \exp (- ξ ω Δ t) {(\bar{ω} \cos \bar{ω} Δ t - ξ ω \sin \bar{ω} Δ t) (\frac{2 ξ^{2} - 1}{ω^{2} \bar{ω} Δ t} + \frac{ξ}{\bar{ω} ω}) \\ - (\bar{ω} \sin \bar{ω} Δ t + ξ ω \cos \bar{ω} Δ t) (\frac{1}{ω^{2}} + \frac{2 ξ}{ω^{3} Δ t})} - \frac{1}{ω^{2} Δ t}, \\ b_{22} = & - \exp (- ξ ω Δ t) {- (\bar{ω} \cos \bar{ω} Δ t - ξ ω \sin \bar{ω} Δ t) \frac{2 ξ^{2} - 1}{\bar{ω} ω^{2} Δ t} + (\bar{ω} \sin \bar{ω} Δ t + ξ ω \cos \bar{ω} Δ t) \frac{2 ξ}{ω^{3} Δ t}} \\ + \frac{1}{ω^{2} Δ t} . \end{matrix}

Damping equal to critical

In this case

\begin{matrix} a_{11} = & \exp (- ξ ω Δ t) (1 + ξ ω Δ t), \\ a_{12} = & \exp (- ξ ω Δ t) Δ t, \\ a_{21} = & - \exp (- ξ ω Δ t) {(ξ ω)}^{2} Δ t, \\ a_{22} = & \exp (- ξ ω Δ t) (1 - ξ ω Δ t), \end{matrix}

\begin{matrix} b_{11} = & - \exp (- ξ ω Δ t) {- (1 + ξ ω Δ t) (\frac{1}{ω^{2}} + \frac{2 ξ}{ω^{3} Δ t}) + \frac{1}{ω^{2}}}, \\ b_{12} = & \exp (- ξ ω Δ t) {(1 + ξ ω Δ t) (\frac{2 ξ}{ω^{3} Δ t} - \frac{1}{ω^{2}}) + \frac{1}{ω^{2}}} - \frac{2 ξ}{ω^{3} Δ t}, \\ b_{21} = & - \exp (- ξ ω Δ t) {{(ξ ω)}^{2} Δ t (\frac{1}{ω^{2}} + \frac{2 ξ}{ω^{3} Δ t}) + \frac{1 - ξ ω Δ t}{ω^{2} Δ t}} - \frac{1}{Δ t ω^{2}}, \\ b_{22} = & - \exp (- ξ ω Δ t) {- {(ξ ω)}^{2} Δ t \frac{2 ξ}{ω^{3} Δ t} - \frac{1 - ξ ω Δ t}{Δ t ω^{2}}} + \frac{1}{ω^{2} Δ t} . \end{matrix}

Damping higher than critical

In this case, with

\bar{ω} = ω \sqrt{ξ^{2} - 1},

we have

\begin{matrix} a_{11} = & \exp (- ξ ω Δ t) (ξ \frac{ω}{\bar{ω}} \sinh \bar{ω} Δ t + \cosh \bar{ω} Δ t), \\ a_{12} = & \exp (- ξ ω Δ t) \frac{1}{\bar{ω}} \sinh \bar{ω} Δ t, \\ a_{21} = & - \exp (- ξ ω Δ t) \frac{ω}{\sqrt{ξ^{2} - 1}} \sinh \bar{ω} Δ t, \\ a_{22} = & \exp (- ξ ω Δ t) (\cosh \bar{ω} Δ t - \frac{ξ ω}{\bar{ω}} \sinh \bar{ω} Δ t), \end{matrix}

\begin{matrix} b_{11} = & - \exp (- ξ ω Δ t) {(\frac{ξ}{ω \bar{ω}} + \frac{2 ξ^{2} - 1}{ω^{2} \bar{ω} Δ t}) \sinh \bar{ω} Δ t + (\frac{1}{ω^{2}} + \frac{2 ξ}{ω^{3} Δ t}) \cosh \bar{ω} Δ t} + \frac{2 ξ}{ω^{3} Δ t}, \\ b_{12} = & \exp (- ξ ω Δ t) {\frac{2 ξ^{2} - 1}{ω^{2} \bar{ω} Δ t} \sinh \bar{ω} Δ t + \frac{2 ξ}{ω^{3} Δ t} \cosh \bar{ω} Δ t} + \frac{1}{ω^{2}} - \frac{2 ξ}{ω^{3} Δ t}, \end{matrix}

\begin{matrix} b_{21} = & - \exp (- ξ ω Δ t) {(\bar{ω} \cosh \bar{ω} Δ t - ξ ω \sinh \bar{ω} Δ t) (\frac{2 ξ^{2} - 1}{ω^{2} \bar{ω} Δ t} + \frac{ξ}{\bar{ω} ω}) \\ - (- \bar{ω} \sinh \bar{ω} Δ t + ξ ω \cosh \bar{ω} Δ t) (\frac{1}{ω^{2}} + \frac{2 ξ}{ω^{3} Δ t})} - \frac{1}{ω^{2} Δ t}, \\ b_{22} = & - \exp (- ξ ω Δ t) {- (\bar{ω} \cosh \bar{ω} Δ t - ξ ω \sinh \bar{ω} Δ t) \frac{2 ξ^{2} - 1}{\bar{ω} ω^{2} Δ t} \\ + (- \bar{ω} \sinh \bar{ω} Δ t + ξ ω \cosh \bar{ω} Δ t) \frac{2 ξ}{ω^{3} Δ t}} + \frac{1}{ω^{2} Δ t} . \end{matrix}

Rigid body mode with damping

If there are rigid body modes in the finite element model, there will be one or several eigenvalues that are zero. The equation of motion (Equation 1) is reduced to

\begin{matrix} (5) & \ddot{q} + α \dot{q} = f_{t - Δ t} + \frac{Δ f}{Δ t} Δ t . \end{matrix}

Only Rayleigh damping can be specified for a rigid body mode, since the critical damping is zero. Furthermore, since it is a rigid body mode, only the mass damping factor, $α$ , appears (stiffness damping requires that there be straining of the body). For this case

\begin{matrix} a_{11} = & 1, \\ a_{12} = & \frac{1}{α} (1 - \exp (- α Δ t)), \\ a_{21} = & 0, \\ a_{22} = & \exp (- α Δ t), \end{matrix}

\begin{matrix} b_{11} = & - \frac{1}{α^{2}} (1 - \exp (- α Δ t)) (1 + \frac{1}{α Δ t}) + Δ t (1 + \frac{1}{α Δ t}) - \frac{Δ t}{2 α}, \\ b_{12} = & \frac{1}{α} (1 - \exp (- α Δ t)) \frac{1}{α^{2} Δ t} - \frac{1}{α^{2}} + \frac{Δ t}{2 α}, \\ b_{21} = & \frac{1}{α} (1 + \frac{1}{α Δ t}) (1 - \exp (- α Δ t)) - \frac{1}{α}, \\ b_{22} = & - \frac{1}{α^{2} Δ t} (1 - \exp (- α Δ t)) + \frac{1}{α} . \end{matrix}

Rigid body mode without damping

For the particular case of a rigid body mode without damping, the equation of motion (Equation 1) is reduced to

\begin{matrix} (6) & \ddot{q} = f_{t - Δ t} + \frac{Δ f}{Δ t} Δ t . \end{matrix}

For this case

\begin{matrix} a_{11} = & 1, \\ a_{12} = & Δ t, \\ a_{21} = & 0, \\ a_{22} = & 1, \\ b_{11} = & \frac{Δ t^{2}}{3}, \\ b_{12} = & \frac{Δ t^{2}}{6}, \\ b_{21} = & \frac{Δ t}{2}, \\ b_{22} = & \frac{Δ t}{2} . \end{matrix}

Structural damping in modal dynamic analysis

Structural damping is a commonly used damping model that represents damping as complex stiffness. This representation causes no difficulty for frequency domain analysis such as steady-state dynamics for which the solution is already complex. However, in the time domain, the solution must remain real-valued. To allow users to apply their structural damping model in the time domain, a method has been developed to convert structural damping to an equivalent viscous damping. This technique was designed so that, in the frequency domain, the viscous damping applied is identical to structural damping if the projected damping matrix is diagonal.

We start with the equation of a single degree-of-freedom oscillator,

m \ddot{x} + c \dot{x} + k (1 + i γ) x = f,

where m is the mass, c is the viscous damping factor, k is the stiffness, x is the response, f is the force, and $γ$ is the structural damping factor. If we normalize by the mass and assume harmonic input, we obtain the relation

(ω^{2} (1 + i γ) - Ω^{2} + i 2 ξ ω Ω) x_{0} = f_{0},

where $ω$ is the natural frequency, $Ω$ is the excitation frequency, $ξ$ is the coefficient of the harmonic response, and $f_{0}$ is the coefficient of the harmonic input. The relation between the viscous damping factor c and the critical damping factor $ξ$ is such that $c = 2 ξ ω m$ .

Based on the two previous equations, if we want the viscous damping factor to have the same effect as the structural damping factor, we must have

ξ = \frac{ω γ}{2 Ω} .

If we further assume that $ω = Ω$ , the relation simplifies to

ξ = \frac{γ}{2} .

Now we consider the finite element system of equations,

M^{N M} {\ddot{x}}^{M} + C^{N M} {\dot{x}}^{M} + (K^{N M} + i K_{s}^{N M}) x^{M} = f^{N},

where $C^{N M}$ is the finite element viscous damping matrix and $K_{s}^{N M}$ is the finite element structural damping matrix. If we project these equations onto the eigenmodes, the damping matrices become

\hat{C} = Φ^{N} C^{N M} Φ^{M} and {\hat{K}}_{s} = Φ^{N} K_{s}^{N M} Φ^{M},

which could be fully populated. If we take the case that $C^{N M}$ and $K_{s}^{N M}$ are diagonal, their diagonal elements can be written

{\hat{c}}_{α} = 2 ξ_{α} ω_{α} and {\hat{k}}_{s α} = γ_{α} ω_{α}^{2},

where $α$ is a particular mode. Maintaining the relation developed for the single degree-of-freedom system $ξ = \frac{γ}{2}$ for this diagonal case would require

{({\hat{c}}_{e q u i v a l e n t})}_{α} = \frac{{\hat{k}}_{s α}}{ω_{α}} .

For nondiagonal matrices the expression for equivalent viscous damping becomes

{\hat{C}}_{e q u i v a l e n t} = D {\hat{K}}_{s} D,

where D is a diagonal matrix whose entries are given by

d_{α} = \frac{1}{\sqrt{ω_{α}}} .

Response of nodal and element variables

The time integration is done in terms of the generalized coordinates, and the response of the physical variables is then immediately available by summation:

u = \sum_{α} ϕ_{α} q_{α},

ε = \sum_{α} ε_{α} q_{α},

σ = \sum_{α} σ_{α} q_{α},

R = \sum_{α} R_{α} q_{α},

where $ϕ_{α}$ are the modes, $ε_{α}$ are the modal strain amplitudes, $σ_{α}$ are the modal stress amplitudes, and $R_{α}$ are the modal reaction force amplitudes corresponding to each eigenvector $α$ .

Initial conditions

At the beginning of the step initial displacements and initial velocities must be converted to equivalent values of the generalized coordinates, which can be done exactly only if the number of eigenvectors equals the number of degrees of freedom. Since this is usually not the case, the initial values of the generalized coordinate displacement and velocity are calculated as

q_{α} = \frac{1}{m_{α}} ϕ_{α}^{M} M^{M N} x_{0}^{N},

where $m_{α}$ is the generalized mass for the eigenvector $α$ , $ϕ_{α}^{M}$ is the eigenvector, $M^{M N}$ is the mass matrix, and $x_{0}^{N}$ are the initial displacements.

Similarly, for the initial velocities

{\dot{q}}_{α} = \frac{1}{m_{α}} ϕ_{α}^{M} M^{M N} {\dot{x}}_{0}^{N} .

For the case where the initial conditions are given by a previous modal dynamic analysis, the generalized displacement, velocity, and acceleration are simply taken from the previous analysis.

Base motion definition

Many linear dynamic problems involve finding the response of a structure to a “base motion”: a time history of displacement, velocity, or acceleration given for the points where the displacements of the structure are prescribed. In all cases these base motions are converted into an acceleration history. If the displacement or velocity history has nonzero values at time $t =$ 0, corrections to the acceleration histories are made at times $t =$ 0 and $t = Δ t$ . If velocities are given, the acceleration at $t =$ 0 is

{\ddot{z}}_{0} = \frac{{\dot{z}}_{Δ t}^{*} - {\dot{z}}_{0}^{*}}{Δ t} .

If displacements are given, the accelerations at $t =$ 0 and $t = Δ t$ are

{\ddot{z}}_{0} = \frac{3 z_{Δ t}^{*}}{Δ t^{2}} and {\ddot{z}}_{Δ t} = \frac{z_{2 Δ t}^{*} - \frac{5}{2} z_{Δ t}^{*}}{Δ t^{2}} .

In the above expressions and in the expressions to follow, a superscript * indicates user-defined amplitude data.

If velocities are given for tabular or equally spaced amplitude curve definitions, the acceleration is defined by the central difference rule

{\ddot{z}}_{t + Δ t} = \frac{{\dot{z}}_{t + 2 Δ t}^{*} - {\dot{z}}_{t}^{*}}{2 Δ t} .

If displacements are given, the acceleration is defined by the central difference rule

{\ddot{z}}_{t + Δ t} = \frac{z_{t + 2 Δ t}^{*} - 2 z_{t + Δ t}^{*} + z_{t}^{*}}{Δ t^{2}} .

The response is calculated relative to the base. If total values of nodal variables are required, the motion at the base is added to the relative values:

\begin{matrix} {\ddot{u}}_{abs} & = \ddot{z} + \ddot{x}, \\ {\dot{u}}_{abs} & = \dot{z} + \dot{x}, \\ u_{abs} & = z + x, \end{matrix}

where

\begin{matrix} {\dot{z}}_{t + Δ t} & = {\dot{z}}_{t} + \frac{Δ t}{2} ({\ddot{z}}_{t} + {\ddot{z}}_{t + Δ t}), \\ z_{t + Δ t} & = z_{t} + {\dot{z}}_{t} Δ t + \frac{Δ t^{2}}{3} ({\ddot{z}}_{t} + \frac{1}{2} {\ddot{z}}_{t + Δ t}) . \end{matrix}