Development of a Blended Time-Domain Program for Predicting the Motions of a Wave Energy Structure

Traditional linear time-domain analysis is used widely for predicting the motions of floating structures. When it comes to a wave energy structure, which usually is subjected to larger relative (to their geometric dimensions) wave and motion amplitudes, the nonlinear effects become significant. This paper presents the development of an in-house blended time-domain program (SIMDYN). SIMDYN’s “blend” option improves the linear option by accounting for the nonlinearity of important external forces (e.g., Froude-Krylov). In addition, nonlinearity due to large body rotations (i.e., inertia forces) is addressed in motion predictions of wave energy structures. Forced motion analysis reveals the significance of these nonlinear effects. Finally, the model test correlations examine the simulation results from SIMDYN under the blended option, which has seldom been done for a wave energy structure. It turns out that the blended time-domain method has significant potential to improve the accuracy of motion predictions for a wave energy structure.


Introduction
Wave energy is distributed extensively in coastal areas and holds huge potential for wider utilization [1]. Many wave energy structures (WESs) were designed to convert wave energy to electricity [2,3]. While some WESs are fixed onto the shore or the seabed, many are floating or submerged structures (see Figure 1). It is worth noting that WESs may be any structure subjected to wave energy. In this paper, the floating power system (FPS) used in the model test correlations is a necessary unit in the wave energy conversion system, but the FPS itself is not a wave energy converter (WEC). That is the reason why we used the more general term "WESs" to refer to it.
In designing a WES, motions, internal loads (represented by local acceleration), and mooring loads are all essential quantities [4]. Predicting the motions accurately and efficiently is important given that the sea states (especially the extreme ones); and other quantities are closely correlated with the motions. A range of methods can model the hydrodynamics of a floating or submerged WES [4]. Generally, the cost and time involved in the application of a modelling method (see Table 1) increase with its fidelity [5]. The methods outlined below in Table 1 have been used comprehensively to take advantage of either their efficient turnaround or high fidelity. In designing a WES, motions, internal loads (represented by local acceleration), and mooring loads are all essential quantities [4]. Predicting the motions accurately and efficiently is important given that the sea states (especially the extreme ones); and other quantities are closely correlated with the motions. A range of methods can model the hydrodynamics of a floating or submerged WES [4]. Generally, the cost and time involved in the application of a modelling method (see Table 1) increase with its fidelity [5]. The methods outlined below in Table 1 have been used comprehensively to take advantage of either their efficient turnaround or high fidelity.
The blended time-domain method usually accounts for the nonlinearity in the Froude-Krylov forces, hydrostatic forces, and the equations of motion (inertia forces), while the remaining forces are computed using either linear or nonlinear modeling [10]. The blended method has been applied to simulate ship motions by many researchers [20][21][22]. The Wave Energy Converter SIMulator (WEC-Sim), developed by the Sandia National Laboratory and the National Renewable Energy Laboratory, integrates the blended method option [23], but it has not been compared with the modeling test results of a WES.
The nonlinear time-domain potential flow method solves the fluid flow problem using the fully nonlinear free-surface boundary conditions instead of the linearized free-surface boundary
The blended time-domain method usually accounts for the nonlinearity in the Froude-Krylov forces, hydrostatic forces, and the equations of motion (inertia forces), while the remaining forces are computed using either linear or nonlinear modeling [10]. The blended method has been applied to simulate ship motions by many researchers [20][21][22]. The Wave Energy Converter SIMulator (WEC-Sim), developed by the Sandia National Laboratory and the National Renewable Energy Laboratory, integrates the blended method option [23], but it has not been compared with the modeling test results of a WES.
The nonlinear time-domain potential flow method solves the fluid flow problem using the fully nonlinear free-surface boundary conditions instead of the linearized free-surface boundary conditions used in the linear and blended time-domain methods [24]. These nonlinear modeling approaches have been reviewed in past research [25].
Computational fluid dynamics (CFD) methods overcome the discrepancies of all the other methods based upon the assumptions of the potential flow theory (by including fluid viscosity); they are also used to simulate WESs [26]. CFD methods include Reynolds-averaged Navier-Stokes (e.g., [14]), smoothed-particle hydrodynamics (e.g., [15]), and large eddy simulation (e.g., [16]). These CFD methods are capable of capturing the full range of phenomena in extreme conditions [4]. Therefore, CFD methods are more accurate than the previously described methods, but the computational time is significantly longer.
Our research group implemented the blended time-domain method using the SYMDYN time-domain program (an in-house program developed in FORTRAN). It was originally developed to simulate the six-degrees-of-freedom motions of ships [19], SIMDYN also has the capability to simulate the linear time domain for ships. In this work, SIMDYN was improved by calculating Froude-Krylov forces, hydrostatic forces, and inertia forces to account for their nonlinear effects. After coupling with the open-source Mooring Analysis Program (MAP++) [27] and implementing other forces when necessary, SIMDYN became capable of simulating WES motions. "Other forces" include slowly varying wave drift forces (when using the linear time-domain option for random waves), Morison element forces, and power take-off forces.
Experimental verification and study have been conducted on typical WESs such as the oscillating water column (e.g., [28]), the point absorber (e.g., [29]), and the oscillating surge wave energy converter (e.g., [30]). The model test correlations in these studies have examined and improved simulation tools in different ways. In 2014, the Hydraulics and Maritime Research Centre (HMRC) (Beaufort Research) conducted model tests on a floating power system (FPS) designed to provide power and connection to the grid [31]. The FPS was used in the model test correlation because the suitable (and accessible) model test data was very limited. The FPS was tested under multiple regular (single-frequency) wave cases and multiple irregular (random) wave cases. The FPS did not generate electricity in the model tests, thus reducing uncertainty. These facts make it a perfect benchmark model for the study herein. Therefore, the motions of the FPS under the sea states tested in the experiments were simulated in SIMDYN to verify the accuracy of the program. The blended time-domain method and our program, SIMDYN, are general and can be applied to different types of WECs (e.g., the oscillating water column, the point absorber). In fact, the geometry and dimensions of the floating power system are similar to a point absorber type WEC. Therefore, we expect the model should work (similarly) well on a real, point absorber type WEC.
The rest of the paper is organized as follows. Section 2 gives the mathematical model of SIMDYN, including formulating the nonlinear equations of motion and the external forces. Section 3 discusses the nonlinear effects of Froude-Krylov and hydrostatic forces through a series of forced motion analyses. Section 4 details the nonlinear inertia forces due to a large-angle of rotations and discusses their significance. Section 5 presents the model test and simulations and, more importantly, the model test correlations for the regular wave cases and irregular wave cases. Section 6 is the discussions while Section 7 is the conclusion.

