Higher-Harmonic Response of a Slender Monopile to Fully Nonlinear Focused Wave Groups

: The “ringing” response of a monopile foundation to focused wave groups was investigated in this study. Such responses are of practical interest in the context of offshore wind turbine foundations. Moderately steep transient focused wave groups were generated in a novel numerical wave tank, based on the high-order spectral method, which was veriﬁed and in good agreement with published laboratory data. The monopile was simpliﬁed into a slender linear elastic cantilever beam to measure the severity of the structural response. The resonant behavior was excited at the triple-wave frequency of the wave loads. The inﬂuence of the damping ratio on the ringing response was considered. Different incident wave models and hydrodynamic models were used to predict the wave loads and induced responses over different wave steepness. The wavelet transform was successfully applied to reveal the local characteristics of the wave loads and ringing response.


Introduction
The rapid expansion of the offshore wind industry has aroused great interest in the behavior of monopile foundations. It is necessary to develop suitable and efficient tools for assessing their reliability and survivability when exposed to steep transient water waves. The peak frequencies of the most common sea states are lower than the first natural frequency of the structure. However, Ridder et al. [1] found that waves with higher harmonics close to the first natural frequency of the structure can still induce a resonant response. This phenomenon, called "ringing," typically occurs when the natural frequency of the monopile is approximately three to five times the wave frequency. The ringing phenomenon was discussed intensively in the 1990s for deep-water tension leg platforms (TLPs). For instance, as reported for model tests by Natvig and Teigen [2], the third-order forcing caused a resonant response for TLP structures. Bredmose et al. [3] investigated the hydro-elastic problem of structural response to fully nonlinear regular wave forcing of a cantilever beam. However, as Atkins et al. [4] mentioned, large transient wave groups are more likely to induce the ringing phenomenon instead of a regular wave train. The response is characterized by a rapid onset of vibration in the natural mode after the passage of a large wave, and it damps out eventually. In this study, an attempt was made to interpret a more realistic and complex case of large focused wave groups, which can be essential in the quest for structural strength design.
In recent decades, considerable efforts have been made to actualize numerical simulations of wave-structure interactions. Analytical theories based on conventional perturbation theory, such as linear and second-order wave diffraction-radiation models, cannot calculate well the higher-harmonic wave loads leading to ringing. These theories, formally assuming a small steepness of the incident wave, are only valid for a moderate wave amplitude. Therefore, to ensure survival in a hostile environment, Rainey [5,6]  (1) On the instantaneous free surface, both the fully nonlinear kinematic and dynamic boundary conditions, derived from the material derivatives and the Bernoulli equation, respectively, must always be fulfilled: where g is the acceleration caused by gravity, and η denotes the free-surface elevation above the still-water level. The kinematic boundary condition on the seabed is: where d denotes the finite water depth. The boundary-value problem can be solved using spectral methods by means of a set of basis functions. The free-surface elevation η(x,t) and the velocity potential ϕ(x,z,t) can be decomposed as follows: where ki are the basis wave numbers, and ψi(x) are the horizontal basis functions (sine and cosine functions or their complex exponential equivalents, satisfying the Laplace equation On the instantaneous free surface, both the fully nonlinear kinematic and dynamic boundary conditions, derived from the material derivatives and the Bernoulli equation, respectively, must always be fulfilled: where g is the acceleration caused by gravity, and η denotes the free-surface elevation above the still-water level. The kinematic boundary condition on the seabed is: where d denotes the finite water depth. The boundary-value problem can be solved using spectral methods by means of a set of basis functions. The free-surface elevation η(x,t) and the velocity potential φ(x,z,t) can be decomposed as follows: where k i are the basis wave numbers, and ψ i (x) are the horizontal basis functions (sine and cosine functions or their complex exponential equivalents, satisfying the Laplace equation and the seabed boundary condition automatically). Here, N is the number of Fourier components, and A η i (t) and A  (2) and (3). Stream function wave theory [18] is adopted for regular waves, and the HOS-NWT for focused wave groups.

Stream Function Theory for Regular Waves
For the progressive regular wave, Rienecker and Fenton [18] provided a fully nonlinear solution. The stream function Ψ is decomposed using the basis functions as follows: where k is the wave number, and c is the phase velocity of each compositing wave. These amplitudes, B j , are solved numerically from the systems of nonlinear equations using Newton's method. More information can be found in the report by Rienecker and Fenton [18].

High-Order Spectral Numerical Wave Tank for Focused Wave Groups
For irregular waves, a two-dimensional wave tank is established with horizontal dimension L x and finite depth d, as shown in Figure 2. An absorbing zone is included near the end wall x = L x through a local modification of the free-surface dynamic boundary condition, which is described using a local modification of pressure at the free surface: where ρ is the fluid density, and the function ν(x) is nonzero, where an absorption zone is required.
where x 0 is the starting point of the absorbing zone, µ = (x − x 0 )/(L x − x 0 ), and α is the strength of the absorbing zone (see previous research [19,20] for more details). Then, the free-surface boundary condition in Equation (3) is rewritten as follows: J. Mar. Sci. Eng. 2021, 9, 286 4 of 25 and the seabed boundary condition automatically). Here, N is the number of Fourier components, and A η i (t) and A ϕ i (t) are the modal amplitudes determined by the free-surface boundary conditions (FSBCs) in Equations (2) and (3). Stream function wave theory [18] is adopted for regular waves, and the HOS-NWT for focused wave groups.

Stream Function Theory for Regular Waves
For the progressive regular wave, Rienecker and Fenton [18] provided a fully nonlinear solution. The stream function Ψ is decomposed using the basis functions as follows: where k is the wave number, and c is the phase velocity of each compositing wave. These amplitudes, Bj, are solved numerically from the systems of nonlinear equations using Newton's method. More information can be found in the report by Rienecker and Fenton [18].

