Modeling wave energy converters (WECs) 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.


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 direction is moving in the +X direction). The positive Z-axis is in the vertical upwards direction, and the positive 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.


WEC-Sim Coordinate System


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).

Boundary Element Method (BEM)

In WEC-Sim, wave forcing components are modeled using linear coefficients obtained from a frequency-domain potential flow Boundary Element Method (BEM) solver (e.g., WAMIT [4], AQWA [5], and Nemoh [6]). 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 [4].

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 scales the hydrodynamic coefficients according to the equations below, where \rho is the water density, \omega is the wave frequency in rad/s, and g is gravity:

|\overline{F_{exc}(\omega)}|=\frac{|F_{exc}(\omega)|}{\rho g}

\overline{A(\omega)} = \frac{A(\omega)}{\rho}

\overline{B(\omega)} = \frac{B(\omega)}{\rho \omega}

\overline{K_{hs}} = \frac{K_{hs}}{\rho g}

where F_{exc} is the wave-excitation force/torque coefficient, A is the radiation added mass coefficient, B is the radiation wave damping coefficient, and K_{hs} is the linear hydrostatic restoring coefficient.

Time-Domain Formulation

A common approach to determining the hydrodynamic forces is to use linear wave theory which assumes the waves are the sum of incident, radiated, and diffracted wave components. The dynamic response of the system is calculated by solving WEC system equations of motion [3][7]. 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 and torque (6-element) vector, F_{rad}(t) is the force and torque vector resulting from wave radiation, F_{pto}(t) is the PTO force and torque vector, F_{v}(t) is the damping force and torque vector, F_{me}(t) is the Morison Element force and torque vector, F_{B}(t) is the net buoyancy restoring force and torque vector, and F_{m}(t) is the force and torque 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 buoyancy term F_{B}(t) depends on the hydrostatic stiffness K_{hs} coefficient, displacement of the body, and its mass.

Numerical 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