Coordinate Systems
To describe rigid body motions conveniently, two right-handed coordinate systems were used. The global coordinate system (GCS) is fixed to the earth (see Figure 1), a point of which is x = (x, y, z). The local coordinate system (LCS) translates and rotates with the rigid body; a point in the LCS is x = (x , y , z ).
Vector ξ = (ξ 1 , ξ 2 , ξ 3 ) represents the translation from the GCS origin (0,0,0) at the calm water line to the LCS origin (ξ 1 , ξ 2 , ξ 3 ). Rotation vector α = (ξ 4 , ξ 5 , ξ 6 ) consists of the Euler angles between the GCS and the LCS. The coordinates in the LCS, x , are related to the coordinates in the GCS, x, through rotation matrix R: In Equation (1), c i represents cos(ξ i ), and s i represents sin(ξ i ). In SIMDYN, the order of rotation used was roll first, then pitch, and then yaw following the convention specified by Ogilvie [32]. The angular velocity in the GCS is ω = R T ω , where ω is the angular velocity in the LCS. If we express the angular velocity in the matrix form: In this paper, . ξ means the time derivative, and ..
ξ means the second time derivative.

Governing Equations
The blended time-domain method is based on the nonlinear equations of motion. According to previous work [12], the exact equations of motion are as follows: where m is the mass, and I is the mass moment of inertia with respect to the center of gravity. The force in Equation (3) is applied at the center of gravity of the rigid body in the GCS. The moment in Equation (3) is with respect to the center of gravity of the rigid body in the LCS. The left hand side of Equation (3) are the inertia forces while the right-hand side are the external forces F (force at the center of gravity in the GCS) and moments M (the moment with respect to the origin of the GCS in the GCS): The Sandia National Laboratory and the National Renewable Energy Laboratory developed WEC-Sim specifically to model a wave energy converter [33]. Compared to WEC-Sim, SIMDYN includes an additional force: slowly varying drift forces F sv (used under the linear time-domain option). The other forces are Froude-Krylov forces (F FK ), diffraction forces (F dif ), radiation forces (F rad ), viscous forces (F vis ), hydrostatic forces (F hyd ), forces from the mooring system (F mor ), Morison element forces (F me ), and power take-off (PTO) forces (F PTO ). Equation (4) includes the power take-off term to describe the program completely. However, in the simulations (linear and blended) performed for this paper, the PTO modeling is not used. This is because the floating power system in the model tests had no PTO, which helps reduce the uncertainty of model test correlations (to examine SIMDYN as a simulation tool).
The viscous forces, radiation forces, diffraction forces, and Morison forces can be calculated in the usual way adopted by any linear time-domain program (e.g., [7]). Table 2 compares several time-domain programs that can simulate WESs on the important forces that can be calculated differently. The blended time-domain algorithm in SIMDYN is one of the tool's most advanced features, including the nonlinear inertia forces (not applied in other programs shown in Table 2) and the drift forces (not applied in most other programs except in the commercial software AQWA).

Mooring Forces/Moments
Mooring in SIMDYN is modeled using the open-source MAP++. MAP++ performs mooring analysis with quadratic modeling [27]. MAP++ ignores the inertia forces of the mooring lines and the fluid drag forces on the mooring lines [34]. Figure 2 shows how MAP++ is coupled with SIMDYN. The blended time-domain algorithm in SIMDYN is one of the tool's most advanced features, including the nonlinear inertia forces (not applied in other programs shown in Table 2) and the drift forces (not applied in most other programs except in the commercial software AQWA).

Mooring Forces/Moments
Mooring in SIMDYN is modeled using the open-source MAP++. MAP++ performs mooring analysis with quadratic modeling [27]. MAP++ ignores the inertia forces of the mooring lines and the fluid drag forces on the mooring lines [34]. Figure 2 shows how MAP++ is coupled with SIMDYN.