High-Order Spectral Numerical Wave Tank for Focused Wave Groups
For irregular waves, a two-dimensional wave tank is established with horizontal dimension Lx and finite depth d, as shown in Figure 2. An absorbing zone is included near the end wall x = Lx through a local modification of the free-surface dynamic boundary condition, which is described using a local modification of pressure at the free surface: where ρ is the fluid density, and the function ν(x) is nonzero, where an absorption zone is required.
where x0 is the starting point of the absorbing zone, μ = (x − x0)/(Lx − x0), and α is the strength of the absorbing zone (see previous research [19,20] for more details). Then, the free-surface boundary condition in Equation (3) is rewritten as follows: At the inflow boundary, a theoretical particle velocity profile is prescribed from the bottom to the free surface: where u is the horizontal velocity of the underlying water particle kinematics. In addition, a nonflow condition is used on other boundaries: At the inflow boundary, a theoretical particle velocity profile is prescribed from the bottom to the free surface: where u is the horizontal velocity of the underlying water particle kinematics. In addition, a nonflow condition is used on other boundaries: Owing to the stochastic nature of the focused wave groups, according to the NewWave model, the amplitude spectrum of discrete wave components can be determined by a practical spectrum as: where A is the input linear sum of the component wave amplitudes, which occurs at the focus location and time, N is the total number of discrete wave components, A i is the amplitude of each wave component, f is the frequency (inversion of wave period), and S(f ) denotes the corresponding wave spectrum. The energy spectrum proposed by the Joint North Sea Wave Project (JONSWAP) [21] with a peak enhancement factor of γ = 3.30 is used to obtain the underlying wave components. Then, the second-order Stokes theory is used to specify the corresponding velocity: where the linear velocity u (1) and the second-order velocity u (2) can be referred to previous work [16] and thus are not elaborated here. The computation domain of the NWT is nonperiodic owing to the presence of the wavemaker. Such cases lead to the idea of the potential decomposition concept [22], where the original potential φ is separated into two components: where φ p is the periodic potential derived by imposing reflecting lateral boundary conditions, and φ a is the additional potential accounting for the wavemaker boundary condition, describing the evanescent modes that decay exponentially away from the boundary. The boundary conditions of φ p can be written as: The additional potential φ a satisfies the wavemaker condition at the left end of the tank: A nonflow condition is used on other boundaries: To formulate the FSBCs with dependence on the horizontal coordinate only following Zakharov's method [23], the periodic surface potential on the free surface is partially defined as: Introducing the derivatives that stem from this definition: Using the separated potential representation Equation (15), the fully nonlinear freesurface conditions can be expressed as: where W = ∂ z φ p (x, z = η, t), which can be computed following the original HOS method [14,15] as follows: where M is the order of nonlinearity, which is chosen as M = 7, whereas the value of N is selected specifically for each simulation. Because φ p has zero flux through the left, right, and bottom boundaries, the spectral expansions chosen for it are: where A (m) n (t) is the modal amplitudes, the natural eigenmodes of the wave tank are k n = nπ/L x , and N x is the number of wave modes that are considered. Furthermore, each φ (m) p is expanded in a Taylor series around z = 0. The following scheme is arranged according to the order m: Following Bonnefoy et al. [19], φ a is solved in a new extended domain of computation, as shown in Figure 3. This new domain is chosen to have the ability of a fully spectral formulation of the additional problem, which can be solved conveniently by using the FFT. The total height of the extended domain is L z = h add + d, and the horizontal velocity along the input boundary is antisymmetric above the central line z c = (h add − d)/2. The h add is selected as 3d following Bonnefoy et al. [19] with central line height z c = d, such that the central line z c is higher than the maximum wave height. Another issue worth noting is that Bonnefoy et al. [19] specified the motion of the wave generator from the bottom of the tank to the static water level. Instead, a more reasonable one is proposed for a steep nonlinear wave to feed the input wave properties from the bottom to the free surface. The third-order polynomial function from z = η(t) to z = h add − d − η(t) is chosen as a matching curve to close smoothly the domain D add , which is defined as: The coefficients considered at an instantaneous time step are given as: ( ) where Bn(t) is the modal amplitudes, and kn = (2n − 1)π/(hadd + d). Subsequently, the additional derivatives ∂xϕa are expanded analytically as: x a n n n 1 n x sinh ( , , ) Substituting the special derivatives into Equation (17), the time-dependent coefficients Bn(t) are determined through the explicit known u by means of the FFT method. Once the time-dependent coefficients Bn(t) are determined, ϕa is solved eventually.
Instead of using backward finite difference to estimate ∂tϕa, the time derivative is directly obtained through the explicit known ut by means of the same scheme for ϕa. According to Equation (17), the time derivatives also satisfy the wavemaker condition: Here, ∂tϕa and ∂xtϕa are expanded as: xt a n n n 1 n x sinh ( , , ) Substituting Equation (33) into Equation (31), the time-dependent coefficients Cn(t) are determined through the explicit known ut by means of the FFT method. Subsequently, ∂tϕadd is solved eventually.
To reduce any impulse-like behavior at the inlet boundary, a cosine ramping function is applied at the start of the computation. The ramping function is given by: where Tm is the ramp time, chosen as 2T, with T the characteristic wave period. To obtain a time-dependent solution, the initial calm water surface condition is applied in the present model as follows: The solution algorithm for the HOS-NWT is schematized in Figure 4. It is divided into two main parts. One addresses the solution of the wave generation problem in Dadd, In this extended domain, the additional potential φ a is expanded as: where B n (t) is the modal amplitudes, and k n = (2n − 1)π/(h add + d). Subsequently, the additional derivatives ∂ x φ a are expanded analytically as: Substituting the special derivatives into Equation (17), the time-dependent coefficients B n (t) are determined through the explicit known u by means of the FFT method. Once the time-dependent coefficients B n (t) are determined, φ a is solved eventually.
Instead of using backward finite difference to estimate ∂ t φ a , the time derivative is directly obtained through the explicit known u t by means of the same scheme for φ a . According to Equation (17), the time derivatives also satisfy the wavemaker condition: Here, ∂ t φ a and ∂ xt φ a are expanded as: Substituting Equation (33) into Equation (31), the time-dependent coefficients C n (t) are determined through the explicit known u t by means of the FFT method. Subsequently, ∂ t φ add is solved eventually.
To reduce any impulse-like behavior at the inlet boundary, a cosine ramping function is applied at the start of the computation. The ramping function is given by: where T m is the ramp time, chosen as 2T, with T the characteristic wave period. To obtain a time-dependent solution, the initial calm water surface condition is applied in the present model as follows: The solution algorithm for the HOS-NWT is schematized in Figure 4. It is divided into two main parts. One addresses the solution of the wave generation problem in D add , and the other addresses the fully nonlinear FSBCs to obtain φ s p (or equivalently φ p ) and η. The fully spectral solution scheme (including the wave generation) makes the FFT-based solution scheme, which exhibits great efficiency, possible. The spatial derivatives of the velocity potentials φ p and φ a are straightforwardly reconstructed from semi-analytical basis expressions. Equations (22) and (23) then provide the time derivatives of the unknowns φ s p and η, which are further used in a time-marching fourth-order Runge-Kutta Cash-Karp scheme with an adaptive step size [24]. J. Mar. Sci. Eng. 2021, 9, 286 8 of 25 and the other addresses the fully nonlinear FSBCs to obtain ϕ s p (or equivalently ϕp) and η. The fully spectral solution scheme (including the wave generation) makes the FFT-based solution scheme, which exhibits great efficiency, possible. The spatial derivatives of the velocity potentials ϕp and ϕa are straightforwardly reconstructed from semi-analytical basis expressions. Equations (22) and (23) then provide the time derivatives of the unknowns ϕ s p and η, which are further used in a time-marching fourth-order Runge-Kutta Cash-Karp scheme with an adaptive step size [24].

Hydro-Elastic Model for Monopile
In this section, the complicated hydro-elastic problem for monopile is solved. The monopile has a diameter D, length l, structural mass per unit length , and lateral bending stiffness EI, where E is the modulus of elasticity and I is the identical cross-sectional moment of inertia along full length. The bottom of the monopile is embedded in the seabed. The natural frequency of the first mode of vibration is denoted as follows: Under the assumption of a small ratio of monopile diameter to the wave length, the wave loads can be calculated by neglecting the disturbed waves, using Rainey's model [5,6]. Instead of a conventional perturbation approach, Rainey's model was derived from energy considerations, allowing the numerical computation of wave loading induced by regular and irregular waves. The loads per unit length of a monopile can be written as: where Cm is the added mass coefficient, CD is the drag coefficient, ρ is the density of water, A = πD 2 /4 is the structural cross-sectional area, u and w are the horizontal and vertical particle velocity components, respectively, and ( , ) are the structural velocity and acceleration, respectively. Furthermore, is the Lagrangian water particle acceleration ( = ut + uux + wuz for the present 2D flow), and the last inertia term in Equation (37) is the axial divergence correction. In addition, a point force [6] is added at the intersection with the instantaneous water level to represent the change in kinetic energy associated with the change in the wetted area:

Hydro-Elastic Model for Monopile
In this section, the complicated hydro-elastic problem for monopile is solved. The monopile has a diameter D, length l, structural mass per unit length m, and lateral bending stiffness EI, where E is the modulus of elasticity and I is the identical cross-sectional moment of inertia along full length. The bottom of the monopile is embedded in the seabed. The natural frequency of the first mode of vibration is denoted as follows: Under the assumption of a small ratio of monopile diameter to the wave length, the wave loads can be calculated by neglecting the disturbed waves, using Rainey's model [5,6]. Instead of a conventional perturbation approach, Rainey's model was derived from energy considerations, allowing the numerical computation of wave loading induced by regular and irregular waves. The loads per unit length of a monopile can be written as: where C m is the added mass coefficient, C D is the drag coefficient, ρ is the density of water, A = πD 2 /4 is the structural cross-sectional area, u and w are the horizontal and vertical particle velocity components, respectively, and ( . X, .. X) are the structural velocity and acceleration, respectively. Furthermore, . u is the Lagrangian water particle acceleration ( . u = u t + uu x + wu z for the present 2D flow), and the last inertia term in Equation (37) is the axial divergence correction. In addition, a point force [6] is added at the intersection with the instantaneous water level to represent the change in kinetic energy associated with the change in the wetted area: Then, the calculation of the overturning moment at the base of the monopile integrated f (t,z) up to the instantaneous water surface, and the expression can be written as: In all calculations, fixed generic values of C D = 1.0 and C m = 1.0 were used. The structural motion is modeled by the linear beam equation: with standard cantilever beam boundary conditions: The coefficient b in Equation (40) is the damping, assessed through the proportional Rayleigh damping method (see previous work [25] for more details). Equation (40) is timeintegrated with the Newmark integral scheme [26]. The top-point motion X T is chosen as one measure of the severity of the structural response.

Wavelet Transform-Based Analyses
The wave loads and the related ringing response are transient phenomena rather than stationary processes. To investigate their local characteristics, the wavelet transform, which can generate localized time-frequency information from time series, is adopted. The continuous wavelet transform is defined by the inner product of a time function χ(t) and a family of continuously translated and dilated wavelets ψ a,τ (t): where the asterisk * denotes the complex conjugate. A family of functions ψ a,τ (t), called "wavelets", can be constructed by translating in time τ and dilation with scale a of a specified mother wavelet ψ(t). The scale a is interpreted as the reciprocal of frequency, f ≈ 1/a. The ψ a,τ (t) expression is defined as: One of the most extensively used mother wavelets in continuous wavelet analysis is the Morlet wavelet [27]. The complex Morlet wavelet is a plane wave modulated by a Gaussian envelope: where ω 0 is the center frequency of the wavelet, usually chosen to be 6.0, to meet the allowable conditions [28].