\frac{1}{2}(1+\cos(\pi+\frac{\pi t}{t_{r}}) & \frac{t}{t_{r}}<1\\
1 & \frac{t}{t_{r}}\geq1

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

Sinusoidal Steady-State Response

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

F_{exc}(t)=\Re\left[ R_{f}(t)\frac{H}{2}F_{exc}(\omega, \theta)e^{i\omega t} \right]

where \Re denotes the real part of the formula, R_{f} is the ramp function, H is the wave height, F_{exc} is the frequency dependent complex wave-excitation amplitude vector, and \theta is the wave direction.

Convolution Integral Formulation

To include the fluid memory effect, 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. This representation also assumes that there is no motion for t<0.

For regular waves, the equation described in the last subsection is used to calculate the wave excitation vector. 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, generally 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

F_{exc}(t)=\Re\left[ R_{f}(t) \sum_{j=1}^{N}F_{exc}(\omega_{j}, \theta)e^{i(\omega_{j}t+\phi_{j})} \sqrt{2S(\omega_{j})d\omega_{j}} \right]

where \phi is the randomized phase angle and N is the number of frequency bands selected to discretize the wave spectrum. For repeatable simulation of an irregular wave field S(\omega), WEC-Sim allows specification of \phi, refer to the following wave features section.

State Space

It is highly desirable to represent the radiation convolution integral described 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:

\dot{X}_{r} \left( t \right) = \mathbf{A_{r}} X_{r} \left( t \right) + \mathbf{B_{r}} \mathbf{u} (t);~~X_{r}\left( 0 \right) = 0~~ \nonumber \\
\int_{0}^{t} \mathbf{K_{r}} \left( t- \tau \right) d\tau \approx \mathbf{C_{r}} X_{r} \left( t \right) + \mathbf{D_{r}} \mathbf{u} \left( t \right)~~

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 u 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

\dot{x} =  \mathbf{A_{r}}x + \mathbf{B_{r}} u~~ \\
y = \mathbf{C_{r}}x~~

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

\dot{x} = \mathbf{A_{r}} x~~, \\
y = \mathbf{C_{r}} x~~

is clearly equivalent to the zero-state impulse response. 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:

x(t) = e^{\mathbf{A_{r}}t} x(0) + \int_{0}^{t} e^{\mathbf{A_{r}}(t-\tau)} \mathbf{B_{r}} u (\tau) d\tau~~

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

K_{r}(t) = \mathbf{C_{r}}e^{\mathbf{A_{r}}t}\mathbf{B_{r}}~~

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:

\mathbf{\tilde{K}_{r}} \left( t_{k} \right) = \mathbf{C_{d}}\mathbf{A_{d}}^{k}\mathbf{B_{d}}~~

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

H = \begin{bmatrix}
        \mathbf{K_{r}}(2) & \mathbf{K_{r}}(3) & \ldots & \mathbf{K_{r}}(n) \\
        \mathbf{K_{r}}(3) & \mathbf{K_{r}}(4) & \ldots & 0 \\
        \vdots & \vdots & \ddots & \vdots \\
        \mathbf{K_{r}}(n) & 0 & \cdots & 0
\end{bmatrix} &\\

can be reproduced exactly by the SVD as

H = \mathbf{U} \Sigma \mathbf{V^{*}}

where \Sigma is a diagonal matrix containing the Hankel singular values 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 here and the reader is referred to [10][12].

Regular Waves

Regular waves are defined as planar sinusoidal waves, where the incident wave is defined as \eta(x,y,t) :

\eta(x,y,t)= \frac{H}{2} \cos( \omega t - k (x\cos \theta + y\sin \theta) + \phi)

where H is the wave height, \omega is the wave frequency (\omega = \frac{2\pi}{T}), k is the wave number (k = \frac{2\pi}{\lambda}), \theta is the wave direction, and \phi is the wave phase.

Irregular Waves

Irregular waves are modeled as the linear superposition of a large number of harmonic waves at different frequencies and angles of incidence, where the incident wave is defined as \eta(x,y,t) :

\eta(x,y,t)= \sum_{i}\frac{H_{i}}{2} \cos( \omega_{i} t - k_{i} (x\cos \theta_{i} + y\sin \theta_{i}) + \phi_{i})

where H is the wave height, \omega is the wave frequency (\omega = \frac{2\pi}{T}), k is the wave number (k = \frac{2\pi}{\lambda}), \theta is the wave direction, and \phi is the wave phase (randomized for irregular waves).

Wave Spectra

The linear 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 wave spectra that are used by the offshore industry are discussed in the following sections. The general form of the wave spectra available in WEC-Sim is given by:

S\left( f , \theta \right)= S\left( f \right)D\left( \theta \right)~~

where S\left( f\right) is the wave power spectrum, f is the wave frequency (in Hertz), D\left( \theta \right) is the directional distribution, and \theta is the wave direction (in Degrees). The formulation of D\left( \theta \right) requires that

\int_{0}^{\infty}\int_{-\pi}^{\pi} S\left( f \right)D\left( \theta \right) d\theta df = \int_{0}^{\infty} S\left( f \right) df

so that the total energy in the directional spectrum must be the same as the total energy in the one-dimensional spectrum.

S\left( f \right) = A f^{-5}\exp\left[-B f^{-4} \right]~~

where A and B are coefficients that vary depending on the wave spectrum and \exp stands for the exponential function. Spectral moments of the wave spectrum, denoted m_{k}~,~k=0, 1, 2,..., are defined as

m_{k} = \int_{0}^{\infty} f^{k} S \left( f \right) df ~~

The spectral moment, m_{0} is the variance of the free surface which allows one to define the mean wave height of the tallest third of waves, significant wave height H_{m0} (in m), as:

H_{m0} = 4 \sqrt{m_{0}}~~

Pierson–Moskowitz (PM)

One of the simplest spectra, the Pierson–Moskowitz spectrum, 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:

S\left( f \right) = \frac{\alpha_{PM}g^{2}}{\left( 2 \pi \right)^{4}}f^{-5}\exp\left[-\frac{5}{4} \left( \frac{f_{p}}{f}\right)^{4} \right]~~ \\

This implies coefficients of the general form:

A = \frac{\alpha_{pm}g^{2}}{\left( 2 \pi \right)^{4}},~~B = \frac{5}{4} {f_{p}}^{4}~~

where the parameter \alpha_{PM} = 0.0081 typically, g=9.81 m/s is gravitational acceleration and f_{p} is the peak frequency of the spectrum.


Pierson-Moskowitz does not use significant wave height to define spectrum

Bretschneider (BS)

The two-parameter Bretschneider 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:

S\left( f \right) = \frac{{H_{m0}}^2}{4}\left(1.057f_{p}\right)^{4}f^{-5}\exp\left[-\frac{5}{4} \left( \frac{f_{p}}{f}\right)^{4} \right]~~ \\

This implies coefficients of the general form:

A =\frac{{H_{m0}}^2}{4}\left(1.057f_{p}\right)^{4} \approx \frac{5}{16} {H_{m0}}^2 {f_{p}}^{4} \approx \frac{B}{4}{H_{m0}}^2~~ \\

B = \left(1.057f_{p}\right)^{4} \approx \frac{5}{4} {f_{p}}^{4}~~ \\


The JONSWAP (Joint North Sea Wave Project) spectrum was purposed by Hasselmann et al. [14], and the original formulation was given as:

& S\left( f \right) = \frac{ \alpha_{js} g^{2} }{ (2\pi)^{4}} f^{-5}\exp\left[-\frac{5}{4} \left( \frac{f_{p}}{f}\right)^{4} \right]\gamma^\Gamma \nonumber  ~~ &\\

&\Gamma = \exp \left[ -\left( \frac{\frac{f}{f_{p}}-1}{\sqrt{2} \sigma}\right)^{2} \right],~~ \sigma = \begin{cases} 0.07 & f \leq f_{p} \\0.09 & f > f_{p} \end{cases} ~~ &\\

with general form coefficients thus defined:

& A =\frac{\alpha_{js} g^{2}}{(2\pi)^{4}} & \\

& B=\frac{5}{4}{f_{p}}^{4} &\\

where \alpha_{js} 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 (default) , \alpha_{js} , and S^{*}\left( f \right) can be calculated by:

\alpha_{js} = \frac{H_{m0}^{2}}{16\int_{0}^{\infty} S^{*} \left( f \right) df}

S^{*}\left( f \right) = \frac{ g^{2} }{ (2\pi)^{4}} f^{-5}\exp\left[-\frac{5}{4} \left( \frac{f_{p}}{f}\right)^{4} \right]\gamma^\Gamma ~~


S\left( f \right) =  S^{*}\left( f \right) \alpha_{js} \\

Power Take-Off (PTO)

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

Linear 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 instantaneous 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)