Slowly Varying Drift Forces/Moments
Slowly varying drift forces affect WES responses under irregular waves. Therefore, when the linear time-domain option is used, SIMDYN calculates the drift forces. In this study, the quadratic transfer functions were output from the Marine Dynamic Laboratory's hydrodynamic analysis program, MDL-HydroD (a frequency domain program developed by our research group). Details on how quadratic transfer functions (QTFs) are evaluated can be found in previous work [35,36]. Figure  3 show sub-harmonic force quadratic transfer functions for surge and pitch.  If the random incident wave is decomposed as

Slowly Varying Drift Forces/Moments
Slowly varying drift forces affect WES responses under irregular waves. Therefore, when the linear time-domain option is used, SIMDYN calculates the drift forces. In this study, the quadratic transfer functions were output from the Marine Dynamic Laboratory's hydrodynamic analysis program, MDL-HydroD (a frequency domain program developed by our research group). Details on how quadratic transfer functions (QTFs) are evaluated can be found in previous work [35,36]. Figure 3 show sub-harmonic force quadratic transfer functions for surge and pitch. The blended time-domain algorithm in SIMDYN is one of the tool's most advanced features, including the nonlinear inertia forces (not applied in other programs shown in Table 2) and the drift forces (not applied in most other programs except in the commercial software AQWA).

Mooring Forces/Moments
Mooring in SIMDYN is modeled using the open-source MAP++. MAP++ performs mooring analysis with quadratic modeling [27]. MAP++ ignores the inertia forces of the mooring lines and the fluid drag forces on the mooring lines [34]. Figure 2 shows how MAP++ is coupled with SIMDYN.

Slowly Varying Drift Forces/Moments
Slowly varying drift forces affect WES responses under irregular waves. Therefore, when the linear time-domain option is used, SIMDYN calculates the drift forces. In this study, the quadratic transfer functions were output from the Marine Dynamic Laboratory's hydrodynamic analysis program, MDL-HydroD (a frequency domain program developed by our research group). Details on how quadratic transfer functions (QTFs) are evaluated can be found in previous work [35,36]. Figure  3 show sub-harmonic force quadratic transfer functions for surge and pitch. If the random incident wave is decomposed as If the random incident wave is decomposed as The corresponding (second-order) slowly varying wave drift force are (see [7]) The sum frequency components in the drift forces are neglected in the equation because their contributions are usually much smaller than the difference frequency components. a i , a j are the amplitude of the wave components with the frequency ω i and ω j and the phase i and j .
NF is the number of frequency used to define the frequency range. P − ij are the in-phase components of the quadratic transfer function for the difference frequency ω − . Q − ij are the out-of-phase components of the quadratic transfer function for difference frequency ω − .

Viscous Forces/Moments
There are two terms related to viscous drag: F vis and F me . Either of them can be used in SIMDYN. F vis : SIMDYN can model the viscous forces using the linearized (equivalent) damping coefficient B eq : Alternatively, the viscous forces can be modeled in the quadratic form of C dv are the quadratic damping coefficients and . ξ j is the jth degree of freedom velocity of the body's center of gravity. The drag coefficients should be determined by model test correlations (e.g., free decay tests, system identification, empirical formula or data). Practically, it depends upon the user which degree of freedom they want to apply the damping to, and what are the appropriate damping values. For example, in SIMDYN: *external_damping 4 1 4E4 → a linear viscous rolling damping of 4 × 10 4 Ns/m is applied. *quadratic_damping 1 2 2E3 → a quadratic viscous surge damping of 2 × 10 3 Ns 2 /m 2 is applied.

Formulation of Nonlinear Froude-Krylov and Hydrostatic Forces
In the blended time-domain method, Froude-Krylov forces are calculated by integrating the nonlinear dynamic pressure over the instantaneous wetted surface of the rigid body [25]. In this way, the method accounts for the effects of the instantaneous body motions and the instantaneous-incident free surface. Linear incident wave potential at point (x, y, z) in the GCS due to a unidirectional irregular sea incident at counterclockwise angle β to the body's longitudinal axis is given by the following: For wave frequency component i, H i is the wave height, ω i is the frequency, k i is the wave number, h is the water depth, and NF is the total number of wave frequencies. The linear incident wave potential is not defined for the points above the mean water line (z = 0). SIMDYN employs Wheeler stretching to extrapolate the incident and hydrostatic pressure profiles to provide an expression for pressure inside the incident wave crest [37,38]. Mathematically, this can be considered as scaling the z-coordinate to compute pressure up to instantaneous free surface elevation η (measured from the calm water plane) due to the incident wave [37,38]. The vertical coordinate z is modified to zw through Wheeler stretching: Dynamic pressure p is as follows: Note that dynamic pressure − ρ 2 ∇φ I (t, x, y, zw) 2 is also nonlinear. The surface panels of the rigid body satisfying condition z ≤ η(t, x, y) form instantaneous wetted panels P B . The Froude-Krylov forces/moments are as follows: where n is the normal vector of the panel, and x is the position vector of the wetted panel centroid (in the GCS). SIMDYN nonlinearly-integrates the results of a pre-processed potential problem solved under linear conditions Like the Froude-Krylov forces, in the blended time-domain method, the hydrostatic forces and moments are also calculated by integrating the hydrostatic pressure over the instantaneous wetted surface area [23]: where −mgk is the rigid body weight vector, and P B are the instantaneous wetted panels satisfying z ≤ η(t, x, y). Note that the hydrostatic pressure, −ρgz, is calculated using the actual vertical coordinate z instead of zw. In this way, Equations (12) and (13) satisfy the dynamic free-surface boundary condition over the incident wave profile.