Validation Tests of Focused Wave Groups
The capability of the proposed formulation for focused wave group generation is verified here and compared with the experimental results from Ning et al. [16]. In the physical experiment, the wave flume is 69 m long and 3 m wide, and the maximum working depth is 1.5 m. The present study only focused on Cases 1 and 2 from the work by Ning et al. The water depth used in the experiments was d = 0.5 m. The details of the chosen focused wave groups are given in Table 1. The peak frequency f p = 0.83 Hz (corresponding to wavelength L p = 2.0 m and wavenumber k p = 3.14 m −1 ). The dimensionless wave steepness is defined as ε = k p A I = 0.10 for Case 1 and ε = k p A I = 0.20 for Case 2. In the numerical simulation, the length of the tank is selected as 10 times the characteristic wavelength (10L p ), of which 4L p is used as a damping layer at the downstream end of the tank to dissipate the outgoing waves, and the depth is 0.5 m. Each focused wave group consisted of 25 wave components, and the focus position and focus time are defined as x f = 1.5L p and t f = 8T p , respectively. At the start of the computation, a ramp function is applied over a time 2T p to prevent any impulse-like behavior at the flux input boundary and to reduce the corresponding unnecessary transient waves. First, for the convergence analysis, Figure 5 shows the time series of free-surface elevations of Case 1 from the work by Ning et al. [16], calculated using two different mesh sizes and two different time steps. For Mesh A, results are performed with N x = 512 (i.e., 52 modes per wavelength), and N z = 128 vertical modes on the wavemaker. The number of modes of Mesh B (N x = 1024, and N z = 256 vertical modes on the wavemaker) is twice that of Mesh A. The time steps were selected as ∆t = T p /60 and ∆t = T p /120. Figure 5 shows that the results are almost identical, indicating that the numerical result is convergent using Mesh A, with ∆t = T p /60. The same grid and time step were adopted in the following simulation for Case 2. chosen focused wave groups are given in Table 1. The peak frequency fp = 0.83 Hz (corresponding to wavelength Lp = 2.0 m and wavenumber kp = 3.14 m −1 ). The dimensionless wave steepness is defined as ε = kpAI = 0.10 for Case 1 and ε = kpAI = 0.20 for Case 2. In the numerical simulation, the length of the tank is selected as 10 times the characteristic wavelength (10Lp), of which 4Lp is used as a damping layer at the downstream end of the tank to dissipate the outgoing waves, and the depth is 0.5 m. Each focused wave group consisted of 25 wave components, and the focus position and focus time are defined as xf = 1.5Lp and tf = 8Tp, respectively. At the start of the computation, a ramp function is applied over a time 2Tp to prevent any impulse-like behavior at the flux input boundary and to reduce the corresponding unnecessary transient waves. First, for the convergence analysis, Figure 5 shows the time series of free-surface elevations of Case 1 from the work by Ning et al. [16], calculated using two different mesh sizes and two different time steps. For Mesh A, results are performed with Nx = 512 (i.e., ≃52 modes per wavelength), and Nz = 128 vertical modes on the wavemaker. The number of modes of Mesh B (Nx = 1024, and Nz = 256 vertical modes on the wavemaker) is twice that of Mesh A. The time steps were selected as Δt = Tp/60 and Δt = Tp/120. Figure 5 shows that the results are almost identical, indicating that the numerical result is convergent using Mesh A, with Δt = Tp/60. The same grid and time step were adopted in the following simulation for Case 2. The simulation results were compared with Cases 1 and 2 laboratory data from Ning et al. [16] and the input theoretical calculations in Figure 6. The present numerical results of Case 1 agree well with the second-order analytical solutions and the experimental measurements. Because the amplitude of the individual wave components was small, the effects of nonlinear wave-wave interactions were also very small. For Case 2, the wave crest at the focal point was larger and narrower, and the wave trough became wider and flatter, i.e., stronger nonlinearity. The present results agree quite well with the experimental measurements in terms of the wave crest, trough, and phase. The second-order analytical solution, whose shape is always perfectly symmetric, underestimates the main The simulation results were compared with Cases 1 and 2 laboratory data from Ning et al. [16] and the input theoretical calculations in Figure 6. The present numerical results of Case 1 agree well with the second-order analytical solutions and the experimental measurements. Because the amplitude of the individual wave components was small, the effects of nonlinear wave-wave interactions were also very small. For Case 2, the wave crest at the focal point was larger and narrower, and the wave trough became wider and flatter, i.e., stronger nonlinearity. The present results agree quite well with the experimental measurements in terms of the wave crest, trough, and phase. The second-order analytical solution, whose shape is always perfectly symmetric, underestimates the main wave crest and troughs. In addition, there is a small phase difference when using it. Thus, the results of HOS-NWT are superior to the theoretical second-order Stokes wave, which is convincing evidence of the generation of focused wave groups in the HOS-NWT. Different types of wave nonlinearity can be easily discerned in corresponding amplitude spectra in Figure 7a,b. Starting from low frequency, the main components of the focused wave groups contain a second-order difference term, a linear term, and a secondorder sum term, where the present numerical results are consistent with the input secondorder analytical solution. The smaller third-order sum terms and even the fourth-order terms can also be discerned with a larger wave steepness.

Validation Tests of Hydro-Elastic Model
The ability of the hydro-elastic model to predict the monopile response to nonlinear wave forcing can be validated by comparing the results with those of Bredmose et al. [3].  Different types of wave nonlinearity can be easily discerned in corresponding amplitude spectra in Figure 7a,b. Starting from low frequency, the main components of the focused wave groups contain a second-order difference term, a linear term, and a second-order sum term, where the present numerical results are consistent with the input second-order analytical solution. The smaller third-order sum terms and even the fourth-order terms can also be discerned with a larger wave steepness. Different types of wave nonlinearity can be easily discerned in corresponding amplitude spectra in Figure 7a,b. Starting from low frequency, the main components of the focused wave groups contain a second-order difference term, a linear term, and a secondorder sum term, where the present numerical results are consistent with the input secondorder analytical solution. The smaller third-order sum terms and even the fourth-order terms can also be discerned with a larger wave steepness.

Validation Tests of Hydro-Elastic Model
The ability of the hydro-elastic model to predict the monopile response to nonlinear wave forcing can be validated by comparing the results with those of Bredmose et al. [3].

Validation Tests of Hydro-Elastic Model
The ability of the hydro-elastic model to predict the monopile response to nonlinear wave forcing can be validated by comparing the results with those of Bredmose et al. [3]. The main simulation parameters were chosen based on a typical offshore wind turbine foundation: a depth of d = 20 m, a diameter of D = 5 m, a length of l = 40.0 m, g = 9.81 m/s 2 , density of water ρ = 1025 kg/m 3 , structural mass per unit length m = 20,000.0 kg/m 3 , and the flexural rigidity of section EI = 10219084252.53 N/m 2 . As a result, the natural frequency of the monopile was set as f 1 = 0.25 Hz. The above parameters were used for all the following simulations. The regular wave parameters were set as wave height H = 5.00 m, wave length L = 155.91 m, wave period T = 12 s, and wave frequency f = 1/12 Hz, along with f 1 /f = 3.0. The dependence of the results on the mesh size and time step was carefully checked by performing a convergence study. Equation (40) was discretized spatially using the finite-element method, with 200 Euler beam elements for Mesh A and 300 Euler beam elements for Mesh B. The time step interval was set as ∆t = T/120 and ∆t = T/160. The time series of the top-point motion X T , normalized with wave height, are presented in Figure 8. One can see that the numerical result is convergent using Mesh A, with ∆t = T/120, which is adopted in the work discussed in Section 4. Overall, a satisfactory agreement with Bredmose et al. was reached, demonstrating the feasibility of the present hydro-elastic model for the monopile.

