The harmonic analysis shows the dynamic behavior of a system under a defined load case. It can be used to assess how a structure vibrates under harmonic loads. Given certain damping mechanisms, the dissipation of vibrational energy can also be taken into account in the harmonic analysis.

Within the model of a harmonic excitation it is assumed that the load has been acting on the system since infinite time. This requirement can be relaxed by requiring the time since the start of the excitation to be much larger than the eigenperiod of the system T , so that \frac{t}{T} \gg 1. Then it can be assumed, that the system is in the steady state vibration. This means, that the transient response, which is present after the beginning of the excitation, has sufficiently decayed. Examples for harmonic excitations are:

  • Machines with rotating components
  • Paper factories (large rolls)
  • Rotating machinery
  • Wind (to certain extent)

Eigenmode of a harmonically excited plate

Eigenmode of a harmonically excited plate

Theory

In general, the equation of motion for a damped linear multi degree of freedom system with n degrees of freedom is given by:

(1) \mathbf{M} \ddot{\mathbf{u}} (t) + \mathbf{C} \dot{\mathbf{u}} (t) + \mathbf{K} \mathbf{u} (t) = \mathbf{f} (t)

Example:

\begin{align}\mathbf{K}= \begin{bmatrix} 1.5 & -0.5\\ -0.5 & 0.5 \\ \end{bmatrix}, \quad \mathbf{M}= \begin{bmatrix} 1 &0\\ 0&0.5 \end{bmatrix},\quad \mathbf{f}=\begin{bmatrix} 1\\0 \end{bmatrix} \end{align}

In order to solve , we apply the Fourier transformation operator to both sides of (1):

\hat{\mathbf{x}} (\omega) = \int\limits_{-\infty}^\infty \mathbf{x} (t) \mathrm{e}^{- \mathrm{i} \omega t} \text{d}t

We now make use of the property of the Fourier transformation that the derivative in time \frac{\text{d}}{\text{d}t} is replaced by a multiplication with \mathrm{i} \omega in the frequency domain.

Thus, the  reads as follows in the transformed space:

(2) - \omega^2 \mathbf{M} \hat{\mathbf{u}} (\omega) + \mathrm{i} \omega \mathbf{C} \hat{\mathbf{u}} (\omega) + \mathbf{K} \hat{\mathbf{u}} (\omega) = \hat{\mathbf{f}} (\omega)

And we finally obtain obtain the equation of motion in frequency domain

\left( - \omega^2 \mathbf{M} + \mathrm{i} \omega \mathbf{C} + \mathbf{K} \right) \hat{\mathbf{u}} (\omega) = \hat{\mathbf{f}} (\omega)

We can define the dynamic stiffness matrix as

\mathbf{D} (\omega) = - \omega^2 \mathbf{M} + \mathrm{i} \omega \mathbf{C} + \mathbf{K}

The solution to the (2) is then given by:

\hat{\mathbf{u}} (\omega) = \mathbf{D}^{-1} (\omega) \hat{\mathbf{f}} (\omega) = \left( - \omega^2 \mathbf{M} + \mathrm{i} \omega \mathbf{C} + \mathbf{K} \right)^{-1} \hat{\mathbf{f}} (\omega)

The solution in time domain can then be found by applying the inverse Fourier transformation operator

\mathbf{x} (t) = \int\limits_{-\infty}^\infty \hat{\mathbf{x}} (\omega) \mathrm{e}^{\mathrm{i} \omega t} \text{d}\omega

We now consider a harmonic load in complex notation

\mathbf{f} (t) = \widehat{\mathbf{f}}_{+} \mathrm{e}^{\mathrm{i} \Omega t} + \widehat{\mathbf{f}}_{-} \mathrm{e}^{- \mathrm{i} \Omega t}

where \Omega is the excitation frequency. The complex notation can be used to represent the harmonic character of the load:

We can show that the load has the following form in frequency domain

\hat{\mathbf{f}} (\omega) = \widehat{\mathbf{f}}_{+} \delta(\omega - \Omega) + \widehat{\mathbf{f}}_{-} \delta(\omega + \Omega),

where \delta(x) is the Dirac-delta function. Note that the Fourier transform of the complex exponential function can only be given in terms of a distribution, here the Dirac-delta function, as it violates the condition of absolute integrability.

For this load, the response reads:

\hat{\mathbf{u}} (\omega) = \mathbf{D}^{-1} (\omega) \left( \widehat{\mathbf{f}}_{+} \delta(\omega - \Omega) + \widehat{\mathbf{f}}_{-} \delta(\omega + \Omega) \right) = \underbrace{\mathbf{D}^{-1} (\omega) \widehat{\mathbf{f}}_{+}}_{\widehat{\mathbf{u}}_{+}} \delta(\omega - \Omega) + \underbrace{\mathbf{D}^{-1} (\omega) \widehat{\mathbf{f}}_{-}}_{\widehat{\mathbf{u}}_{-}} \delta(\omega + \Omega)

Thus, \widehat{\mathbf{u}}_{+} will only contribute at \omega = \Omega, and \widehat{\mathbf{u}}_{-} at \omega = - \Omega, i.e., at the negative excitation frequency, due to the multiplication with the Dirac-delta function.

Therefore the inverse of the dynamic stiffness matrix only needs to be evaluated at \omega = \Omega and \omega = - \Omega.

Furthermore, as \widehat{\mathbf{f}}_{+} and \widehat{\mathbf{f}}_{-} as well as \mathbf{D}^{-1} (\Omega) and \mathbf{D}^{-1} (- \Omega) are complex conjugate to each other, also \widehat{\mathbf{u}}_{+} and \widehat{\mathbf{u}}_{-} will be complex conjugate to each other.

For this reason, only one of the two parts of the solution needs to be computed. Usually the solution is obtained for \widehat{\mathbf{u}}_{+}.

Finally, the response in the frequency domain reads:

\hat{\mathbf{u}} (\omega) = \widehat{\mathbf{u}}_{+} \delta(\omega - \Omega) + \widehat{\mathbf{u}}_{-} \delta(\omega + \Omega)

The solution in the time domain can then be found by the inverse Fourier transformation and is given by

\mathbf{u} (t) = \widehat{\mathbf{u}}_{+} \mathrm{e}^{\mathrm{i} \Omega t} + \widehat{\mathbf{u}}_{-} \mathrm{e}^{- \mathrm{i} \Omega t}

Hence, the response of the system to harmonic excitation will also be harmonic with the same frequency \Omega .

SDOF system

We consider a damped single degree of freedom system with stiffness k, viscous damping constant c and mass m.

The natural frequency is given by

\omega_n=\sqrt{\frac{k}{m}}

and its dynamic stiffness as

d (\omega) = - \omega^2 m + \mathrm{i} \omega c + k = k \left( 1 - r^2 + \mathrm{i} 2 \zeta r \right)

where r = \frac{\omega}{\omega_n} is the frequency ratio, and \zeta = \frac{c}{2 \sqrt{km}} is the damping ratio.

The system is excited by a harmonic force with the complex pointer

\widehat{f}_{+} = 0.5 + 0.25 \mathrm{i}

and an excitation frequency \Omega = 2 \omega_n that is twice the eigenfrequency of the system, i.e., r = 2.

To calculate the response of the system, we need to evaluate the inverse of the dynamic stiffness for r = 2. We assume that the system is lightly damped with a damping ratio of \zeta = 5\%:

d^{-1} (\omega = 2 \omega_n) = \frac{1}{k \left(1 - 2^2 + \mathrm{i} \cdot 2 \cdot 0.05 \cdot 2\right)} = \frac{1}{k \left( -3 + \mathrm{i} \cdot 0.2 \right)}

In order to obtain the displacement, we can now multiply the foregoing result with the complex pointer \widehat{f}_{+}. We do not write the Dirac-delta function explicitly, nevertheless we inherently consider its sampling property by only evaluating the response for \omega = \Omega (The load in the frequency domain \hat{\mathbf{f}} (\omega) is zero for all \omega \neq \Omega).

\widehat{u}_+ = d^{-1} (\omega = \Omega) \widehat{f}_{+} = \frac{0.5 + 0.25 \mathrm{i}}{k \left( -3 + \mathrm{i} \cdot 0.2 \right)} = \frac{1}{k} (-0.16 - 0.094 \mathrm{i})

The corresponding pointer \widehat{u}_- can be found by complex conjugation and thus reads:

\widehat{u}_- = \widehat{u}_+^{\ast} = \frac{1}{k} (-0.16 - 0.094 \mathrm{i})^{\ast} = \frac{1}{k} (-0.16 + 0.094 \mathrm{i})

The full solution in the frequency domain then reads:

\hat u (\omega) = \widehat{u}_+ \delta (\omega - \Omega) +\widehat{u}_- \delta (\omega + \Omega) = \frac{1}{k} \left[ (-0.16 - 0.094 \mathrm{i}) \delta (\omega - \Omega) + (-0.16 + 0.094 \mathrm{i}) \delta (\omega + \Omega) \right]

The solution can also be visualized as pointers on the complex plane: