This paper presents a numerical investigation on the harmonic structure of hydrodynamic forces on a fixed and simplified representative floating production, storage and offloading (FPSO) vessel hull under dispersive phase-focused wave groups. The high-fidelity numerical model utilizes the two-phase flow solver in the open-source toolbox OpenFOAM. A series of cases were computed using the numerical model, where the effects of wave steepness, bow diameter, and length of the FPSO are investigated. It is found that given an FPSO under different wave steepness, the non-dimensional inline force exhibits remarkable similarity in terms of the temporal development. The harmonic structure of the inline force is only weakly dependent on the steepness of the incident wave group and the bow diameter, but strongly dependent on the FPSO length. When k p L = 2.27, where L is the length of the FPSO and kp is the wave number at peak frequency, the incident wave group is diffracted significantly by the FPSO. The entire wave–structure interaction process is largely linear, where transfer between different harmonics is rarely seen. However, when kpL is further reduced to 0.57, globally the disturbance of the FPSO on the far field incident wave group is reduced, but locally a strongly nonlinear flow occurs at the rear of the FPSO, where severe run-up occurs at the downstream stagnation point. Higher-order harmonics of inline forces are excited, and the interaction process becomes much more nonlinear.

## I. INTRODUCTION

Floating production storage and offloading (FPSO) vessels have been used extensively in the offshore oil and gas industry for production and processing of hydrocarbons, and for the storage of oil. With the development of the industry, they are being designed to work in increasingly deep water. In contrast to ships or naval vessels, an FPSO usually operates at specific locations and do not sail away even during adverse weather conditions. Therefore, it is of vital importance to examine the survivability of FPSOs under extreme sea states.

It has been shown that, for a given sea state with a specific spectrum, the average shape of the largest and steepest non-breaking wave crests can be represented by a theoretical wave form, which is the normalized auto-correlation function of a random ocean surface based on the underlying spectrum scaled by the crest amplitude.^{1,2} When the distribution of the sea surface elevation follows a Gaussian process, the average shape of the large waves corresponds to the NewWave model, which is typically generated by adjusting the phase of each component and making all the components of a wave group come into phase.^{3,4} The applicability of the NewWave model has been supported by the field measurement from the North Sea, which shows that it is a reasonable model for large waves.^{5}

The NewWave-type focused wave groups have been successfully generated both experimentally and numerically. Laboratory studies of focusing wave train generated by both linear and second-order wavemaker theories were reported, where the spurious free wave components generated by the linear wave-maker theory were discussed.^{6–8} The nonlinear evolution of unidirectional focused wave groups propagating on a flat bottom was investigated both numerically and experimentally, from which it has been recognized that during the formation of focused waves, the nonlinear wave–wave interaction process alters both the amplitude of the wave components and their dispersive properties.^{9–11} Significant energy is transferred into both the higher and lower harmonics, and a steeper wave envelope is produced, in which the central wave crest is higher and narrower, while the adjacent wave troughs are broader and less deep.^{12} Owing to their highly transient characteristics, their broad spectral frequencies and high nonlinearity, it is generally not straightforward to examine the properties of such wave groups. A phase manipulation method was proposed to extract their harmonic structure.^{13} Underlying this method is the assumption that the hydrodynamic force in focused waves possesses a Stokes-like structure, which is the prerequisite for the first four sum harmonics to be separated by phase control and linear combinations of the resultant time-histories.

The evolution of focused wave groups can be further complicated when wave breaking occurs. In the presence of wave breaking, the spectral bandwidth is significantly reduced immediately following breaking, and eventually becomes much smaller than its initial level.^{14,15} The energy dissipation can be quantified based on the breaking event.^{16,17} In addition, directionality has a profound effect upon the nonlinearity of focused wave groups. For increasing the directional spread of the wave field, the energy tends to be redistributed in the frequency domain, leading to less nonlinearity.^{18–20} The evolution of nonlinear focused wave group is also affected by bottom topography.^{21,22} Analysis of focused wave groups on a sloping bed indicates that the crest of the focused wave decreases with the increase in bed sloping. When propagating on the slope, high-frequency short waves are reflected, and the wave energy concentrated at high frequencies is lower for the case of a sloping bed than those for a flat bottom.^{23}

Focused wave groups have been widely used as the representative wave conditions for design of coastal and offshore structures.^{24–29} It is believed that the use of this kind of dispersive phase-focused wave groups can better represent the spectral broadband properties of ocean waves, which is preferred compared to regular waves with corresponding wave height and period.^{30} The resonant fluid response in the gap between two identical fixed rectangular boxes was investigated experimentally in a wave basin under focused wave groups, where the spatial and temporal structure of the resonant fluid response in the narrow gap were investigated. It is shown that for an incident group with appropriate frequency content, the linear gap response may be substantially smaller than the second harmonic component.^{31–33} A fully nonlinear Boussinesq model was adopted to simulate focused transient wave groups and their interactions with the harbor, which reveals the capability of focused transient wave groups to trigger harbor resonance.^{34} Extreme wave run-up in transient focused wave groups around a floating FPSO was studied using computational fluid dynamics. The strong nonlinear local wave field was found to be the main reason for nonlinear extreme run-up, as both low- and high-frequency second harmonic components can lead to wave run-up at significantly higher levels than predicted by a linear analysis. However, the vessel motions are very close to linear.^{35,36}

