|
Dynamic analysis: Frequency domain solution |
|
|
Dynamic analysis: Frequency domain solution |
The frequency domain solver is aimed at solving for the dynamic response of the system at either wave frequency or low frequency, as determined by the user-specified solution frequencies.
The response at wave frequency is defined to be that of the system subject to first order dynamic loading associated with the wave elevation stochastic process. As in the time domain, the wave elevation process is described by user-specified wave elevation spectra. Unlike the time domain, the wave elevation spectra are not used to synthesise a particular realisation of the sea state. Instead, frequency domain analysis solves for the wave frequency response process of the system by applying a series of linear mappings to the wave elevation process.
The response at low frequency is defined to be that of the system subject to both second order wave drift dynamic loading, associated with the wave elevation stochastic process, and wind dynamic loading, associated with the wind velocity stochastic process. The solution at wind component frequencies is very similar to the wave frequency solution; the response process at wind component frequencies is found by applying a series of linear mappings to the wind velocity process. For wave drift frequencies, the calculation is more complex due to the second order nature of the wave drift load, and the response at these frequencies cannot simply be determined by applying a series of linear mappings to the wave elevation process. Instead, wave drift load processes must be calculated from the wave spectra and the vessel QTF data (vessels are the only OrcaFlex objects which respond to wave drift frequencies), and a series of linear mappings applied to find the response process.
The input and output of the frequency domain solver are stochastic; the result process describes the statistics of the result over the whole ensemble of realisations, i.e. over all possible synthesised wave and wind realisations.
In more detail then, the frequency domain solution is conducted by applying a series of linear mappings to some underlying stochastic environmental or loading process. In the case of the wave frequency solution, the underlying process is the wave elevation and the solution procedure is as follows.
A discretised scalar-valued wave elevation process, $\eta(f_n)$, is calculated from the wave elevation spectrum $S_{\eta\eta}(f_n)$, using the same discretisation methods as are available in the time domain, such that the wave elevation process associated with the $n$th wave component, over the frequency band $\dfn$ centred at $f_n$, is given by \begin{equation} \eta(f_n) = \left( S_{\eta\eta}(f_n)\,\dfn \right)^{1/2} \end{equation} The response process is calculated in OrcaFlex by deriving, in stages, a set of linear transfer functions that map the wave elevation process through to the response process.
First we derive a complex vector-valued linear transfer function, $\vec{\lambda}(f_n)$, mapping $\eta(f_n)$ to the load process on the free degrees of freedom of the system, $\vec{l}(f_n)$, such that \begin{equation} \vec{l}(f_n) = \vec{\lambda}(f_n)\,\eta(f_n) \end{equation} To calculate $\vec{\lambda}(f_n)$, we apply linear Airy wave theory to map the wave elevation process to the wave velocity and wave acceleration processes at the model object positions, in combination with the model objects' properties which result in a load on the object due to the presence of these wave elevations / velocities / accelerations. Such properties are, typically, added mass, linear damping, viscous drag, RAOs etc.
Next we derive a complex matrix-valued linear transfer function, $\mat{X}(f_n)$ that maps the load process, $\vec{l}(f_n)$, to the response process, $\vec{x}(f_n)$ such that \begin{equation} \vec{x}(f_n) = \mat{X}(f_n)\,\vec{l}(f_n) \end{equation} where \begin{equation} \mat{X}(f_n) = \left( -\left(2\pi f_n\right)^2\mat{M} + \mathrm{i}2\pi f_n \mat{C} + \mat{K} \right)^{-1} \end{equation} Here
$\mat{M}$ is the system inertia matrix
$\mat{C}$ is the system damping matrix
$\mat{K}$ is the system stiffness matrix.
Finally the response process $\vec{x}$ is determined, for each wave component $f_n$, by applying these transfer functions \begin{equation} \vec{x}(f_n) = \mat{X}(f_n)\,\vec{\lambda}(f_n)\,\eta(f_n) \end{equation}
Unfortunately, in general neither the wave elevation-loading relation $\vec{\lambda}$ nor the loading-response relation X is linear: they must be linearised in order to derive linear transfer functions. There are many potential sources of nonlinearity in an OrcaFlex model; barring the exceptions discussed below, their linearised form is obtained from the Jacobian matrix of the relationship, evaluated with the model in the static state. Essentially the Jacobian is the first order derivative, or gradient, of the function at the point of evaluation. For example, the system's stiffness is linearised by using the tangent stiffness matrix for the system calculated after statics has been run. The exceptions are friction and drag.
Linearisation of friction divides into two distinct cases: friction between the seabed and nodes, and all other forms of friction. The reason for the distinction is that only for seabed friction on nodes is it possible for the deflection from the friction target position to exceed the critical deflection $D_\mathrm{crit}$.
Seabed friction for nodes is, if possible, linearised using the tangent stiffness as above. However, if after the static analysis a node is slipping, i.e. the deflection from its friction target position exceeds $D_\mathrm{crit}$, then the tangent stiffness is not well-defined. This is because in dynamics the slipped node is always on the transition between the pre-slip part of the friction curve (which has non-zero linear stiffness) and the post-slip part (which has zero stiffness). When performing frequency domain analysis, OrcaFlex linearises assuming that nodes on the seabed are always on the pre-slip part of the curve, so essentially that they are restrained by a purely linear stiffness effect determined by the seabed's shear stiffness, $k_s$ and the node's contact area, $a$. This corresponds to a linear spring of stiffness $k_s a$, acting in the plane of the seabed, connecting the node and its target position.
For other friction relationships, i.e. between objects other than nodes and the seabed, friction is not active in statics. The friction target positions are defined by the objects' static positions and so the tangent stiffness is well defined and is given by the linear pre-slip part of the friction curve.
Hydrodynamic drag models are not linear in relative velocity $\vec{v}$, and therefore must be linearised before the system's response can be solved for. Because the drag loads often contribute significantly to the loading, and the damping, in an OrcaFlex model, and because the Jacobian will generally yield a poor linearisation in this case, this linearisation is handled in a special way known as equivalent linearisation. Equivalent linearisation is employed to linearise the line drag, lumped buoy drag, spar buoy and towed fish drag, 3D buoy drag, Morison elements, vessel other damping, and vessel yaw rate drag models.
In equivalent linearisation, rather than simply neglecting the nonlinear terms, a linear transfer function is instead calculated such that some measure of the difference between the linear and nonlinear solutions is minimised. The measure chosen here is the mean square error, and the method employed is the so called minimum mean square error (MMSE) method as detailed by Langley and Leira. Essentially the aim is to replace the nonlinear model for the drag force $\vec{f}$ \begin{equation} \vec{f} \propto \vec{v} \vert\vec{v}\vert \end{equation} with a linear model of the form \begin{equation} \vec{f} \propto \mat{L}_1 \vec{v} \end{equation} by calculating a linear matrix-valued transfer function $\mat{L}_1$ which minimises the mean square of the error $\vec{\epsilon}$ between the linear and nonlinear models \begin{equation} \PD{\E\left[\vec{\epsilon}^T\vec{\epsilon}\right]}{\mat{L}_1} = 0 \end{equation} where \begin{equation} \vec{\epsilon} = \mat{L}_1 \vec{v} - \vec{v} \vert\vec{v}\vert \end{equation} and $\E$ is the expectation operator. Solving this equation for $\mat{L}_1$ yields a linear transfer function \begin{equation} \label{L1} \mat{L}_1 = \mat{\Sigma}_\mathrm{v}^{-1} \E\left[\vert\vec{v}\vert \vec{v} \vec{v}^T \right] \end{equation} where $\mat{\Sigma}_\mathrm{v}$ is the relative velocity covariance matrix. The expectations in equation (\ref{L1}) are of functions of the relative velocity: they can be calculated, assuming the components of the relative velocity are distributed with the multivariate normal distribution, by using Gauss-Hermite quadrature to evaluate the integral \begin{equation} \E\left[f(\vec{v})\right] = \int_{-\infty}^\infty f(\vec{v})\, \P(\mat{\Sigma}_\mathrm{v})\, \ud\vec{v} \end{equation} where $\P$ is the multivariate normal distribution. This results in a linear transfer function that is a function of the covariance matrix of the relative velocity vector; unfortunately, this is in turn a function of the response velocity, and so an iterative approach becomes necessary. In the case of pure wave frequency or low frequency solution frequencies, the procedure is as follows:
This process continues until the convergence criteria are met or the maximum number of iterations is reached.
In the special case of combined linearisation, an initial wave frequency solve is conducted (using the procedure outlined above), followed by a low frequency solve. In this latter low frequency solve, the procedure is as follows:
Again, this second process continues until the convergence criteria are met or the maximum number of iterations is reached. This combined linearisation method is used to allow for the coupling that occurs across the response frequencies introduced by the quadratic nature of the drag models. Because it is broken into a distinct initial wave frequency solve and then a subsequent low frequency solve, which is dependent on the initial wave frequency solve, there is an implicit assumption that the wave frequency motion can increases the drag at low frequency but not vice versa, i.e. it is assumed the low frequency motion does not significantly influence the wave frequencies response; generally the response velocities associated with the wave frequency motion can be expected to significantly larger than those associated with the low frequency motion.
The MMSE method has the advantage of providing a frame-invariant linearisation scheme: the linearisation is not dependent on the arbitrary choice of the object's local axes directions.
In the presence of a steady current, the mean drag load is not necessarily equal to the static drag load. At present, this effect is not included by the frequency domain solver, i.e. the mean and static drag loads are assumed to be equal.
The frequency domain solution is deemed to have converged if the relative error in the linear drag transfer function's matrix norm, measured between the current iteration and the previous iteration, is less than the user-specified tolerance $\textit{tol}$ \begin{equation} \label{Convergence} \abs\big( \lVert\mat{L}_1\rVert - \lVert\mat{L}_1^\textrm{prev}\rVert \big) \lt \max\big\{ \lVert\mat{L}_1\rVert,\ \lVert\mat{L}_1^\textrm{prev}\rVert \big\} . \textit{tol} \end{equation} All model objects which support drag data have at least one linear drag transfer function associated with them. Inequality (\ref{Convergence}) must be satisfied for all linear drag transfer functions, across all model objects, for convergence to be achieved.
Once a converged response process has been found, the simulation is complete and frequency domain results are available through the results form.
To produce any results in the frequency domain, OrcaFlex first calculates the result process, $r(f_n)$, for the particular result of interest. To do this, a complex linear transfer function, $\vec{\rho}(f_n)$, which maps the response process $\vec{x}(f_n)$ to the result process is calculated and applied \begin{equation} r(f_n) = \vec{\rho}(f_n)^T\,\vec{x}(f_n) \end{equation} The results offered are limited to those for which a good linear mapping exists. Some nonlinear results, such as vector magnitudes, which often do not have a good linear transform and which in the time domain would result in non-Gaussian processes, are not available in the frequency domain.
As described above, the results process is linearly mapped from some underlying process, such as the wave elevation in the case of a wave frequency solve. This means that, when viewed about its static value, the result process can also be assumed to be a zero-mean stationary Gaussian process. Consequently, its statistics are well known.
The spectral density of the results process, at the frequency associated with the $n$th wave component, $f_n$, is calculated as \begin{equation} S_{rr}(f_n) = \frac{r(f_n)\, r(f_n)^*}{\dfn} \end{equation} where $^*$ denotes complex conjugate.
The $i$th moment of the continuous result spectral density is defined as \begin{equation} m_i = \int_0^\infty f^i\, S_{rr}(f)\, \ud f \end{equation} However, in OrcaFlex the moments are calculated from the discretised result process, so \begin{equation} m_i = \sum_n f_n^i\, r(f_n)\, r(f_n)^* \end{equation} The standard deviation of the result process is given by the square root of the zeroth moment of the result spectral density \begin{equation} \sigma = \sqrt{m_0} \end{equation} The mean period of zero up-crossings, $T_z$, is calculated as \begin{equation} T_z = \sqrt{\frac{m_0}{m_2}} \end{equation} The mean period of crests, $T_c$, is calculated as \begin{equation} T_c = \sqrt{\frac{m_2}{m_4}} \end{equation} The spectral bandwidth $\epsilon$ is calculated as \begin{equation} \epsilon = \sqrt{1 - \frac{m_2^2}{m_0m_4}} = \sqrt{1 - \left(\frac{T_c}{T_z}\right)^2} \end{equation} The most probable maximum (MPM) of the dynamic part of the result process, occurring in duration $T$, is \begin{equation} \textrm{MPM} = \sigma \sqrt{2 \ln \left(\frac{T}{T_z}\right)} \end{equation} The MPM is calculated under the assumption of the Rayleigh distribution, which is appropriate for a stationary Gaussian process. When solving at low frequency for a vessel's slow drift response process, which can be strongly non-Gaussian, the assumption that the maxima follow a Rayleigh distribution might be inadequate and yield poor estimations of the extremes.
Having conducted a wave frequency solve, if the wave elevation is defined by a single wave spectrum, without spreading, then the amplitude of the total transfer function, $\lvert\textrm{RAO}\rvert$, that maps the wave elevation process to the result process is given by the square root of the ratio of the spectra \begin{equation} \lvert\textrm{RAO}\rvert = \sqrt{\frac{S_{rr}(f)}{S_{\eta\eta}(f)}} \end{equation}
The above theory pertains to stochastic, or random wave, analysis, i.e. where the wave elevation is defined by wave spectra. Such stochastic analyses are the primary intended application of the frequency domain solver. However, for wave frequency solution frequencies the solver can alternatively be used with a single regular linear Airy wave as the input.
The regular wave frequency domain solver has been included in the software principally to conduct in-house testing and so has not been documented in full detail. A lot of the above theory is still relevant to the frequency domain solver when run with a regular wave. However, rather than solving for the stochastic response process, a single harmonic wave component is mapped, using the same linear transfer functions described above, to a single harmonic response component. Drag linearisation is still conducted by minimising the mean square error between the nonlinear and equivalent linear drag models. Because the wave elevation is deterministic, the expectations required for the drag linearisation can be evaluated by integrating in time over a single wave cycle.
Furthermore, again because the regular wave analysis is deterministic, result types associated with the stochastic solution (e.g. spectral density, spectral moments, MPM, $T_z$, $T_c$, and spectral response) are not available.