Wave-frequency and low-frequency motions of a deep-draft spar buoy in irregular waves based on a consistent second-order theory

This study investigates the effect of horizontal low-frequency (LF) displacements and velocities on the responses of floating structures in irregular waves, with a focus on a deep-draft spar buoy. It employs a consistent second-order hydrodynamic model in the time domain and its two variations to assess the influence of LF motions. To accurately estimate LF velocities from total velocities, an improved wavelet low-pass filter is proposed and applied, overcoming challenges posed by one-sided data from simulations. For the examined spar buoy, the incorporation of LF displacements and velocities in the seakeeping analysis is deemed essential, as they are found to contribute to the reduction of slowly varying, and consequently, total surge and pitch responses. The study finds that the standard deviations of LF surge and pitch motions scale with significant wave height 𝐻 𝑠 as ( 𝐻 𝑠 ) 𝑏 , with 𝑏 close to 4∕3 for surge motion, highlighting viscous damping as the dominating damping mechanism. For the storm seastate, third-order viscous drag loads may also have made a non-negligible contribution. Meanwhile, pitch motion exhibits 𝑏 values around 1.1 for both low and storm seastates, indicating the significance of both wave-drift and viscous damping in LF pitch motion prediction. The effects of mooring stiffness are also discussed.


Introduction
The analyses of slowly-varying LF loads and motions have received great interest in recent years, especially after mooring failures suffered by some moored structures [e.g., see [1][2][3].Since the leading-order nonlinearity that governs the LF responses is second-order in wave amplitude, the second-order potential-flow weakly-nonlinear theory is of preferred interest.This theory can be enhanced by empirical correction to account for viscous effects based on the drag term of Morison's equation, making it well-suited for addressing the analyses of slowly-varying LF loads and motions in engineering applications.
To address the modeling of LF responses, various second-order weakly-nonlinear methods have been developed.These methods primarily fall into three categories: frequency-domain method [4], indirect time-domain method based on Volterra series expansion [5] using inputs from the frequency-domain method, and the direct time-domain method which directly solves the first-and second-order initial-boundary-value problems in the time domain [e.g., see 6,7].The first two methods require pre-calculated difference-frequency Quadratic Transfer Functions (QTFs) in the frequency domain.In order to acquire the full difference-frequency QTFs, the second-order velocity potential needs to be solved, which is intricate owing to the in-homogeneous forcing terms in the second-order free-surface conditions.Consequently, numerous approximations of full QTFs have been put forward.Some wellknown examples are: Newman approximation [8], Pinkster approximation [9] and full QTFs analysis neglecting free-surface integral [10].These methods can yield accurate evaluation of QTFs for deep-water conditions, provided the natural frequencies of the LF motions are small.Otherwise, the evaluation of QTFs from these approximations become less reliable when the second-order velocity potential has a significant contribution in the circumstances of, for instance, shallow water and increasing difference frequency [e.g., see [11][12][13].Thus, exact solutions to complete second-order problem, [e.g., see [14][15][16] are still pursued in engineering practice.Significant progresses have been made over the last years, and it is fair to say that the full QTFs can be readily obtained from the state-of-the-art radiation/diffraction solvers, e.g.WAMIT, WADAM and HydroStar.However, challenges still exist, which will be summarized as follows.
Firstly, the LF horizontal motions can exhibit comparable or larger magnitude than the corresponding wave-frequency (WF) motions for moored structures in storm seas.This violates the implied assumption of perturbation in the traditional weakly-nonlinear models formulated in the Inertial Coordinate System (IneCS).On the other hand, with these large LF horizontal motions, it becomes theoretically inconsistent to apply Taylor expansion with respect to the equilibrium of a structure.These inconsistencies inevitably diminish the accuracy of both first-and second-order wave loads and motions of floating structures.
Secondly, it is challenging to adequately account for the effects of LF horizontal motions on the QTFs in a frequency-domain or an indirect time-domain method as they solve the QTFs and LF motions separately.Whereas, these effects may be non-negligible.For instance, large LF horizontal motions and LF velocities can alter WF responses and thus the QTFs [e.g., see 17,18].The wave-current interaction can have substantial impact on both equilibrium and dynamic motions of the floating units under certain directional wave and current loading cases [19].Analogous effects of wave-current interaction are expected when one considers a LF velocity as an opposite current and ignores the wave-current coupling.Besides, Kinoshita and Takaiwa [20],Kinoshita et al. [21,22] found that the wave-drift added mass, which is similar to the wave-drift damping to be quadratically dependent on wave amplitudes, has non-negligible effects on the natural frequencies of the horizontal motion of floating structures.Further, the nonlinear stiffness of the mooring system varies depending on the horizontal excursion of the floater, and thus can also alter the resonant frequencies of the system [e.g., see 6,23].So far, rather limited improvement has been made among these mainstream engineering tools to consistently address the above effects of LF displacements, velocity and acceleration on the WF wave loads and in return on the second-order LF wave loads, except for the inclusion of wave-drift damping to reduce the predicted LF motions and mooring tension [24].
Thirdly, it is not straightforward to include other nonlinear effects in the frequency domain, for instance, viscous drag loads, though they may contribute to the second-order LF loads by the relevant quadratic terms (e.g., see; [25][26][27] and cubic terms arising from the splash zone of the structures [e.g., see 23,28].Further, viscous drag loads can influence the WF motions and relative wave elevation, and consequently the LF loads [1,29]. The direct time-domain methods are more powerful to handle these challenges in the frequency-domain and indirect time-domain methods based on the traditional IneCS models.Teng et al. [6] was perhaps the first in the literature to introduce the time-varying LF displacements into the IneCS model by Taylor-expanding the boundary conditions at an instantaneous mean position of the structure.The cost of this approach is to approximate the variables (velocity potential and wave elevation) at the instantaneous position by a Taylor expansion, which however, involves additional numerical approximation of second-derivative terms.Alternatively, Shao et al. [7] proposed a consistent second-order model formulated in the body-fixed coordinate system (BFCS) in the time domain, which can consistently include the time-variant LF displacements and velocities in the seakeeping analysis.This model theoretically eliminates the difficulties associated with large horizontal LF motions without requiring the utilization of twice Taylor expansion in the horizontal plane.This model was an extension of the previous work by Shao and Faltinsen [30,31,32], where the BFCS formulation was proposed and applied to theoretically avoid higher derivatives on the body surface, e.g. the so-called m-terms.
To date, Teng et al. [6] and Shao et al. [7] have only focused on moored monopiles in regular or bi-chromatic waves.In this paper, the consistent second-order model by Shao et al. [7] is extended to study a moored deep-draft spar buoy under realistic irregular waves.Since this model requires instantaneous LF displacements and velocities as inputs to the seakeeping analysis, an imminent challenge needs to be resolved before further implementation of the irregular wave simulation in the time domain, that is how to accurately estimate the LF components from the total motions of broadband for moored structures exposed to irregular waves.In general, standard low-pass filters are designed for central-type stencils where data is available on both sides of the to-be-filtered point.For one-sided stencils, as they exist in our application that only the current and historical information is available from simulations, standard filters may produce unreliable results, for instance, errors of amplitude and phase shift in the time domain due to end effects.In this paper, an improved low-pass filter based on the wavelet transform, a re-sampling strategy using larger time steps, and a symmetric data padding, will be presented to address this challenge.Verification of this filter will be performed by considering manufactured signals before it is applied in the simulation of the moored spar buoy in irregular waves.
Using the consistent models, Teng et al. [6] and Shao et al. [7] have demonstrated that the LF motions have substantial influences on wave loads and motions responses for a spar/monopile in bi-chromatic waves.A very relevant question is: how does the traditional IneCS model perform if more realistic seastate conditions are considered, even though it is now commonly understood that the traditional models have inconsistencies in dealing with large horizontal motions?On the other hand, the individual effects of the LF displacements and LF velocities are not yet fully understood for moored floating structures with large horizontal motions.To address these questions, we employ the consistent second-order model by Shao et al. [7] and two variations for comparative studies in which the LF displacement and velocity effects are either completely ignored or partially included.To simplify the analysis and focus on the hydrodynamic aspects, a quasi-static approach is adopted to model the mooring system as linear springs, whose stiffness is applied to provide realistic natural periods for surge/sway and roll/pitch motions.Besides, the present study accounts for the viscous loads empirically by the drag term in the Morison's equation.The role of drag loads and natural periods (alternatively the mooring stiffness) on the floater responses are also investigated.Furthermore, we will examine the statistical dependency of motions on the significant wave height   , focusing on the mean and standard deviation of surge and pitch responses.The rest of the paper is organized as follows.Section 2 presents the first-and second-order Boundary Value Problems (BVPs) in the consistent model formulated in BFCS, followed by the description of numerical methods in Section 3. Section 4 presents preliminary verification of the numerical model using a bottom-mounted vertical cylinder with or without forward speeds in monochromatic and bi-chromatic waves, where comparisons are made with other numerical and semi-analytical results.In Section 5, a deep-draft spar buoy is introduced and a further verification is made.Afterwards, WF and LF responses of the spar buoy in irregular waves are investigated in Section 6, where three hydrodynamic models, two mooring systems, numerous seastates and two empirical drag coefficients are considered.Section 7 summarizes the present work.

Governing equations and coordinate systems
The velocity potential of the fluid particle  satisfies the Laplace equation with the assumptions of inviscid and incompressible fluid, and irrotational flow: In the complete and consistent second-order model formulated in BFCS, three Cartesian coordinate references are introduced, as displayed in Fig. 1.An Earth-fixed reference frame         is used to describe the incident wave fields.  is a reference frame moving with steady or slowly-varying LF motions in the horizontal plane, and  is the BFCS moving with rigid-body motions in six Degrees of Freedom (6-DOFs) relative to  .These reference frames initially coincide with each other at the initial equilibrium of the structure, with their third axes pointing vertically upward and origins at the still-water plane.        and  are always inertial and non-inertial, respectively.  system is non-inertial with the presence of slowly-varying LF velocities.Unless otherwise proclaimed, all variables are defined in BFCS in this paper.

BVP with steady or slowly varying velocities
The velocity potential on the instantaneous free surface (, , ) satisfies the fully nonlinear kinematic and dynamic boundary conditions in BFCS [33]: =  + W denotes the forward speed on the free surface, with  and W as the steady and LF velocities, respectively.
are the translatory and angular rigid-body velocities, respectively. is the position vector and   = − ⋅  its corresponding gravitational potential, with  as the gravitational acceleration.∇ and ∇ are the three-dimensional (3D) and two-dimensional (2D) gradient operators, respectively.The subscripts , ,  and  denote the partial derivatives of the variables.
In the consistent BFCS model,  in Eqs. ( 2) and ( 3) are decomposed into the following forms: Here η denotes the summation of the incident and scattered waves observed in BFCS.η is the relative wave elevation induced by the vertical motions. =  is wave steepness, indicating the smallness of variables, e.g., the free surface elevation and the rigid-body motions.
The velocity potential on the instantaneous wetted body surface   satisfies the no-flow condition: where  is the positive outward normal vector of   , which is time-independent in BFCS.
Fig. 2 shows a flowchart of the time-domain solution based on the consistent BFCS model.The numerical models in the brackets will be demonstrated in the following sections.

Time-integration of free-surface conditions
Eqs. ( 6) and (7) need to be time-integrated to provide the required velocity potential on the free surface for the BIEs in Eq. (15).Whereas, in the presence of forward speeds such that the convective terms  (0) − ∇ (0) exist on the left-hand sides of Eqs. ( 6) and (7), it becomes less stable for the time-integration of free-surface conditions.Thus, a recently proposed explicit time integration scheme approximated from the implicit schemes for the convective-type equations [7,37] will be applied.When the forward speeds are involved, Eqs. ( 6) and ( 7) can be treated as the convective-type equations in the following generic form: where  = ∇ (0) −  (0) denotes the convective velocity with an amplitude of   . is the unknown wave elevation or potential, and  (, , ) is the forcing term in conformity with that in Eqs. ( 6) and (7).By applying an weighting factor for the convective terms, Eq. ( 16) can be discrsetized as follows: where the superscripts  and  + 1 are the time step indices. denotes the time interval. is a weight factor, with  = 0, 0.5, 1.0 representing the explicit Euler, implicit Crank-Nicolson and implicit Euler schemes, respectively. = ⋅∇   is the directional derivative operator (matrix).Without losing generality, Eq. ( 17) can be re-organized in the form of: where  =  +    ℎ ℎ, and  =  − (1 − ), with  as an identity matrix, ℎ the characteristic mesh size.The factor  =   ∕ℎ is the same as the CFL number.Considering the CFL number (or ) as a small parameter, the inverse of  (denoted as  −1 ) satisfying  −1 =  can be expanded into a series expansion in  as follows: After truncating the above infinite summation in Eq. ( 19) with only the first   terms retained and applying it into Eq.( 17), the final expression of the variable in Eq. ( 17) can be explicitly given as: In this paper,   = 6 is used for all the time-domain analyses in Sections 4-6.Details of the methodology and verification of the scheme are referred to Shao et al. [7].

Low-pass filter close to the waterline
To suppress the short-wave instabilities on the free surface and the waterline in the time-domain analyses, particularly in the presence of steady and LF velocities that rapid variation of wave elevation close to the waterline exists, an optimized low-pass filter based on weighted least squares (WLS) recently proposed by Shao et al. [7] is applied.
The commonly used WLS filter [38] can be written as follows: + (21) where  denotes the th stencil center and   the stencil point number for the WLS filter.Here   with  =  − 1 2 (  − 1), … ,  + 1 2 (  − 1) are the stencil values and   their convolution coefficients. +  is considered as the order of the polynomial.A Gaussian function is introduced to determine the weighting coefficients of the WLS, which is expressed as follows: where  0 indicates the stencil center for filtering. 0 is twice the local element size, which is used to define the shape of the Gaussian function together with . 0 is estimated as: Thus, the weighting coefficients for optimized WLS are given as follows The trait of the optimized WLS filter is that the Gaussian function in Eq. ( 22) has been elaborately designed such that the transfer function, or in other words the remaining energy at  0 is sufficiently small.In this case, the short-wave instabilities can be minimized.The optimized WLS filters have been implemented for both 1D and 2D stencil with satisfactory results [7].In this paper, a 2D optimized WLS filter with 45 stencil nodes and 5th order polynomial will be applied for the velocity potential and wave elevation on the free surface.For each free-surface node, three surrounding layers of nodes are included as stencil nodes.The 1D filters are only applied for waterline points.Besides, the following mild filter modified from Eq. ( 21) are applied at each time step: + (25) where  ∈ [0, 1] is a constant typically much smaller than 1.0. approximately from 0.05 to 0.1 for the first 2 or 3 layers of points adjacent to the waterline and  from 0.1 to 0.3 for points away from the waterline, will be applied for the time-domain simulations in the following sections.

Low-pass filter based on wavelet transform
In this section, an improved low-pass filter based on wavelet transform, which involves a symmetric padding and a re-sampling of the raw signal for the improvement of both filtering accuracy and efficiency has been proposed to instantaneously extract the LF velocities.The commonly used low-pass filter based on wavelet transform will be briefly introduced in the first place.Interested readers are referred to Torrence and Compo [39] for more details.
For a structure under large horizontal LF motions, i.e. surge, sway and yaw motions, the low-pass filtering can be implemented separately for each motion.Let us take surge motion or velocity as an example.By choosing a discrete set of equally spaced time series: the continuous wavelet transform is expressed as the convolution of  with scaled and translated wavelet basis function : Here   ,  = 0, … ,   , denotes the wavelet scales.  , which will be defined in Eq. ( 34), is an integer determining the largest wavelet scale. is the time interval of .The superscript * denotes the complex conjugate.  (  ) represents the wavelet coefficient at the th wavelet scale for   in the data set .In convenience of directly comparing the results from wavelet transform at each scale, the wavelet function is normalized by a mother wavelet  0 who has unit energy in the Fourier space.The normalization can be written in the following forms: where ,  = 0, … ,  − 1 and  = 0, … ,   .Alternatively, the wavelet transform for  can also be done in the Fourier space by inverse Fourier transform of the products of x and ψ (     ) according to the convolution theorem, Z. Zheng et al.
Note that here and hereafter i denotes a unit pure complex.x and ψ (     ) are the DFT of time series  and wavelet function , respectively, which are given as follows: where  = 0, 1, … ,  − 1.The angular frequency   is defined as: ,  = 0, … ,   , needs to be properly chosen to accomplish the continuous wavelet transform.For convenience,   for  = 0, 1, … ,   is generally written as fractional power of two: Here  0 is the smallest scale, with the corresponding Fourier period of approximately 2. should be sufficiently small and depends on the width in spectral space of the wavelet function.  can be evaluated as follows: Once the wavelet coefficients at each wavelet scale are obtained, it is straightforward to reconstruct the original time series using either deconvolution or the inverse filter.For the continuous wavelet transform, it also allows to reconstruct the time series using a different wavelet function.Herein, a  function [40] is employed for this purpose, which gives: where   is a constant emerged from the reconstruction of a  function that has utilized the wavelet function  0 () for wavelet transform.Re{⋅} is an operator to take the real part in the case that the original time series contain only the real part.
For the present application of moored offshore structure with LF motions in irregular waves, it is unnecessary to do wavelet transform at all scales, as the only interested scale levels are those close to the natural periods of horizontal motions, which are generally in the magnitude of 1 to 2 min.Hence, we only need to perform wavelet transform and reconstruction for  at larger scales.By defining a proper truncation scale  cut and the largest scale  max , among which the filtered wavelet scales are distributed based on Eq. ( 33), a low-pass filter based on wavelet transform can be designed.
Appropriate wavelets should be chosen for the low-pass filter based on wavelet transform, which significantly depends on the characteristics of the time series.In this case, a nonorthogonal complex wavelets of Morlet wavelet basis has been adopted to properly preserve the smooth and continuous variations of the time series of LF velocities in both amplitude and phase.The applied Morlet wavelets is expressed as: and its expression at   ,  = 1, … ,   , in Fourier space: where  0 is given as  0 = 6.0 rad/s for Morlet wavelets.() is Heaviside step function, with () = 1 for  > 0 and () = 0 elsewise.For Morlet wavelets,   = 0.776 and  0 =  − 1 4 are defined for reconstruction in Eq. ( 35).

Reducing the end effects of the low-pass filter
A standard low-pass filter based on wavelet transform can be constructed with the prescribed methodology and have been successfully implemented to filter the LF displacements from the total motions in the time-domain analyses of, for instance, bichromatic waves as presented in Teng et al. [6].In general, wavelet-based filters, or other low-pass filters based on, for instance, Fourier Transform or weight least square [38], exhibit optimal filtering performance for two-sided stencils.However, for one-sided stencils, as they behave in the present application of extracting the instantaneous LF velocities that only the current and historical information is available from simulations, these standard filters may yield unreliable results, for instance, errors of amplitude and phase shift for points near the two ends where the filtering results we acquire.This issue becomes particularly pronounced when dealing with finite-length non-periodic time series [39].
A strategy is proposed in this paper to effectively reduce the end effects of the wavelet low-pass filter, which involves two steps: (1) re-sampling the raw data at a larger time interval, and (2) performing symmetric padding on both ends of the re-sampled data set before wavelet transform and remove them afterward, and modification of the formulations of wavelet transform and reconstruction.
As the first step, the temporal results of horizontal motions are re-sampled using a much larger time interval  multiple of the time step  used in the numerical simulations.The re-sampling of a data set can be easily implemented by an interpolation of Fig. 3. Re-sampling for a manufactured signal  = 0.5 cos(1.0)+ 2.0 cos(0.1).The total number of nodes are reduced from 800 with 40 nodes per WF period to 200 after re-sampling.100 nodes per LF period is used for re-sampling.the time series using, for instance, a higher order B-spline fitting. is normally very small as it is selected in the simulation such that the WF and sum-frequency oscillations are necessarily modeled in the case of irregular waves.In contrast,  is only required to be able to capture the LF band information, in spite of the WF and high-frequency responses.In fact, re-sampling with a larger time interval, serving as a low-frequency bandpass, has already filtered out some unnecessary parts of the WF components without undermining the energy of LF components.Another benefit of the re-sampling is that the resulting number of points in the data set is significantly reduced, thus improving the efficiency of the subsequent filtering based on wavelet transform.An example of re-sampling for a manufactured signal  () = 0.5 cos(1.0) + 2.0 cos(0.1 ) is given in Fig. 3. Within two LF periods, the original data set contains 20 WF periods, and the total number of raw data is  = 800.A re-sampling based on 50 points per LF period reduces the total number of data points to  re = 100.
As the second step, if the re-sampled data set is: where  re denotes the number of re-sampled data set, the symmetric padding performed for the re-sampled data can be expressed as: Modification of the formulations of wavelet-transform in Eq. ( 27) and reconstruction in Eq. ( 35) is also proposed here: Here   is a modification coefficient.  = 1.4 is used in later verification of the proposed improved wavelet low-pass filter and timedomain simulations.Besides, one should note that once re-sampling is applied, which means the data to be filtered has changed, the parameters in Eqs. ( 40) and ( 41), for instance, , , etc., should be adjusted accordingly.Two manufactured signals are considered to verify the performance of the improved wavelet low-pass filter before it is further applied in the time-domain hydrodynamic analyses in the next sections.The first signal is defined as  1 () = 0.2 cos() + 0.3 cos(1.1 ) + 1.5 sin(0.1 ), comprising two WF components and one LF component.The frequency of LF component is equivalent to the difference of the two WF frequencies.The second signal  2 () is generated from a JONSWAP wave spectrum with a peak period   = 12 s and a significant wave height   = 10 m, superposed to a sinusoidal signal with a period of 60 s.
Both the low-pass filters based on standard wavelet transform in Section 3.3.1 and the improved strategy as discussed in Section 3.3.2are adopted.In order to be aligned with the requirement of the time-domain simulations in the next sections, only a one-sided stencil participates for filtering at each time step.Besides, the choice of the time window length of the stencil and thus the number of data points is arbitrary.However, based on our experiences, a time window of at least one LF period is needed to give satisfactory filtered results.Moreover, a proper truncation scale  cut should be specified in the reconstruction of the LF signal using Eq. ( 35) for both filters.We have chosen  cut with the corresponding Fourier periods 30 s and 40 s for the two tested signals  1 () and  2 (), respectively.In particular, for the improved wavelet low-pass filter, 100 points per LF period are used to re-sample the raw data so that the LF component of the signal is practically not influenced by the re-sampling.
Comparisons of the time series and Power Spectrum Density (PSD) for the raw signal and filtered results are displayed in Figs.   with the improved wavelet low-pass filter, the majority of the high-frequency components are eliminated while the LF components are preserved.In contrast, the performance of the standard wavelet low-pass filter is far from satisfactory, as a noticeable decrease in amplitude and phase shift of LF components is observed from the filtered results compared to the original LF components.
Performance of the standard and the improved wavelet low-pass filters can also seen from Figs. 4(b) and 5(b) by comparing the PSDs of the raw signal and the filtered results for each time series.The energy is mostly preserved in the LF range and effectively eliminated in the high-frequency range by the improved wavelet filter, while the standard wavelet filter has removed a significant amount of energy in the LF range.In general, the improved low-pass filter based on wavelet transform performs well in extracting the instantaneous LF components in both amplitude and phase, even though only a one-sided stencil participated, demonstrating its potential to approximate the LF motions instantaneously in the time domain.

Preliminary verification of the numerical model
The higher-order BEM solver adopted in this study has been extensively verified and validated in previous studies [7,32].In this section, we will only briefly present two verification studies, both using a bottom-mounted cylinder with a radius of  and water depth ℎ = .The first verification concerns the second-order mean and sum-frequency loads on the cylinder with and without a forward speed.The second case examines the second-order difference-frequency and sum-frequency wave loads in bi-chromatic waves.
Figs. 6(a) and 6(b) show the comparisons of the mean-drift forces and sum-frequency forces based on the present second-order model and the corresponding results from Cheung et al. [41] and Skourup et al. [42], respectively.Three different Froude numbers   =  ∕ √  = −0.1,0.0, 0.1 have been studied.Accurate solutions of mean-drift and sum-frequency forces are obtained using the present second-order model.
The same bottom-mounted cylinder, fixed in bi-chromatic waves whose parameters are given in Table 1, is also studied.The sum-and difference-frequency QTFs of the wave loads in surge are compared to the semi-analytical solutions by Cong et al. [43] in Fig. 7, with very good agreement.In the figure, the nondimensional  is defined as  =  2 ∕ and the subscripts  and  denote the two wave components of the bi-chromatic waves.7. The difference-frequency and sum-frequency QTFs of the second-order wave forces in surge for a bottom-mounted cylinder.The draft of the cylinder ℎ is equal to its radius . and  denote the two wave components of the bi-chromatic incident waves.The wave parameters are defined in Table 1.

A floating deep-draft spar: description and a further verification
In this section, the floating deep-draft spar buoy and its loading condition that will be studied in the rest of the paper are presented.To gain further confidence of the numerical model, a further verification of the second-order difference-frequency wave excitation loads on the spar buoy, either fixed or freely floating in bi-chromatic waves, will be carried out.In principle and within the context of the second-order theory, a more complicated scenario involving irregular waves can be viewed as an amalgamation of multiple pairs of bichromatic waves.The same floating spar buoy exposed to irregular waves will be further studied in Section 6.

Geometry and loading condition
A deep-draft spar with dimensions and loading conditions given in Drake [44] is considered.It has a radius  = 12.5 m and a draught  = 200 m.The considered water depth is 1000 m, representing a rather deep-water condition.The gravity and buoyancy centers are 105 m and 100 m below the still-water level, respectively.The main particulars of the structure and the loading condition are summarized in Table 2.The gyration radius for roll and pitch is defined with respect to the center of gravity.

A further verification in bi-chromatic waves
The spar, either fixed or undergoing 6-DOF rigid-body motions in bi-chromatic waves, which is the simplest form of irregular waves, will be studied as a further verification of the present model.The difference-frequency full QTFs, including contributions  The mesh model for the considered spar is so designed that, the overall computational domain of the free surface is a circular disk truncated at a radius of 600 m, corresponding to at least 2.4 times the longest linear wave component among the investigated waves.An illustration of the meshes is displayed in Fig. 8.The free surface mesh consists of two annular domains, namely, an inner domain where meshes are refined to better represent the sum-frequency waves with shorter wavelengths, and an outer domain having relatively coarser meshes with larger stretch rate along the radial direction.For all of the investigated bi-chromatic waves, the same mesh model with an inner domain of a radius of 90 m is used while the damping zone may be different.A larger wave generation zone will be applied for longer bi-chromatic waves.Overall, the damping zone starts from one time the average wavelength of the bi-chromatic waves, i.e. 0.5 ( 1 +  2 ), and ends at the free surface boundary.Here  1 and  2 are the linear wavelengths of the two components of the bi-chromatic waves.The coefficients  and  in Eq. ( 9) are chosen as  = 1.1 and  = 2 for all simulations.In total, 1856 cubic elements (9184 nodes) on the free surface and 1047 cubic elements (5188 nodes) on the mean wetted body surface are used for spatial discretization.In the near field of the body surface, at least 32 elements per linear wavelength are used to ensure a sufficiently fine resolution.A time step of  =  2 ∕400 is employed to give convergent QTFs for both sum-and difference-frequency results, where  2 represents the smaller wave period of the two wave components.Given that the primary focus of this study is on the hydrodynamics of moored structures, with particular interest in the LF components, only these results of difference-frequency QTFs will be presented.
Fig. 9 shows the comparison of the difference-frequency full QTFs between the present model and a frequency-domain analytical model [24].Results for the fixed and floating conditions are presented in Figs.9(a) and 9(b), respectively.As seen from the excellent concurrence in the comparisons, the present consistent second-order model can provide very accurate second-order difference-frequency wave loads for both fixed and floating deep-draft spar.

Numerical analysis of the floating deep-draft spar in irregular waves
In this section, we will continue to investigate the same spar buoy under the loading condition in Table 2.However, to better represent the real ocean environment and wave loads, irregular waves will be considered.The drag loads will also be empirically  accounted for as part of the hydrodynamic loads.Three different hydrodynamic models are tailor-made and compared with each other to investigate the individual effects of LF displacements and velocities.In total, our numerical analysis encompasses 30 simulations, involving a combination of 8 seastates, 2 mooring stiffness, 2 drag coefficients, 3 hydrodynamic models.Each simulation has a duration of 5400 s, while the first 1800 s are discarded to eliminate the transient effect resulting from initial conditions.Consequently, only the results from the final 1 h are utilized in subsequent statistical analyses.These numerical results provide insight into the effects of LF displacements and velocities, viscous drag loads, mooring stiffness, wave steepness and peak periods of the seastates, which will be discussed in the later subsections.

Seastates
The JONSWAP wave spectrum with a peak-enhancement factor =2.4 is applied in the present study.Eight seastates, categorized based on two different peak periods   and four different wave steepness are chosen to represent the irregular waves the spar buoy may encounter under operational or survival conditions.Here  is defined as  =     in the case of irregular waves and used as an indicator of the steepness of the waves, with   as the significant wave height.The parameters for these seastates are detailed in Table 4.
For the incident wave spectrum discretization, an equal energy method is employed as it enables to effectively model the wave spectrum with a relatively small number of frequency components, in contrast to arithmetic progression or geometric progression [45].In this method, the frequency components are selected so that each frequency component represents an equal amount of spectral energy.Additionally, appropriate lower and upper cut-off frequencies are determined to ensure negligible energy outside this range.For each seastate, the number of wave frequency components are chosen as 100, and the phase of each wave frequency component is selected by a random number in a range of [0, 2].Taking Seastate 2, the largest seastate to be investigated as an example, the theoretical description, denoted as 'Target' and the corresponding discretization based on an equal energy method, denoted as 'Modeled' are displayed in Fig. 10.

Two mooring stiffness
To effectively model the dynamic responses of wave-frequency and low-frequency motions of the spar buoy while preventing excessive drifting, we employ soft linear springs to approximate the quasi-static restoring forces of the mooring system.These springs are utilized for both first-and second-order motions of the spar.The choice of two distinct mooring systems, each leading to different natural surge and pitch periods, is a crucial aspect of our study.These natural periods play a pivotal role in governing the system's  responses to external forces.Table 5 provides a summary of the respective natural periods of the surge and pitch motions for the two mooring systems.These values have been estimated through single-DOF free-decay tests conducted in still water using our time-domain numerical model.Importantly, the first mooring system results in both larger surge and pitch natural periods compared to the second mooring system.This deliberate variation in natural periods allows us to investigate how the effects of large-amplitude LF motions and the influence of drag loads change under different dynamic conditions.
In a later section of our analysis, we will delve into these variations to gain insights into the behavior of the spar buoy under different mooring configurations.By doing so, we aim to uncover the intricate interplay between LF and WF motions, the influence of mooring stiffness, and the impact of drag loads on the system's response.

Three hydrodynamic models
To assess the significance of incorporating LF displacements and velocities in the first-and second-order seakeeping models, we will apply and compare three hydrodynamic models.The key characteristics of these models are summarized in Table 6.
• 'Model 1' does not consider either LF displacements or LF velocities.
By examining these three models, we aim to gain insights into the impact of LF motion components on the seakeeping behavior of the system.
Here,  =  −  0 represents the motion vector of a point fixed on the -plane of BFCS due to rigid-body motions of the structure. () (where  = 1, 2) is the th order component of .If we approximate  (1) ( 0 +  (1) , 0, ) −  (1) ( 0 , 0, ) by  (1) ⋅ ∇ (1) ( 0 , 0, ) using Taylor expansion, 'Model 1' is theoretically equivalent to the traditional model formulated in IneCS.This traditional model is commonly employed in existing state-of-the-art hydrodynamic tools.However, it is worth noting that for a moored structure undergoing large-amplitude motions, theoretical inconsistencies may arise when applying the above two equations or a traditional IneCS model.Firstly, the truncation of higher-order terms in  assumes that || is asymptotically small compared to the characteristic dimension of the structure or the wavelengths.Furthermore, || 2 is assumed to be one order of magnitude smaller than ||.These two assumptions can be easily violated for structures moored by a compliant mooring system, as LF motion responses in the horizontal plane can become comparable to the dimensions of the structures and much larger than the corresponding wave-frequency (WF) components.
Different from 'Model 1' or 'Model 2', 'Model 3' incorporates the LF velocities ẋLF as input for both first-and second-order seakeeping analyses.In this study, ẋLF are obtained by applying the improved wavelet low-pass filters in Section 3.3 to the total velocities, which contain frequency components   ,   +   and |  −   |.Here   and   represent the frequencies of the th and th wave components respectively in irregular waves.In practical application, the improved wavelet low-pass filters are employed for both surge and pitch motions.Subsequently, the LF motions at any location, induced by the surge and pitch motions can be determined through geometrical relationship.

Estimation of drag loads 6.4.1. 'Independent flow field' form of morison equation
For a structure under both WF and LF responses, it is a common practice to consider the viscous effects as empirical corrections to potential-flow solutions.This correction is typically implemented using what is known as the 'independent flow field' form of the Morison equation, as detailed in references such as DNV [46], which separates the drag loads arising from WF and LF motions.It is still debatable how the Morison equation should be adopted to model the viscous-drag induced slowly varying hydrodynamic loads.However, as suggested in Molin [47], the 'independent flow field' formulation might be preferred from engineering point view, because it provides less viscous damping and thus is relatively more conservative.
In addition, the drag loads in the splash zone, which are third order in the wave amplitude, are included.This consideration is essential as their contribution may not be negligible, particularly in storm seas.As a result, the viscous drag forces and moments on the spar buoy are estimated accordingly: and +  (  4 ) . ( Here,  is the diameter of the spar. ()   (, ) and  ()   (, ) with  = 1, 2, denote the th order ambient-flow speed and rigid-body velocities at a cross-section located at the -coordinate.  () signifies the relative wave elevation observed at still waterline in BFCS. ()   is the viscous drag coefficient at th order.This derivation is based on perturbation scheme and assumes the cross-flow velocity in the splash zone remains constant.For more detailed information, please refer to Faltinsen [28].

Determination of drag coefficient 𝐶 𝑑
In general, the drag coefficient depends on Reynolds number (), Keulegan-Carpenter () values (or a combination of these two factors as  number, defined as  =   ) and surface roughness , etc. Empirical formula for drag coefficients as function of  has been suggested by, for instance [46], at high  values for circular cylinders.For the studied deep-draft spar and seastates, the  numbers are practically very high, e.g.typically in order of 10 6 , and the  numbers remain below 5.0 for both WF and LF motions.This makes it, as discussed in Molin [47], very difficult to justify the selection of drag coefficients, which may take very different values, depending on the actual  and .
Furthermore, both  and  vary along the spar buoy, due to combined surge and pitch motions, making it even more challenging to accurately estimate the time-and location-dependent drag coefficients based on empirical formulas.
Last and not least, the presence of the free surface can further complicate the selection of   coefficients in the splash zone and the estimation of the viscous drag loads.
Therefore, as a very rough estimation, we will in this study simply use the 'independent flow field' form of Morison equation, and apply a constant drag coefficient   along the spar buoy.To examine the sensitivity of motion responses to   and the uncertainties, two different drag coefficients   = 0.2 and 1.0 will be considered in this study.Those two   values are close to the lower and upper bounds, estimated using the empirical formula in DNV [46] and the ranges of actual  and  numbers in the considered simulations.
More accurate and efficient computations of the viscous drag loads are pursued in the future towards the slow-drift motions of moored floating structures in irregular waves.This can be achieved by calibration towards dedicated CFD analysis with proper turbulence modeling, model tests or full-scale measurements, which, however, falls outside the scope of this study.

Setup of numerical model and the low-pass filter
A similar meshing strategy, as applied in the bi-chromatic wave cases in Section 5.2, is employed in simulations in irregular waves.A slight difference is the use of a larger free surface domain with a radius of 1000 m, corresponding to about 3.4 to 7.9 times the linear wavelength   .Here   is estimated by the peak period   of the wave spectrum.To improve the wave resolution near the structure, we employ an annual refinement zone with a radius of 90 m, approximately 0.3   to 0.7   .A damping zone with a width of at least 2.2  is exploited to absorb the outgoing scattered waves.Besides, 40 elements per   are utilized in the near field of the free surface.A time step of  =   ∕400 is used to achieve convergent first-and second-order LF results.In terms of the improved wavelet low-pass filters in time-domain simulation, the truncation scale  cut and the largest scale  max representing the lower and upper limit for band-pass filtering are necessitated, which primarily rely on the wave spectrum and the natural periods of the spar for irregular waves. cut and  max are chosen such that their corresponding Fourier periods are 40 s and 300 s, respectively.Eighty scales are evaluated by Eq. ( 33) while only these from  cut to  max are selected for the applied filters.Besides, a historical time window of total velocities with a length of 15   , roughly equivalent to at least one surge natural period, is chosen, and 100 nodes per surge natural period is applied for re-sampling of the raw data in the improved wavelet low-pass filters.Ultimately, the number of the stencil for the improved wavelet low-pass filters is reduced by at least 90% at each time step, significantly enhancing the efficiency of extracting the LF velocities in the time domain.

Results and discussion
In this section, the WF and LF responses in surge and pitch in irregular waves will be presented and discussed.The significance of incorporating the LF displacements and velocities in the first-and second-order seakeeping analyses is investigated for eight seastates defined in Section 6.3, two mooring systems described in Section 6.2 and two drag coefficients specified in Section 6.4.

Performance of the wavelet low-pass filter
The performance of the improved wavelet low-pass filter is illustrated for both the low seastate (Seastate 1) and the storm seastate (Seastate 2) in Table 4. Figs.11 and 12 show the raw and filtered LF horizontal velocities at the mean waterline (MWL) of the structure using improved wavelet low-pass filters and 'Model 3' for these two seastates.Two different drag coefficients   = 0.2 and 1.0 are considered.The estimated LF horizontal velocities at MWL will be used as input to the seakeeping model as described in Section 2.
In these figures, two types of results are presented, denoted by 'LF Vel.' and 'Total Vel.'.'LF Vel.' represents LF velocity results obtained using the improved low-pass filter, which exclusively employs one-sided data stencils for calculations.'Total Vel.' results are derived through a post-processing step conducted after the entire simulation has been completed, using two-sided data stencils.The 'Total Vel.' results are entirely free of end-effect errors.The difference between the two results represents the error of the improved wavelet low-pass filter.
After a visual inspection of the time series in Figs.11(a) and 11(c) for different   values, it is evident that the majority of the high-frequency components have been effectively filtered out.The corresponding power spectral densities (PSDs) of the time series , which represents the square of the standard deviation (STD), is defined as the total area below a PSD curve, and   ẋ1 represents the error in the LF surge velocity STD,  ẋ1 , before and after applying the filter.For Seastate 1,   ẋ1 is approximately 5.2% for both cases with   = 0.2 and   = 1.0, respectively.The improved wavelet low-pass filter demonstrates good filter performance for Seastate 2, as indicated by the time series and PSDs of the total and filtered LF velocities in Figs.12(a)-12(b) and Figs.12(c)-12(d).In these cases,   ẋ1 is 3.9% and 6.2% for   = 0.2 and   = 1.0, respectively.
Overall, it is concluded that the improved wavelet low-pass filter is capable of extracting the LF velocities from the total velocities with reasonably good accuracy.Nonetheless, it is important to note that filtering one-sided data remains a challenging issue, as previously highlighted in this paper and in BV-NI691 [48] (Section 3.4.3).To the best of the authors' knowledge, there is currently no perfect solution to address this matter.Interested readers can also refer to relevant discussions in BV-NI691 [48].

The effects of LF motions and drag coefficients
Taking Seastate 1 (a low seastate) as defined in Table 4 and employing an empirical drag coefficient of   = 0.2, we present the time series of surge and pitch motions at the center of gravity (COG) of the spar in Figs.13(a)-13(b) for the first-order components and in Figs.13(c)-13(d) for the second-order components.
Comparing the results from the three hydrodynamic models prescribed in Section 6.3, we observe significant LF resonant surge and pitch motions in all three models.Non-negligible LF resonant surge and pitch motions are also evident in the first-order solutions in Figs.13(a)-13(b), particularly in 'Model 3', demonstrating that incorporating LF displacements and velocities in the seakeeping analysis can modify LF responses through changes in the first-order responses.
The total surge and pitch motions at COG are presented in Figs. 14 and 15 for   = 0.2 and   = 1.0, respectively.Total motions represent the summation of the first-and second-order components or, within the context of second-order theory, the combination of WF and LF motions.The statistical mean and STD values for WF, LF, and total (WF+LF) responses using different   and hydrodynamic models are presented in Fig. 16.Based on the comparison of the temporal and statistical results in Figs.13-16, the following observations and discussions can be made: • For the considered spar buoy exposed to Seastate 1 (a low seastate), the STDs of surge and pitch WF responses are significantly smaller than the corresponding LF responses, regardless of the choice of drag coefficient   .In contrast to LF responses, which are significantly affected by the choice of drag coefficient (  ), the WF responses exhibit lower sensitivity to changes in   .• These three hydrodynamic models exhibit more diverse motion responses for the studied spar when a smaller   value of 0.2 is used.For a low seastate, with   = 0.2 applied, 'Model 1,' which neglects both LF displacements and velocities ( LF = 0, ẋLF = 0), predicts larger LF surge and pitch responses compared to 'Model 2' ( LF ≠ 0 and ẋLF = 0).The over-prediction of LF responses becomes even more pronounced when compared to 'Model 3' ( LF ≠ 0 and ẋLF ≠ 0).This can be partly explained by the presence of  wave-drift damping loads induced by the slowly-varying LF velocities of the spar buoy, which are consistently included in 'Model 3' but ignored by 'Model 1' and 'Model 2'.For the considered deep-draft spar buoy, the present numerical results seem to suggest that, the LF displacements and velocities tend to individually reduce the surge and pitch responses of the spar if they are incorporated in the seakeeping analysis.In general, ignoring the LF displacements and velocities may lead to moderately or significantly conservative results, depending on the applied drag coefficients, or more broadly, the relative significance of drag loads compared to the potential-flow hydrodynamic loads.
• The viscous drag loads have played a significant role in the slow-drift LF motions of the spar in the considered seastates.
Specifically, a smaller drag coefficient   = 0.2 results in larger LF responses in both surge and pitch, compared to that for   = 1.0.When using the relative velocity (  −   ) for drag load calculation in Eqs. ( 48) and ( 49), both wave excitation and viscous damping are considered.As the fluid velocities induced by the wave field decay exponentially along the vertical direction, the drag-induced wave excitation on the spar appears primarily in a relatively small layer close to the waterline.On the other hand, the entire submerged portion of the spar buoy can contribute significantly to the viscous damping loads, owing to the rigid-body motions of the spar.Given that the studied spar buoy has a deep draft (see definition in Table 2), it appears that the viscous damping effect has outweighed the viscous-drag excitation.
The above observations and discussions were made based on a specific seastate (Seastate 1 in Table 4) and a particular mooring stiffness (Mooring 1 in Table 5).It is then natural to inquiry if those observations still hold when considering different seastates and mooring stiffness.In the following subsections, the cases combining different seastates in Table 4 and mooring stiffness in Table 5 will be investigated for further analyses of the effects of seastate and mooring stiffness.

The effect of seastates
Since the relevant physical mechanism may be different between low and higher seastates, the subsequent investigation into the seastate effects will be conducted in two stages, encompassing eight seastates defined in Table 4.In the first stage, the results  for Seastates 1 and 2, which share the same largest wave steepness given as  =     = 0.048 and different peak period   will be compared.Note that Seastates 1 and 2 represent the low and storm seastates, respectively.
Figs. 17-18 and 19 show the time series of total motion responses and the statistical results (mean and STDs) for surge and pitch motions at COG under Seastates 2. Results of two drag coefficients and three hydrodynamic models are presented.Through the comparison of results from seastate 1 and 2, the subsequent conclusions can be drawn: • The LF responses are the dominating contribution to the total responses in storm seastate.
• Like the lower seastate, the surge and pitch responses, in particular their LF components, are very sensitive to the drag coefficient   that a larger value (  = 1.0) results in smaller responses.Again, this might be attributed to the fact that the viscous damping has out-weighted the viscous-drag excitation.• Similar to the low seastate, the higher seastate also exhibits discrepancies of STDs for LF and total pitch motions between the three hydrodynamic models under the same drag coefficient, particularly when   = 0.2 is selected.Specifically, 'Model 1', neglecting both the LF displacements and velocities notably overestimates the LF pitch responses than 'Model 2' and 'Model 3'.
It is worth noting that all results in Figs.[13][14][15][16][17][18][19] are presented for the COG of the spar buoy, which is located 105 m below the mean water level.The pitch motions can induce significant surge motions for a deep-draft spar.Fig. 20 shows the statistical values of the surge motion at MWL under the low and storm seastates.With a smaller drag coefficient (  = 0.2), the numerical results show much more significant differences between the three models.In summary, our results suggest that the effects of LF displacements and velocities on both simulated LF and total motion responses are significant and should be carefully incorporated into the seakeeping analysis for the deep-draft spar buoy in both low and storm seastates.
In    .The values of  1 ,  4∕3 and  2 are determined by curve fitting.Based on Newman's approximation for the second-order wave excitation force and the assumption of narrow-band surge LF resonant motions of moored structures, it has been demonstrated by Faltinsen [28] and Molin [47]  As seen from the comparison between the simulated results and the fitted curves, the STDs of LF surge and pitch motions under different   do not consistently follow any of the three curves, except for the STDs of LF surge motions at COG under the lower seastates of   = 9.0 s that appears to be a reasonable approximation, even though neither Newman's approximation for the second-order LF excitation loads nor an assumption of narrow-band dynamic responses at resonance is used.The power of the fitted function for LF surge motion STDs suggests that quadratic damping dominates over other effects, such as wave-drift damping and higher-order viscous damping, under lower seastates with   = 9.0 s.For longer waves with   = 13.8 s,    1 ∼  4∕3  still seems to hold except for the highest seastate with   = 14.4 s, which might be due to increased importance of the third-order viscous drag terms in Eqs. ( 48)- (49).In contrast, the STDs of LF pitch motions, which consistently scale approximately with  1.1   for both   values, are primarily influenced by both wave-drift damping and viscous damping.Similar fitting has been implemented for STDs of WF surge and pitch motions, with comparison of results and linear curve presented in Fig. 22.The numerical results indicate that the WF responses clearly exhibit a linear relationship with   for both surge and pitch motions under a smaller   of 9 s.However, when the peak period increases to   = 13.8 s, the WF responses seem to follow a higher power of   than 1, attributed to the increased significance of the viscous drag loads for longer waves and larger wave heights.

The effect of mooring stiffness
Figs. 23 and 24 compare PSDs and STDs of the total motions, respectively, under storm seastate (Seastate 2 in Table 4).The figures include results from all three hydrodynamic models and two different mooring systems.Note that the PSDs in the WF range are not included in Fig. 23 as the WF motions are much smaller.In general, the mooring system has significant influence on the LF responses that the stiffer mooring (Mooring 2) leads to surge motions at COG with much larger amplitudes and oscillatory frequencies, and thus significantly larger velocities than those by the softer mooring (Mooring 1).
From the surge and pitch PSDs in Figs.23(a) and 23(b) for the soft mooring, the individual influence of LF displacement and velocities can be identified by comparing the three numerical models.Through comparing the pitch PSDs by the three models in Fig. 23(b), the presence of LF velocity are found to lead to a significant reduction of LF pitch responses, while the LF displacement contributed to a much smaller reduction.The pitch induced surge at COG shows a similar reduction effects, as seen from the peaks at larger frequency in Fig. 23(a).Interestingly, 'Model 3' shows a large surge response at COG than the other two numerical models, which might be due to an increase of second-order wave excitation force due to the presence of LF speeds.The surge PSD predicted by 'Model 3' also exhibits a slightly smaller surge natural frequency, which might be due to the presence of wave-drift added mass.However, this effect seems to be very small for the considered deep-draft spar buoy.
Unlike the softer mooring for which the PSD results are shown in Figs.23(a) and 23(b), as one can see from Figs. 23(c) and 23(d) for the stiff mooring that, the three numerical models show very little differences in the predicted LF surge and pitch motions.This is probably due to the dominating role of viscous damping over other damping sources.
With the stiffer mooring, stronger surge-pitch coupling is also observed, evident from the double peaks in the pitch PSD in Fig. 23(d).Conversely, the softer mooring results in pitch PSDs with a relatively smaller contribution from surge coupling, as seen from the corresponding pitch PSD in Fig. 23(b).It is worth noting that the apparent coupled natural frequencies estimated in Fig. 23 differ from the uncoupled ones obtained from single-DOF free-decay tests.
In Fig. 25, the statistical results of the two different mooring systems under the low seastate (Seastate 1) are compared.Again, a drag coefficient   = 0.2 has been applied.The three hydrodynamic models yield very different STDs for both surge and pitch motions, emphasizing the importance of incorporating LF displacements and velocities in the seakeeping analysis under low seastate and softer mooring conditions.

Conclusions and future perspective
This study utilizes a consistent second-order model to analyze a moored deep-draft spar buoy in irregular waves.The main focus is on assessing the importance of including LF displacements and velocities in seakeeping analysis.To enhance accuracy when estimating LF motions from one-sided data, a low-pass filter is proposed.The study also explores the impact of factors like seastates, mooring stiffness, and drag coefficient (  ).results highlight key insights into the numerical modeling of WF and LF responses in this context: • The LF motions dominate WF motions across various seastates, mooring stiffness, and drag coefficients (  ).
• The inclusion of LF displacements and velocities in the second-order analysis proves crucial, particularly when viscous drag loads are not dominant.Neglecting LF displacements and velocities in slow drift motion analysis can lead to conservative results, potentially underestimating motion responses.• It appears that for shorter waves (smaller   ), the standard deviation of WF surge and pitch motions follows a linear relationship with significant wave height   .However, for longer waves, this relationship deviates upward due to the increasing significance of quadratic damping and even higher-order effects.• The standard deviation for LF surge and pitch motions appears to increase with (  )  .The exponent  is close to 1.1 for LF pitch motions under different   values.For LF surge motions, it  increases from 1.30 for shorter waves (  = 9 s) to 1.55 for longer waves (  = 13.8 s).This suggests that LF pitch motions are influenced by both wave-drift damping and quadratic damping, while LF surge motions are primarily affected by quadratic viscous damping and even higher-order wave loads.These observations are based on certain assumptions and simplifications related to narrow-band LF resonant motions and wave excitation at resonance.• A larger drag coefficient results in smaller LF surge and pitch responses, suggesting that viscous damping dominates over viscous excitation.This observation can be explained by the fact that wave excitation decays exponentially with depth, whereas viscous damping can be contributed by all parts of the structure.The conclusions drawn in this study are specific to the considered deep-draft spar buoy, and they may not necessarily apply to other types of structures.The relative importance of viscous effects can vary for different types of floating structures.Future studies should investigate various types of floating structures before generalizing these conclusions.In the case of slender structures like the spar buoy studied here, it is essential to develop accurate and efficient empirical formulations for LF drag loads in irregular waves to meet practical engineering needs.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Since we still linearize the free-surface conditions using a Taylor expansion with respect to -plane in BFCS, the hydrodynamic loads due to the wetted-surface fluctuation, contributed both by wave elevation and rigid-body motions at the waterline, have to be approximated.The approximation is similar in the IneCS and BFCS formulations, represented by the waterline integrals in Eqs.(60)-(63).
After obtaining the total forces/moments, the motion equations can be formulated in BFCS based on momentum conservation and conservation of angular momentum.A Newmark-Beta method is employed to solve the motion equations in the time domain.To mitigate potential instabilities associated with the ∕ term when solving the motion equations, a stable approach by introducing an infinite-frequency added mass term on both sides of the motion equations, as presented in Kim et al. [50], is applied.For further details, please refer to the works of Shao [36] and Shao and Faltinsen [31].

Fig. 2 .
Fig. 2. Time-domain process based on the consistent second-order BFCS model.
4 and 5 for the two tested signals, respectively.It can be seen from the time series of the filtered results in Figs.4(a) and 5(a) that Z. Zheng et al.

Fig. 5 .
Fig. 5. Comparison of the raw and filtered signals by the improved low-pass filter.The raw data is generated from a JONSWAP wave spectrum with a peak period   = 12 s and significant wave height   = 10 m, superposed to sinusoidal signal with a period of 60 s.

Fig. 6 .
Fig.6.Second-order wave loads acting on a bottom-mounted cylinder with/without forward speed   = 0, −0.1, 0.1.The draft of the cylinder ℎ is equal to its radius .

Fig. 9 .
Fig. 9. Difference-frequency QTFs of the surge excitation force of the spar buoy.Results are presented for both the fixed and freely floating conditions.  is the frequency of the first wave component in the bi-chromatic waves.

Fig. 10 .
Fig. 10.Discretization of irregular waves based on JONSWAP wave spectrum of   =13.8 s,   =14.4 m and  = 2.4 by an equal energy method.One hundred wave frequency components with random phase in a range of [0, 2] are used for discretization.

Fig. 11 .
Fig. 11.The total and LF velocities at MWL of the spar platform in No. 1 seastate (a low seastate) obtained by the improved wavelet low-pass filter.

Fig. 12 .
Fig. 12.The total and LF velocities at MWL of the spar platform at mean waterline in No. 2 seastate (a storm seastate) obtained by improved wavelet low-pass filter.

Fig. 13 .
Fig. 13.Time series of the first-order and second-order solution of the motion responses at COG of the spar platform in No. 1 seastate (a low seastate).Drag coefficient   = 0.2 is used in Morison's equation to empirically account for the viscous loads.

Fig. 14 .
Fig. 14.Time series of the total motion responses (first-and second-order) of the spar platform in No. 1 seastate (a low seastate).Drag coefficient   = 0.2 is used in Morison's equation to empirically account for the viscous darg loads.

Fig. 15 .
Fig. 15.Time series of the total motion responses (first-and second-order) at COG of the spar platform in No. 1 seastate (a low seastate).Drag coefficient   = 1.0 is used in Morison's equation to empirically account for the viscous loads.

Fig. 16 .
Fig. 16.Mean and STD of the surge and pitch motions at COG of the spar platform in No. 1 seastate (a low seastate) using   = 0.2 and 1.0.

Fig. 17 .
Fig. 17.Time series of the total motion responses (first-and second-order) at COG of the spar platform in No. 2 seastate (a storm seastate).Drag coefficient   = 0.2 is used in Morison's equation to empirically account for the viscous loads.

Fig. 18 .
Fig. 18.Time series of the total motion responses (first-and second-order) at COG of the spar platform in No. 2 seastate (a storm seastate).Drag coefficient   = 1.0 is used in Morison's equation to empirically account for the viscous loads.

Fig. 19 .
Fig. 19.Mean and STD of the surge and pitch motions at COG of the spar platform in No. 2 seastate (a storm seastate) using   = 0.2 and 1.0.

Fig. 20 .
Fig. 20.Mean and STDs of the surge motions at MWL of the spar buoy based on different numerical models and   coefficients.(a) Seastate 1 (a low seastate); (b) Seastate 2 (a storm seastate).

Fig. 21 .
Fig. 21.The LF STDs for surge and pitch at COG under different significant wave heights, and the fitted curves.The curves representing the linear and quadratic relationship to   are also shown.(a) Surge at COG; (b) Pitch.

Fig. 22 . 1 ∼ 5 ∼ 1 ∼ 5 ∼
Fig. 22.The WF STDs for different significant wave heights, and the fitted curves.The curves representing the linear and quadratic relationship to   are also shown.

Fig. 23 .
Fig. 23.PSDs of the total motion (first-+ second-order) at COG of the spar platform in No. 2 seastate (a storm seastate).Two different mooring systems are compared.Drag coefficient   = 0.2 is used in Morison's equation to empirically account for the viscous loads.

Fig. 24 .
Fig. 24.Mean and STD of the surge and pitch motions at COG of the spar platform in No. 2 seastate (a storm seastate).Two different mooring systems are compared.Drag coefficient   = 0.2 is used in Morison's equation to empirically account for the viscous loads.

Fig. 25 .
Fig. 25.Mean and STD of the surge and pitch motions at COG of the spar buoy in No. 1 seastate (a low seastate).Two different mooring systems are compared.Drag coefficient   = 0.2 is used in Morison's equation to empirically account for the viscous loads.

Table 1
Parameters of the bi-chromatic waves for a bottom-mounted cylinder. =  2 ∕.

Table 2
Main particulars and loading condition of the spar buoy.

Table 4
Parameters for the four irregular wave conditions.

Table 5
Surge and pitch natural periods of the spar buoy for two different mooring system estimated from single-DOF decay tests.(unit: second)

Table 6
Characterization of incident and scattered wave fields based on three models.
the second stage, we focus on seastates with the same   but different   values, with varying wave steepness  =  = 0.032 to 0.048.Specifically, we compare Seastate 1 with Seastates 3, 5, and 7, and Seastate 2 with Seastates 4, 6, and 8, as defined in Table4.The STDs for these seastates are shown as functions of   in Figs.21(a) and 21(b) for LF surge and pitch motions at COG, respectively.To model the relationship between STDs and   for seastates with the same   , we use the power function represent the STDs of LF surge and pitch motions, respectively.The fitted curves for LF surge STDs based on