Numerical Results and Discussion
The ringing response induced by the third-order harmonic wave loads of focused wave groups generated in the HOS-NWT was investigated, and the results are discussed in this section. The simulation parameters of the monopile are the same as those in Section 3.2. The peak spectral frequency of the incident focused wave group was set to one-third of the natural frequency of the monopile, i.e., fp = f1/3 = 1/12 Hz. The focused wave group consisted of 25 continuous sinusoidal components whose frequencies were equally spaced within the range 0.5fp < f < 1.5fp. The wave number kp = 0.0412 m −1 is the linear dispersion wave number corresponding to fp, along with the wavelength Lp = 152.33 m. Following Ning et al. [16], the focus position, where the monopile is placed, and focus time are defined as xf = 1.50Lp and tf = 8Tp, respectively. The length of the wave tank was L = 10Lp = 1523.30 m. The length of the damping zones at the end of the computation domain was 4Lp. The simulations covered seven different input amplitudes A. Detailed information is listed in Table 2. The time sequence of the surface elevation at the focus position, with the highest input amplitude (i.e., Case 7), and the corresponding wavelet spectrum are shown in Figures 9 and 10. In Figure 9, the results are compared with the linear, second-order solutions. The numerically calculated crest of focused wave groups is higher and narrower, while troughs become wider and shallower than those of theoretical solutions, which is attributed to the nonlinear interaction between wave components. In Figure 10, the color bar indicates the modulus of the wavelet coefficients, and the image on the right is a partial enlarged view of the left image for higher-order components in the normalized fre-

Numerical Results and Discussion
The ringing response induced by the third-order harmonic wave loads of focused wave groups generated in the HOS-NWT was investigated, and the results are discussed in this section. The simulation parameters of the monopile are the same as those in Section 3.2. The peak spectral frequency of the incident focused wave group was set to one-third of the natural frequency of the monopile, i.e., f p = f 1 /3 = 1/12 Hz. The focused wave group consisted of 25 continuous sinusoidal components whose frequencies were equally spaced within the range 0.5f p < f < 1.5f p . The wave number k p = 0.0412 m −1 is the linear dispersion wave number corresponding to f p , along with the wavelength L p = 152.33 m. Following Ning et al. [16], the focus position, where the monopile is placed, and focus time are defined as x f = 1.50L p and t f = 8T p , respectively. The length of the wave tank was L = 10L p = 1523.30 m. The length of the damping zones at the end of the computation domain was 4L p . The simulations covered seven different input amplitudes A. Detailed information is listed in Table 2. The time sequence of the surface elevation at the focus position, with the highest input amplitude (i.e., Case 7), and the corresponding wavelet spectrum are shown in Figures 9 and 10. In Figure 9, the results are compared with the linear, second-order solutions. The numerically calculated crest of focused wave groups is higher and narrower, while troughs become wider and shallower than those of theoretical solutions, which is attributed to the nonlinear interaction between wave components. In Figure 10, the color bar indicates the modulus of the wavelet coefficients, and the image on the right is a partial enlarged view of the left image for higher-order components in the normalized frequency range of 1.5~4.5. The wavelet analysis results show that only the first and second-harmonic components of focused wave groups are evident for theoretical solutions. However, the third and fourth harmonics appear using the focused wave generated in the HOS-NWT. J. Mar. Sci. Eng. 2021, 9,286 13 of 25 harmonic components of focused wave groups are evident for theoretical solutions. However, the third and fourth harmonics appear using the focused wave generated in the HOS-NWT.  The dimensionless maximum crest elevations versus the wave steepness kpA are shown in Figure 11. The dimensionless maximum crest elevations are further enhanced by increasing the input amplitude, indicating that the nonlinearity of the freak wave increases with the wave steepness. harmonic components of focused wave groups are evident for theoretical solutions. However, the third and fourth harmonics appear using the focused wave generated in the HOS-NWT.  The dimensionless maximum crest elevations versus the wave steepness kpA are shown in Figure 11. The dimensionless maximum crest elevations are further enhanced by increasing the input amplitude, indicating that the nonlinearity of the freak wave increases with the wave steepness. The dimensionless maximum crest elevations versus the wave steepness k p A are shown in Figure 11. The dimensionless maximum crest elevations are further enhanced by increasing the input amplitude, indicating that the nonlinearity of the freak wave increases with the wave steepness.

Damping Ratio Study
Although the damping of the monopile is typically low [29], it has a significant effect on the resonant response. Different damping ratios ξ = 0.00, 0.01, 0.02, and 0.05 for Case 7 were investigated. The time histories of the sectional moment at the base, calculated using Rainey's model with focused wave groups generated in the HOS-NWT, and the corresponding wavelet spectrum are shown in Figures 12 and 13. In Figure 13, the image on the right is a partial enlarged view of the left image for higher-order components in the normalized frequency range of 2.5~4.5. Figure 12 shows that the magnitude of the damping ratio has a slight influence on the peak of the wave loads. From the wavelet spectrum in Figure 13, one can see that the wave load is dominated by the first-order component, and the higher-order components of wave loads are much smaller, as expected. Their wavelet amplitudes in the freak wave region are more concentrated, larger in amplitude, and wider in frequency band, which indicates that considerable high-order harmonic overturning moments tend to appear in waves with transient characteristics. In Figure 13, another larger peak of the third-order moment after the first peak can be found with ξ = 0.00.

Damping Ratio Study
Although the damping of the monopile is typically low [29], it has a significant effect on the resonant response. Different damping ratios ξ = 0.00, 0.01, 0.02, and 0.05 for Case 7 were investigated. The time histories of the sectional moment at the base, calculated using Rainey's model with focused wave groups generated in the HOS-NWT, and the corresponding wavelet spectrum are shown in Figures 12 and 13. In Figure 13, the image on the right is a partial enlarged view of the left image for higher-order components in the normalized frequency range of 2.5~4.5. Figure 12 shows that the magnitude of the damping ratio has a slight influence on the peak of the wave loads. From the wavelet spectrum in Figure 13, one can see that the wave load is dominated by the first-order component, and the higher-order components of wave loads are much smaller, as expected. Their wavelet amplitudes in the freak wave region are more concentrated, larger in amplitude, and wider in frequency band, which indicates that considerable high-order harmonic overturning moments tend to appear in waves with transient characteristics. In Figure 13, another larger peak of the third-order moment after the first peak can be found with ξ = 0.00.