Nonlinear Effects of Froude=Krylov and Hydrostatic Forces
This section presents the nonlinear effects through the forced motion tests. The forced motion tests are important because: (1) The forced motion tests show the different forcing corresponding to the same (large) motions.
That can indicate the motion (as the final result) differences in an implicit way; (2) Any simulation tool comes with limitations. The forced motion test avoids the problem by disabling modules not robust enough and not very relevant (for example, the mooring module is not the focus of this study); (3) The forced motion test is a control-variable test. It eliminates the effects of other forces, which makes the effect from each force component clearer.
The geometry used in these tests is a floating power system (its details are presented in Section 5) as plotted in Figure 4. The forced motion tests refer to a series of SIMDYN simulations that use the specified (forced) motion time series to study a force component (e.g., Froude-Krylov). In these tests, Note that the motions in the forced motion tests (except 1) were not simulated but were specified (as input time series). This is quite a harsh condition (large wave and motions) for the FPS; therefore, nonlinear effects could be observed clearly. Figure 5 compares the Froude-Krylov and hydrostatic heave forces. The hydrostatic heave forces from the two simulations are very close, while the Froude-Krylov heave forces show slight differences.   Note that the motions in the forced motion tests (except M1) were not simulated but were specified (as input time series). This is quite a harsh condition (large wave and motions) for the FPS; therefore, nonlinear effects could be observed clearly.  Note that the motions in the forced motion tests (except 1) were not simulated but were specified (as input time series). This is quite a harsh condition (large wave and motions) for the FPS; therefore, nonlinear effects could be observed clearly.       3) The translation motions of the structure are not captured in the linear method, so the wetted surfaces and the corresponding pressures are different. Note that the surge and sway influence the relative phase of the incident wave, so they also contribute to the differences.

Derivations of the Nonlinear Inertia Forces
For an irregular wave, the radiation forces are as follows [39]:   3) The translation motions of the structure are not captured in the linear method, so the wetted surfaces and the corresponding pressures are different. Note that the surge and sway influence the relative phase of the incident wave, so they also contribute to the differences.

Derivations of the Nonlinear Inertia Forces
For an irregular wave, the radiation forces are as follows [39]: Note that the surge and sway influence the relative phase of the incident wave, so they also contribute to the differences.

Derivations of the Nonlinear Inertia Forces
For an irregular wave, the radiation forces are as follows [39]: where [A(∞)] is the 6 × 6 added mass matrix at infinite frequency. We can write the following:  (14) should be moved to the left side of Equation (3). The equations of motion become the following: Angular velocities ω and ω defined in Equation (2) can be expressed in matrix form as shown below: Differentiating Equation (17) .
Substituting Equations (17) and (18) into Equation (3) and rearranging to keep only the terms containing acceleration on the left side results in the following: Substituting Equations (17) to (19) into Equation (16) yields the following: Equations (19) and (20) are rearranged into the matrix form. Let v denote a generic vector. Cross multiplication x G × v can be written in the matrix form: Similarly, cross multiplication v × (x G − ξ) can be written in the matrix form: Using Equation (22), (19) can be written as follows: Substituting Equations (21) and (22) into Equation (20) and rearranging to keep only the terms containing acceleration on the left side results in the following: In the matrix form, Equations (23) and (24) can be written as follows: Equation (25) is the nonlinear equation of motion used in SIMDYN; it can be solved using the fourth-order Runge-Kutta method [12].