Transient focused wave groups have also been used as the design wave conditions to examine survivability of marine renewable energy devices in extreme waves. The focused wave trains based on the 100 year wave were numerically reproduced in OpenFOAM, in which a fixed truncated circular cylinder and a floating hemispherical-bottomed buoy were subject to these wave events. The numerical wave tank was demonstrated to have strong capability in predicting the fully nonlinear interaction between extreme waves and wave energy converters.^{37} The NewWave-type wave groups in current were used to test a fully instrumented 1:15 scale tidal turbine, where variations in rotor-based loads, power, and blade root bending moments were reported.^{38} Within the Smoothed Particle Hydrodynamics (SPH) framework, numerical simulations were conducted to model a taut-moored point-absorber wave energy converter together with a linear power takeoff device system in focused wave groups, which seamlessly exploited its functions of energy harvesting and load bearing in such wave conditions.^{39} The response of a multi-body wave energy converter system known as M4 in focused wave groups, which consists of three cylindrical floats with rounded bases, was investigated mainly by experiments. It was shown that the WEC system response is predominantly linear, with weak nonlinearity in beam bending moment.^{40} A series of numerical simulations on focused wave groups interaction with semi-submersible floating offshore wind turbine were performed, and the higher-harmonic wave loads and low frequency resonance response were extracted and analyzed.^{41}

The dispersive phase-focused wave groups with large wave steepness may induce highly nonlinear wave loads on marine structures, comprising a linear harmonic component around the incident spectral peak frequency, and high-order harmonics which are at multiples of the linear peak frequencies, owing to nonlinear wave–wave and wave–structure interactions.^{42,43} For instance, it was reported that under certain conditions, the linear loading on a bottom-mounted circular cylinder is actually less than 40% of the total wave loading and the high-order harmonics contribute more than 60% of the loading, which reveals the importance of high-order nonlinear wave loading on offshore structures.^{30} These higher-order harmonics are non-negligible under such extreme wave conditions and are known to cause highly intense nonlinear structural behaviors, i.e., springing (at double frequency) and ringing (at triple frequency).^{44–46} Therefore, it is of concern for structural design, and improving the understanding of those higher harmonics of wave loading is of importance.

The nonlinearity of the hydrodynamic forces is found to be associated with local diffraction effects, where violent free surface motion has been observed due to steep incident waves. Paulsen *et al.*^{47} numerically computed regular wave forces on a bottom-mounted cylinder in finite water depth ( $ k h \u2248 O ( 1 )$ where *k* is the wave number and *h* is the water depth). Specifically, the second load cycle was discussed, which was shown to be induced by a strong return flow due to overturning of the free surface behind the cylinder. Kristiansen and Faltinsen^{48} also observed the local rear run-up in physical laboratory tests, but no second load cycle was observed, probably due to the relatively small steepness of the incident waves. Actually, two types of nonlinear scattered wave field have been defined.^{49} The first is driven by the run-up and wash-down on the surface of the column in the vicinity of the upstream and downstream stagnation points. The second concerns the circulation of fluid around the column, leading to the scattering of a pair of non-concentric wave fronts. Type-2 scattered wave field introduces an additional timescale that is dependent upon the time taken for the fluid to move around the body. Therefore, this timescale cannot be described by existing perturbation based diffraction solutions. In Paulsen *et al.*,^{47} the force signal was reconstructed using the first sixth harmonics. However, the second load cycle was not captured. Hereby, they conjectured that the second load cycle is more like an indicator of a highly nonlinear local flow around the cylinder, which in turn implies a strong harmonic force content at the lower harmonics.^{47}

Most of the research works as seen above are focused on the forcing of a bottom-mounted circular cylinder rather than an FPSO. For FPSOs, Mai *et al.*^{50} carried out a set of physical experiments to examine the scattered wave field. The nonlinearity of the surface elevations at selected wave gauges around the FPSO-shaped body was evaluated under different parameter sets, but the forces were not measured. The Blind Test Series 1, which was part of the Collaborative Computational Project in Wave–Structure Interaction (CCP-WSI), simulated blindly the interaction between a fixed FPSO-shaped structure and focused waves ranging in steepness and direction. Numerical results from models with different fidelities were compared against corresponding physical data, and the predictive capability of each method was assessed based on pressure and run-up measurements.^{51}

The major difference between an FPSO and a surface-piercing cylinder is that the FPSO-shaped body does not penetrate to the sea bottom. Hence, flow is allowed to develop beneath the body. Furthermore, an important new parameter introduced in this problem is the global length of the FPSO. It may significantly affect the flow in the rear part and alter the harmonic structure of the inline force on the FPSO, since the time taken for the fluid to move around the body is changed. However, this is yet addressed in previous works.