Damping Ratio Study
Although the damping of the monopile is typically low [29], it has a significant effect on the resonant response. Different damping ratios ξ = 0.00, 0.01, 0.02, and 0.05 for Case 7 were investigated. The time histories of the sectional moment at the base, calculated using Rainey's model with focused wave groups generated in the HOS-NWT, and the corresponding wavelet spectrum are shown in Figures 12 and 13. In Figure 13, the image on the right is a partial enlarged view of the left image for higher-order components in the normalized frequency range of 2.5~4.5. Figure 12 shows that the magnitude of the damping ratio has a slight influence on the peak of the wave loads. From the wavelet spectrum in Figure 13, one can see that the wave load is dominated by the first-order component, and the higher-order components of wave loads are much smaller, as expected. Their wavelet amplitudes in the freak wave region are more concentrated, larger in amplitude, and wider in frequency band, which indicates that considerable high-order harmonic overturning moments tend to appear in waves with transient characteristics. In Figure 13, another larger peak of the third-order moment after the first peak can be found with ξ = 0.00.    Figure 14 shows that the ringing phenomenon is obvious. The oscillation of the beam bursts a little after the focus time and lasts for several peak wave periods, which indicates that the ringing response occurs immediately after the incident wave crest passes the location of the monopile. As the damping ratio increases, the burst motions of the monopile are almost constant; however, the duration of the subsequent resonant response is much shorter compared with that of the pattern obtained without considering the damping. Figure 15 shows that, although the first-and second-harmonic components of wave loads are predominant, the response induced by them is almost negligible, while the third one dominates the resonant response because of the dynamic amplification. After the focus time, the third-harmonic components appear to be quite significant and last for several peak periods. When a larger damping ratio is considered, the amplitudes of the third-order response are significantly reduced.  The top-point motion and corresponding wavelet spectrum are shown in Figures 14 and 15. Figure 14 shows that the ringing phenomenon is obvious. The oscillation of the beam bursts a little after the focus time and lasts for several peak wave periods, which indicates that the ringing response occurs immediately after the incident wave crest passes the location of the monopile. As the damping ratio increases, the burst motions of the monopile are almost constant; however, the duration of the subsequent resonant response is much shorter compared with that of the pattern obtained without considering the damping. Figure 15 shows that, although the first-and second-harmonic components of wave loads are predominant, the response induced by them is almost negligible, while the third one dominates the resonant response because of the dynamic amplification. After the focus time, the third-harmonic components appear to be quite significant and last for several peak periods. When a larger damping ratio is considered, the amplitudes of the third-order response are significantly reduced.  Figure 14 shows that the ringing phenomenon is obvious. The oscillation of the beam bursts a little after the focus time and lasts for several peak wave periods, which indicates that the ringing response occurs immediately after the incident wave crest passes the location of the monopile. As the damping ratio increases, the burst motions of the monopile are almost constant; however, the duration of the subsequent resonant response is much shorter compared with that of the pattern obtained without considering the damping. Figure 15 shows that, although the first-and second-harmonic components of wave loads are predominant, the response induced by them is almost negligible, while the third one dominates the resonant response because of the dynamic amplification. After the focus time, the third-harmonic components appear to be quite significant and last for several peak periods. When a larger damping ratio is considered, the amplitudes of the third-order response are significantly reduced.   In addition, Figure 16 shows the top-point horizontal acceleration, which is further amplified owing to the ω 2 effect. The maximum horizontal acceleration can reach 0.53 g with ξ = 0.00. Therefore, the ringing response induced by the high-order wave loads will be noticed in the design and operation of monopiles. In summary, the wavelet transform method is promising for providing reliable and rich detailed information for the characteristics of the higher-harmonic wave loads and ringing response. In the present study, following Bredmose et al. [3], a smaller value of damping ratio ξ = 0.01 is used in the following simulations.

Different Wave Models for Wave Loads and Ringing Response
For the ultimate limit state, it is interesting to see how large the contribution from the nonlinearity of the incident wave for wave loads. The overturning moments for Case 7 are presented as an example, using Rainey's model. In Figure 17, the differences between the In addition, Figure 16 shows the top-point horizontal acceleration, which is further amplified owing to the ω 2 effect. The maximum horizontal acceleration can reach 0.53 g with ξ = 0.00. Therefore, the ringing response induced by the high-order wave loads will be noticed in the design and operation of monopiles. In addition, Figure 16 shows the top-point horizontal acceleration, which is further amplified owing to the ω 2 effect. The maximum horizontal acceleration can reach 0.53 g with ξ = 0.00. Therefore, the ringing response induced by the high-order wave loads will be noticed in the design and operation of monopiles. In summary, the wavelet transform method is promising for providing reliable and rich detailed information for the characteristics of the higher-harmonic wave loads and ringing response. In the present study, following Bredmose et al. [3], a smaller value of damping ratio ξ = 0.01 is used in the following simulations.

Different Wave Models for Wave Loads and Ringing Response
For the ultimate limit state, it is interesting to see how large the contribution from the nonlinearity of the incident wave for wave loads. The overturning moments for Case 7 are presented as an example, using Rainey's model. In Figure 17, the differences between the In summary, the wavelet transform method is promising for providing reliable and rich detailed information for the characteristics of the higher-harmonic wave loads and ringing response. In the present study, following Bredmose et al. [3], a smaller value of damping ratio ξ = 0.01 is used in the following simulations.