Nonlinear Effects of the Inertia Force
This subsection presents the nonlinear effects of the inertia force (see also [40][41][42]) through forced motion tests on the FPS geometry shown in Figure 4. In Equation (16), all the nonlinear terms are related to rotations. Note that Equation (16) includes the inertia forces due to the added masses on the left sides. Equation (16) can be rewritten as follows: Using Equation (26), nonlinear inertia forces F and moments M can be calculated. Note that in the linear time-domain method, the motions and rotations are assumed small: where x G are the coordinates of the center of gravity in the LCS, which are constant. Taking these assumptions to remove the nonlinear terms in Equation (25), we obtain the following: Note that the terms containing L 1 x G and L 2 x G convert the momentum of inertia from the center of gravity to the origin of the LCS in the linear model, (0, 0, 0). This is the typical linear equation of motion. From Equation (28), the linear inertia forces and moments can be calculated.
Using Equations (26) and (28), the inertia forces in the forced motion tests (i.e., the time series of the motions, velocities, and accelerations are given) can be obtained. Note that the moments are given in the GCS from Equations (26) and (28); therefore, the inertia moments from the blended time-domain method should be transformed to the instantaneous origin of the body coordinate system, (ξ 1 , ξ 2 , ξ 3 ), to be compared with those from the linear time-domain method.
In the forced motion tests, the translational motions were set to zero to reveal the nonlinear effects due to the rotation. The first set of forced motions of the FPS are quite mild, referred to as M2 (see Figure 8). M2 is specified to show the inertia force differences (between the two options), so the relative phases of the rotations are not necessarily realistic. system, ( 1 , 2 , 3 ), to be compared with those from the linear time-domain method.
In the forced motion tests, the translational motions were set to zero to reveal the nonlinear effects due to the rotation. The first set of forced motions of the FPS are quite mild, referred to as 2 (see Figure 8). 2 is specified to show the inertia force differences (between the two options), so the relative phases of the rotations are not necessarily realistic.  The inertia forces corresponding to 10 times the reference motions, 10 × M2, are plotted in Figure 9. The surge and sway forces from the two methods are slightly different, indicating that the nonlinearity due to rotation is not significant at these rotation amplitudes. The inertia forces corresponding to 10 times the reference motions, 10 × 2, are plotted in Figure 9. The surge and sway forces from the two methods are slightly different, indicating that the nonlinearity due to rotation is not significant at these rotation amplitudes.
The heave force from the blended model looks very different from the linear model but bear in mind that when the actual heave motion is in place, the contribution from the heave motion itself is much greater (> 1,000 times) than that from the rotations. Under the forced rotations, 10 × 2, the differences in roll and pitch moments between the linear model and the blended model are noticeable (see Figure 10). The significant difference of the yaw moments indicates that for yaw, in this case, the linear model is no longer valid. The heave force from the blended model looks very different from the linear model but bear in mind that when the actual heave motion is in place, the contribution from the heave motion itself is much greater (>1000 times) than that from the rotations.
Under the forced rotations, 10 × M2, the differences in roll and pitch moments between the linear model and the blended model are noticeable (see Figure 10). The significant difference of the yaw moments indicates that for yaw, in this case, the linear model is no longer valid.
The blended results contain higher (mostly double) frequency components in sway, heave, yaw, and roll. They are the superharmonic components contributed by the nonlinear terms in the equation of motion. For example, in a case with zero yaw accelerations (an axisymmetric buoy), the combination of roll and pitch rotations will provide a nonlinear contribution to the yaw inertial forces using the nonlinear model. Under the forced rotations, 10 × 2, the differences in roll and pitch moments between the linear model and the blended model are noticeable (see Figure 10). The significant difference of the yaw moments indicates that for yaw, in this case, the linear model is no longer valid. The blended results contain higher (mostly double) frequency components in sway, heave, yaw, and roll. They are the superharmonic components contributed by the nonlinear terms in the equation of motion. For example, in a case with zero yaw accelerations (an axisymmetric buoy), the

Model Test Setup and Input Information
The blended time-domain program SIMDYN (manufacturer, city and country) was examined using correlations with a series of 1:25-scale model tests [31] performed by Beaufort Research/HMRC in Ireland. In these model tests, the incident wave headings included 0 • , 30 • , and 60 • . The 0 • wave heading test was selected as it consists of more sea states. The 0 • wave test configuration is shown in Figure 11. combination of roll and pitch rotations will provide a nonlinear contribution to the yaw inertial forces using the nonlinear model.