In the present paper, a series of numerical simulations are performed using the two-phase flow solver in OpenFOAM. Force history and scattered wave field are extracted from the numerical simulations to evaluate the effects of several parameters on the harmonic structure of the inline forces under phase-focused wave groups. This paper is organized as follows. The setup of the numerical model is described in Sec. II, and its validation against experimental data is given in Sec. III. The analysis on the harmonic structure of the inline force is given in Sec. IV, where a large number of cases are performed to examine the effect of wave steepness, diameter, and length of FPSO on the harmonic structure of the inline force. In Sec. V, the scattered wave field from FPSO with different lengths is compared, and the difference between them is discussed. Finally, Sec. VI contains the main conclusions.

## II. NUMERICAL MODELS

### A. Governing equations

^{52}

**is the velocity,**

*u***is the gravitational acceleration,**

*g**ρ*is the density,

*μ*is the viscosity, and $ p *$ is the excess pressure, where the hydrostatic pressure has been subtracted. A scalar field of volume fraction field

*α*is used to track the interface between two phases, where

*α*= 0 represents pure air and

*α*= 1 gives pure water, with any intermediate value representing a mixture. $ u r$ here is referred to as the compression velocity,

^{53}which aids in retaining a sharp interface. The term $ \alpha ( 1 \u2212 \alpha )$ vanishes everywhere except at the interface.

Equations (1) are solved within the open-source CFD environment OpenFOAM (ESI version 1706). The time derivatives are discretized using a first-order implicit Euler scheme. Convection terms are discretized using a blend of central difference and upwind schemes, depending on the ratio of the successive gradients. The remaining terms are discretized with second-order central difference schemes. In all the forthcoming cases, an adjustable time step is applied such that a maximum Courant number of 0.25 is maintained at all times.

### B. Setup of numerical model

The numerical wave tank was set up with a dimension of $ 12 \xd7 3 \xd7 3.43 \u2009 m 3$. The coordinate system in the numerical wave tank was defined that the x axis was aligned with the wave propagation direction, the y axis was in the transverse direction, and the z axis was in the vertical direction. The gravity force acted in the negative *z* direction. The origin of the coordinate was located at the theoretical focal position of the wave group based on the linear wave theory.

^{54}OlaFlow applies fixed wave-making boundary conditions, where the flux is directly introduced into the computational domain based on the user-prescribed surface elevation and velocity profile. We applied the second-order irregular wave theory

^{55}on the wave-making boundary to drive the numerical wave tank generating phase-focused wave groups, of which the surface elevation is expressed as

*a*,

_{i}*ω*,

_{i}*k*, and

_{i}*ε*are the amplitude, circular frequency, wave number, and initial phase of each wave component.

_{i}*x*and

_{f}*t*are the

_{f}*x*coordinate of the focal position and focal time, and

*x*is the

_{P}*x*coordinate of the center of the wave paddle. The second-order irregular wave theory employed is actually an extension of the second-order bi-chromatic wave theory, which contains the wave–wave interaction between two wave components. This results in generation of super and sub-harmonics, referred to as bound waves. Here,

*B*

^{+}and

*B*are the transfer functions for super and sub-harmonics, indicating the strength of the interactions between wave components.

^{−}The incident wave groups were focused at 12 s at the bow of the FPSO. At the other side, the active absorption boundary in OlaFlow was applied to absorb the diffracted waves. The geometry of the FPSO was simplified as a rectangular box with a half circular cylinder at the bow and stern. The rectangular box was set to be 0.9 m long and 0.3 m wide. The connecting cylinder had a radius of 0.15 m, and the height of the FPSO was 0.3 m. It was half immersed at a distance of 6 m from the bow of the FPSO to the wave-maker side.

A snapshot of the mesh structure is given in Fig. 1. The utility snappyHexMesh was applied to generate the mesh. It first generated a background hex mesh covering the whole domain. The background mesh was refined near the free surface area in the vertical direction. In the transverse and horizontal directions, the mesh was stretched to provide locally high resolution at the free surface and around the FPSO. Then, the stereolithography surface of the FPSO was utilized by snappyHexMesh which removed all the cells inside the FPSO, and moved the cell vertex points onto the surface geometry to remove the jagged castellated surface from the mesh.

### C. Separation of harmonic structures

Phase-focused wave groups are highly transient events. Thus, for such cases, Fourier transformation may not be suitable to extract the spectrum structure of the hydrodynamic signals. Instead, either the phase-inversion method^{12,13,30,50} or wavelet transformation^{56} is used to separate the higher-order terms. In the present paper, the former method was adopted. We compared the results using different bandwidth for the digital filter and different combinations of phase separation in the Appendix. Eventually, the two-phase based harmonic separation method with a digital filter centered in *nf _{p}* and a bandwidth of

*f*was selected to extract the higher-order harmonics of the hydrodynamic signals in this paper, where

_{p}*f*is the peak frequency of the spectrum and

_{p}*n*represents the integer of harmonics ranging from one to four.

## III. MODEL VERIFICATION AND VALIDATION

### A. Experiments

