# Theory¶

The section provides an overview of the underlying theory behind the WEC-Sim code.

## Introduction¶

Modeling a wave energy converter (WEC) involves the interaction between the incident waves, device motion, power-take-off (PTO mechanism), and mooring. WEC-Sim uses a radiation and diffraction method [2][3] to predict power performance and design optimization. The radiation and diffraction method generally obtains the hydrodynamic forces from a frequency-domain boundary element method (BEM) solver using linear coefficients to solve the system dynamics in the time domain.

*WEC-Sim Methodology*

## Coordinate System¶

The WEC-Sim Coordinate System figure illustrates a 3-D floating point absorber subject to incoming waves in water. The figure also defines the coordinates and the 6 degree of freedom (DOF) in WEC-Sim. The WEC-Sim coordinate system assumes that the positive X-axis defines a wave angle heading of zero (e.g., a wave propagating with along a zero-degree heading is moving in the +X direction). The Z-axis is in the vertical upwards direction, and the Y-axis direction is defined by the right-hand rule. In the vectors and matrices used in the code, Surge (x), Sway (y), and Heave (z) correspond to the first, second and third position respectively. Roll (Rx), Pitch (Ry), and Yaw (Rz) correspond to the fourth, fifth, and sixth position respectively.

## Units¶

All units within WEC-Sim are in the MKS (meters-kilograms-seconds system) and angular measurements are specified in radians (except for wave directionality which is defined in degrees).

## Time-Domain Numerical Method¶

A common approach to determining the hydrodynamic forces is to use the linear wave theory assumption that waves are the sum of incident, radiated, and diffracted wave components. The dynamic response of the system is calculated by solving the equation of motion for WEC systems [3][4]. The equation of motion for a floating body about its center of gravity can be given as:

where \(\ddot{X}\) is the (translational and rotational) acceleration vector of the device, \(m\) is the mass matrix, \(F_{exc}(t)\) is the wave excitation force (6-element) vector, \(F_{rad}(t)\) is the force vector resulting from wave radiation, \(F_{PTO}(t)\) is the PTO force vector, \(F_{v}(t)\) is the damping force vector, \(F_{ME}(t)\) is the Morison Element Force vector, \(F_{B}(t)\) is the net buoyancy restoring force vector, and \(F_{m}(t)\) is the force vector resulting from the mooring connection.

\(F_{exc}(t)\) , \(F_{rad}(t)\) , and \(F_{B}(t)\) are calculated using hydrodynamic coefficients provided by the frequency-domain BEM solver. The radiation term includes an added-mass term, matrix \(A(\omega)\), and wave damping term, matrix \(B(\omega)\), associated with the acceleration and velocity of the floating body, respectively, and given as functions of radian frequency (\(\omega\)) by the BEM solver. The wave excitation term \(F_{exc}(\omega)\) includes a Froude-Krylov force component generated by the undisturbed incident waves and a diffraction component that results from the presence of the floating body. The bouyancy term \(F_{B}(t)\) depends on the hydrostatic stiffness \(K_{hs}\) coefficient, displacement of the body, and its mass.

## Boundary Element Method¶

In WEC-Sim, wave forcing components are modeled using linear coefficients obtained from a frequency-domain potential flow BEM solver (e.g., WAMIT [5], AQWA [6], and Nemoh [7]). The BEM solutions are obtained by solving the Laplace equation for the velocity potential, which assumes the flow is inviscid, incompressible, and irrotational. More details on the theory for the frequency-domain BEM can be found in [5].

WEC-Sim imports non-dimensionalized hydrodynamic coefficients from an `*.h5`

data structure generated by BEMIO for the WAMIT, AQWA or Nemoh BEM solvers.
Alternatively, the `*.h5`

data structure can be manually defined by the user.
The WEC-Sim code then scales the hydrodynamic coefficients according to the scaling defined below, where \(\rho\) is the water density, \(\omega\) is the wave frequency in rad/s, and g is gravity:

## Solver Methods¶

WEC-Sim can be used for regular and irregular wave simulations, but note that \(F_{exc}(t)\) and \(F_{rad}(t)\) are calculated differently for sinusoidal steady-state response scenarios and random sea simulations. The sinusoidal steady-state response is often used for simple WEC designs with regular incoming waves. However, for random sea simulations or any simulations where fluid memory effects of the system are essential, the convolution integral method is recommended to represent the fluid memory retardation force on the floating body. To speed computation of the convolution integral, the state space representation method can be specified to approximate this calculation as a system of linear ordinary differential equations.