Model Test Setup and Input Information
The blended time-domain program SIMDYN (manufacturer, city and country)was examined using correlations with a series of 1:25-scale model tests [31] performed by Beaufort Research/HMRC in Ireland. In these model tests, the incident wave headings included 0°, 30°, and 60°. The 0° wave heading test was selected as it consists of more sea states. The 0° wave test configuration is shown in Figure 11. The 0° wave heading tests consisted of 43 regular wave cases. There were 20 random (Bretschneider spectrum) wave cases with 0° wave heading, but some cases were skipped due to indications of incomplete measurements (over the time span of 1,748 s) or unexpectedly large roll and yaw measurements (which should not be the case, as the configuration is symmetric about the xaxis). The unexpectedly large roll and yaw measurements may be partially attributed to parametric instabilities, which many wave-activated WECs suffer from (see [42][43][44][45][46][47][48][49][50]). These dynamic instabilities deserve systematic studies in the future. Consequently, 12 irregular wave cases were compared. The correlated sea states are plotted in Figure 12. The 0 • wave heading tests consisted of 43 regular wave cases. There were 20 random (Bretschneider spectrum) wave cases with 0 • wave heading, but some cases were skipped due to indications of incomplete measurements (over the time span of 1748 s) or unexpectedly large roll and yaw measurements (which should not be the case, as the configuration is symmetric about the x-axis).
The unexpectedly large roll and yaw measurements may be partially attributed to parametric instabilities, which many wave-activated WECs suffer from (see [42][43][44][45][46][47][48][49][50]). These dynamic instabilities deserve systematic studies in the future. Consequently, 12 irregular wave cases were compared. The correlated sea states are plotted in Figure 12. The 0° wave heading tests consisted of 43 regular wave cases. There were 20 random (Bretschneider spectrum) wave cases with 0° wave heading, but some cases were skipped due to indications of incomplete measurements (over the time span of 1,748 s) or unexpectedly large roll and yaw measurements (which should not be the case, as the configuration is symmetric about the xaxis). The unexpectedly large roll and yaw measurements may be partially attributed to parametric instabilities, which many wave-activated WECs suffer from (see [42][43][44][45][46][47][48][49][50]). These dynamic instabilities deserve systematic studies in the future. Consequently, 12 irregular wave cases were compared. The correlated sea states are plotted in Figure 12. MDL-HydroD was used to perform the frequency domain analysis. The added masses, the radiation dampings, the diffraction forces, and the QTFs necessary for the blended time-domain analysis were obtained. More details about MDL-HydroD can be found in previous work [51][52][53][54]. Figure 13a shows the 1:25-scale FPS model used in the wave basin. The FPS panel model in Figure 13b was used in SIMDYN (for the integration of instantaneous pressure). The FPS has an octagon cross-section with a decreasing cross-section area below the calm water line. It consisted of 8,919 panels (see Figure 7), which proved fine enough to capture the nonlinear effects. The wetted surface was not re-meshed during the simulations. A panel is either fully submerged or fully emerged (by determining whether the instantaneous centroid of the panel is above/below the free surface). The FPS is positioned by a 3-leg catenary mooring system with a 120-degree azimuth difference between each mooring leg. MDL-HydroD was used to perform the frequency domain analysis. The added masses, the radiation dampings, the diffraction forces, and the QTFs necessary for the blended time-domain analysis were obtained. More details about MDL-HydroD can be found in previous work [51][52][53][54]. Figure 13a shows the 1:25-scale FPS model used in the wave basin. The FPS panel model in Figure 13b was used in SIMDYN (for the integration of instantaneous pressure). The FPS has an octagon cross-section with a decreasing cross-section area below the calm water line. It consisted of 8,919 panels (see Figure 7), which proved fine enough to capture the nonlinear effects. The wetted surface was not re-meshed during the simulations. A panel is either fully submerged or fully emerged (by determining whether the instantaneous centroid of the panel is above/below the free surface). The FPS is positioned by a 3-leg catenary mooring system with a 120-degree azimuth difference between each mooring leg.  Table 3 lists the inputs used in the model test correlations. VCG is the vertical center of gravity and it is measured from the calm water plane (instead of from the bottom of the body). Kxx, Kyy, and Kzz are the gyration radius, around the center of gravity. EA is the element axial stiffness.   Table 3 lists the inputs used in the model test correlations. VCG is the vertical center of gravity and it is measured from the calm water plane (instead of from the bottom of the body). Kxx, Kyy, and Kzz are the gyration radius, around the center of gravity. EA is the element axial stiffness.

