Modal analysis: Theory

A modal analysis calculates the undamped natural modes of a system, characterised by their modal frequency and mode shape. These modes are numbered, from 1, in order of increasing frequency. Modal analysis is a standard technique, well-documented in the literature: we give here a brief description of the underlying theory.

Outline theory

The theory is illustrated by a single degree of freedom system consisting of a mass $m$ attached to a linear spring of stiffness $k$. The undamped equation of motion is \begin{equation} \label{eqOfMotion} m x''(t) = -k x(t) \end{equation} where $x(t)$ is the offset at time $t$ from the mean position and $x''(t)$ the acceleration. This analysis neglects any damping, which is what gives rise to the term undamped modes.

The solution of this equation is known to be simple harmonic, i.e. of the form \begin{equation} x(t) = a\sin(\omega t) \end{equation} where $a$ and $\omega$ are unknowns to be determined. Differentiating $x(t)$ twice gives \begin{equation} x''(t) = -\omega^2 a\sin(\omega t) = -\omega^2 x(t) \end{equation} and substituting into the equation of motion (\ref{eqOfMotion}) we obtain \begin{equation} \label{eigenProblem} -m \omega^2 x(t) = -k x(t) \end{equation} which can be rearranged to give the familiar result \begin{equation} \omega = \sqrt{\frac{k}{m}} \end{equation} This is the angular natural frequency of the oscillation: the natural period $T$ is the equally familiar \begin{equation} T = 2\pi\sqrt{\frac{m}{k}} \end{equation}

Discretisation of continuous systems

For the simple harmonic oscillator above, there is just a single undamped natural mode, corresponding to the single degree of freedom $x$. For a continuous riser, say, there are an infinite number of degrees of freedom, and hence an infinite number of undamped natural modes – but computers must work with discrete models with finite numbers of degrees of freedom.

Consider, then, a discretised line in OrcaFlex with $n$ degrees of freedom. The above equations still apply, but they now have to be interpreted as matrix equations: while $\omega$ and $T$ remain scalars, we now have vectors $\vec{a}$, $\vec{x}$ and $\vec{x''}$, each with $n$ elements (one for each degree of freedom), and $n{\times}n$ matrices $\mat{M}$ and $\mat{K}$. Equation (\ref{eigenProblem}) becomes \begin{equation} \label{discreteEigenProblem} \omega^2 \mat{M} \vec{x} = \mat{K} \vec{x} \end{equation} This is an eigen-problem with $n$ solutions, the $i$th solution having frequency $\omega = \omega_i$ and associated amplitude $\vec{x}=\vec{a}_i$, say: $\omega_i$ is a scalar, $\vec{a}_i$ an $n$-vector. This $i$th solution is called the $i$th natural mode, and represents a standing wave oscillation of the line in which the $n$ degrees of freedom oscillate at the single angular frequency $\omega_i$, but with different amplitudes given by the components of $\vec{a}_i$. For the $i$th mode, the components of $\vec{a}_i$ represent the maximum offset from the mean position of each degree of freedom.

But what is the effect of calculating the natural modes of the discretised model, rather than those of the real continuous system? Well, the discretised modes are, broadly speaking, close to the continuous ones, and for a given mode number they become closer as the discretisation becomes finer – as more and more elements are used to model the system. For any given level of discretisation, the accuracy is found to be better for the lower modes and progressively worse for higher modes. The highest numbered modes represent oscillations whose wavelengths are of the same order as the segment length: it is unlikely that these will be physically realistic.

Note that the eigenvectors $\vec{a}_i$ are only unique up to a scalar multiple: it is clear from equation (\ref{discreteEigenProblem}) that, if $\vec{x}=\vec{a}_i$ is a solution, then so is $\lambda \vec{a}_i$ for any real value of $\lambda$. It is the relative differences between the components of $\vec{a}_i$ which constitutes the $i$th mode shape.

Seabed friction

The theory above requires that the mass and stiffness matrices are symmetric which is not always the case in an OrcaFlex model. The most important example of this is the case of friction. Friction is a non-conservative effect, leading directly to non-symmetric terms in the stiffness matrix. Clearly this presents a problem.

This non-conservatism occurs when a node is slipping, that is when the deflection from its friction target position exceeds $D_\textrm{crit}$. When performing a modal analysis OrcaFlex assumes that nodes on the seabed are restrained by a linear stiffness calculated as the product of the seabed's shear stiffness $k_\mathrm{s}$ and the node's contact area $a$. This stiffness term corresponds to that of a linear spring acting in the plane of the seabed, connecting the node and its target position with a stiffness of $k_\mathrm{s}a$.