### Ramp Function¶

A ramp function (\(R_{f}\)), necessary to avoid strong transient flows at the earlier time steps of the simulation, is used to calculate the wave excitation force. The ramp function is given by

where \(t\) is the simulation time and \(t_{r}\) is the ramp time.

### Sinusoidal Steady-State Response Scenario¶

This approach assumes that the system response is in sinusoidal steady-state form; therefore, it is only valid for regular wave simulations. The radiation term can be calculated using the added mass and the wave radiation damping term for a given wave frequency, which is obtained from

where \(\dot{X}\) is the velocity vector of the floating body.

The free surface profile is based on linear wave theory for a given wave height, wave frequency, and water depth. The regular wave excitation force is obtained from

where \(\Re\) denotes the real part of the formula, \(R_{f}\) is the ramp function, \(H\) is the wave height, and \(F_{exc}\) is the excitation vector, including the magnitude and phase of the force.

### Convolution Integral Formulation¶

To include the fluid memory effect on the system the convolution integral formulation based upon the Cummins equation [8] is used. The radiation term can be calculated by

where \(A_{\infty}\) is the added mass matrix at infinite frequency and \(K_{r}\) is the radiation impulse response function.

For regular waves, the equation described in the last subsection is used to calculate the wave excitation force. For irregular waves, the free surface elevation is constructed from a linear superposition of a number of regular wave components. Each regular wave component is extracted from a wave spectrum, \(S(\omega)\), describing the wave energy distribution over a range,of wave frequencies, characterized by a significant wave height and peak wave period. The irregular excitation force can be calculated as the real part of an integral term across all wave frequencies as follows

where \(\phi\) is a random phase angle. For repeatable simulation of an irregular wave field \(S(\omega)\), WEC-Sim allows specification of \(\phi\).

*Example irregular wave surface elevation generated by WEC-Sim*

### State Space Representation of \(K_{r}\)¶

It is highly desirable to represent the radiation convolution integral discribed in the last subsection in state space (SS) form [9]. This has been shown to dramatically increase computational speeds [10] and allow utilization of conventional control methods that rely on linear state space models. An approximation will need to be made as \(K_{r}\) is solved from a set of partial differential equations where as a linear state space is constructed from a set of ordinary differential equations. In general, a linear system is desired such that:

with \(\mathbf{A_{r}},~\mathbf{B_{r}},~\mathbf{C_{r}},~\mathbf{D_{r}}\) being the time-invariant state, input, output, and feed through matrices, while \(\dot{\zeta}\) is the input to the system and \(X_{r}\) is the state vector describing the convolution kernal as time progresses.

#### Calculation of \(K_{r}\) from State Space Matrices¶

The impulse response of a single-input zero-state state-space model is represented by

where \(u\) is an impulse. If the initial state is set to \(x(0)= \mathbf{B_{r}} u\) the response of the unforced (\(u=0\)) system

is clearly equivalent to the zero-state impulse reponse. The impulse response of a continuous system with a nonzero \(\mathbf{D_r}\) matrix is infinite at \(t=0\); therefore, the lower continuity value \(\mathbf{C_{r}}\mathbf{B_{r}}\) is reported at \(t=0\). The general solution to a linear time invariant (LTI) system is given by:

where \(e^{\mathbf{A_{r}}}\) is the matrix exponential and the calculation of \(K_{r}\) follows:

#### Realization Theory¶

The state space realization of the hydrodynamic radiation coefficients can be pursued in the time domain (TD). This consists of finding the minimal order of the system and the discrete time state matrices (\(\mathbf{A_{d}},~\mathbf{B_{d}},~\mathbf{C_{d}},~\mathbf{D_{d}}\)) from samples of the impulse response function. This problem is easier to handle for a discrete-time system than for continuous-time. The reason being is that the impulse response function of a discrete-time system is given by the Markov parameters of the system:

where \(t_{k}=k\Delta t\) for \(k=0,~1,~2,~\ldots\) with \(\Delta t\) being the sampling period. The feedthrough matrix \(\mathbf{D_d}\) is assumed to be zero in order to maintain causality of the system, as a non-zero mathbf{D_d} results in an infinite value at \(t=0\).