The experimental work was carried out in the Ocean Basin at the COAST Lab in Plymouth University, UK.^{50} The ocean basin has a length of 35 m and a width of 15.5 m. At the wavemaker, the water depth was set to 4 m. There was a slope connecting the basin bottom from the wavemaker to the working area, in which the water depth was set to 2.93 m. At the far end of the basin, a parabolic absorbing beach was applied to absorb waves. The geometry of the FPSO was exactly the same as in the numerical model. It was half immersed at a distance of 13.886 m from the bow of the FPSO to the wave-maker side.

The wave condition for validation of the model is given in Table I. The wave was focused at 13.886 m downstream from the wave paddle. Each wave group was created using linear superposition of 144 wave fronts with frequencies evenly spaced between 0.1 and 2 Hz. The amplitudes of the frequency components were derived by applying the NewWave theory to a JONSWAP spectrum (peak enhancement factor $ \gamma = 3.3$). A number of wave gauges were installed at different places to measure the incident and diffracted waves. Their positions are presented in Fig. 2, where Gauge 05 and Gauge 19 were deliberately neglected in the experiments. Meanwhile, the local pressure at the bow of the FPSO was also measured by a number of pressure sensors whose positions are depicted in Fig. 3.

### B. Incident phase-focused waves

#### 1. Convergence analysis

Convergence analysis was performed to examine the influence of grid resolution on the surface elevation profile for the trough-focused wave group. Totally four geometrically similar grid resolutions were tested, and the relevant parameters are given in Table II.

Grid no. . | Mesh size . | Mesh no. (×10^{6})
. | Deviation (%) . | ||
---|---|---|---|---|---|

$ \Delta x$ (mm) . | $ \Delta y$ (mm) . | $ \Delta z$ (mm) . | |||

1 | 25 | 57 | 10 | 1.8 | 3.04 |

2 | 22 | 50 | 8.5 | 2.6 | 1.79 |

3 | 20 | 44 | 7 | 3.4 | 0.89 |

4 | 17 | 36 | 5.9 | 4.5 | ⋯ |

Grid no. . | Mesh size . | Mesh no. (×10^{6})
. | Deviation (%) . | ||
---|---|---|---|---|---|

$ \Delta x$ (mm) . | $ \Delta y$ (mm) . | $ \Delta z$ (mm) . | |||

1 | 25 | 57 | 10 | 1.8 | 3.04 |

2 | 22 | 50 | 8.5 | 2.6 | 1.79 |

3 | 20 | 44 | 7 | 3.4 | 0.89 |

4 | 17 | 36 | 5.9 | 4.5 | ⋯ |

The surface elevation at the focal position is given in Fig. 4 using different grid resolutions. It is shown that from grid 1 to grid 3, a simulation with a finer grid is likely to produce a deeper trough at the focal position. However, when the mesh resolution further increases from grid 3 to grid 4, the magnitude of the trough is decreased slightly. The deviation of the trough between different grids is also presented in Table II, which decreases consistently with the increase in the mesh resolution. Between grid 3 and grid 4, the deviation is reduced to less than 1%, indicating that the numerical results are convergent and independent of mesh resolution. In the present work, we decide to adopt grid 3 as the standard mesh in the following simulations, which is a good compromise between the computational cost and accuracy.

#### 2. Comparison of surface elevation

Figure 5 directly compares the surface elevation between the numerical and experimental data. In general, the incident wave field is accurately reproduced by the numerical model. The root mean square deviation of the trough at focal time for the five gauges is 0.0037 m, effectively 5.36% of the crest of the wave group *A*. The spectrum structures for the odd and even harmonics of the surface elevation are compared in Fig. 6. For surface elevation at all the wave gauges, we observe a good agreement for the linear harmonic, but a consistent overestimation of second-order harmonic from the numerical model, and consistent underestimation of the fourth-order harmonic from the numerical model. However, their magnitudes are significantly smaller than the linear harmonic, and the overall agreement between the numerical and experimental results is good.

### C. Numerical wave tank with FPSO in place

In the same way as in Sec. III B 2, we compare the surface elevation and its spectrum structure at the same locations for the numerical wave tank with the FPSO in place. These are given in Figs. 7 and 8, respectively. Generally, the agreement is good before and during the focusing process, though at the focal time, the trough is slightly underestimated by the numerical model at all the gauges. The root mean square deviation of the trough at focal time for the five gauges is 0.0042 m, corresponding to 6.1% of the crest of the wave group *A*. However, at defocusing process, we observe a larger deviation, where the numerical model mostly overestimates the surface elevation, particularly at WG09, WG10, and WG16. The root mean square deviation of the peak at the last two cycles of all the gauges reaches 9.23%. Figure 8 shows the spectrum structure of the surface elevation from the numerical model and the experiments. Similar to Fig. 6, it is found that overestimation of the surface elevation is primarily on the higher order harmonics, rather than the linear harmonic.