Different Wave Models for Wave Loads and Ringing Response
For the ultimate limit state, it is interesting to see how large the contribution from the nonlinearity of the incident wave for wave loads. The overturning moments for Case 7 are presented as an example, using Rainey's model. In Figure 17, the differences between the different focused wave models are therefore easy to identify. There is an obvious secondary load cycle [30], which only occurs using focused waves generated in the HOS-NWT, suggesting the severest kinematics among the incident wave models and higher harmonics appearance. Its main peak moment is the largest among all three incident wave models. The wavelet spectra overturning moments were examined, as shown in Figure 18. The image on the right is a partial enlarged view of the left image for higherorder components in the normalized frequency range of 2.5~4.5. For the linear incident wave model, only the first-harmonic components are evident. For the second-order Stokes theory, except for the linear component, the second harmonic is also obvious. The third and fourth harmonics appear using the focused wave generated in the HOS-NWT.
J. Mar. Sci. Eng. 2021, 9,286 17 of 25 different focused wave models are therefore easy to identify. There is an obvious secondary load cycle [30], which only occurs using focused waves generated in the HOS-NWT, suggesting the severest kinematics among the incident wave models and higher harmonics appearance. Its main peak moment is the largest among all three incident wave models. The wavelet spectra overturning moments were examined, as shown in Figure 18. The image on the right is a partial enlarged view of the left image for higher-order components in the normalized frequency range of 2.5~4.5. For the linear incident wave model, only the first-harmonic components are evident. For the second-order Stokes theory, except for the linear component, the second harmonic is also obvious. The third and fourth harmonics appear using the focused wave generated in the HOS-NWT.   different focused wave models are therefore easy to identify. There is an obvious secondary load cycle [30], which only occurs using focused waves generated in the HOS-NWT, suggesting the severest kinematics among the incident wave models and higher harmonics appearance. Its main peak moment is the largest among all three incident wave models. The wavelet spectra overturning moments were examined, as shown in Figure 18. The image on the right is a partial enlarged view of the left image for higher-order components in the normalized frequency range of 2.5~4.5. For the linear incident wave model, only the first-harmonic components are evident. For the second-order Stokes theory, except for the linear component, the second harmonic is also obvious. The third and fourth harmonics appear using the focused wave generated in the HOS-NWT.    Figure 19 presents the total moment and its components based on the focused wave groups generated in the HOS-NWT. It shows that the peak value of the total moment is dominated by the inertial part, and the peak of the inertial moment and the axial divergence moment appear in the rising phase of the focused wave groups, while the drag moment peak and the surface intersection moment trough are synchronized with the occurrence of the focused wave groups. Because of the reverse phase, although the peak value of the drag moment is considerable in this case, it contributes disproportionately to the peak of the total moment, which is still synchronized with the inertial part.
. 2021, 9,286 18 of 25 Figure 19 presents the total moment and its components based on the focused wave groups generated in the HOS-NWT. It shows that the peak value of the total moment is dominated by the inertial part, and the peak of the inertial moment and the axial divergence moment appear in the rising phase of the focused wave groups, while the drag moment peak and the surface intersection moment trough are synchronized with the occurrence of the focused wave groups. Because of the reverse phase, although the peak value of the drag moment is considerable in this case, it contributes disproportionately to the peak of the total moment, which is still synchronized with the inertial part. Comparisons of four components of the Rainey's model using different incident wave models are shown in Figure 20a-d. For inertia-dominated structures, as here, the wave load depends more on the particle acceleration and, therefore, becomes larger the steeper the waves are. The inertial part is the cause of the secondary loading cycle. For the remaining three moments, the wave generated in the HOS-NWT gives the largest amplitude, owing to the severest kinematics, as explained above. Comparisons of four components of the Rainey's model using different incident wave models are shown in Figure 20a-d. For inertia-dominated structures, as here, the wave load depends more on the particle acceleration and, therefore, becomes larger the steeper the waves are. The inertial part is the cause of the secondary loading cycle. For the remaining three moments, the wave generated in the HOS-NWT gives the largest amplitude, owing to the severest kinematics, as explained above. Figure 21 shows the time series of the overturning moment, using the focused wave groups generated in the HOS-NWT with different wave steepness. As the wave steepness increases, the peak moment (i.e., the maximum moment) tends to appear simultaneously with the occurrence of the focused wave groups, i.e., t = 8T p . As the wave steepness decreases, the secondary loading cycle disappears. In Figure 22, the dimensionless peak moments are presented as a function of wave steepness. The dimensionless peak moment increases with an increase in the wave steepness. Based on the linear theory, the results should be independent of the wave steepness. Thus, the relationship between the peak nondimensional moment and wave steepness is strongly nonlinear rather than linear. In general, the ultimate loads calculated by the two different nonlinear wave theories agree with each other well. However, the discrepancies between the nonlinear and linear wave models increased monotonically with increasing wave steepness.  Figure 21 shows the time series of the overturning moment, using the focused wave groups generated in the HOS-NWT with different wave steepness. As the wave steepness increases, the peak moment (i.e., the maximum moment) tends to appear simultaneously with the occurrence of the focused wave groups, i.e., t = 8Tp. As the wave steepness decreases, the secondary loading cycle disappears. In Figure 22, the dimensionless peak moments are presented as a function of wave steepness. The dimensionless peak moment increases with an increase in the wave steepness. Based on the linear theory, the results should be independent of the wave steepness. Thus, the relationship between the peak nondimensional moment and wave steepness is strongly nonlinear rather than linear. In general, the ultimate loads calculated by the two different nonlinear wave theories agree with each other well. However, the discrepancies between the nonlinear and linear wave models increased monotonically with increasing wave steepness.     Figure 23 shows that the peak values o the dimensionless motion, using the focused wave groups generated in the HOS-NWT are larger than those of the remaining two wave models. The linear wave model is in rather poor agreement with other nonlinear wave models. In Figure 24, the wavelet anal ysis shows that differences are obvious among different incident wave models. Second order wave theory performs well for the prediction of peak wave loading, but it has poorer prediction of the ringing response because the higher harmonics of the moment are no captured, as shown in Figure 18. Therefore, the predictions of wave loads and induced responses are closely correlated with the incident wave model, and the fully nonlinear focused wave generated in the HOS-NWT are preferred as the incident focused wave model.  Figure 23 shows that the peak values of the dimensionless motion, using the focused wave groups generated in the HOS-NWT, are larger than those of the remaining two wave models. The linear wave model is in rather poor agreement with other nonlinear wave models. In Figure 24, the wavelet analysis shows that differences are obvious among different incident wave models. Second-order wave theory performs well for the prediction of peak wave loading, but it has poorer prediction of the ringing response because the higher harmonics of the moment are not captured, as shown in Figure 18. Therefore, the predictions of wave loads and induced responses are closely correlated with the incident wave model, and the fully nonlinear focused wave generated in the HOS-NWT are preferred as the incident focused wave model.

