Time-Domain Implementation and Analyses of Multi-Motion Modes of Floating Structures

<jats:p>The study of wave-structure interactions involving nonlinear forces would often make use of the popular hybrid frequency–time domain method. In the hybrid method, the frequency-domain analysis could firstly provide the reliable and accurate dynamic parameters and responses; then these parameters and responses are transformed to the parameters to establishing the basic time-domain equation. Additionally, with the addition of the required linear and nonlinear forces, the time-domain analysis can be used to solve for the practical problems. However, the transformation from the frequency domain to the time domain is not straightforward, and the implementation of the time-domain equation could become increasingly complicated when different modes of motion are coupled. This research presents a systematic introduction on how to implement the time-domain analysis for floating structures, including the parameter transformations from the frequency domain to the time domain, and the methods for calculating and approximating the impulse functions and the fluid-memory effects, with special attention being paid to the coupling terms among the different motion modes, and the correctness of the time-domain-equation implementation. The main purpose of this article is to provide relevant information for those who wish to build their own time-domain analyses with the open-source hydrodynamic analysis packages, although commercial packages are available for time-domain analyses.</jats:p>

section_content: Frequency-domain analysis is a widely used method for wave-structure interactions and, in the method, the dynamic system is assumed to be fully linear: the governing equation is the linear Laplace equation (based on the incompressible potential flows); the body and seabed boundary conditions are linear; and the free-surface condition can be linearized for practical applications. Although the assumptions are strict, the frequencydomain analysis could provide reliable and accurate assessments for the hydrodynamic parameters and responses, such as the added mass, the radiation-damping coefficients, the wave-excitation forces, and the response amplitude operators (RAOs) [1] [2] [3] [4] [5], and, in some cases, it can be even extended to hydro-elastic analysis [6] and the elastic wave-energy converters (WECs) [7, 8]. For its applications in wave-energy converters, the conventional frequency-domain analysis may provide an accurate calculation for the resonance frequency/period of the device, which is generally regarded as the most important parameter for wave-energy converters, since most wave-energy converters would be ideally optimized to have resonance with the wave for efficient energy extraction from the waves [9, 10]. In addition, the frequency-domain analysis is also a powerful tool for optimizing WECs and power take-off (PTO) to maximize wave-energy production [11, 12]. 