Figure 9 shows the hydrodynamic pressure at PS 03 and PS 06. The local pressure at the focal time $ t = 12 \u2009 s$ is found to be zero, due to the fact that the free surface is below the pressure sensor, since the wave group is trough-focused. Adjacent to the trough, the neighboring peaks are accurately captured by the numerical model. However, further away from the focal time, the numerical model slightly underestimates the peaks at the focusing process but overestimates them at the defocusing process. The root mean square deviation of the peak at the last two cycles of the pressure at these two gauges is 79.4 Pa. Taking 968 Pa as a reference, which is the maximum pressure appearing in the series, the deviation is 8.2%, in the same order of magnitude as the surface elevation.

### D. Brief summary

Up to now, we have successfully validated our numerical model against the available experimental data for both cases, i.e., with an empty wave tank and with the FPSO in place. Good agreement has been achieved for the surface elevation in empty wave tank, where the numerical model can reproduce it at different locations accurately up to the fourth order. For the case with the FPSO in place, the model can still reproduce the local pressure and surface elevation adequately, though overestimation of the elevation is noticed at defocusing process. In the following section, this model will be applied to simulate the phase-focused wave groups interaction with the FPSO under a variety of different conditions.

## IV. INTEGRATED INLINE HYDRODYNAMIC FORCES

### A. Test conditions for the present study

*A*is the wave crest at the focal position,

*h*is the water depth,

*λ*is the wavelength at the peak frequency,

_{p}*d*is the draft of the FPSO,

*D*and

*L*are the diameter and the length of the FPSO, and

*α*is the wave propagation direction. Then, by dimension analysis, the non-dimensional hydrodynamic forces are written as a function of the following parameters:

*k*is the wave number at the peak frequency, and

_{p}*Re*is the Reynolds number.

As introduced in Sec. I, FPSOs are usually operated in deep water. Hereby, we neglect the water depth effects. In addition, the effect of the draft and the wave propagation angle is also outside the scope of the present paper. Only normal wave groups are considered, and the FPSO is half immersed in still water. The wave groups are always made to focus at the bow of the FPSO, corresponding to the most severe condition. Due to the wave–wave interaction, the actual focal position will be downshifted from the theoretical focal position, which is found by examining the surface elevation near the theoretical focal position. The location that produces the largest surface elevation represents the actual focal position for that wave train, and it is also the location at which the FPSO will be placed. The main focus in the present paper is to investigate the effects of wave steepness and size of the FPSO. Totally five sets of cases are computed, and a complete list of all the cases is given in Table III. Within each set there are four cases, where the non-dimensional wave steepness *k _{pA}* ranges between 0.13 and 0.22. Across sets A–C, the effects on the diameter of the bow of FPSO can be revealed. Meanwhile, the effect of the length of the FPSO is depicted by comparisons between sets A, D, and E.

ID . | A (m)
. | T (s)
. _{p} | h (m)
. | k (-)
. _{pA} | D (m)
. | L (m)
. | $ k p D \u2009 ( m )$ . | k (-)
. _{pL} | KC (-)
. | Re (10^{5})
. |
---|---|---|---|---|---|---|---|---|---|---|

A-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.3 | 1.2 | 0.57 | 2.27 | 0.56 | 5.54 |

A-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.3 | 1.2 | 0.57 | 2.27 | 0.77 | 6.33 |

A-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.3 | 1.2 | 0.57 | 2.27 | 0.88 | 8.64 |

A-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.3 | 1.2 | 0.57 | 2.27 | 0.97 | 9.60 |

B-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.2 | 1.2 | 0.38 | 2.27 | 0.56 | 5.54 |

B-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.2 | 1.2 | 0.38 | 2.27 | 0.77 | 6.33 |

B-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.2 | 1.2 | 0.38 | 2.27 | 0.88 | 8.64 |

B-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.2 | 1.2 | 0.38 | 2.27 | 0.97 | 9.60 |

C-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.1 | 1.2 | 0.19 | 2.27 | 0.56 | 5.54 |

C-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.1 | 1.2 | 0.19 | 2.27 | 0.77 | 6.33 |

C-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.1 | 1.2 | 0.19 | 2.27 | 0.88 | 8.64 |

C-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.1 | 1.2 | 0.19 | 2.27 | 0.97 | 9.60 |

D-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.3 | 0.6 | 0.57 | 1.14 | 1.12 | 2.77 |

D-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.3 | 0.6 | 0.57 | 1.14 | 1.53 | 3.80 |

D-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.3 | 0.6 | 0.57 | 1.14 | 1.76 | 4.36 |

D-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.3 | 0.6 | 0.57 | 1.14 | 1.94 | 4.81 |

E-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.3 | 0.3 | 0.57 | 0.57 | 2.24 | 1.39 |

E-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.3 | 0.3 | 0.57 | 0.57 | 3.07 | 1.90 |

E-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.3 | 0.3 | 0.57 | 0.57 | 3.53 | 2.20 |

E-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.3 | 0.3 | 0.57 | 0.57 | 3.89 | 2.41 |

ID . | A (m)
. | T (s)
. _{p} | h (m)
. | k (-)
. _{pA} | D (m)
. | L (m)
. | $ k p D \u2009 ( m )$ . | k (-)
. _{pL} | KC (-)
. | Re (10^{5})
. |
---|---|---|---|---|---|---|---|---|---|---|