This has the effect of restraining movement of the nodes on the seabed in the plane of the seabed, which is desirable for a modal analysis of a system with seabed contact. This modification to the seabed friction model, for the purposes of modal analysis only, results in a symmetric, conservative system which is amenable to modal analysis.

Stiffness terms due to fluid loading

As we have seen, undamped modal analysis excludes damping terms (those dependent on $\vec{x'}$). So, for example, the effects of drag loading are not included.

Fluid loads do, however, also contribute stiffness terms, despite their dependence on velocity, because perturbations of position and orientation can result in a change of direction of the fluid load vector. These are, strictly speaking, stiffnesses rather than damping, because they arise from translational and rotational displacement; these fluid load terms are not, however, included by OrcaFlex in the modal analysis.

Nonlinear line type stiffness

A fundamental assumption of modal analysis, clear from the equations above, is linearity of the system under analysis. Axial and bend stiffness at their simplest, are indeed linear in OrcaFlex, and present no problems to modal analysis. For nonlinear elastic stiffness we use the value of the local tangent stiffness. For small oscillations about the static configuration, this should result in modal analysis being adequately accurate.

For hysteretic stiffness, the local tangent stiffness is in general ill-defined: it is multi-valued and depends on whether the perturbation increases or decreases strain or curvature.

For hysteretic bending, the stiffness used for modal analysis depends on the statics model. If depressurised behaviour is selected, modal analysis uses a linear bend stiffness determined by the slope of the final two rows of the bend stiffness table. For pressurised behaviour modal analysis takes an average of the possible stiffness values.

Axial or bend stiffness may also be externally calculated. The external function interface provides no mechanism for specifying the local tangent bend stiffness; instead, OrcaFlex uses the nominal stiffness value which the external function is required to provide.

Clearly, the modal analysis is likely to be less accurate for hysteretic and externally calculated stiffness cases. It is often found that tension, rather than bend stiffness, is the dominant contribution to lateral stiffness; in the case of hysteretic and externally calculated bend stiffness this inaccuracy is therefore often not significant.

Vessel added mass

In a whole system modal analysis which includes vessel degrees of freedom, the vessel added mass is usually a significant factor in the analysis. However, if the vessel added mass and damping data are frequency dependent, there is no easy way to account for this dependence on frequency; added mass is therefore neglected in modal analysis if the data are frequency dependent.

To account for added mass in the modal analysis, you can instead use a single constant added mass matrix, which is included in the system-wide mass matrix. Because you can only specify a single matrix, you must first assess what modes are of interest and choose appropriate added mass values based on their associated frequency. If the added mass varies significantly with frequency over the range of modal frequencies under consideration, then you may need to perform multiple modal analyses with different added mass matrices.

Line contact with containment

In the case of line contact with containment, the modal analysis excludes most of the fluid inertia contributions between contained line nodes and the contents fluid of the containing line. In particular,

These effects are excluded because the associated inertia terms introduce inertial couplings between the inner and outer lines. These couplings may significantly increase the computational cost of the modal analysis and for reasons of expediency we choose to neglect them.

Mode offset distribution

For VIV analysis of a line, OrcaFlex classifies modes according to the distribution of mode offset between inline, transverse, axial and rotational components.

Recall that, for a given mode $m$, we have a mode offset vector $\vec{a}_m$ whose components are the degrees of freedom of the line. Then, for a given node $n$ on the line, we can obtain a node offset vector, $\vec{a}_{mn}$, the components of which are the appropriate degrees of freedom components of $\vec{a}_m$. Recall too, that $\vec{a}_m$ is unique only up to a scalar multiple so, in order to do this classification, we must choose an appropriate scalar value for each mode. We do this by normalising each mode offset $\vec{a}_m$.

Mode offset normalisation

To normalise a mode offset vector $\vec{a}_m$, we first partition each of its node offset vectors, $\vec{a}_{mn}$, into translational ($\vec{a}_{mn}^\mathrm{T}$, in length units) and rotational ($\vec{a}_{mn}^\mathrm{R}$, in radians) parts. We then determine the largest magnitude of these vectors over all nodes, $\max\limits_n \{ \vec{a}_{mn}^\mathrm{T}, \vec{a}_{mn}^\mathrm{R} \}$, and scale the entire mode offset vector $\vec{a}_m$ by this value such that the magnitude of the largest node offset vector, $\vec{a}_{mn}^\mathrm{T}$ or $\vec{a}_{mn}^\mathrm{R}$, for that mode, is $1$.

Distribution

Writing the individual components of $\vec{a}_{mn}^\mathrm{T}$ after this normalising process, in the inline, transverse, and axial directions, as $t_{mn\mathrm{i}}$, $t_{mn\mathrm{t}}$ and $t_{mn\mathrm{a}}$ respectively, OrcaFlex then calculates the inline, transverse and axial components of the overall translation in the mode shape $m$ as \begin{equation} \begin{aligned} M_{m\mathrm{i}} &= \left( \sum_n l_n\, t_{mn\mathrm{i}}^2 \right)^{\small{1/2}} \\ M_{m\mathrm{t}} &= \left( \sum_n l_n\, t_{mn\mathrm{t}}^2 \right)^{\small{1/2}} \\ M_{m\mathrm{a}} &= \left( \sum_n l_n\, t_{mn\mathrm{a}}^2 \right)^{\small{1/2}} \end{aligned} \label{MTranslation} \end{equation} where the summations are over all nodes $n$ in the line and $l_n$ is the length of line represented by that node.

This root sum of squares formulation is the multi-dimensional equivalent of the usual 3D formula, $|v| = (v\urm{i}^2 + v\urm{t}^2 + v\urm{a}^2)^{\small{1/2}}$. The scaling by $l_n$ is introduced for independence of the level of discretisation.

If the line includes torsion, then a further component $M_{m\mathrm{r}}$ is required to represent the overall rotational content of the mode shape. Analogous to (\ref{MTranslation}) for translation, components $r_{mn\mathrm{i}}$, $r_{mn\mathrm{t}}$ and $r_{mn\mathrm{a}}$ are obtained from the rotational partition $\vec{a}_{mn}^\mathrm{R}$ of the mode offset vector $\vec{a}_m$, after normalisation, for each node $n$; we multiply these by the node's radius of gyration about the line axis, $g_n$, to give them dimensions of length and allow an appropriate comparison between the magnitudes of translation and rotation. Again we sum over all the line nodes, but this time we also sum the three components because all we require is a single measure of how much rotational content is present in the mode offset vector. Hence \begin{equation} M_{m\mathrm{r}} = \left( \sum_n l_n\, (g_n r_{mn\mathrm{i}})^2 \right)^{\small{1/2}} + \left( \sum_n l_n\, (g_n r_{mn\mathrm{t}})^2 \right)^{\small{1/2}} + \left( \sum_n l_n\, (g_n r_{mn\mathrm{a}})^2 \right)^{\small{1/2}} \end{equation} If the line does not include torsion then $M_{m\mathrm{r}} = 0$.

Finally, OrcaFlex expresses the components $M_{m\mathrm{i}}$, $M_{m\mathrm{t}}$, $M_{m\mathrm{a}}$, and $M_{m\mathrm{r}}$ of the mode shape $m$ as percentages of their sum $M_m = M_{m\mathrm{i}}{+}M_{m\mathrm{t}}{+}M_{m\mathrm{a}}{+}M_{m\mathrm{r}}$, and reports these percentages as offset distribution values for each mode $m$.

Modal mass and stiffness

The modal mass, associated with mode $m$, is calculated as \begin{equation} m_m = \vec{a}_m^\mathrm{T} \mat{M} \vec{a}_m \end{equation} where $\vec{a}_m$ is the normalised mode shape vector, $\vec{a}_m^\mathrm{T}$ is its transpose (row vector) and $\mat{M}$ is the system's mass matrix.

The modal stiffness is calculated as \begin{equation} k_m = \omega_m^2 m_m \end{equation} where $\omega_m$ is the angular frequency of the mode.

This is equivalent to \begin{equation} k_m = \vec{a}_m^\mathrm{T} \mat{K} \vec{a}_m \end{equation} where $\mat{K}$ is the system's stiffness matrix.

Eigen-solvers

A further technical detail is the choice of solution method. OrcaFlex uses two different eigen-solvers to perform modal analysis. The choice of which to use is based on the number of modes extracted, $m$, and the number of degrees of freedom, $n$.

If $m\leq n/3$ and $m\leq 1000$ then an iterative Lanczos algorithm is used, otherwise we use a direct method based on tridiagonal matrix diagonalisation. For large problems the iterative Lanczos algorithm is much faster and requires much less memory, and so we use it in preference to the direct method wherever possible.

One final subtlety concerns the precise definition of $m$, the number of modes extracted. The Lanczos algorithm works by finding the largest (or smallest) eigenvalue first, then the next largest (or smallest) and so on. Consequently, if you ask for (say) modes 5 to 10, the solver has to first find modes 1 to 4, and so $m$ is effectively 10 in this case.