The assumption of fully linear dynamic systems in the frequency-domain analysis could limit its direct applications in WECs since, for some practical WECs, nonlinear PTO systems may be employed, such as hydraulic PTOs [13] [14] [15] [16] or impulse turbines [17] [18] [19]. Additionally, when a control is applied to improve wave-energy conversion [20] [21] [22] [23] [24], the corresponding dynamic system in wave-energy conversion would be identically nonlinear, even if the PTO is linear. For a nonlinear dynamic system, the frequency-domain analysis can not be used directly; instead, the most applied method is the so-called hybrid frequency-time domain method [9, [25] [26] [27]. In the hybrid method, the frequencydomain analysis is firstly used to obtain the relevant hydrodynamic parameters (that is, the frequency-dependent parameters), and these parameters and the frequency-domain equation are then transformed to the corresponding time-domain parameters and the Cummins' time-domain equation [28]. To realize the transformation, the Ogilvie relation [29] is used for transforming the relevant frequency-dependent parameters to the parameters in the time domain, such as the impulse function, the fluid-memory effect, and the added mass at infinite frequency. For real problems, the transformed Cummins' time-domain equation can be modified by adding the required forces (linear and/or nonlinear) for the actual problems. 

However, such a transformation from the frequency domain to the time domain is not straightforward. It involves the transformation from the frequency-dependent radiationdamping coefficients to the impulse functions and the subsequent fluid-memory effects to represent the fluid-damping effects on the floating structures, and the transformation from the frequency-dependent added mass to the added mass at infinite frequency. Additionally, to calculate the memory effects more effectively, the impulse functions are frequently approximated, so are the fluid memory effects [30] [31] [32]. 

The transformations and approximations would become increasingly complicated when different modes of motion are coupled together [12] (compared to those time-domain analyses of the single-motion mode, such as on heave [9], or on pitch [26, 27] ). However, in reality, the multiple motion modes must be analyzed for some floating structures and WECs, for instance, the multi-axis WECs [16]. The current research systematically presents and explains how the correct transformations can be made, how the memory effects are approximated, and how the multiple coupled motions are solved. Moreover, a validation method is introduced, which would allow the researchers to check whether the numerical modeling implementation for the time-domain modeling is correct. When this important step is completed, the researchers can be more confident to add the required forces to the dynamic equation for the practical problems. 

The rest of the paper is arranged as follows: Section 2 presents a brief introduction of the multi-axis TALOS WECs; Section 3 introduces the frequency-domain equation and responses, paying attention to the coupled-motion modes of the TALOS WECs; Section 4 establishes the time-domain analysis, introducing how to transform the frequency-domain equation into the time-domain equation; Section 5 presents the methods for approximating the impulse function and the memory effect of the fluid on the structure; Section 6 discusses the implementation and validation method of the time-domain analysis; and, finally, in Section 7, the conclusions are presented regarding the implementation of the time-domain analyses for floating structures.
section_name: Introduction

section_content: TALOS is a point absorber-type WEC, which was proposed and initially developed at Lancaster University, UK [16]. The general shape for TALOS I is an octagonal-shaped floating structure (see the panels of the wet surfaces of the TALOS in Figure 1a ), and it was the target floating structure in this research for studying the wave-structure interaction and the implementation of the time-domain analysis for the multiple-motion modes of a floating structure. 

between the structure and the heavy ball to drive the PTOs (including the springs dampers) to convert the wave energy into useful mechanical energy. As a result of an arrangement, the TALOS WEC had the advantage of being a multi-axis wave-en converter (the other advantages of the TALOS WEC include the fully enclosed PTO in structure, thus affording a good survivability for the WEC device). Recently, the TALOS multi-axis wave-energy-converter technology (the proj short name: 'NHP-WEC') was selected as one of eight projects supported by UK Rese and Innovation (UKRI) [33], which aim to unlock the potentials of marine wave ene The main focuses of the NHP-WEC project are to optimize the TALOS WEC and the vant PTO system, as well as to advance WEC technology by developing the essentia vice control and monitoring systems that are integrated with high-fidelity sea-state casting using the SmartWave toolset [34]. The current work is part of the ongoing rese of the NHP-WEC project, an important step towards the implementation of the nume model for the multi-axis TALOS wave-energy technology.
section_name: TALOS WEC

section_name: Frequency-Domain Governing Equation and Responses

section_content: Frequency-domain analysis is a very popular method for the wave-structure int tions [35] [36] [37], in which the governing equation for wave-structure interactions is so in the frequency domain, with the relevant hydrodynamic parameters (added mass, ation-damping coefficient, and wave excitation) being assessed using the panel met such as the open source, HAMS (Hydrodynamic Analysis of Marine Structure [38] ), the commercial package, WAMIT [39]. The panel method for wave-structure interac is based on the potential flow theory and solves for the problems of wave-structure i actions under the assumption of the dynamic system being fully linear. These assu tions seem very strict, but the panel method is in fact very reliable and accurate for a l range of practical problems, while, for many other problems of wave-structure int tions, it could still provide the required hydrodynamic parameters reliably for fur The TALOS WEC is a wave-energy converter with a multi-axis PTO system, which consists of a heavy ball enclosed in the wave-energy converter, supported with the springs and dampers, see the PTO system in TALOS II from the cut-off in Figure 1b, and the multiaxis PTO test rig (Figure 1c ). This multi-axis wave-energy converter converts wave energy from the multiple-motion modes of the structure in waves for a more efficient wave-energy extraction. In converting the wave energy, the heavy ball remained relatively stationary while the structure motions under the wave excitation induced the relative motion between the structure and the heavy ball to drive the PTOs (including the springs and dampers) to convert the wave energy into useful mechanical energy. As a result of such an arrangement, the TALOS WEC had the advantage of being a multi-axis wave-energy converter (the other advantages of the TALOS WEC include the fully enclosed PTO in the structure, thus affording a good survivability for the WEC device). 

Recently, the TALOS multi-axis wave-energy-converter technology (the project's short name: 'NHP-WEC') was selected as one of eight projects supported by UK Research and Innovation (UKRI) [33], which aim to unlock the potentials of marine wave energy. The main focuses of the NHP-WEC project are to optimize the TALOS WEC and the relevant PTO system, as well as to advance WEC technology by developing the essential device control and monitoring systems that are integrated with high-fidelity sea-state forecasting using the SmartWave toolset [34]. The current work is part of the ongoing research of the NHP-WEC project, an important step towards the implementation of the numerical model for the multi-axis TALOS wave-energy technology.
section_name: Frequency-Domain Governing Equation

section_content: Frequency-domain analysis is a very popular method for the wave-structure interactions [35] [36] [37], in which the governing equation for wave-structure interactions is solved in the frequency domain, with the relevant hydrodynamic parameters (added mass, radiationdamping coefficient, and wave excitation) being assessed using the panel method, such as the open source, HAMS (Hydrodynamic Analysis of Marine Structure [38] ), and the commercial package, WAMIT [39]. The panel method for wave-structure interaction is based on the potential flow theory and solves for the problems of wave-structure interactions under the assumption of the dynamic system being fully linear. These assumptions seem very strict, but the panel method is in fact very reliable and accurate for a large range of practical problems, while, for many other problems of wave-structure interactions, it could still provide the required hydrodynamic parameters reliably for further analyses, including the time-domain analyses as the time-domain implementation in this research. 

The main purpose here is to implement time-domain modeling based on the predictions of the open-source panel method, HAMS. However, for the purpose of basic validations, the commercial code, WAMIT (widely regarded as an industrial standard for studying wave-structure interactions), was employed for the purpose of comparisons (more detailed comparisons between HAMS and WAMIT can be found in [40] ). Here, the comparisons between HAMS and WAMIT were made to simply check the correctness of the relevant frequency-dependent parameters and of the corresponding impulse functions. 

Following WAMIT [39], the frequency-domain governing equation of 6-degrees-offreedom (DOFs) motions of a rigid structure is presented in a form of a mass-spring-damper system, as

where ω is the circular frequency of the wave excitation, and the parameters with the variable ω represent their frequency dependency. 

M jk, M E jk, A jk (ω) (j, k = 1, 2,. . ., 6) are the structure and external and added-mass matrices, and the former two must be specified for the numerical modeling, while the latter can be assessed using the panel method. In some wave-energy converters, the external-mass matrix M E jk may be from the generic PTO, see Babarit [41]. 

) are the radiation and external-damping coefficients, with the former being assessed using the panel method, while the latter must be specified in the numerical modeling; the external-damping coefficient B E jk may be from the generic PTO as above or from other linear dampings. 

C jk, C E jk (j, k = 1, 2,. . ., 6) are the structure hydrostatic and external restoring coefficients (both must be specified or calculated). The definition of the hydrostatic restoring coefficients C jk can be found in the WAMIT manual [39] ; the external restoring coefficient C E jk may be from the PTO or from the linearized mooring system. F j (ω) (j = 1, 2,. . ., 6) is the frequency-dependent complex amplitude of the wave excitation, which can be calculated using the panel method. 

ξ k (ω) (k = 1, 2,. . ., 6) corresponds to the motions of surge, sway, heave, roll, pitch, and yaw, respectively, as the frequency-dependent complex amplitudes of motions for the floating structure. Additionally, these are solved from the above governing equation and, as a convention, the response amplitude operator (RAO) is a more useful expression, defined as

where A is the wave amplitude (here, the wave amplitude A is without a subscript or superscript). Obviously, in the wave of a unit amplitude, the frequency-dependent ξ k itself is the RAO. 

In the conventional plots, the module of the RAO may be seen more often, which is defined as

From the terms regarding the added mass, the radiation damping and the restoring coefficients in Equation (1), we can observe that the motions of a free-floating structure may couple together through the cross-coupling coefficients. The motion couplings can be either through the wave radiation, that is, through the cross-coupling added-mass and radiation-damping coefficients, such as A 15 (as the surge-pitch coupling), B 24 (the sway-roll coupling), if these coupling coefficients are not zero, nor very small when compared to other terms, say, A 11 and B 22 ; or, through the coupling hydrostatic restoring coefficients, for instance, the coupled restoring coefficients between heave and roll (pitch), as C 34 or C 35. It should be noted that some of the coupled motions are inherently existing, e. g., surge-pitch coupling and sway-roll coupling, while some other coupled motions may exist or not, very much depending on the shape of the floating structure, for instance, a structure with a symmetry about x axis may lead to the de-coupling of the heave and pitch motions. 

In the following subsections, the hydrodynamic parameters and responses are presented and discussed for the TALOS floating structure, and the comparisons between HAMS and WAMIT aim to confirm that the open-source panel method, HAMS, is reliable and accurate.
section_name: Frequency-Domain Governing Equation

section_content: The conventional added mass of the TALOS wave-energy converter is presented in Figure 2, including the added mass for surge (A 11 ), sway (A 22 ), heave (A 33 ), roll (A 44 ), pitch (A 55 ), and yaw (A 66 ). We can observe from the comparisons between HAMS and WAMIT that they are almost identical, and some differences (spikes) can be observed in the so-called irregular frequencies in the short waves, when the wave period, T, is less than 3 s. The reason for the irregular frequencies may relate to the panel size: as a rule of thumb, the maximal panel size should be less than 1/8 of the wave length (see [37] ). cients, for instance, the coupled restoring coefficients between heave and roll (pitch), as 34 or 35. It should be noted that some of the coupled motions are inherently existing, e. g., surge-pitch coupling and sway-roll coupling, while some other coupled motions may exist or not, very much depending on the shape of the floating structure, for instance, a structure with a symmetry about x axis may lead to the de-coupling of the heave and pitch motions. 

In the following subsections, the hydrodynamic parameters and responses are presented and discussed for the TALOS floating structure, and the comparisons between HAMS and WAMIT aim to confirm that the open-source panel method, HAMS, is reliable and accurate.
section_name: Added Mass

section_content: The conventional added mass of the TALOS wave-energy converter is presented in Figure 2, including the added mass for surge ( 11 ), sway ( 22 ), heave ( 33 ), roll ( 44 ), pitch ( 55 ), and yaw ( 66 ). We can observe from the comparisons between HAMS and WAMIT that they are almost identical, and some differences (spikes) can be observed in the socalled irregular frequencies in the short waves, when the wave period,, is less than 3 s. The reason for the irregular frequencies may relate to the panel size: as a rule of thumb, the maximal panel size should be less than 1/8 of the wave length (see [37] ). 

However, the spiky responses of the short waves (of high frequencies) would have no significant influence on the performance of the floating structure (see the smooth-motion responses, the RAOs in Section 3. 5), or in the interest for wave-energy extraction due to the corresponding low-wave-energy density of short waves. For the free-floating structure, TALOS, we can observe that there are some strong couplings between the motions in terms of the added mass, such as the surge pitch (see 15, Figure 3a ) and sway roll (see 24, Figure 3b ). It can also be noted that, due to the However, the spiky responses of the short waves (of high frequencies) would have no significant influence on the performance of the floating structure (see the smooth-motion responses, the RAOs in Section 3. 5), or in the interest for wave-energy extraction due to the corresponding low-wave-energy density of short waves. 

For the free-floating structure, TALOS, we can observe that there are some strong couplings between the motions in terms of the added mass, such as the surge pitch (see A 15, Figure 3a ) and sway roll (see A 24, Figure 3b ). It can also be noted that, due to the symmetry of the couplings, A 51 = A 15, and A 42 = A 24, and also the corresponding added mass of the strong coupling terms, which have a similar order of magnitude as the conventional added mass (compared to Figure 2a, b ). Here, we can observe that A 24 is a negative added mass, which means that the motion of the roll would produce an added mass on the sway in an opposite phase. Obviously, in the applications, these strong coupling terms must be included in both the frequency-domain and time-domain analyses. 

and WAMIT predictions, probably due to the slightly different methods employed to approximate the free-surface Green function. However, the magnitude of the added mass for the weak couplings would be several orders smaller than the conventional added mass (~300 vs. 3 × 10 6 ). Therefore, the weak coupling terms can be fully omitted in the numerical modeling since these terms would not have any influence on the overall performance of the structure motions and responses, although such a weak coupling presented between the heave and pitch of TALOS may be an indication that the heave and pitch could be strongly coupled if a different structure, for instance, a floating structure with no symmetry for the y axis was applied (here, TALOS has 2 symmetries for the x and y axes, respectively).
section_name: Added Mass

section_content: The conventional radiation-damping coefficients of the TALOS wave-energy converter are presented in Figure 4. These include the radiation-damping coefficients, 11 for surge, 22 for sway, 33 for heave, 44 for roll, 55 for pitch, and 66 for yaw. Again, we can observe that the results from HAMS and WAMIT are almost identical, with differences at the irregular frequencies in the short waves (the wave periods,, less than 4 s). Again, these irregular frequencies are due to the limited panel size, similar to the cases in the added mass. For a comparison, we also present the weak couplings in the heave pitch (see A 35, Figure 3c ), and, here, we can observe that they are quite different in relation to the HAMS and WAMIT predictions, probably due to the slightly different methods employed to approximate the free-surface Green function. However, the magnitude of the added mass for the weak couplings would be several orders smaller than the conventional added mass (~300 vs. 3 × 10 6 ). Therefore, the weak coupling terms can be fully omitted in the numerical modeling since these terms would not have any influence on the overall performance of the structure motions and responses, although such a weak coupling presented between the heave and pitch of TALOS may be an indication that the heave and pitch could be strongly coupled if a different structure, for instance, a floating structure with no symmetry for the y axis was applied (here, TALOS has 2 symmetries for the x and y axes, respectively).
section_name: Radiation-Damping Coefficients

section_content: The conventional radiation-damping coefficients of the TALOS wave-energy converter are presented in Figure 4. These include the radiation-damping coefficients, B 11 for surge, B 22 for sway, B 33 for heave, B 44 for roll, B 55 for pitch, and B 66 for yaw. Again, we can observe that the results from HAMS and WAMIT are almost identical, with differences at the irregular frequencies in the short waves (the wave periods, T, less than 4 s). Again, these irregular frequencies are due to the limited panel size, similar to the cases in the added mass. 

Similarly, we can observe that there are some strong couplings between the motions, such as surge pitch (see B 15, Figure 5a ) and sway roll (see B 24, Figure 5b ), and these coupling terms would correspond to those of the strong coupling terms on added mass. These strongly coupled radiation-damping coefficients would have similar orders of magnitude as those conventional radiation-damping coefficients (refer to Figure 4a, b ). One interesting factor is that the coupled radiation-damping coefficient can be negative, for instance, the coupled radiation-damping coefficients for sway and roll, B 24 (Figure 5b ), which is very different from the positive conventional radiation-damping coefficient. This can be understood as follows: in the conventional radiation-damping coefficients, the selfradiation would simply radiate the wave away due to the structure motion, and dissipate the energy from the moving structure, hence the radiation-damping coefficient must be positive; while in the coupled damping coefficients, for instance, B 24, which is the damping coefficient due to the radiated wave from the roll motion (the subscript '4') contributing to the damping effect on the sway motion (the subscript '2'), the coupling effect is in an opposite phase, hence it is physically negative (for a comparison, B 15 is positive, see Figure 5a ). 

verter are presented in Figure 4. These include the radiation-damping coefficients, 11 for surge, 22 for sway, 33 for heave, 44 for roll, 55 for pitch, and 66 for yaw. Again, we can observe that the results from HAMS and WAMIT are almost identical, with differences at the irregular frequencies in the short waves (the wave periods,, less than 4 s). Again, these irregular frequencies are due to the limited panel size, similar to the cases in the added mass. Similarly, we can observe that there are some strong couplings between the motions, such as surge pitch (see 15, Figure 5a ) and sway roll (see 24, Figure 5b ), and these coupling terms would correspond to those of the strong coupling terms on added mass. These strongly coupled radiation-damping coefficients would have similar orders of magnitude as those conventional radiation-damping coefficients (refer to Figure 4a, b ). One interesting factor is that the coupled radiation-damping coefficient can be negative, for instance, the coupled radiation-damping coefficients for sway and roll, 24 (Figure 5b ), which is very different from the positive conventional radiation-damping coefficient. This can be understood as follows: in the conventional radiation-damping coefficients, the self-radiation would simply radiate the wave away due to the structure motion, and dissipate the energy from the moving structure, hence the radiation-damping coefficient must be positive; while in the coupled damping coefficients, for instance, 24, which is the damping coefficient due to the radiated wave from the roll motion (the subscript '4') contributing to the damping effect on the sway motion (the subscript '2'), the coupling effect is in an opposite phase, hence it is physically negative (for a comparison, 15 is positive, see Figure 5a ). 

Similarly, we can observe the weak coupling presented in the heave pitch (see 35, Figure 5c ). Again, the magnitude of the radiation-damping coefficient for the weak coupling is a few orders smaller than the conventional radiation-damping coefficient (~200 vs. 2 × 10 6 ). Therefore, the weak coupling terms can be omitted from the numerical modeling. In fact, in the corresponding time-domain analysis, such weak coupling terms must be dropped; otherwise, it may cause divergence problems in the numerical modeling due to the possibly corresponding divergent-impulse functions (the details are discussed later in the research). Similarly, we can observe that there are some strong couplings between the motions, such as surge pitch (see 15, Figure 5a ) and sway roll (see 24, Figure 5b ), and these coupling terms would correspond to those of the strong coupling terms on added mass. These strongly coupled radiation-damping coefficients would have similar orders of magnitude as those conventional radiation-damping coefficients (refer to Figure 4a, b ). One interesting factor is that the coupled radiation-damping coefficient can be negative, for instance, the coupled radiation-damping coefficients for sway and roll, 24 (Figure 5b ), which is very different from the positive conventional radiation-damping coefficient. This can be understood as follows: in the conventional radiation-damping coefficients, the self-radiation would simply radiate the wave away due to the structure motion, and dissipate the energy from the moving structure, hence the radiation-damping coefficient must be positive; while in the coupled damping coefficients, for instance, 24, which is the damping coefficient due to the radiated wave from the roll motion (the subscript '4') contributing to the damping effect on the sway motion (the subscript '2'), the coupling effect is in an opposite phase, hence it is physically negative (for a comparison, 15 is positive, see Figure 5a ). 

Similarly, we can observe the weak coupling presented in the heave pitch (see 35, Figure 5c ). Again, the magnitude of the radiation-damping coefficient for the weak coupling is a few orders smaller than the conventional radiation-damping coefficient (~200 vs. 2 × 10 6 ). Therefore, the weak coupling terms can be omitted from the numerical modeling. In fact, in the corresponding time-domain analysis, such weak coupling terms must be dropped; otherwise, it may cause divergence problems in the numerical modeling due to the possibly corresponding divergent-impulse functions (the details are discussed later in the research). Similarly, we can observe the weak coupling presented in the heave pitch (see B 35, Figure 5c ). Again, the magnitude of the radiation-damping coefficient for the weak coupling is a few orders smaller than the conventional radiation-damping coefficient (~200 vs. 2 × 10 6 ). Therefore, the weak coupling terms can be omitted from the numerical modeling. In fact, in the corresponding time-domain analysis, such weak coupling terms must be dropped; otherwise, it may cause divergence problems in the numerical modeling due to the possibly corresponding divergent-impulse functions (the details are discussed later in the research).
section_name: Radiation-Damping Coefficients

section_content: Wave-excitation forces are dependent on both the floating structure itself and the direction of the incoming waves. Here, an angle of 45 • for the incoming wave was used for calculating the wave excitations acting on the structure. From Figure 6, we can observe the wave-excitation responses (under the waves of unit amplitude) for different motion modes (surge, F ex1 ; sway, F ex2 ; heave, F ex3 ; roll, F ex4 ; pitch, F ex5 ; and yaw, F ex6 ). Additionally, the wave excitations for surge and sway are the same (Figure 6a, b ) in such a wave condition (45 • ) and due to the symmetry of the TALOS structure. Additionally, for the same reason, the excitations for roll and pitch are the same (see Figure 6d, e ). 

tionally, the wave excitations for surge and sway are the same (Figure 6a, b ) in such a wave condition (45°) and due to the symmetry of the TALOS structure. Additionally, for the same reason, the excitations for roll and pitch are the same (see Figure 6d, e ). 

The excitations for the surge, sway, roll, and pitch all have maximums at certain wave periods (see Figure 6a, b, d, e ), while the excitation for the heave monotonically increases in short waves (Figure 6c ), and its increase is in long waves, but in very long waves, the increase would be small. The reason for these is because the excitations for the surge, sway, roll, and pitch are all dependent on the slope of the wave, while the excitation for heave is dependent on the dynamic pressure acting on the bottom of the structure, thus it will be greater in longer waves. Therefore, the excitation for heave increases with the longer waves.
section_name: Wave Excitation

section_content: For the free-floating marine structures, only three motion modes (out of six DOFs) had inherent restoring forces: heave, roll, and pitch, and the latter two motions would also be very dependent on the position of the center of gravity of the structure. In the current modeling, the center of gravity was set at = -5. 0 (that is, 5. 0 m below the water's surface). 

For the other three motion modes (surge, sway, and yaw), there were no inherent restoring forces. Hence, for most floating structures, the required restoring forces/moments must be provided by the mooring system to the structure or other means. Here, for simplifying the analysis, we assumed the linear restoring coefficients for the surge, sway, The excitations for the surge, sway, roll, and pitch all have maximums at certain wave periods (see Figure 6a, b, d, e ), while the excitation for the heave monotonically increases in short waves (Figure 6c ), and its increase is in long waves, but in very long waves, the increase would be small. The reason for these is because the excitations for the surge, sway, roll, and pitch are all dependent on the slope of the wave, while the excitation for heave is dependent on the dynamic pressure acting on the bottom of the structure, thus it will be greater in longer waves. Therefore, the excitation for heave increases with the longer waves.
section_name: Response Amplitude Operators (RAOs)

section_content: For the free-floating marine structures, only three motion modes (out of six DOFs) had inherent restoring forces: heave, roll, and pitch, and the latter two motions would also be very dependent on the position of the center of gravity of the structure. In the current modeling, the center of gravity was set at z g = -5. 0 m (that is, 5. 0 m below the water's surface). 

For the other three motion modes (surge, sway, and yaw), there were no inherent restoring forces. Hence, for most floating structures, the required restoring forces/moments must be provided by the mooring system to the structure or other means. Here, for simplifying the analysis, we assumed the linear restoring coefficients for the surge, sway, and yaw motions to be: C 11 = C 22 = 2. 5 × 10 6 N/m and C 66 = 5 × 10 6 Nm (and these added restoring coefficients mean a relatively stiff mooring system). It should be noted that these linear restoring coefficients can be obtained by linearizing the mooring forces, and the linearization of the mooring forces can be valid if the structure motions are small in magnitude. 

Under the above mooring setting, the surge and sway motion responses are the same, as can be observed in Figure 7a,b, due to the symmetrical TALOS and the same added restoring-force coefficients for the surge and sway. The large responses presented in Figure 7a,b at T = 10. 83 s correspond to the resonance period of the surge and sway, respectively, while the very small responses at T = 9. 24 s are due to the respective responses of the pitch and roll motions. 

restoring-force coefficients for the surge and sway. The large responses presented in Figure 7a,b at T = 10. 83 s correspond to the resonance period of the surge and sway, respectively, while the very small responses at T = 9. 24 s are due to the respective responses of the pitch and roll motions. 

From Figure 7c, it can be observed that the heave motion response would only be dependent on the structure itself, and the resonance period is T = 7. 22 s, while for the roll and pitch (Figure 7d,e), their responses are identical, simply due to the symmetrical structure of the TALOS and incoming wave (45°). The roll and pitch have large responses at the resonance period of T = 10. 83 s due to the couplings of the surge and sway, and the small responses at period T = 9. 24 s, which would correspond to the resonances of the roll and pitch motions. 

For the yaw response, due to the symmetries of the TALOS structure (a nearly axisymmetrical structure), it is very small, so small as an effective zero (see Figure 7f ).
section_name: Response Amplitude Operators (RAOs)

section_name: Time-Domain Dynamic Equation and Analysis

section_content: The time-domain dynamic equation is in fact the dynamic equation based on Newton's second law of motion for multiple motion modes of a rigid body, and, generally, this dynamic equation is written in the form of mass-spring-damper, with the damping force and the restoring force acting on the structure being moved to the left-hand side. With the From Figure 7c, it can be observed that the heave motion response would only be dependent on the structure itself, and the resonance period is T = 7. 22 s, while for the roll and pitch (Figure 7d,e), their responses are identical, simply due to the symmetrical structure of the TALOS and incoming wave ( 45 • ). The roll and pitch have large responses at the resonance period of T = 10. 83 s due to the couplings of the surge and sway, and the small responses at period T = 9. 24 s, which would correspond to the resonances of the roll and pitch motions. 

For the yaw response, due to the symmetries of the TALOS structure (a nearly axisymmetrical structure), it is very small, so small as an effective zero (see Figure 7f ).
section_name: Dynamic Equation

section_content: The time-domain dynamic equation is in fact the dynamic equation based on Newton's second law of motion for multiple motion modes of a rigid body, and, generally, this dynamic equation is written in the form of mass-spring-damper, with the damping force and the restoring force acting on the structure being moved to the left-hand side. With the inclusions of the wave-radiation forces, wave excitation as well as the additional external forces, the dynamic equation would be presented as:

where

x k (t) (k = 1, 2,. . ., 6) are the 6-DOF motions of the structure in the time domain. f j (t) (j = 1, 2,. . ., 6) is the time-domain wave-excitation force, which can be calculated based on the wave-excitation response from the panel method analysis. 

For a linear dynamic equation, that means, the mass matrices, the damping-coefficient matrices and the restoring-coefficient matrices are all independent of the motions. Additionally, under a sinusoidal excitation, for instance, f j (t) = F j e iωt (here, F j is a complex amplitude of the jth wave excitation), the motion would be sinusoidal too, then we shall have the following relations:

Motion: x k (t) = ξ k e iωt Velocity:

. 

x k (t) = iωξ k e iωt Acceleration:

. . 

x k (t) = -ω 2 ξ k e iωt where ξ k is the complex amplitude of the motion. Substituting the above relations into Equation (4), and dropping the time factors in all terms, we can obtain the frequency-domain equation, which would be exactly same as in Equation ( 1 ). However, it should be noted that, when transforming the frequencydomain Equation (1) back to the time-domain Equation ( 4 ), we should be very careful in determining the frequency-dependent added-mass and radiation-damping coefficients. 

The direct transformation from Equation (1) to Equation ( 4 ) can only be correct if the dynamic system is linear and under a sinusoidal excitation, because in such a dynamic system, the excitation and motions are all sinusoidal (i. e., of single frequency), hence the frequency-dependent parameters, A jk, B jk and F j, can be determined and all clearly have a physical significance, as does the time-domain Equation (4). 

However, such a fully linear dynamic system is not the purpose of employing the time-domain approach, and this is especially true for most wave-energy converters. For instance, to examine the wave-energy converter's performance in irregular waves, the excitation and the motions in irregular waves would all have multi-frequencies. As such, the frequency-dependent terms, A jk and B jk, become undecided, although the time-domain excitation is calculated, see Equation ( 6 ) below. A more general case is that the wave-energy converter may include some nonlinear forces in the dynamic system, such as the force from the nonlinear PTO or the force from the PTO control for improving wave-energy extraction, or the nonlinear forces from the mooring system; the motions would surely be have multiple frequencies, even if the excitation is sinusoidal. As such, A jk and B jk become undetermined too. 

The most used time-domain equation for the wave-structure interaction would be the Cummins' time-domain equation [28] with the Ogilvie relation [29], namely, the hybrid frequency-time domain method, according to [25], and the time-domain equation has the following form:

where

x k (t) (k = 1, 2,. . ., 6) are the structure motions (the motions of 6 DOFs of the structure, which would be solved in the dynamic equation);

A jk (∞) (j, k = 1, 2,. . ., 6) is the added-mass matrix at infinite frequency, which can be assessed using the frequency-domain analysis data (see the next section); K jk (t) (j, k = 1, 2,. . ., 6) is the impulse function, which would be calculated based on the frequency-domain results (see the next section);

The parameters with superscript ' E ' are the parameters/forces externally added to the dynamic system, such as f E j (t) (j = 1, 2,. . ., 6) the external force, for instance, the force from power take-off (PTO), or the control force, or the force from the mooring system, and these forces can be both linear and nonlinear. 

f j (t) (j = 1, 2,. . ., 6) is the wave-excitation force in an irregular wave of a given wave spectrum. The corresponding wave-excitation force can be calculated using the frequency-dependent wave-excitation responses, F j. The time-domain wave excitation is calculated as

where F j (ω) is a complex wave-excitation response (assessed using the panel method); ε n is the random phase for the wave components in an irregular wave;

A n is the wave amplitude of the wave component corresponding to the frequency ω n (of a bandwith ∆ω n ), calculated as

with S(ω) being the wave spectral density, such as the popular wave spectra, Bretschneider or JONSWAP spectrum [42].
section_name: Dynamic Equation

section_content: In the transformation from frequency domain to time domain, the relevant frequencydependent parameters must be transformed accordingly, including the added mass and the radiation damping coefficients, while the wave excitation in irregular waves can be assessed in Equation (6). 

From the Cummins' time-domain equation, the frequency-dependent damping terms in the frequency-domain equation are all replaced using the corresponding convolution terms, or, more specifically, the memory effects, due to the relevant fluid motions. The memory effect is presented in a convolution term as

Accordingly, in the Cummins' time-domain equation, the added mass at infinite frequency, A jk (∞), replaces the frequency-dependent added mass in the frequency-domain equation, see Equation (5).
section_name: Memory Effect and the Impulse Functions

section_content: Based on the Ogilvie relation [29], the impulse function for the wave-structure interaction can be calculated from the frequency-dependent radiation-damping coefficients, K jk, given as

In the numerical modeling, the sampled impulse function can be obtained from the sampled radiation-damping coefficient (here, the subscripts are dropped for a simplicity of the expression), calculated as

where N is the total number of the frequencies in the frequency-domain analysis, and m = 1, 2,. . ., M is the number of the sampled impulse function, with a time interval ∆t, and the time series is presented as t m = (m -1)∆t. One can observe that the calculation of the impulse function is from a low frequency to a reasonably high frequency, not necessarily from zero to infinite frequency (note: the radiation-damping coefficient is asymptotic to zero on both sides of the wave frequencies). A MATLAB function for calculating the impulse function can be found in Appendix A. 

In the numerical modeling, based on the assessment of the frequency-domain radiationdamping coefficients, we can use Equation (10) to calculate the impulse function and, for the purpose of validation, WAMIT F2T utility was used to obtain the impulse function [39]. Figure 8 shows the comparisons of the results calculated from the HAMS radiation-damping coefficients, as described above, and the results obtained from WAMIT F2T for the conventional impulse functions, K 11, K 22, K 33, K 44, K 55, K 66. It can be observed that they are all almost identical between the results obtained for HAMS and WAMIT. 

the time series is presented as = ( -1). One can observe that the calculation of the impulse function is from a low frequency to a reasonably high frequency, not necessarily from zero to infinite frequency (note: the radiation-damping coefficient is asymptotic to zero on both sides of the wave frequencies). A MATLAB function for calculating the impulse function can be found in Appendix A. 

In the numerical modeling, based on the assessment of the frequency-domain radiation-damping coefficients, we can use Equation (10) to calculate the impulse function and, for the purpose of validation, WAMIT F2T utility was used to obtain the impulse function [39]. Figure 8 shows the comparisons of the results calculated from the HAMS radiationdamping coefficients, as described above, and the results obtained from WAMIT F2T for the conventional impulse functions, 11, 22, 33, 44, 55, 66. It can be observed that they are all almost identical between the results obtained for HAMS and WAMIT. When we examined the dynamically coupling terms between the different motion modes, we observed the coupling terms between the surge pitch, 15, and sway roll, 24 (see Figure 9a, b ). These were strong coupling terms, since their magnitudes would be in orders similar to those of the conventional impulse functions. In comparison, the impulse functions for the strong coupling terms from HAMS and WAMIT were almost identical. 

Figure 9c shows a weak coupling between the heave and pitch, 35. Here, we can observe a great difference between HAMS and WAMIT, and also the weak coupling impulse functions are very different from those of the conventional impulse functions: the weak coupling impulse function can be divergent. 

Having observed the great differences, their magnitudes are essentially much smaller than the conventional impulse functions (~10. 10 6 ). In fact, for a practical purpose, these weak coupling terms must be dropped in the numerical modeling, otherwise it may When we examined the dynamically coupling terms between the different motion modes, we observed the coupling terms between the surge pitch, K 15, and sway roll, K 24 (see Figure 9a, b ). These were strong coupling terms, since their magnitudes would be in orders similar to those of the conventional impulse functions. In comparison, the impulse functions for the strong coupling terms from HAMS and WAMIT were almost identical.
section_name: Impulse Function

section_content: Once the impulse function has been calculated, the added mass at infinite frequency can the be assessed as

In a discretized form, for a given frequency,, the added mass at infinite frequency is calculated as Figure 9c shows a weak coupling between the heave and pitch, K 35. Here, we can observe a great difference between HAMS and WAMIT, and also the weak coupling impulse functions are very different from those of the conventional impulse functions: the weak coupling impulse function can be divergent. 

Having observed the great differences, their magnitudes are essentially much smaller than the conventional impulse functions (∼ 10 vs. 10 6 ). In fact, for a practical purpose, these weak coupling terms must be dropped in the numerical modeling, otherwise it may cause a divergence problem in the time-domain calculations due to the divergent impulse function.
section_name: Added Mass at Infinite Frequency

section_content: Once the impulse function has been calculated, the added mass at infinite frequency can the be assessed as

In a discretized form, for a given frequency, ω n, the added mass at infinite frequency is calculated as

In principle, the added mass at infinite frequency can be calculated based on a given frequency, say ω n. However, due to the possible errors in the numerical modeling, such as irregular frequencies, a more reliable calculation for the added mass at infinite frequency is to average the added mass at infinite frequency based on the corresponding multiple frequencies, ω 1, ω 2, ω 3,. . ., ω N, as

Using the above method, we can calculate the added mass at infinite frequency based on the frequency-dependent added mass and the calculated impulse function (see also Appendix A). Table 1 shows the comparison of the added mass at infinite frequency from the HAMS calculation and the assessment using WAMIT. Here, we can observe that they are very close, especially for the conventional added mass at infinite frequency, with an error less than 0. 5%, while a slightly large difference occurs for the coupled added mass A 15, with an error of about 2. 5%. Hence, we can conclude that all the errors should be within the acceptable range in the engineering calculations.
section_name: Added Mass at Infinite Frequency

section_name: Approximations of Impulse Function and Memory Effect

section_content: Generally, the normal impulse function K(t) would be a decayed function, which would converge to zero with time. This feature of the impulse function corresponds to the fluid-memory effect in physics: the events closer to the current time would have a greater influence on the current motion, and with time lapsing further, the influence becomes less and less, and eventually fades to zero (that is, K(t) = 0 when t is large). 

There are generally two different ways to approximate such impulse functions, namely, the Prony approximation method [31, 32, 43] and the state-space method [30]. Here, the Prony function approximation was employed, because the Prony function has a more straightforward physical meaning for approximating the impulse function. 

Following Duclos et al. [31], the impulse function can be approximated using the Prony function in the following form:

where N 0 is the order of the Prony function, and the complex coefficients α k and β k can be obtained from the Prony method. For a convergent impulse function, the real part of the parameter, β k, must be negative. If it is not, the corresponding term must be dropped since it is divergent in nature (a MATLAB function for obtaining α k and β k from the impulse function, K(t), can be found in Appendix B). 

Figure 10 shows the approximations of the impulse functions for the self-radiations; here, N 0 = 10 (that is, the 10th order Prony function was used). It can be observed that all the impulse functions have been well reproduced. Figure 11 shows the impulse functions for the strong coupled terms, such as K 15 (note: K 51 is symmetric to K 15 ) and K 24. Again, these impulse functions have been well reproduced by the Prony approximations.
section_name: Approximation of Impulse Function

section_content: Generally, the normal impulse function () would be a decayed function, which would converge to zero with time. This feature of the impulse function corresponds to the fluid-memory effect in physics: the events closer to the current time would have a greater influence on the current motion, and with time lapsing further, the influence becomes less and less, and eventually fades to zero (that is, () = 0 when is large). 

There are generally two different ways to approximate such impulse functions, namely, the Prony approximation method [31, 32, 43] and the state-space method [30]. Here, the Prony function approximation was employed, because the Prony function has a more straightforward physical meaning for approximating the impulse function. 

Following Duclos et al. [31], the impulse function can be approximated using the Prony function in the following form:

where 0 is the order of the Prony function, and the complex coefficients and can be obtained from the Prony method. For a convergent impulse function, the real part of the parameter,, must be negative. If it is not, the corresponding term must be dropped since it is divergent in nature (a MATLAB function for obtaining and from the impulse function, (), can be found in Appendix B). 

Figure 10 shows the approximations of the impulse functions for the self-radiations; here, 0 = 10 (that is, the 10th order Prony function was used). It can be observed that all the impulse functions have been well reproduced. Figure 11 shows the impulse functions for the strong coupled terms, such as 15 (note: 51 is symmetric to 15 ) and 24. Again, these impulse functions have been well reproduced by the Prony approximations. For the weak coupling terms, such as 13, 35, and 46 (Figure 12 ), it seems more difficult to reproduce the weak impulse functions, though 35 seems good (Figure 12b ). Since these weakly coupled impulse functions are unlike the well-convergent normal impulse functions, they may not decay with time (Figure 12a, c ), or simply diverge with time For the weak coupling terms, such as K 13, K 35, and K 46 (Figure 12 ), it seems more difficult to reproduce the weak impulse functions, though K 35 seems good (Figure 12b ). Since these weakly coupled impulse functions are unlike the well-convergent normal impulse functions, they may not decay with time (Figure 12a, c ), or simply diverge with time (Figure 12b ). For such impulse functions, the Prony method would not correctly reproduce the impulse function due to the omission of the terms when the real part of β k is positive (see Figure 12b ). For the weak coupling terms, such as 13, 35, and 46 (Figure 12 ), it seems more difficult to reproduce the weak impulse functions, though 35 seems good (Figure 12b ). Since these weakly coupled impulse functions are unlike the well-convergent normal impulse functions, they may not decay with time (Figure 12a, c ), or simply diverge with time (Figure 12b ). For such impulse functions, the Prony method would not correctly reproduce the impulse function due to the omission of the terms when the real part of is positive (see Figure 12b ). 

However, if these weak impulse functions are compared to the strong coupled impulse functions, their magnitudes could be many orders smaller than the strong coupling terms, for instance, 35 ~10 vs. 15 ~10 6. Thus, in numerical modeling, these weak coupled impulse functions must be dropped to ensure the convergent solution in the timedomain analysis, while these small amplitude impulse functions have no influence on the overall motions of the structure. The order of the Prony function, 0, would be an important factor to determine how accurately the impulse functions can be approximated. Taking 33 as an example, we can determine the differences using different orders of the Prony function to approximate the impulse function (Figure 13 ): 0 = 6 is quite good, while 0 = 5 seems much worse; 0 = 8 could provide a much closer approximation, while for 0 = 10, the approximation would be almost identical to the original impulse function. However, if these weak impulse functions are compared to the strong coupled impulse functions, their magnitudes could be many orders smaller than the strong coupling terms, for instance, K 35 ~10 vs. K 15 ~10 6. Thus, in numerical modeling, these weak coupled impulse functions must be dropped to ensure the convergent solution in the time-domain analysis, while these small amplitude impulse functions have no influence on the overall motions of the structure. 

The order of the Prony function, N 0, would be an important factor to determine how accurately the impulse functions can be approximated. Taking K 33 as an example, we can determine the differences using different orders of the Prony function to approximate the impulse function (Figure 13 ): N 0 = 6 is quite good, while N 0 = 5 seems much worse; N 0 = 8 could provide a much closer approximation, while for N 0 = 10, the approximation would be almost identical to the original impulse function.
section_name: Approximation of Impulse Function

section_content: Substituting the Prony approximation, Equation ( 14 ), to the formula of the memory effect, Equation ( 8
section_name: Calculation of the Memory Effect

section_content: Substituting the Prony approximation, Equation (14), to the formula of the memory effect, Equation (8), yields

x(τ)dτ =

in which the order of integral and summation has been swapped, and

Up to here, the calculation of the memory effect would be concentrated on how well and how fast we can calculate or solve for I k (t) (k = 1, 2,. . ., N 0 ). Subsequently, two different ways to solve I k (t) would be introduced. 

(1) Method '1': solving additional differential equations [31] Taking a differentiation of I k (t) with regard to t, we obtain. 

x(τ)dτ + α k e β k t e -β k t. 

where the expression in the square brackets is a differentiation of the integral in Equation ( 16 ) with regard to time, t. 

Obviously, Equation ( 17 ) is an actual first-order differential equation for I k. In the classical Prony method [31], these first-order differential equations are simply added to the dynamic equations for solving I k, when solving for the motions of the structure. Once I k is solved, the memory effect can be calculated using Equation (15). 

In the numerical modeling using the classical Prony method, if multiple motion modes are involved, the number of the additional differential equations would be large. For example, taking 2 coupled motions (say, surge-pitch coupling), we would have 4 memory-effect terms (K 11, K 15, K 51, and K 55 ). Now, if we use a 10th-order Prony function, then we would have a total of 40 additional first-order differential equations, added to 2 s-order dynamic equations for the motions (i. e., surge and pitch). If we write the second-order differential equation in the form of first-order differential equations, we have 4 first-order differential equations for the motions. Hence, here, for 2 s-order differential equations (equivalent to 4 first-order differential equations) for motions, we need to add 40 first-order equations to solve the memory effects. Thus, the computation time would be greatly increased. 

(2) Method '2': a recursive method for calculating I k [32] Based on Equation ( 16 ), we can easily calculate

x(τ)dτ = α k e β k (t+∆t) t 0 e -β k τ. 

x(τ)dτ + α k e β k (t+∆t) t+∆t t e -β k τ. 

x(τ)dτ

x(τ)dτ

where the first term of the last expression can be easily obtained, and the integral from t to t + ∆t in the second term is formulated using the mid-point value of the expression e -β k τ ; that is, at τ = t + ∆t/2 (from the 2nd to 3rd rows in the square brackets). In a sampled system, similar to the numerical modeling, a sampled form of I k can be written in a form as

It can be observed that I k (n + 1) is a renewal based on the previous result I k (n) by timing a time factor e β k ∆t to modify the previous term, and add a modification term due to the motion increment ∆x(n) (the 2nd term of RHS in Equation ( 19 )). By employing such a scheme, the memory effect can be calculated in a very straightforward manner to avoid adding any differential equations to the dynamic system. As shown in [32], using the recursive method for calculating the memory effect, the computational time would not increase, even if a higher-order Prony function is used for many coupling terms. In this research, this approach was applied to calculate the fluid-memory effects.
section_name: Calculation of the Memory Effect

section_name: Implementation and Validation of Time-Domain Analysis

section_content: In solving the dynamic equation numerically, the time-domain equation can be simply written in a sampled form as

where {M} is the mass matrix, which includes the mass matrix of the structure, the external- mass and the added-mass matrices at infinite frequency (hereafter, the parameter in the curly brackets represents a matrix);

. . 

x (t + ∆t) is the acceleration vector at the time step t + ∆t (the parameter in the square brackets is a vector);

[F(t)] is the total-force vector at the time step t. This force vector includes all the forces, wave excitation, the damping force (and memory effect), restoring force, and any other added forces, defined as

Then, the acceleration vector can be obtained as. . 

where {M} -1 is the inverse mass matrix of {M}. Such a differential equation can be easily solved using a numerical method, for instance, using the popular MATLAB scheme 'ode45'. For a given irregular wave of significant wave height H s = 2. 0 m and spectral peak period T p = 8. 0 s (Bretschneider spectrum), the wave incoming angle is 45 •, and we can solve the motions of the structure for the waves. Figure 14 shows the time histories of the surge, heave, and pitch motions of the TALOS structure. Additionally, in the following section, a validation method is presented for examining whether the time-domain numerical scheme is correctly made or not.
section_name: Implementation of Time-Domain Analysis

section_content: In building the hybrid frequency-time domain analysis, there are a few transformations in the process, including:
section_name: Validations of Time-Domain Modeling

section_content: The transformation from the radiation-damping coefficients to impulse functions, as well as how to efficiently approximate the impulse functions;
section_name: •

section_content: The transformation of the frequency-dependent added mass to the added mass at the infinite frequency;
section_name: •

section_content: The transformation from the radiation-damping effects to the memory effects, together with the method for how we can reliably and rapidly calculate the memory effect;
section_name: •

section_content: The inclusions of the coupling terms between different motion modes, especially the calculation of the coupled memory-effect terms;
section_name: •

section_content: The transformation of the excitation responses to the forces in the time domain for a given wave spectrum. 

In implementing the time-domain numerical scheme, things can go wrong quite easily, especially with the multiple coupled terms involved in the coupled memory effects. For instance, if the program can run as we hoped, and generate the corresponding time histories, as shown in Figure 14, how we can make sure the numerical modeling is correctly implemented? Or, do we have a validation method for the numerical result?

where {} -1 is the inverse mass matrix of {}. Such a differential equation can be easily solved using a numerical method, for in stance, using the popular MATLAB scheme 'ode45′. For a given irregular wave of signifi cant wave height = 2. 0 m and spectral peak period = 8. 0 s (Bretschneider spectrum the wave incoming angle is 45°, and we can solve the motions of the structure for th waves. Figure 14 shows the time histories of the surge, heave, and pitch motions of th TALOS structure. Additionally, in the following section, a validation method is presented for examining whether the time-domain numerical scheme is correctly made or not.
section_name: •

section_content: In building the hybrid frequency-time domain analysis, there are a few transfor mations in the process, including:
section_name: Validations of Time-Domain Modeling

section_content: The transformation from the radiation-damping coefficients to impulse functions, a well as how to efficiently approximate the impulse functions;
section_name: •

section_content: The transformation of the frequency-dependent added mass to the added mass at th infinite frequency;
section_name: •

section_content: The transformation from the radiation-damping effects to the memory effects, to gether with the method for how we can reliably and rapidly calculate the memor effect; Generally, in the hybrid frequency-time domain method, the approximations of the impulse functions and of the memory-effect terms can be complicated in the numerical implementation, especially for those complicated coupled multiple-motion modes, in which the matrices and vectors are both involved in the calculation (in comparison, the addition of the external forces to the dynamic system would be relatively easy, since only vectors are involved), and mistakes may be easily made. Hence, a validation method is required to check the correctness of the implementation of the time-domain analysis, and it is possible to examine the transformation of the equation from the frequency domain to the time domain, if the dynamic system is linear, similar to the Cummins' time-domain equation, since both the frequency-domain and time-domain analyses can be conducted in exactly the same conditions, so they are comparable:

(1) The time-domain analysis can be obtained by solving Equation (22), in which only the linear forces are applied. For the purpose of comparison, all these linear forces must be presented in the frequency-domain analysis, too, and the corresponding motion responses in the frequency domain can be obtained under same conditions as those in the time-domain analyses. (2) For a comparison, a time history can be generated directly based on the RAOs in the frequency-domain analysis, in which a transformation is made in a formula as

x FD j (t) = ∞ ∑ n=1 χ j (ω n )e iω n t A n e iε n = ∞ ∑ n=1 χ j (ω n )A n e i(ω n t+ε n )

where χ j (ω) is the complex RAO of jth motion mode, and for 6-DOF motions, j = 1, 2,. . ., 6. ε n the random phase must be the same as that for generating the wave-excitation force in Equation (6). 

In comparison, Figure 15 shows the time histories of the 6-DOF motions of TALOS in a wave of H s = 2. 0 s, T p = 8. 0 s, an angle of incoming wave of 45 •, with a linear mooring, which is the same as in Section 3. 5. Here, we can observe that the results of the time-domain ('TD' in the figure, directly solving the time-domain equation, Equation ( 22 )) and the frequency-domain analyses ('FD' in the figure, from Equation ( 23 )) are very close after about 15 wave cycles, except for the yaw motion. In the first 10-15 wave cycles, the timedomain solutions may be different from those transformed from the RAOs subject to the preset initial motions, and, generally, after a few wave cycles, the effects of the preset initial conditions are finally dissipated; thus, the very close results can be observed for both the time-domain (solving Equation ( 22 )) and frequency-domain solutions (from Equation ( 23 )). Such comparisons can confirm that the numerical scheme for the transformation from the frequency domain to the time domain is correct. This would be a solid first step towards a more complicated time-domain analysis, for instance, the inclusions of the nonlinear forces from PTO from the mooring system and other sources. Moreover, the addition of the external forces would be much easier when compared to the calculation of the memory effects of all the coupled terms. For the yaw motion (x6), the results from the time-domain and frequency-domain analyses are very different. However, the yaw motions are so small in magnitude (~10 -9 ), and they can be fully ignored in the analysis.
section_name: •

section_content: In this research, we presented a hybrid frequency-time domain approach for wavestructure interactions, in which all the transformations to and implementations in the time-domain scheme were discussed, concentrating on the multiple coupled-motion modes. Such transformations from the frequency domain to the time domain, and the motion coupling would make the implementation of the time-domain scheme much more difficult. As a result, to ensure the correctness of the time-domain scheme, a validation method was presented. In this research, the following conclusions can be drawn:
section_name: Conclusions

section_content: A discussion of how the transformation from the frequency domain to the time domain can be made, and a direct transformation is possible, but not very useful due to its inherent limitations;
section_name: •

section_content: The method for calculating the impulse function and the added mass at infinite frequency based on the frequency-domain prediction was presented, and a MATLAB function was presented for reference;
section_name: •

section_content: The approximation of the impulse function using the Prony approximation method was introduced, and a comparison with the results from WAMIT F2T were made to ensure that the calculation method was reliable; • A simple recursive method for calculating the memory-effect based on the Prony approximation was introduced and the results were validated for its accuracy; • A validation method for the time-domain implementation was explained, which can For the yaw motion (x 6 ), the results from the time-domain and frequency-domain analyses are very different. However, the yaw motions are so small in magnitude (~10 -9 ), and they can be fully ignored in the analysis.
section_name: •

section_content: The method for calculating the impulse function and the added mass at infinite frequency based on the frequency-domain prediction was presented, and a MATLAB function was presented for reference;

• The approximation of the impulse function using the Prony approximation method was introduced, and a comparison with the results from WAMIT F2T were made to ensure that the calculation method was reliable; • A simple recursive method for calculating the memory-effect based on the Prony approximation was introduced and the results were validated for its accuracy; • A validation method for the time-domain implementation was explained, which can be used for ensuring the correctness of the time-domain analysis;

Having examined the correctness of the time-domain scheme, as shown in this research, the future work would implement the PTO systems in the time-domain method for the TALOS wave-energy converter.
section_name: •

section_content: Funding: This research was funded by Engineering and Physics Science Research Council (EPSRC), grant No. EP/V040561/1. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable.

section_name: Data Availability Statement: Not applicable.

section_content: Acknowledgments: This work was supported by the EPSRC (grant number EP/V040561/1) for the project Novel High Performance Wave Energy Converters with advanced control, reliability and survivability systems through machine-learning forecasting (NHP-WEC).
section_name: Data Availability Statement: Not applicable.

section_name: Conflicts of Interest:

section_content: Author Contributions: Conceptualization and methodology, W. S. and G. A. ; formal analysis and investigation, W. S. ; writing-original draft preparation, W. S. ; writing-review and editing, E. T. ; X. M. ; C. J. T. ; R. D. ; D. R. P. and G. A. ; funding acquisition, G. A. All authors have read and agreed to the published version of the manuscript.

section_content: There are no conflict of interest.
section_name: Conflicts of Interest:

section_content: A function for function and added mass at infinite frequency function [K,A_inf] = F2T(w, A, B, t) % inputs: % A, B: the added-mass and damping coefficients % w: the frequency in the frequency-domain analysis % t: the time vector for impulse function % Outputs: % K: impulse function % A_inf: the added mass at infinite frequency %====================================================== n = length(t); dt = t(2)t(1); nw = length(w); dw = w(2)w(1); K = zeros(n,1); for j = 1:n for k = 1:nw K(j,1) = K(j,1) + 2. 0/pi()*B(k)*cos(w(k)*t(j))*dw; end; end; IT = zeros(nw,1); Ma = 0; for k = 1:nw for j = 1:n IT(k,1) = IT(k,1) + K(j)*sin(w(k)*t(j))*dt; end; IT(k,1) = A(k) + IT(k,1)/w(k); Ma = Ma + IT(k,1); end;
section_name: Appendix A