A-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.3 | 1.2 | 0.57 | 2.27 | 0.56 | 5.54 |

A-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.3 | 1.2 | 0.57 | 2.27 | 0.77 | 6.33 |

A-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.3 | 1.2 | 0.57 | 2.27 | 0.88 | 8.64 |

A-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.3 | 1.2 | 0.57 | 2.27 | 0.97 | 9.60 |

B-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.2 | 1.2 | 0.38 | 2.27 | 0.56 | 5.54 |

B-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.2 | 1.2 | 0.38 | 2.27 | 0.77 | 6.33 |

B-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.2 | 1.2 | 0.38 | 2.27 | 0.88 | 8.64 |

B-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.2 | 1.2 | 0.38 | 2.27 | 0.97 | 9.60 |

C-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.1 | 1.2 | 0.19 | 2.27 | 0.56 | 5.54 |

C-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.1 | 1.2 | 0.19 | 2.27 | 0.77 | 6.33 |

C-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.1 | 1.2 | 0.19 | 2.27 | 0.88 | 8.64 |

C-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.1 | 1.2 | 0.19 | 2.27 | 0.97 | 9.60 |

D-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.3 | 0.6 | 0.57 | 1.14 | 1.12 | 2.77 |

D-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.3 | 0.6 | 0.57 | 1.14 | 1.53 | 3.80 |

D-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.3 | 0.6 | 0.57 | 1.14 | 1.76 | 4.36 |

D-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.3 | 0.6 | 0.57 | 1.14 | 1.94 | 4.81 |

E-1 | 0.069 | 1.456 | 2.93 | 0.13 | 0.3 | 0.3 | 0.57 | 0.57 | 2.24 | 1.39 |

E-2 | 0.090 | 1.456 | 2.93 | 0.17 | 0.3 | 0.3 | 0.57 | 0.57 | 3.07 | 1.90 |

E-3 | 0.105 | 1.456 | 2.93 | 0.20 | 0.3 | 0.3 | 0.57 | 0.57 | 3.53 | 2.20 |

E-4 | 0.120 | 1.456 | 2.93 | 0.22 | 0.3 | 0.3 | 0.57 | 0.57 | 3.89 | 2.41 |

In the present paper, considering that an FPSO usually operates in the offshore area with rather large water depth, *k _{ph}* is selected as 5.52 for all the cases, i.e., the water depth effect has been neglected as mentioned above. The non-dimensional wave steepness $ k p A$ ranges between $ 0.13 \u2009 and \u2009 0.22$, which is in general the same as in other publications. The size of the FPSO is chosen so that it crosses the linear diffraction regime where $ \lambda p / L < 5$, and the inertia-drag regime where $ \lambda p / L > 5$. The

*KC*number in our case is generally lower than in other works, due to the difference of the FPSO shape from a circular cylinder (note that the FPSO length is used as the characteristic length for calculating the

*KC*number).

### B. Temporal development of the force curve

The inline hydrodynamic forces induced by the phase-focused wave groups are depicted in Fig. 10 for the cases in sets A–E. For the cases with the same FPSO but under different wave steepness, the non-dimensional forces almost collapse perfectly into a single line, indicating their remarkable similarity in the time domain. By comparing sets A–C, it is found that the force signals exhibit quite different amplitudes near the focal time, due to the different bow diameter of the FPSO. The crest of the force is increasing with the growing diameter. However, the general shape of the force curves is rather similar, which behave like the same force curve just with a different scaling factor, indicating that the harmonic structure of the forces should be quite similar to each other.

However, regarding sets A, D, and E, the temporal development of the forces is completely different from each other, indicating distinct harmonic structures of the force with each FPSO length. With the decrease in the FPSO length, the signal is becoming more flat in the crest as shown in $ ( t \u2212 t f ) g k p \u223c$ [−4 −2] and [2 4], where the asymmetry of the signal is increasing. This indicates that the nonlinearity is likely to increase with decreasing length of the FPSO, which will be further analyzed in detail in Sec. IV C.

### C. Harmonic structure of the force

The harmonic structures of the hydrodynamic force signals are extracted using the two-phase based separation method, as introduced in Sec. II C. The relative contribution of each harmonic is presented in Fig. 11. Here, we have several comments: (1) The linear harmonic is generally dominant over the higher-order harmonics for all cases. The relative contribution ranges from 70% to 90% depending on the variation of the parameters. The nonlinearity of the hydrodynamic force in the present paper is generally weaker than for the circular bottom-mounted cylinder with similar diameter and wave conditions.^{30} (2) The third- and fourth-order harmonics are of minor importance, and the contribution is mostly negligible. For the cases in set E, which shows the strongest nonlinearity for the hydrodynamic forces, the third- and fourth-order harmonics are around 1%–3% of the total forces. (3) Note that the wave steepness is within the range of 0.13–0.24. Within such a range, we notice a weak dependence of the harmonic structure of the hydrodynamic forces on the wave steepness. The linear component drops within 5%. (4) The harmonic structure of the hydrodynamic force is also weakly dependent on the bow diameter. When $ k p D$ ranges between 0.19 and 0.57, the linear harmonic varies within 8%. (5) The harmonic structure of the hydrodynamic force is strongly dependent on the length of the FPSO. By comparing between cases C-1 and E-1, it is found that the linear harmonic decreases dramatically from 90% to about 70%, given that the FPSO length $ k p L$ is reduced from 2.27 to 0.57.