Different Hydrodynamic Models for Wave Loads and Ringing Response
The Morison equation and Rainey's model, which contains two additional potentia flow loads, i.e., the axial divergence part and surface intersection part, to improve accu racy, were applied in assessing the wave loads for comparison purposes. The nonlinear freak waves generated in the HOS-NWT were used as incident waves under Case 7 as an example. As shown in Figure 25, the results from the Morison equation and Rainey's model agree well, except for the main peak and trough value. The Morison equation may underestimate the peak wave loads and trough wave loads. When the axial divergence

Different Hydrodynamic Models for Wave Loads and Ringing Response
The Morison equation and Rainey's model, which contains two additional potential flow loads, i.e., the axial divergence part and surface intersection part, to improve accuracy, were applied in assessing the wave loads for comparison purposes. The nonlinear freak waves generated in the HOS-NWT were used as incident waves under Case 7 as an example. As shown in Figure 25, the results from the Morison equation and Rainey's model agree well, except for the main peak and trough value. The Morison equation may underestimate the peak wave loads and trough wave loads. When the axial divergence part is considered, the peak overturning moment agrees better with Rainey's model. In other words, the axial divergence part plays a more important role in prediction for peak wave loads, compared with the surface intersection part. Particular attention should be given to the secondary loading cycle, which appears even if only the Morison equation is considered. In Figure 26, the image on the right is a partial enlarged view of the left image for higher-order components in the normalized frequency range of 2.5~4.5. As shown in Figure 26, the first-and second-order components of the wavelet transform resulting from the Morison equation agree well with Rainey's model, but the results around the higher frequency using Rainey's model appear to be slightly larger because of the inclusion of slender-body terms.
given to the secondary loading cycle, which appears even if only the Morison equation is considered. In Figure 26, the image on the right is a partial enlarged view of the left image for higher-order components in the normalized frequency range of 2.5~4.5. As shown in Figure 26, the first-and second-order components of the wavelet transform resulting from the Morison equation agree well with Rainey's model, but the results around the higher frequency using Rainey's model appear to be slightly larger because of the inclusion of slender-body terms.  given to the secondary loading cycle, which appears even if only the Morison equation is considered. In Figure 26, the image on the right is a partial enlarged view of the left image for higher-order components in the normalized frequency range of 2.5~4.5. As shown in Figure 26, the first-and second-order components of the wavelet transform resulting from the Morison equation agree well with Rainey's model, but the results around the higher frequency using Rainey's model appear to be slightly larger because of the inclusion of slender-body terms.   Figure 27 shows the relationship between the peak nondimensional moment and the incident wave steepness using different hydrodynamic models. The predictions of Rainey's model exceed those of the Morison equation by a margin that increases rapidly with wave steepness. The agreement is better when the axial divergence part is included, as explained above. Figure 27 shows the relationship between the peak nondimensional moment and the incident wave steepness using different hydrodynamic models. The predictions of Rainey's model exceed those of the Morison equation by a margin that increases rapidly with wave steepness. The agreement is better when the axial divergence part is included, as explained above. Figure 27. Peak normalized moments as a function of wave steepness, computed using different hydrodynamic models with focused wave groups generated in the HOS-NWT. Figure 28 shows the time series of the dimensionless top-point motion. No matter which hydrodynamic model is used, the monopile is excited at its natural frequency, yielding a large ringing response. Similar to the wave loads in Figure 25, a higher peak is observed when using Rainey's model compared with the Morison predictions after the focus time.

Conclusions
The ringing response of a monopile to higher-harmonic wave loads induced by fully nonlinear focused wave groups was investigated. This study reveals that for the steepest Figure 27. Peak normalized moments as a function of wave steepness, computed using different hydrodynamic models with focused wave groups generated in the HOS-NWT. Figure 28 shows the time series of the dimensionless top-point motion. No matter which hydrodynamic model is used, the monopile is excited at its natural frequency, yielding a large ringing response. Similar to the wave loads in Figure 25, a higher peak is observed when using Rainey's model compared with the Morison predictions after the focus time.
Rainey's model. Figure 27 shows the relationship between the peak nondimensional moment and the incident wave steepness using different hydrodynamic models. The predictions of Rainey's model exceed those of the Morison equation by a margin that increases rapidly with wave steepness. The agreement is better when the axial divergence part is included, as explained above.  Figure 28 shows the time series of the dimensionless top-point motion. No matter which hydrodynamic model is used, the monopile is excited at its natural frequency, yielding a large ringing response. Similar to the wave loads in Figure 25, a higher peak is observed when using Rainey's model compared with the Morison predictions after the focus time.

Conclusions
The ringing response of a monopile to higher-harmonic wave loads induced by fully nonlinear focused wave groups was investigated. This study reveals that for the steepest

Conclusions
The ringing response of a monopile to higher-harmonic wave loads induced by fully nonlinear focused wave groups was investigated. This study reveals that for the steepest case, quantitatively, the largest absolute values of the normalized wave loads M y (i.e., the main crest) using HOS-NWT are approximately 1.36 times that using the linear wave theory, and the largest absolute values of the normalized ringing response X T (i.e., the main trough) using HOS-NWT are 2.68 times that using the linear wave theory. Wavelet-transform-based analyses were performed to investigate the local characteristics of the sectional moment at the base of a monopile and induced ringing response. The main conclusions drawn are as follows.
(1) The magnitude of the damping ratio has a slight influence on the peak wave loads and induced response. However, a larger damping ratio leads to a faster decay of the subsequent resonance response. (2) Nonlinear wave models are important for the prediction of ultimate wave loads and induced responses. The prediction of peak wave loads using nonlinear wave models exceeds that using a linear wave model by a margin that increases rapidly with wave steepness. Although second-order wave theory performs well for ultimate loads, it performs worse for the prediction of the ringing response because the higher harmonics of the moment cannot be captured. (3) Rainey's model gives a larger peak sectional moment than the Morison equation because of slender-body corrections, which are essential for accurate prediction of the wave loads. As a result, a higher peak of the dimensionless top-point motion was observed. The discrepancies in the dimensionless peak moments between the two hydrodynamic models tend to increase with wave steepness. Institutional Review Board Statement: Not applicable.