Hydraulic PTO

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

F_{pto}=\Delta{} p_{piston}A_{piston}

where \Delta{} p_{piston} is the differential pressure of the hydraulic piston and A_{piston} is the piston area. The instantaneous hydraulic 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), \lambda_{fd} is the flux linkage of the stator d-axis winding due to flux produced by the rotor magnets, and i_{sq} is the stator q-axis current. The instantaneous mechanical power absorbed by the PTO is given by:


For more information about application of pto systems in WEC-Sim, refer to Constraint and PTO Features.


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.

Mooring Matrix

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].


MoorDyn mooring model elements

For more information about application of mooring systems in WEC-Sim, refer to Mooring Features .

Nonlinear Buoynancy and Froude-Krylov Wave Excitation

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 using the BEM calculated linear wave-excitation and hydrostatic 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 [17], which applies a correction to the instantaneous wave elevation that forces its height to be 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.


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

For more information about application of nonlinear hydrodynamics in WEC-Sim, refer to Nonlinear Buoyancy and Froude-Krylov Wave Excitation.

Viscous Damping and Morison Elements

Additional damping and added-mass can be added to the WEC system. This facilitates experimental validation of the WEC-Sim code, particularly in the event that the BEM hydrodynamic outputs are not sufficiently representative of the physical system.

Viscous Damping

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

& F_{v}=-C_{v}\dot{X}-\frac{C_{d} \rho A_{d}}{2}\dot{X}|\dot{X}| & \\

&  =-C_{v}\dot{X}-C_{D}\dot{X}|\dot{X}| & \\

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

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, \lambda, which is generally met when D/\lambda < 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 [18]. The formulation for each element on the body can be given as

F_{me}=\rho\forall\dot{v}+\rho\forall C_{a}(\dot{v}-\ddot{X})+\frac{C_{d}\rho A_{d}}{2}(v-\dot{X})|v-\dot{X}|

where v is the fluid particle velocity due to wave and current, C_{a} is the coefficient of added mass, and \forall is the displaced volume.


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

For more information about application of Morison Elements in WEC-Sim, refer to Morison Elements.



Matthew Hall. Moordyn user’s guide. 2015.


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.


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.


CH Lee and JN Newman. WAMIT User Manual. 2006.




Nemoh a Open source BEM. URL:


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.


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


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


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.


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.


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.


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.


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.


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.


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.


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


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.