The most common algorithm to obtain the realization is to perform a Singular Value Decomposition (SVD) on the Hankel matrix of the impulse response function, as proposed by Kung [11]. The order of the system and state-space parameters are determined from the number of significant singular values and the factors of the SVD. The Hankel matrix (\(H\)) of the impulse response function

can be reproduced exactly by the SVD as

where \(\Sigma\) is a diagonal matrix containing the Hankel singular vales in descending order. Examination of the Hankel singular values reveals there are only a small number of significant states and that the rank of \(H\) can be greatly reduced without a significant loss in accuracy [10][12]. Further detail into the SVD method and calculation of the state space parameters will not be discussed and the reader is referred to [10][12].

## Wave Spectrum¶

A superposition of regular waves of distinct amplitudes and periods is characterized in the frequency domain by a wave spectrum. Through statistical analysis, spectra are characterized by specific parameters such as significant wave height, peak period, wind speed, fetch length, and others. Common types of spectra that are used by the offshore industry are discussed in the following sections. The general form of the sea spectrums available in WEC-Sim is given by:

where \(f\) is the wave frequency (in Hertz), \(A\) is wave amplitude (m), and \(\exp\) stands for the exponential function. Spectral moments of the wave spectrum, denoted \(m_{k}~,~k=0, 1, 2,...\), are defined as

The spectral moment, \(m_{0}\) is the variance of the free surface which allows one to define

where \(H_{m0}\) is a definition of the significant wave height (m), the mean wave height of the tallest third of waves.

### Pierson–Moskowitz¶

One of the simplest spectra was proposed by [13]. It assumed that after the wind blew steadily for a long time over a large area, the waves would come into equilibrium with the wind. This is the concept of a fully developed sea where a “long time” is roughly 10,000 wave periods and a “large area” is roughly 5,000 wave-lengths on a side. The spectrum is calculated from:

This implies coefficients of the general form:

where where parameter \(\alpha_{PM}\) = 0.0081 typically, \(g=9.81\) m/s is gravitational acceleration and \(f_{p}\) is the peak frequency of the spectrum. However, this spectrum representation does not allow the user to define the significant wave height: WEC-Sim calculates the value of \(\alpha_{PM}\) such that a power matrix can be developed as a function of the desired significant wave height. The \(\alpha_{PM}\) parameter was calculated as follows:

Where:

### Bretschneider Spectrum¶

This two-parameter spectrum is based on significant wave height and peak wave frequency. For a given significant wave height, the peak frequency can be varied to cover a range of conditions including developing and decaying seas. In general, the parameters depend on strongly on wind speed, and also wind direction, fetch, and locations of storm fronts. The spectrum is given as:

This implies coefficients of the general form:

### JONSWAP (Joint North Sea Wave Project) Spectrum¶

The spectrum was purposed by Hasselmann et al. [14], and the original formulation was given as:

with general form coefficients thus defined:

where \(\alpha_{j}\) is a nondimensional variable that is a function of the wind speed and fetch length. Empirical fits were applied in an attempt to find a mean value that would capture the spectral shape of most measured sea states. For a given significant wave height, setting \(\gamma = 3.3\) , \(\alpha_{j}\) , and \(S^{*}\left( f \right)\) can be calculated by:

Where:

## Power Take-off Forces¶

Throughout the following sections, unless specification is made between linear and rotary PTOs, units are not explicitly stated.

### Linear Spring-damper PTO¶

The PTO mechanism is represented as a linear spring-damper system where the reaction force is given by:

where \(K_{PTO}\) is the stiffness of the PTO, \(C_{PTO}\) is the damping of the PTO, and \(X_{rel}\) and \(\dot{X}_{rel}\) are the relative motion and velocity between two bodies. The power absorbed by the PTO is given by:

\[P_{PTO} = -F_{PTO}\dot{X}_{rel}=\left(K_{PTO}X_{rel}\dot{X}_{rel}+C_{PTO}\dot{X}^{2}_{rel}\right)\]

However, the relative motion and velocity between two bodies is out of phase by \(\pi/2\), resulting in a time-averaged product of 0. This allows the absorbed power to be written as

### Hydraulic PTO¶

The PTO mechanism is modeled as a hydraulic system [15], where the reaction force is given by:

where \(\Delta{} p_{piston}\) is the differential pressure of hydraulic piston and \(A_{piston}\) is the piston area. The power absorbed by the PTO is given by:

### Mechanical PTO¶

The PTO mechanism is modeled as a direct-drive linear generator system [15], where the reaction force is given by:

where \(\tau_{pm}\) is the magnet pole pitch (the center-to-center distance of adjacent magnetic poles, m), \(\lambda_{fd}\) is the flux linkage of the stator \(d\)-axis winding due to flux produced by the rotor magnets (Volt-second), and \(i_{sq}\) is the stator \(q\)-axis current (Amps). The power absorbed by the PTO is given by:

## Mooring Forces¶

The mooring load is represented using a linear quasi-static mooring stiffness or by using the mooring forces calculated from MoorDyn [1], which is an open-source lumped-mass mooring dynamics model.

When linear quasi-static mooring stiffness is used, the mooring load can be calculated by

where \(K_{m}\) and \(C_{m}\) are the stiffness and damping matrices for the mooring system, and \(X\) and \(\dot{X}\) are the displacement and velocity of the body, respectively.

MoorDyn discretizes each mooring line in a mooring system into evenly-sized line segments connected by node points (see MoorDyn figure). The line mass is lumped at these node points along with gravitational and buoyancy forces, hydrodynamic loads, and reactions from contact with the seabed. Hydrodynamic drag and added mass are calculated based on Morison’s equation. A mooring line’s axial stiffness is modeled by applying a linear stiffness to each line segment in tension only. A damping term is also applied in each segment to dampen non-physical resonances caused by the lumped-mass discretization. Bending and torsional stiffnesses are neglected. Bottom contact is represented by vertical stiffness and damping forces applied at the nodes when a node is located below the seabed. [16].

## Additional Added-mass & Damping Forces¶

Additional added-mass and damping forces can be added to force definitions. This facilitates experimental validation of the WEC-Sim code, particularly in the even that the BEM hydrodynamic outputs are not sufficiently representative of actual device properties.

### Linear & Quadratic Damping Forces¶

Linear and quadratic damping forces add flexibility to the definition of viscous forcing

\[F_{v}=-C_{ld}\dot{X}-\frac{1}{2}C_{d}\rho A_{D}\dot{X}|\dot{X}|\]

where \(C_{ld}\) is the linear damping coefficient, \(C_{d}\) is the (quadratic) viscous drag coefficient, \(\rho\) is the fluid density, and \(A_{D}\) is the characteristic area for drag calculation.

Because BEM codes are potential flow solvers and neglect the effects of viscosity, \(F_{v}\) generally must be included to accurately model device performance. However, it can be difficult to select representative drag coefficients, as they depend on device geometry, scale, and relative velocity between the body and the flow around it. Empirical data on the drag coefficient can be found in various literature and standards, but is generally limited to simple geometries evaluated at a limited number of scales and flow conditions. For realistic device geometries, the use of computational fluid dynamic simulations or experimental data is encouraged.

### Morison Elements¶

The Morison Equation assumes that the fluid forces in an oscillating flow on a structure of slender cylinders or other similar geometries arise partly from pressure effects from potential flow and partly from viscous effects. A slender cylinder implies that the diameter, D, is small relative to the wave length, \(λ_w\), which is generally met when \(D/λ_w < 0.1 − 0.2\). If this condition is not met, wave diffraction effects must be taken into account. Assuming that the geometries are slender, the resulting force can be approximated by a modified Morison formulation [17]. The formulation for each element on the body can be given as

\[F_{ME}=\rho∀\dot{v}+\rho\forall C_{a}(\dot{v}-\ddot{X})+\frac{1}{2}C_{d}\rho A_{D}(v-\dot{X})|v-\dot{X}|\]

where \(v\) is the fluid particle velocity, \(C_{a}\) is the coefficient of added mass, and \(\forall\) is the displaced volume.

Note

WEC-Sim does not consider buoyancy effects when calculating the forces from Morison elements.

## Nonlinear Hydrodynamic Forces¶

The linear model assumes that the body motion and the waves consist of small amplitudes in comparison to the wavelengths. A weakly nonlinear approach is applied to account for the nonlinear hydrodynamic forces induced by the instantaneous water surface elevation and body position. Rather than the BEM calculated linear hydrodynamic force coefficients, the nonlinear buoyancy and the Froude-Krylov force components can be obtained by integrating the static and dynamic pressures over each panel along the wetted body surface at each time step. Because linear wave theory is used to determine the flow velocity and pressure field, the values become unrealistically large for wetted panels that are above the mean water level. To correct this, the Wheeler stretching method is applied [18], which applies a correction to the instantaneous wave elevation that forces its height to equal to the water depth when calculating the flow velocity and pressure,