Figure 12 presents a detailed example on the effect of FPSO length on the harmonic structure of the force, where the first four harmonics of the forces in cases C-1 and E-1, and the incident wave group are plotted. The incident wave group is actually weakly nonlinear, where the second-order harmonic takes around 5%. From the middle column of Fig. 12, it can be seen that the harmonic structure of the force in case C-1 is quite similar to the incident wave group, where the FPSO with a global length of 1.2 m is placed in the wave tank. The second-order harmonic of the hydrodynamic force is also about 5% of the total force. This indicates that the linear wave–structure interaction process is dominant in this case, where the wave force at each harmonic is directly associated with the incident wave group at that harmonic. However, for case E-1, the harmonic structure of the hydrodynamic force is significantly altered, when the global length of the FPSO is reduced from 1.2 to 0.3 m, i.e., the FPSO is shortened to a circular cylinder. In this case, the linear harmonic is reduced to 74% of the total hydrodynamic force, while the second-order harmonics are increased to 17%, and third and fourth harmonics are also slightly increased by 4%, respectively. This shows that the higher-order harmonics of the hydrodynamic force are generated not only due to the incident wave group at that harmonic, but also due to transfer of energy from the fundamental harmonic.

In Fig. 13, the *n*th order harmonics are scaled by $ \rho g A n D 2 \u2212 n L$. We observe an almost constant slope over the wave steepness across all the sets, except the fourth-order harmonic in set E. Specifically, for the linear force component, we notice that the non-dimensional forces are well collapsed into one line between sets A–C. This actually indicates that the Morison equation is capable of predicting the linear force component with a suitable but single force coefficient for the FPSO with different diameters. However, when the length is varying, the coefficient should also be changed. The magnitude of the second-order harmonics well exhibit the nonlinearity of the wave forces in each set. Set E, of which the forces are associated with the strongest nonlinearity, produces the largest non-dimensional second-order force.

## V. SCATTERED WAVE FIELD

The nonlinearity of the hydrodynamic force is likely to be associated with the nonlinear free surface motions, as introduced in Sec. I. In this section, we will carefully investigate the effects of the length of FPSO on the local free surface motion by comparing cases A-2, D-2, and E-2, where the same incident wave group is applied. Special attention will be given to the run-up at the bow and the stern, since the pressure on these parts significantly contributes to the inline force. As introduced in Sec. I, in the laboratory tests, two types of high-frequency wave scattering phenomenon were observed.^{49} Both of which are found in the present simulations, and a detailed analysis is given below on these two types of wave scattering and its association with the harmonic structure of the inline force on the FPSO.

### A. Upstream wave scattering

Type-1 wave scattering is observed in the upstream direction when the wave crest arrives the cylinder as shown in Fig. 14. The incident wave group creates a wave run-up upon arrival. Then, it washes back down and induce a disturbance on the waves. Essentially the disturbance radiates out in the upstream direction as a concentric wave field.

By comparison between the three different cases, it is found that the length has a minor effect on the upstream wave scattering field. The three models produce quite similar concentric wave fields. This can further be demonstrated in Fig. 15, which plots the surface elevation at the upstream stagnation point. It is shown that the maximum run-up reaches $ 1.4 A$ for all the three cases. Actually, the upstream scattering should not be affected by the length of the FPSO significantly, since the wave condition and the front part of the three FPSOs are exactly the same. Therefore, it is inferred that the distinct harmonic structure of the force for each case should be implied by the differences in the downstream scattering field.

### B. Downstream wave scattering

Type-2 wave scattering field is observed in the rear part of the FPSO as shown in Fig. 16. When the crest of the incident wave group arrives at the FPSO, along with the Type-1 wave scattering which mainly occurs near the upstream stagnation point, a low pressure cavity is formed behind the cylinder. Therefore, the local diffracted wave is driven to travel around the FPSO circumference in both clockwise and counterclockwise directions. Due to the disturbance of the FPSO, the wave is likely to accelerate along the two sides of the FPSO. When they eventually arrive on the back face of the FPSO, they collide and then merge together, creating significant wave run-up at that point. As the incident wave advances, the disturbance is further swept away, forming a symmetric but non-concentric pattern of scattered waves.

The surface elevations around the FPSO at different instants are presented in Fig. 17, where *θ* is the angle representing the position of the elevation. The rear run-up at the stern is in between $ \theta = [ 180 \xb0 \u2013 270 \xb0 ]$. Substantial differences in the strength of the Type-2 scattered wave field have been noticed, where particularly strong scattering is observed for case E-2, which eventually form a sharp point at $ ( t \u2212 t f ) g k p = 0.43$, whereas the run-up is quite smooth even it reaches the highest at $ ( t \u2212 t f ) g k p = 3$ for case A-2. This can further be demonstrated by the run-up at the downstream stagnation point, as shown in Fig. 18. We notice that this rear run-up is well correlated with the nonlinearity of the inline force on the FPSO, where a larger rear run-up implies a more nonlinear harmonic structure of the inline force.