Correlations of the Regular Wave Cases
The FPS geometry and mooring system configuration was symmetric about the x-axis, so relative small sway, roll, and yaw motions were expected. The dominating surge, heave, and pitch motions from the time-domain (the linear and the blended) analysis were compared with the model tests.
Of the total 43 cases, 38 were successfully analyzed with reasonable time-series results. Typical cases are of the form similar to sinusoidal time history (as plotted in Figures 14 and 15). From the plots, we can find both the linear model and the blended model to yield reasonable comparisons with the model tests (for surge and pitch). Motion amplitudes from the blended option look closer (than the linear model) as compared to the model test. The heave plots were omitted because the simulated heave motions were very close to the model test (no noticeable difference).
For the motion time series (longer than 10 min), their standard deviations represent the motion magnitudes. Therefore, the simulation errors are as follows: Figure 16 compares the errors for surge, heave, and pitch. The range [−10%, 10%] bounds 95% of the error dots, meaning that both the linear model and the blended model were found to yield an acceptable error level. While the heave motion errors from the two models are similar, the blended model was found to yield lower surge and pitch motion errors compared with the linear model. The statistics in Table 4 also reflect this. As listed in Table 4, the mean errors of the 38 regular wave cases from the blended method are consistently lower than the mean errors from the linear method for surge, heave, and pitch.
In the case that the large positive errors canceled out the large negative errors (meaning the mean errors did not represent the actual error levels), mean absolute errors are listed in Table 4 to reflect the level of error in another way. The mean absolute errors from the blended model also were shown to be lower than the mean absolute errors from the linear model, except for heave (which was very close). Table 4 indicates that the blended time-domain method has a considerable improvement in accuracy. The advantage of the blended method can be credited to its capability to account for the nonlinearity in Froude-Krylov, hydrostatic, and inertia forces.
Depending on the computer's CPU capability, the absolute calculation time of the two methods is subject to change. However, the relative speed is meaningful. For a simulation of 3 min in real-time, using a step of 0.1 s, there should be 1800 steps (number of steps should be the criteria here as the step itself is flexible). Under this setting, the blended method spends 19 times the simulation time of that spent by the linear method. The FPS geometry and mooring system configuration was symmetric about the x-axis, so relative small sway, roll, and yaw motions were expected. The dominating surge, heave, and pitch motions from the time-domain (the linear and the blended) analysis were compared with the model tests.
Of the total 43 cases, 38 were successfully analyzed with reasonable time-series results. Typical cases are of the form similar to sinusoidal time history (as plotted in Figures 14 and 15). From the plots, we can find both the linear model and the blended model to yield reasonable comparisons with the model tests (for surge and pitch). Motion amplitudes from the blended option look closer (than the linear model) as compared to the model test. The heave plots were omitted because the simulated heave motions were very close to the model test (no noticeable difference). For the motion time series (longer than 10 min), their standard deviations represent the motion magnitudes. Therefore, the simulation errors are as follows: Figure 16 compares the errors for surge, heave, and pitch. The range [−10%, 10%] bounds 95% of the error dots, meaning that both the linear model and the blended model were found to yield an acceptable error level. While the heave motion errors from the two models are similar, the blended model was found to yield lower surge and pitch motion errors compared with the linear model. The statistics in Table 4 also reflect this. The FPS geometry and mooring system configuration was symmetric about the x-axis, so relative small sway, roll, and yaw motions were expected. The dominating surge, heave, and pitch motions from the time-domain (the linear and the blended) analysis were compared with the model tests.
Of the total 43 cases, 38 were successfully analyzed with reasonable time-series results. Typical cases are of the form similar to sinusoidal time history (as plotted in Figures 14 and 15). From the plots, we can find both the linear model and the blended model to yield reasonable comparisons with the model tests (for surge and pitch). Motion amplitudes from the blended option look closer (than the linear model) as compared to the model test. The heave plots were omitted because the simulated heave motions were very close to the model test (no noticeable difference). For the motion time series (longer than 10 min), their standard deviations represent the motion magnitudes. Therefore, the simulation errors are as follows: Figure 16 compares the errors for surge, heave, and pitch. The range [−10%, 10%] bounds 95% of the error dots, meaning that both the linear model and the blended model were found to yield an acceptable error level. While the heave motion errors from the two models are similar, the blended model was found to yield lower surge and pitch motion errors compared with the linear model. The statistics in Table 4 also reflect this. Some factors will influence the relative speed (e.g., the number of panels, the inclusion of the drift forces, ramp time setting). While the algorithm of the linear SIMDYN is quite mature, the algorithm of the blended SIMDYN has the room for more efficient optimization. In general, a fair estimation is that the blended time-domain method is about 10 times (one order of magnitude) slower than the linear time-domain method. Please remember that the next fidelity model (nonlinear time-domain method) is about 10 2 times (two orders of magnitude) slower than the linear time-domain method (see [5]).
Although the model test correlations in this paper show insignificant differences using the two modeling approaches. However, this happens when the linear method results are already fairly close to the model test results, leaving small room for improvements. We believe that under more severe sea states, or for better-tuned devices, the differences will be more significant.  As listed in Table 4, the mean errors of the 38 regular wave cases from the blended method are consistently lower than the mean errors from the linear method for surge, heave, and pitch.
In the case that the large positive errors canceled out the large negative errors (meaning the mean errors did not represent the actual error levels), mean absolute errors are listed in Table 4 to reflect the level of error in another way. The mean absolute errors from the blended model also were shown to be lower than the mean absolute errors from the linear model, except for heave (which was very close). Table 4 indicates that the blended time-domain method has a considerable improvement in accuracy. The advantage of the blended method can be credited to its capability to account for the nonlinearity in Froude-Krylov, hydrostatic, and inertia forces.
Depending on the computer's CPU capability, the absolute calculation time of the two methods is subject to change. However, the relative speed is meaningful. For a simulation of 3 min in realtime, using a step of 0.1s, there should be 1800 steps (number of steps should be the criteria here as the step itself is flexible). Under this setting, the blended method spends 19 times the simulation time of that spent by the linear method. For all the regular wave cases, the general trend of the response amplitude ratio (blended/linear) (including surge, heave, and pitch) is plotted in Figure 17. Some factors will influence the relative speed (e.g., the number of panels, the inclusion of the drift forces, ramp time setting). While the algorithm of the linear SIMDYN is quite mature, the algorithm of the blended SIMDYN has the room for more efficient optimization. In general, a fair estimation is that the blended time-domain method is about 10 times (one order of magnitude) slower than the linear time-domain method. Please remember that the next fidelity model (nonlinear timedomain method) is about 10 2 times (two orders of magnitude) slower than the linear time-domain method (see [5]).
Although the model test correlations in this paper show insignificant differences using the two modeling approaches. However, this happens when the linear method results are already fairly close to the model test results, leaving small room for improvements. We believe that under more severe sea states, or for better-tuned devices, the differences will be more significant.
For all the regular wave cases, the general trend of the response amplitude ratio (blended/linear) (including surge, heave, and pitch) is plotted in Figure 17. From Figure 17a, we can find that (1) when the wave amplitude and motion amplitude are small, the blended method and the linear method are statistically very close and (2) as the wave and motion amplitude increase, the results from the linear method exceed the results from the blended method. The gradual loss of credibility of the linear assumptions with the increase of wave and motion amplitude leads to overprediction of the motions from the linear method. In addition, Figure 17b indicates that the accuracy of the blend method is quite consistent as the wave height changes.