\[z^* = \frac{D(D+z)}{(D+\eta)} - D\]

where \(D\) is the mean water depth, and \(\eta\) is the z-value on the instantaneous water surface.

Note

The nonlinear WEC-Sim method is not intended to model highly nonlinear hydrodynamic events, such as wave slamming and wave breaking.

## References¶

[1] | Matthew Hall. Moordyn user’s guide. 2015. |

[2] | Ye Li and Yi-Hsiang Yu. A Synthesis of Numerical Methods for Modeling Wave Energy Converter-Point Absorbers. Renewable and Sustainable Energy Reviews, 16(6):4352–4364, 2012. |

[3] | (1, 2) Aurélien Babarit, J. Hals, M.J. Muliawan, A. Kurniawan, T. Moan, and J. Krokstad. Numerical Benchmarking Study of a Selection of Wave Energy Converters. Renewable Energy, 41:44–63, 2012. |

[4] | Jerica D. Nolte and R. C. Ertekin. Wave power calculations for a wave energy conversion device connected to a drogue. Journal of Renewable and Sustainable Energy, 2014. |

[5] | (1, 2) CH Lee and JN Newman. WAMIT User Manual. 2006. |

[6] | ANSYS Aqwa. URL: http://www.ansys.com/Products/Other+Products/ANSYS+AQWA. |

[7] | Nemoh a Open source BEM. URL: http://openore.org/tag/nemoh/. |

[8] | WE Cummins. The Impulse Response Function and Ship Motions. Technical Report, David Taylor Model Dasin-DTNSRDC, 1962. |

[9] | Z. Yu and J. Falnes. State-space Modelling of a Vertical Cylinder in Heave. Applied Ocean Research, 5:265–275, 1996. |

[10] | (1, 2, 3) R. Taghipour, T. Perez, and T. Moan. Hybrid Frequency-time Domain Models for Dynamic Response Analysis of Marine Structures. Ocean Engineering, 7:685–705, 2008. |

[11] | Sun-Yuan Kung. A new identification and model reduction algorithm via singular value decomposition. In Proc. 12th Asilomar Conf. Circuits, Syst. Comput., Pacific Grove, CA, 705–714. 1978. |

[12] | (1, 2) Erlend Kristiansen, Åsmund Hjulstad, and Olav Egeland. State-space representation of radiation forces in time-domain vessel models. Ocean Engineering, 32(17):2195–2216, 2005. |

[13] | W. J. Pierson and L. A. Moskowitz. Proposed Spectral Form for Fully Developed Wind Seas Based on the Similarity Theory of S. A. Kitaigorodskii. Geophysical Research, 69:5181–5190, 1964. |

[14] | K. Hasselman, T.P. Barnett, E. Bouws, H. Carlson, D.E. Cartwright, K. Enke, J.A. Ewing, H. Gienapp, D.E. Hasselmann, P. Kruseman, A. Meerburg, P. Mller, D.J. Olbers, K. Richter, W. Sell, and H. Walden. Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP). Technical Report 12, German Hydrographic Institute, 1973. |

[15] | (1, 2) R. So, A. Simmons, T. Brekken, K. Ruehl, and C. Michelen. Development of PTO-SIM: A Power Performance Module for the Open-Souse Wave Energy Converter Code WEC-SIM. In Proceedings of OMAE 2015. 2015. |

[16] | Matthew Hall and Andrew Goupee. Validation of a lumped-mass mooring line model with deepcwind semisubmersible model test data. Ocean Engineering, 104:590–603, 2015. |

[17] | J R Morison, M P O’brien, J W Johnson, and S A Schaaf. The force exerted by surface waves on piles. Petroleum Transactions, AIME, 189:149–154, 1950. |

[18] | J D Wheeler and others. Methods for calculating forces produced by irregular waves. In Offshore Technology Conference. 1969. |

[1] | K. Ruehl, C. Michelen, S. Kanner, M. Lawson, and Y. Yu. Preliminary Verification and Validation of WEC-Sim, an Open-Source Wave Energy Converter Design Tool. In Proceedings of OMAE 2014. San Francisco, CA, 2014. |