## VI. CONCLUSIONS

A numerical model has been set up and applied to investigate the interaction of nonlinear phase-focused wave groups with a fixed FPSO-shaped body under various wave conditions and FPSO sizes. The numerical model has been successfully verified and validated against available experimental data. For a wave group propagating in an empty wave tank, the surface elevation is captured accurately up to the fourth order. For the numerical wave tank with the FPSO in place, slight overestimation is observed from the second-order harmonic onward. In addition, the local pressure on the FPSO is compared with the experimental data, and the deviation is in the same order of magnitude as it is for the surface elevation.

Due to the highly transient characteristics, a phase-inversion method is employed to extract the harmonic structure of the hydrodynamic forces on the FPSO. A series of cases have been computed, and the effects of wave steepness, the diameter of the bow of the FPSO, and the length of the FPSO are carefully examined. For the same FPSO but with different wave steepness, the temporal development of the inline forces is quite similar. The harmonic structure of the force is only weakly dependent on the wave steepness and the bow diameter. Meanwhile, under the same wave conditions, the harmonic structure of the inline force is strongly dependent on the length of the FPSO. With decreasing length, the higher order harmonics are shown to grow rapidly. By looking into the scattered wave field and the free surface motion, it is found that very strong run-up occurs at the rear part of the FPSO when the length is reduced. This rear run-up well implies the nonlinearity of the inline force on the FPSO.

## ACKNOWLEDGMENTS

The physical data presented here were generated as part of a UK Engineering and Physical Sciences Research Council (EPSRC) funded research project (No. EP/J012866/1) and are distributed via the Collaborative Computational Project in Wave Structure Interactions (CCP-WSI) (No. EP/M022382/1) through the CCP-WSI project website (http://www.ccp-wsi.ac.uk/). This work is partially funded by the EPSRC (UK) projects “A Zonal CFD Approach for Fully Nonlinear Simulations of Two Vessels in Launch and Recovery Operations” (No. EP/N008839/1) and “Extreme Loading on FOWT under Complex Environmental Conditions” (No. EP/T004150/1).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Hao Chen:** Conceptualization (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Validation (lead); Writing – original draft (lead). **Ling Qian:** Conceptualization (lead); Funding acquisition (lead); Project administration (lead); Resources (lead); Writing – review & editing (equal). **Deping Cao:** Resources (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: PHASE-INVERSION METHOD

*x*,

_{p}*t*). Therefore, significant nonlinearity only appears around the focal time in the neighborhood of the focus location. Stokes-type structures are used to approximate the focused wave groups, where the idea is to express the signal envelope as a slowly varying function of time

_{p}*A*(

*t*), assuming that the wave energy spectrum is sufficiently narrow-banded. For instance, the hydrodynamic force at the focal location around the focal time can be simplified to

*f*is a wave to force transfer function.

*F*

_{0}and

*F*

_{180}. The odd and even harmonics are given as

*H*is the Hilbert transform used to shift by $ 90 \xb0$ every Fourier component within a signal.

The two-phase based method only requires two time histories, while four time series for the same case must be available in order to employ the four-phase based separation method. Hence, the two-phase based method significantly reduces the computational/experimental work load when a large number of cases are investigated. However, a digital filtering process is needed in order to separate each harmonic from the odd and even harmonics, which may introduce energy leakage between harmonics owing to the smearing of the spectrum. Although in the four-phase based harmonic separation method, the fourth-order harmonic is also filtered from *F*_{24}, minimum energy leakage is expected since there is little overlapping between these two harmonics.^{13}

Figure 19 presents an example of four time series of hydrodynamic forces that are $ 90 \xb0$ out of phase, which are the basis to separate each harmonic. The spectrum of the odd and even harmonics used in the two-phase based method, and the spectrum of the first three harmonics used in the four-phase based method are shown Fig. 20. In this case, the peaks for the first- and second-order harmonics are clearly visible. The third-order harmonic is on the order of 0.5 N, which is quite small compared with the linear harmonic.

Figure 21 compares the separated first four harmonics using the two-phase and four-phase separation methods. It is noticed that actually the results from two-phase based separation are sensitive to the bandwidth, especially for the third and fourth order terms. An almost identical linear harmonic can be separated from both methods, irrespective of the bandwidth of the digital filter. However, the two-phase based separation method significantly overestimates the third and fourth harmonics, which includes contributions from the odd and even harmonic spectrum over an excessive frequency range including probably some of the tail of the first- and second-order peak. A filter with a range of *f _{p}* can produce relatively reasonable results, the envelope of which agrees better with the four-phase based separation method. Considering the large number of cases carried out in Sec. IV, the two-phase based method together with the digital filter centered in

*nf*with a width of

_{p}*f*is used to extract the higher order harmonics.

_{p}## REFERENCES

^{®}