Correlations of Irregular Wave Cases
Similarly, the time-domain results of the selected irregular cases were compared with the model From Figure 17a, we can find that (1) when the wave amplitude and motion amplitude are small, the blended method and the linear method are statistically very close and (2) as the wave and motion amplitude increase, the results from the linear method exceed the results from the blended method. The gradual loss of credibility of the linear assumptions with the increase of wave and motion amplitude leads to overprediction of the motions from the linear method. In addition, Figure 17b indicates that the accuracy of the blend method is quite consistent as the wave height changes.

Correlations of Irregular Wave Cases
Similarly, the time-domain results of the selected irregular cases were compared with the model tests. Typical time series are plotted in Figures 18 and 19.  Figure 17a, we can find that (1) when the wave amplitude and motion amplitude are small, the blended method and the linear method are statistically very close and (2) as the wave and motion amplitude increase, the results from the linear method exceed the results from the blended method. The gradual loss of credibility of the linear assumptions with the increase of wave and motion amplitude leads to overprediction of the motions from the linear method. In addition, Figure 17b indicates that the accuracy of the blend method is quite consistent as the wave height changes.

Correlations of Irregular Wave Cases
Similarly, the time-domain results of the selected irregular cases were compared with the model tests. Typical time series are plotted in Figures 18 and 19. From Figures 18 and 19, we can observe reasonable correlations between the simulations and the model tests. However, as the irregular wave cases were more complex, the close matches of the time series (similar to what we obtained in Figures 14 and 15 in the regular wave cases) were not available in the irregular wave cases. However, the typical heave motions (see Figure 20) correlated better (with the model tests) than the typical surge and pitch motions shown in Figures 18 and 19. From Figures 18 and 19, we can observe reasonable correlations between the simulations and the model tests. However, as the irregular wave cases were more complex, the close matches of the time series (similar to what we obtained in Figures 14 and 15 in the regular wave cases) were not available in the irregular wave cases. However, the typical heave motions (see Figure 20) correlated better (with the model tests) than the typical surge and pitch motions shown in Figures 18 and 19.
The motion errors from the two methods are listed in Table 5. The errors from the blended time-domain model are also lower.

Discussion
The blended, weakly nonlinear time-domain method was implemented to predict WES motions. This method calculates important external forces by directly integrating the instantaneous pressure on the wetted panels of the floating structures, while the remaining external forces are processed in ways adopted by traditional linear time-domain programs (e.g., OrcaFlex). The method also accounts for nonlinearity in the equations of motion (i.e., inertia forces) due to rotations of the rigid body. The simulation results from the blended time-domain method were compared with the experimental results measured under the same incident waves. Correlations with the model test results indicate that the blended time-domain method is more accurate than the widely used linear time-domain method as it reduces the error of the motion responses. Moreover, the method is significantly faster than possibly more exact methods such as fully nonlinear time-domain methods or CFD.
This study also discussed the nonlinear effects of Froude-Krylov, hydrostatic, and inertia forces. The (forced) motion time series were input to both the linear model and the blended model to compare the corresponding Froude-Krylov, hydrostatic, and inertia forces. As in the forced motion tests, the method used (linear or blended) was the only variable, so the differences between the two methods under certain motion and wave amplitude could be compared. In general, the gap between the blended method and the linear method widens as the motion amplitude and wave amplitude increase. While the model test results generally compare well with both the linear and blended time-domain methods under the tested sea states, we reasonably deduce from the forced motion test results that with larger wave amplitudes, the advantages of the blended method may increase substantially.
In this work, a floating power system was used in the model test correlations. This geometry (selected from a very limited suitable model test data accessible to us), though not a real wave energy converter, is of dimensions and geometry similar to the typical point absorber type wave energy converter. At the stage of motion predictions through numerical modeling, especially under relatively large sea states, the WEC is usually in survival mode (when its PTO is deactivated). Therefore, for this study, the floating power system is almost equivalent to a point absorber type WEC. From this study, the blended method was found to be more accurate for the point absorber type WEC that usually undergoes large (relative to their dimensions) wave and motion amplitudes.
At this stage, the focus is to make sure that the motion predictions from the numerical modeling under different sea states are accurate. Even with the real wave energy converter, for this study (examining the modeling method under different sea states), it is preferable that the power take-off is deactivated. That makes the floating power system similar to the point absorber in essence. Actually, the improvements of SIMDYN and the confidences/benchmarks gained in this study are the foundation for any further study. In the future study, the power take-off forces will be activated in our model to simulate the point absorber in their working mode. That should produce better results for predicting the WEC's electrical outputs (e.g., voltage, current and power). The refined motion predictions will not only improve the power take-off predictions for the WEC's working mode but will also improve (more greatly) the load predictions for the point absorber type WEC's survival mode. The more accurate motion predictions from our improved model can be directly used to determine the local accelerations (will determine the loads and pressures) as well as the mooring line tensions and offsets, which are the most critical quantities for designing WECs. Loads under extreme sea states are the key to overcoming risks and reducing costs of the point absorber type WECs, especially before their practical applications become more extensive.

Conclusions
(1) A simulation program, SIMDYN, was developed using both the linear and the blended time-domain methods to predict the motions (six degrees of freedom) of the wave energy converters. (2) The nonlinear Froude-Krylov and hydrostatic forces were implemented in SIMDYN.
The (dominant) external forces can be calculated more accurately without requiring significantly more calculation time. Funding: This research received no external funding.