Uncertainty Assessment of CFD Investigation of the Nonlinear Difference-Frequency Wave Loads on a Semisubmersible FOWT Platform

Abstract

Current mid-fidelity modeling approaches for floating offshore wind turbines (FOWTs) have been found to underpredict the nonlinear, low-frequency wave excitation and the response of semisubmersible FOWTs. To examine the cause of this underprediction, the OC6 project is using computational fluid dynamics (CFD) tools to investigate the wave loads on the OC5-DeepCwind semisubmersible, with a focus on the nonlinear difference-frequency excitation. This paper focuses on assessing the uncertainty of the CFD predictions from simulations of the semisubmersible in a fixed condition under bichromatic wave loading and on establishing confidence in the results for use in improving mid-fidelity models. The uncertainty for the nonlinear wave excitation is found to be acceptable but larger than that for the wave-frequency excitation, with the spatial discretization error being the dominant contributor. Further, unwanted free waves at the difference frequency have been identified in the CFD solution. A wave-splitting and wave load-correction procedure are presented to remove the contamination from the free waves in the results. A preliminary comparison to second-order potential-flow theory shows that the CFD model predicted significantly higher difference-frequency wave excitations, especially in surge, suggesting that the CFD results can be used to better calibrate the mid-fidelity tools.

1. Introduction

Recent years have seen a rapid advancement in offshore wind technologies. Although fixed-bottom systems are still dominant, floating offshore wind turbines (FOWTs) are also being trialed, with several precommercial projects recently deployed or under development. For the most part, the development of such systems relies heavily on computer simulations with simplified engineering models rather than on high-fidelity simulations, especially during the early design/sizing phase of the project. It is still not feasible to perform extensive parametric studies and optimization exclusively with high-fidelity simulations, such as computational fluid dynamics (CFD), because of their prohibitive computational cost; therefore, the success of the offshore wind industry relies on researchers identifying and addressing the limitations of mid-fidelity, engineering-level tools, such as OpenFAST, developed by the National Renewable Energy Laboratory (NREL) [1]. This is the explicit goal of the past Offshore Code Comparison, Collaboration, Continued, with Correlation (OC5) project as well as the ongoing OC6 (OC5 with unCertainty) project under the International Energy Agency (IEA) Wind Task 30 framework.

One major limitation of engineering-level tools commonly used by FOWT designers identified in OC5 and OC6 is the significant underprediction of the low-frequency, nonlinear wave excitation and response of semisubmersible FOWTs [2,3]. Although the nonlinear wave loads are generally much smaller than those from direct (linear) wave excitation, they can potentially trigger large surge and pitch motion in a semisubmersible FOWT, because the surge and pitch resonance frequencies of this type of system are specifically designed to fall in the low-frequency range to avoid direct wave excitation [4,5,6]; therefore, the underprediction of nonlinear wave excitation and the resulting response can have major consequences for the safety and reliability of the FOWT design. As part of OC6 Phase I, a recent paper by Robertson et al. [3] documented the extensive effort made by many institutions worldwide in tuning various engineering-level models to address this particular issue. Some model features, such as the incorporation of full Quadratic Transfer Functions (QTFs) for wave excitation from second-order potential-flow theory and the use of fully nonlinear wave kinematics for the Morison drag computation, have been found to improve results. Nevertheless, it was discovered that varying degrees of underprediction at low frequencies are ubiquitous across engineering-level tools and cannot be addressed simply by further tuning the model coefficients; rather, improvements to the formulation of the engineering models are required. Similar underpredictions have also been encountered by other researchers (see, e.g., [7,8]).

The underprediction of the low-frequency wave-exciting force in surge (in the direction of wave propagation) and moment in pitch (about a transverse axis perpendicular to the wave propagation) on a fixed OC5-DeepCwind semisubmersible is shown in Figure 1. The full-scale simulation results were obtained using an OpenFAST model without any aerodynamic/tower effects to be consistent with the experiment. (For more details on the model and experiment, see [2,3,9].) The low-frequency excitation predicted by the OpenFAST model primarily came from the QTFs of second-order potential-flow theory and the Morison drag. Despite the inclusion of these components, the model still significantly underpredicted the wave excitation for $f<0.05$ Hz. The corresponding motion response of the same floater under a freely floating condition at the surge and pitch resonance frequencies—approximately 0.01 Hz and 0.03 Hz, respectively—were also severely underpredicted [3]. The underprediction of the low-frequency wave excitation on a fixed floater is expected to be at least partially responsible for the underprediction of the resulting motion, with overestimation of the damping of the resonant modes being the other.

To examine the cause of the underprediction and to inform the ongoing effort in improving the engineering models, the OC6 project initiated a study of the problem using CFD tools. Instead of using full irregular waves, the CFD simulations are performed with bichromatic incident waves. Compared with full irregular waves, bichromatic waves have several advantages. First, bichromatic incident waves can be designed to have shorter repeat periods, and the CFD simulations need to run for only a few repeat periods (on the order of 10 minutes in full-scale physical time) rather than a three-hour time window typically required for irregular-wave simulations. The computing time is, therefore, significantly reduced, allowing more wave cases to be investigated. The shorter computing time also makes systematic uncertainty analysis more feasible, which is necessary to ensure that the CFD solutions have adequately converged. Second, the results from CFD simulations with bichromatic waves can be directly compared to the QTFs from second-order potential-flow theory. Such comparisons can help us better understand the limitations of the potential-flow QTFs. In [10], a set of bichromatic-wave pairs were selected for CFD investigation, and the present paper focuses on one key wave case. The two wave frequencies of the bichromatic wave pair investigated in this paper, $f_1$ and $f_2$, are marked in Figure 1 using vertical lines. Both wave frequencies are close to the peak frequency of the wave spectrum. Further, the difference frequency, $f_d=f_2-f_1$, also labeled, is near the surge resonance frequency at 0.01 Hz. We would like to determine whether CFD can better predict the nonlinear wave excitation at this frequency.

Extensive CFD investigations of FOWTs have been reported in the literature. One major motivation to study FOWTs using CFD is to better quantify the viscous damping characteristics of the floater to help tune engineering-level tools. For example, Burmester et al. extracted surge-damping coefficients from the CFD simulations of the forced surge motion of a semisubmersible FOWT [11]. A similar investigation with forced heave and pitch motion was done by Bozonnet and Emery [12]. Alternatively, the damping coefficients can also be extracted from simulations of free-decay motions. For example, the pitch decay of the OC5-DeepCwind platform was simulated by Wang et al. using CFD [13]. Floater response to waves—mostly regular waves—is another focus of CFD investigation. Tran and Kim investigated both free decay and the regular wave-induced motion of the OC5-DeepCwind semisubmersible with comparisons between mid-fidelity, engineering-level tools and higher fidelity CFD simulations [14]. Wang et al. computed the surge, heave, and pitch response amplitude operators for the same floater geometry by simulating the motion of the floater with various incident wave frequencies [15]. Del Águila Ferrandis et al. computed the response amplitude operators of an unconventional FOWT platform in regular waves using CFD and compared the results to both frequency- and time-domain boundary element methods [16]. More recent works incorporated fully coupled aero- and hydrodynamics to investigate the effect of wave-induced motion on the performance and operation of the mounted wind turbines [17,18,19,20].

The present investigation, however, differs from prior CFD efforts by placing a special focus on the nonlinear difference-frequency wave excitation. This leads to a unique set of challenges resulting from the smallness of the nonlinear loads, which can be two orders of magnitude lower than the direct (linear) wave excitation. As a result, it is important to obtain reasonable uncertainty estimates for the CFD solution. Without proper uncertainty quantification, it is difficult to perform meaningful cross-verification among CFD simulations and validation against experimental measurements for the small second-order quantities, which are the two principal objectives of the OC6 project. For the present paper, extensive convergence studies have been carried out for the selected bichromatic wave case to obtain estimates of the various types of numerical uncertainties, which include iterative uncertainty coming from, among other things, the iterative and segregated solution of the nonlinear field equations and pressure-velocity coupling, as well as temporal and spatial discretization uncertainties. For the present problem, spatial discretization error is expected to play an important role in the overall numerical uncertainty because of the potential effect of fluid viscosity on the nonlinear, difference-frequency wave loads. It is often challenging to accurately capture the vortical flow structures generated by fluid viscosity numerically, especially for practical problems where it is generally not feasible to resolve all flow scales. In such situations, the solution of the vortical flow can be sensitive to the numerical algorithm [21]. For instance, Drikakis and Smolarkiewicz investigated the formation of spurious vortical structures on underresolved grids where the truncation error from the discretization of the advective terms acts as a spurious source of vorticity [22]. To quantify the various sources of numerical uncertainties, recommended uncertainty quantification practices from the literature have been adapted for the present application (see, e.g., [23,24,25]). The uncertainty estimates obtained will continue to be used in ongoing and future studies performed as part of OC6. The uncertainties will be especially important in the future when the present CFD solution will be validated against the measurements from an upcoming experimental validation campaign. Finally, a preliminary comparison based on the bichromatic-wave case studied shows that the nonlinear, difference-frequency wave excitation predicted by the present CFD model is significantly higher than that from the second-order potential-flow theory, especially in surge.

2. Problem Description

To investigate the nonlinear low-frequency wave excitation, participants of the OC6 project are performing CFD simulations to study the wave loads on a fixed semisubmersible FOWT platform in bichromatic incident waves. The present paper presents the CFD study by NREL along with the relevant uncertainty analyses. The obtained uncertainty band will continue to be used in the next stage of the OC6 CFD project involving validation of the CFD results against experiment. All simulations were carried out at 1:50 scale to be consistent with an upcoming experimental validation campaign; however, all dimensional values presented in the paper are at full scale based on Froude scaling.

For the floater geometry, a simplified version of the OC5-DeepCwind FOWT platform [2] is used. The geometry of the semisubmersible and the adopted coordinate system are both illustrated in Figure 2. Compared to the original OC5 semisubmersible, the central main column as well as the cross members (pontoons and braces) between the columns have been omitted. This simplified floater geometry is consistent with the physical model to be used in the upcoming experimental validation campaign.

The three columns of the semisubmersible are identical to each other. The upper part of each column has a diameter of $D_U=12$ m and a draft of 14 m. The heave plates (base columns) have a diameter of $D_B=24$ m and a height of 6 m, leading to a total draft of 20 m for each column. The three columns are arranged to form an equilateral triangle with a center-to-center distance between columns of $L=50$ m. The origin of the coordinate system is on the calm water surface at the geometric center of the equilateral triangle. The positive $x$-axis is in the direction of the incident wave propagation, whereas the $+y$-axis and $+z$-axis point starboard and vertically upward, respectively. The mean water depth is 250 m.

The first-order (linear) properties of the unidirectional bichromatic incident wave pair are listed in

Table 1. Linear-wave properties of the bichromatic incident waves.

Wave	Period, T	Frequency, f	Wavelength, λ	Wave Height, H	Repeat Period, TR
1	11.9000 s	0.08403 Hz	221.1 m	3.569 m	190.4 s
2	10.5778 s	0.09454 Hz	174.7 m	3.703 m

. In the rest of the paper, the subscripts 1 and 2 are used when referring to the wave properties to denote the first (lower frequency) and second (higher frequency) components of the bichromatic wave pair. The two wave components selected yield a difference frequency of $f_d=f_2-f_1=0.0105$ Hz, which is very close to the surge natural frequency of the OC5-DeepCwind semisubmersible [2]. Of course, the present paper considers only a fixed floater; therefore, the natural frequencies of floater motion have no direct effect. The repeat period of the bichromatic waves under linear wave assumption is $T_R=190.4$ s. Both wave components are well within the deepwater limit, with the water depth greater than the wavelengths.

The Keulegan–Carpenter (KC) number of the fixed floater in the bichromatic waves can be estimated as:

$$\mathrm{KC}=2\pi \frac{\left(A_1+A_2\right)}{D_U}=1.9,$$

(1)

where $A_1$ and $A_2$ are the wave amplitudes of the first and second components of the bichromatic wave pair. The contributions from the two wave components are combined conservatively through direct summation. At this KC number, inertial effects dominate (see, e.g., [26]), and we do not expect the formation of viscosity-induced asymmetric wakes from the columns.

3. Numerical Method

The commercial software tool STAR-CCM+ [27] was used to perform the CFD simulations. The incompressible Navier–Stokes equation was solved numerically, along with the continuity equation, using the finite-volume method:

$$\frac{\partial u}{\partial t}+\left(u\cdot \nabla \right)u=-\frac{\nabla p}{\rho }+g+\nu {\nabla }^{2}u,$$

(2)

$$\nabla \cdot u=0,$$

(3)

where $u$ is the flow velocity, $p$ is the total pressure, $g$ is the gravitational acceleration, $\rho$ is the density of the fluid, and $\nu$ is the kinematic viscosity of the fluid. The variables in bold are vectors. The field equations were solved in a segregated fashion using the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) for pressure-velocity coupling. During our investigation, we found the common turbulence models currently available in STAR-CCM+ to be too dissipative, with the incident waves and, by extension, the wave loads on the floater decaying continuously with time without ever reaching a steady value. Interestingly, this issue appears to be more pronounced with bichromatic waves than with regular waves. Similar issues with turbulence models in free-surface flows have also been encountered by other researchers (see, e.g., [28]). As a result, the simulations were performed without any turbulence model or wall functions.

The two-phase flow problem (water and air) was modeled using the volume-of-fluid formulation [29], where the phase fraction, ${\alpha }_w$, represents the volume fraction of water in each grid cell. The local fluid properties are given by a mix of air and water properties based on the volume fraction:

$$\rho ={\rho }_a\left(1-{\alpha }_w\right)+{\rho }_w{\alpha }_w,$$

(4)

$$\mu ={\mu }_a\left(1-{\alpha }_w\right)+{\mu }_w{\alpha }_w,$$

(5)

where $\mu$ is the dynamic viscosity. The subscripts $a$ and $w$ denote air and water, respectively. The evolution of ${\alpha }_w$ is given by the following scalar transport equation:

$$\frac{\partial {\alpha }_w}{\partial t}+\left(u\cdot \nabla \right){\alpha }_w=0.$$

(6)

As discussed in Section 2, the CFD simulations were performed at 1:50 scale with ${\rho }_w=998.6$ kg/m³ and ${\mu }_w=8.89\cdot {10}^{-4}$ Pa∙s; however, all dimensional quantities in this paper have been scaled up to full scale based on Froude scaling, with ${\rho }_w=1025$ kg/m³.

To reduce computing time, only a half domain was used for the simulation. This choice is justified by the low KC number encountered. The numerical domain is shown in Figure 3. The boundary conditions used are summarized in

Table 2. Boundary conditions.

Boundary	Type	Velocity	Pressure	Phase Fraction
Floater	Wall	No slip	N/A	N/A
Upstream	Velocity inlet	Linear incident wave field	N/A	Incident wave
Downstream *	Pressure outlet	Extrapolated backflow dir.	Hydrostatic	Calm water
Top	Pressure outlet	Extrapolated backflow dir.	Hydrostatic	Air only
Bottom	Velocity inlet	Linear incident wave field	N/A	Water only
Sides	Symmetry planes	N/A	N/A	N/A

* A wave damping zone of length $2{\lambda }_1$ was placed next to the downstream boundary only.

.

On the upstream and bottom boundaries, the flow velocity was prescribed according to the linear wave theory [16]; however, because of the large water depth relative to the wavelengths, the bottom boundary was essentially a no-slip boundary. In fact, effectively identical results for the present problem down to the small nonlinear wave loads have also been obtained with no-slip condition on the bottom. To avoid wave reflection, a wave-damping zone [30] of length $2{\lambda }_1$ was placed next to the downstream boundary. The coefficients for wave damping were determined with the help of [31]. No wave-forcing/relaxation zone was employed upstream or to the side of the domain. This would, of course, result in diffracted waves being reflected from the upstream and side boundaries. This is a deliberate choice made to ensure that the setup of the numerical simulation is fully consistent with the experimental setup of the upcoming validation campaign. The numerical domain was initialized with calm water.

Formally second-order discretization schemes in both time and space were used. The details are summarized in

Table 3. Spatial and temporal discretization schemes.

Terms	Methods
Momentum advection	Second-order upwinding
Gradient	Second-order hybrid Gauss-LSQ with Venkatakrishnan’s limiter
Phase fraction advection	High-resolution interface-capturing scheme
Time integration	Second-order implicit scheme based on backward differentiation

. The convection of the phase fraction was done with the high-resolution interface-capturing scheme to avoid excessive smearing of the water–air interface. It should be noted that, while we have chosen second-order methods for the present application for their robustness and availability, higher-order methods for finite-volume CFD simulations have also been proposed in the literature (see, e.g., [32]) and successfully applied to practical engineering problems (see, e.g., [33]).

4. Baseline Computational Setup

The baseline computational setup described in this section aims to strike a balance between numerical accuracy and computing time. A systematic numerical uncertainty analysis was performed and is described in Section 5 for the baseline setup by repeating the simulation with different numbers of outer iterations, time steps, and grid sizes.

The overall dimensions of the computational domain are 1950 m long, 225 m wide (half domain), and 350 m tall (100 m between the top boundary and the calm water level). The origin (center of the floater) is exactly 442 m ($2{\lambda }_1$) from the upstream boundary. The depth, width, and distance to the upstream boundary from the floater are matched exactly between the CFD simulations and the upcoming experimental validation campaign.

The trimmer meshing tool of STAR-CCM+, which generates predominantly hexahedral cells, was used to mesh the numerical domain. A reference cell size, $h_{ref}$, can be used to concisely describe the mesh resolution. For the baseline grid, $h_{ref}=0.45$ m or $D_B/h_{ref}=53.3$, where $D_B$ is the diameter of the heave plates. We refined the computational mesh near the wetted portion of the floater as well as near and right below the free surface. The mesh refinement was specified using several mesh-refinement zones in the shape of rectangular boxes. The extent of each mesh refinement box along with the corresponding target cell sizes in the $x$, $y$, and $z$ directions—$\mathsf{\Delta }x$, $\mathsf{\Delta }y$, and $\mathsf{\Delta }z$—as multiples of $h_{ref}$ are listed in

Table 4. Mesh refinement zones for the baseline setup.

Box No.	x [m]	y [m]	z [m]	$\mathbf{\Delta }\mathit{x}/{\mathit{h}}_{\mathit{ref}}$	$\mathbf{\Delta }\mathit{y}/{\mathit{h}}_{\mathit{ref}}$	$\mathbf{\Delta }\mathit{z}/{\mathit{h}}_{\mathit{ref}}$
1	Full length *	[0, 100]	[−4, 4.5]	2	8	0.5
2	Full length *	[0, 100]	[−10, −4]	4	8	1
3	Full length *	[0, 100]	[−25, −10]	4	8	2
4	Full length *	[100, 225]	[−4, 4.5]	2	16	0.5
5	Full length *	[100, 225]	[−10, −4]	4	16	1
6	Full length *	[100, 225]	[−25, −10]	4	16	2
7	[−47.5, 35]	[0, 45]	[−4, 4.5]	2	2	0.5
8	[−45.0, 30]	[0, 40]	[−25, −4]	2	2	2

* Up to the downstream wave-damping zone.

. A portion of the computational domain near the floater is shown in Figure 4 to highlight each mesh refinement box.

Refinement boxes 1 through 6 improve the mesh resolution near and right below the free surface. They cover the full extent of the numerical domain in the $x$-direction, except for the downstream wave-damping zone, where the mesh was suitably coarsened. Boxes 1 through 3 cover the inner (in the $y$-direction) section of the domain containing the floater, whereas boxes 4 through 6 cover the outer flank of the domain with doubled $\mathsf{\Delta }y$ to reduce the number of cells. Near the free surface, we have $\mathsf{\Delta }x={\lambda }_2/194$ and $\mathsf{\Delta }y=\left(H_1+H_2\right)/32$. Both longitudinal and vertical mesh resolutions exceed what would typically be used in regular-wave simulations (see, e.g., [34]). Note that the cell dimensions do not change with $x$ in boxes 1 through 6 to avoid any spurious wave reflection caused by an abrupt change in mesh resolution. Refinement boxes 7 and 8 cover the space near the wetted portion of the floater. The values of $\mathsf{\Delta }x$ and $\mathsf{\Delta }z$ on the free surface near the floater (Box 7) are kept consistent with those in the rest of the domain (Box 1) to ensure uninterrupted wave propagation; however, $\mathsf{\Delta }y$ is reduced in Box 7 to be equal to $\mathsf{\Delta }x$ to better resolve the diffracted waves near the floater. Box 8 refines the mesh near the submerged portion of the floater with isotropic cell sizes to better resolve flow separation from the heave plates. The surface mesh of the floater consists of square patches approximately $h_{ref}$ wide. A prism mesh zone extruded from the floater surface helps better resolve the viscous boundary layer. The total thickness of the prism mesh is 0.4 m divided into 10 layers with a constant expansion ratio and a first layer thickness of 1 mm. Note that because we did not use any turbulence model or wall function, it is likely that the boundary-layer behavior was not fully captured. Because of the low KC number of the present problem, however, viscous effects primarily manifested as flow separation from the sharp corners of the heave plates, which is insensitive to the structure of the boundary layer; therefore, we still expect to obtain reasonable results with the present setup. Outside the refinement zones, the cell size gradually expands to 32$h_{ref}$ in all three directions in the far field. Overall, the baseline grid consists of 6 million cells for a half domain.

The baseline setup has a fixed time step of $\mathsf{\Delta }t=T_2/1030$ with 40 SIMPLE iterations per time step. This choice of time step, together with the value of $\mathsf{\Delta }x$ near the free surface, results in wave celerity-based Courant numbers of 0.21 and 0.19 for the first and second wave components, respectively.

5. Numerical Results

Using the numerical method and baseline setup described in Section 3 and Section 4, CFD simulations were performed both with and without the floater. Each simulation was run for $3.5T_R$ (666.4-second physical time at full scale), and a time window of length $2T_R$ at the very end of the simulation was used for data processing.

The steady-state amplitudes of the nonlinear wave-exciting surge force in the $x$-direction and pitch moment about the $y$-axis at the difference frequency are the primary quantities of interest. The difference-frequency heave force in the $z$-direction is less important, because heave resonance of this semisubmersible occurs near the wave-frequency range, and it is controlled by the much stronger direct linear wave excitation. In the ensuing analysis, the wave-load amplitudes at the two wave frequencies are also included for completeness.

5.1. Wave-Only Simulation

To better understand the characteristics of the bichromatic incident waves, we first carried out the CFD simulation with only the incident waves and without the floater. For this purpose, a 2D computational mesh with equivalent domain size and mesh resolution in the $x$-$z$ plane was used. Wave elevation, $\eta \left(t,x\right)$, was recorded at 71 longitudinal positions evenly spaced between $x=$ −440 m and $x=$ 1325 m at 25-m intervals.

As an example, the time series of $\eta$ at $x=0$ is shown in Figure 5. The initial transient phase can be observed during the first half of the first repeat period. During the last two repeat periods, $t\in \left[1.5T_R,3.5T_R\right)$, $\eta$ becomes approximately periodic; therefore, using a fast Fourier transform (FFT), we can extract the complex wave amplitudes, $A_1$ and $A_2$, at the two wave frequencies, $f_1$ and $f_2$, as well as the wave amplitude, $A_d$, at the difference frequency, $f_d$, from the time series of $\eta$ over the time window $t\in \left[1.5T_R,3.5T_R\right)$. This was done for each $x$-position where the wave elevation was recorded, and the results are shown in Figure 6. Overall, the wave-frequency amplitudes remain relatively constant during the entire length of the numerical domain—except in the downstream wave damping zone, where a rapid decrease can be observed. The amplitude of the higher frequency wave component, $\left|A_2\right|$, shows only slight wave dissipation downstream, which is likely caused by minor numerical diffusion.

The magnitude of the difference-frequency wave amplitude, $\left|A_d\right|$, is significantly lower than those at the wave frequencies, but it varies strongly with $x$. This is the same behavior identified in [35], where the modulation of $\left|A_d\right|$ with $x$ was attributed to the superposition of three different wave components at frequency $f_d$: a forward-propagating (incident) free-wave component with complex amplitude ${\zeta }_{if}$ in the $+x$ direction, a backward-propagating (reflected) free-wave component with amplitude ${\zeta }_{rf}$ in the $-x$ direction, and an incident bound-wave component in the $+x$ direction with amplitude ${\zeta }_{ib}$. Of the three wave components, only the bound waves can be considered an intrinsic part of the nonlinear bichromatic-wave system and can be directly obtained from second-order potential-flow theory. The two free-wave components in the solution were likely generated by imperfect boundary conditions and can be considered “spurious”, potentially leading to contamination of the difference-frequency wave loads on the floater.

Note here that the terms incident and reflected are simply used as shorthand to denote the direction of propagation. In fact, the wavelength of the free-wave components at $f_d$ obtained from the linear dispersion relation—4600 m at full scale—is longer than the numerical domain; the downstream wave-damping zone cannot absorb such a long wave, leading to back-and-forth reflection between the upstream and downstream boundaries. The bound waves, on the other hand, come from the nonlinear interaction between the two primary wave components. Because both primary wave components were effectively absorbed by the downstream damping zone, there are no measurable reflected bound waves. Unlike the free waves, the bound waves do not obey the linear dispersion relation; instead, the wavenumber of the bound waves is given by $k_b=k_2-k_1$, where $k_1$ and $k_2$ are the wavenumbers of the two primary wave components.

By exploiting their different wavelengths and directions, we can reliably separate the three wave components at $f_d$ using a wave-splitting procedure [35,36]. We can assume $A_{d,j}=A_d\left(x_j\right)$ to have the form:

$$A_{d,j}=A_d\left(x_j\right)={\zeta }_{if}{e}^{-i{k}_d{x}_j}+{\zeta }_{rf}{e}^{i{k}_d{x}_j}+{\zeta }_{ib}{e}^{-i{k}_b{x}_j},j=1,2,\dots ,n,$$

(7)

where $k_d$ is the wavenumber of the free waves obtained from the linear dispersion relation. If $A_{d,j}$ is known at $n\ge 3$ different $x$-positions, a complex linear system of equations can be formed with the help of Equation (7) to solve for the three unknown complex amplitudes—${\zeta }_{if}$, ${\zeta }_{rf}$, and ${\zeta }_{ib}$—assuming that the choices of $x_j$ do not result in a singular system. In the present study, $A_{d,j}$ at $n=28$ different $x$-positions evenly spaced between $x=-275$ m and $400$ m were used to solve for the unknown amplitudes in the least-squares sense. The resulting values are $\left|{\zeta }_{if}\right|=0.0146$ m, $\left|{\zeta }_{rf}\right|=0.0317$ m, and $\left|{\zeta }_{ib}\right|=0.0132$ m.

As a check, the values of $A_d\left(x\right)$ obtained directly from the CFD solution are compared to the reconstructed values obtained using Equation (7) with ${\zeta }_{if}$, ${\zeta }_{rf}$, and ${\zeta }_{ib}$ in Figure 7. The agreement between the original wave amplitudes and the reconstructed ones is good for both the real and imaginary parts, lending confidence to the wave-splitting procedure. Further, the bound-wave amplitude $\left|{\zeta }_{ib}\right|=0.0132$ m obtained from the wave-splitting procedure is in reasonable agreement with the prediction by second-order potential-flow theory, $0.5\left|A_1\right|\left|A_2\right|\left(k_2-k_1\right)=0.0125$ m [37], with only a 6% difference.

The amplitudes of the free waves obtained from this wave-splitting procedure can be used to estimate their contribution to the difference-frequency wave loads on the floater. This contribution/contamination can then be subtracted from the results as a correction.

5.2. Simulation with the FOWT Platform

Using the baseline setup, we carried out simulations with the FOWT platform present. The amplitude spectra of the wave-exciting forces/moments are obtained by performing an FFT on the corresponding time series during the time window $t\in \left[1.5T_R,3.5T_R\right)$. The resulting amplitude spectra for the surge force and pitch moment are shown in Figure 8. The difference-frequency wave excitation is two orders of magnitude smaller than the wave-frequency excitation, posing significant numerical challenges.

Further, the amplitudes of the wave excitation (the frequency components highlighted in red in Figure 8) fluctuate slightly depending on the time window used for the FFT analysis. This fluctuation obtained using a sliding-window analysis is illustrated in Figure 9, which shows the wave-excitation amplitudes at the difference frequency and the two wave frequencies as functions of the start of the time window, $t_s$, normalized by the repeat period, $T_R$. The FFT time window is $2T_R$ long, covering the interval $t\in \left[t_s,t_s+2T_R\right)$. Strong transient effects can be observed for $t_s<0.5T_R$. For $t_s>0.5T_R$, the fluctuation becomes negligible at the two wave frequencies, indicating that periodic steady state has been reached. The nonlinear wave excitation at the difference frequency still shows some fluctuation. This fluctuation can be quantified using the standard deviations, ${\sigma }_f$, of the wave-load amplitudes from $t_s=0.5T_R$ to $1.5T_R$. The grey bands in Figure 9 show the $\pm 2{\sigma }_f$ (95% confidence) intervals about the means. To minimize the initial transient effect, the wave-load amplitudes computed from the last time window, $t_s\in \left[1.5T_R,3.5T_R\right)$, are taken as the final results.

Section 5.1 discusses the presence and identification of the “spurious” difference-frequency free waves. The free waves, though small, have nonnegligible contribution to the wave loads at $f_d$, which can be considered a source of contamination. Fortunately, this contamination can be systematically subtracted from the CFD solution using the procedure developed in [35]:

$${\tilde{A}}_{f,j}\left(f_d\right)=A_{f,j}\left(f_d\right)-{\zeta }_{if}{X}_j\left(f_d,\theta =0\right)-{\zeta }_{rf}{X}_j\left(f_d,\theta =\pi \right),$$

(8)

where ${\tilde{A}}_{f,j}\left(f_d\right)$ and $A_{f,j}\left(f_d\right)$ are the corrected and uncorrected wave-excitation amplitudes in direction $j$ ($j=1,\dots ,6$) at frequency $f_d$. The unit-amplitude wave excitation in direction $j$ is given by $X_j$, which can be obtained from linear potential-flow theory, and $\theta$ is the angle between the direction of the wave propagation and the $x$-axis. The complex amplitudes of the incident and reflected free waves, ${\zeta }_{if}$ and ${\zeta }_{rf}$, can be obtained from the wave-only simulation using the procedure described in Section 5.1. Of course, Equation (8) addresses only the leading order contribution from the free waves; however, the amplitudes of the free waves are generally very small—any higher order contribution to the difference-frequency wave loads is likely negligible. The correction procedure given by Equation (8) was found to be important when performing the cross-verification of the CFD solutions in the OC6 project [35].

To facilitate comparison with potential-flow theory, the wave-load amplitudes obtained from the CFD simulations are nondimensionalized using the normalization factors listed in

Table 5. Normalization factors for the wave-load amplitudes at the frequencies of interest.

	$\mathbf{Diff}.\mathbf{Freq}.,{\mathit{f}}_{\mathit{d}}={\mathit{f}}_2-{\mathit{f}}_1$	$\mathbf{First}\mathbf{Wave}\mathbf{Freq}.,{\mathit{f}}_1$	$\mathbf{Sec}\mathbf{ond}\mathbf{Wave}\mathbf{Freq}.,{\mathit{f}}_2$
$F_x$	${\rho }_{w}g{A}_{wp}{A}_1^{*}{A}_2\left(k_2-k_1\right)$	$2{\rho }_{w}g{A}_{wp}{A}_1$	$2{\rho }_{w}g{A}_{wp}{A}_2$
$F_z$	$0.5{\rho }_{w}g{A}_{wp}{A}_1^{*}{A}_2\left(k_2-k_1\right)$	${\rho }_{w}g{A}_{wp}{A}_1$	${\rho }_{w}g{A}_{wp}{A}_2$
$M_y$	$0.5{\rho }_{w}gL{A}_{wp}{A}_1^{*}{A}_2\left(k_2-k_1\right)$	${\rho }_{w}gL{A}_{wp}{A}_1$	${\rho }_{w}gL{A}_{wp}{A}_2$

. Note that the force/moment amplitudes and the wave amplitudes are all complex, containing both magnitude and phase information. When normalizing the difference-frequency loads, $A_1^{*}$, the complex conjugate amplitude of the first low-frequency wave component is used. In Table 5, $A_{wp}$ and $L$ are the waterplane area and the center-to-center distance between columns, respectively. The normalization factors for the wave-frequency surge and heave forces are based on the linear wave excitation on a deep drafted vertical cylinder at the long-wave limit. The theoretical nonlinear bound wave amplitude, $0.5A_1^{*}{A}_2\left(k_2-k_1\right)$, is used to normalize the difference-frequency wave loads.

The normalized wave-excitation amplitudes computed using the baseline CFD setup are listed in

Table 6. Normalized amplitudes (magnitude) of wave excitation at the frequencies of interest.

	$\mathbf{Difference}\mathbf{Frequency},{\mathit{f}}_{\mathit{d}}={\mathit{f}}_2-{\mathit{f}}_1$			$\mathbf{First}\mathbf{Wave}\mathbf{Freq}.,{\mathit{f}}_1$	$\mathbf{Second}\mathbf{Wave}\mathbf{Freq}.,{\mathit{f}}_2$
	Before Correction	After Correction	% Change
$F_x$	1.65	1.57	−4.7%	0.544	0.559
$F_z$	6.51	4.63	−29%	0.447	0.533
$M_y$	1.20	1.18	−1.8%	0.389	0.450

. The difference-frequency amplitudes before and after the correction of Equation (8) are both included. The corrections are small for the amplitudes of $F_x$ and $M_y$ but significant for $F_z$. This is expected, because the wavelength of the difference-frequency free waves is much greater than the dimension of the floater.

6. Uncertainty Estimation

As discussed earlier, the primary goals of the ongoing OC6 CFD effort are the cross-verification of the CFD solutions and the validation of the CFD solutions against experimental measurements regarding the nonlinear difference-frequency wave excitation. Meaningful validation requires not only knowledge of experimental uncertainties, but also the uncertainties of the numerical CFD solutions. With the present uncertainty analysis, ideally, we would like to aim for a 95% confidence interval for our CFD results. The numerical uncertainty estimate is especially important for the present problem because of the low magnitudes of the difference-frequency wave excitation. Small changes in the numerical setup could potentially result in nontrivial changes in the nonlinear wave loads. The following analysis should provide a sense of the variability of the numerical solutions.

The three primary sources of uncertainty in the CFD solution are numerical uncertainty, modeling uncertainty, and statistical uncertainty. The three sources are separately discussed in this section. For consistency, the uncertainty quantification for the difference-frequency wave-load amplitudes are based on the values before the correction of Equation (8). All uncertainties in Section 6 are expressed as percentages of the uncorrected baseline solution in Table 6.

6.1. Numerical Uncertainty

The primary sources of numerical uncertainty in the CFD solution are the round-off errors, iterative errors, and discretization errors. With the use of double-precision calculations, the round-off errors are considered negligible in the present study, and we consider only iterative and discretization uncertainties.

6.1.1. Iterative Uncertainty

In the present study, the iterative error in the numerical solution comes from solving the nonlinear field equations in an iterative and segregated fashion with the SIMPLE algorithm for pressure-velocity coupling as well as from solving the linearized system of equations during each (outer) iteration iteratively. Generally, the iterative error and discretization error are not fully independent. To minimize contamination from the iterative error to the estimation of the discretization error, it is beneficial to first reduce the iterative error as much as possible. To estimate the iterative error, the numerical simulation was repeated using the baseline setup but with 10, 20, 30, and 40 (baseline) outer iterations per time step. As an example, the maximum residual of the $x$-momentum equation at the end of the outer iterations, as computed by STAR-CCM+, during the time window $t\in \left[1.5T_R,3.5T_R\right)$ is shown in Figure 10. Doubling the number of iterations per time step from 20 to 40 resulted in only a 30% reduction in final residual; therefore, significant (an order of magnitude) reduction in residual beyond that of 40 iterations is no longer practical, and 40 iterations per time step was chosen for the baseline setup.

Because it is impractical to reduce the residual to machine precision, an extrapolation procedure is needed to provide an estimate, ${\delta }_{it}$, of the iterative error ${\epsilon }_{it}={\varphi }_{it}-{\varphi }_{it,0}$ where ${\varphi }_{it}$ is a scalar quantity of interest. The hypothetical value of the same quantity obtained at the limit of no residual is given by ${\varphi }_{it,0}$. The iterative error can be assumed to have the form [25]:

$${\epsilon }_{it}\approx {\delta }_{it}={\varphi }_{it}-{\tilde{\varphi }}_{it,0}={\alpha }_{it}{e}^{H\left(r_{it}\right){\beta }_{it}},$$

(9)

where $r_{it}$ is the residual of the solution associated with ${\varphi }_{it}$, and ${\tilde{\varphi }}_{it,0}$ is the estimated value of ${\varphi }_{it,0}$. More specifically, in the current analysis, $r_{it}$ is taken as the maximum residual of the $x$-momentum equation at the end of the outer iterations during the time window for analysis, $t\in \left[1.5T_R,3.5T_R\right)$. As noted in [25], it is important to use the residual as the independent variable for the error estimation rather than the number of iterations, because the convergence might have already stalled.

The function $H\left(r_{it}\right)$ should go to $-\infty$ as $r_{it}\to 0$. Two alternatives have been proposed for $H\left(r_{it}\right)$ [25]: $H\left(r_{it}\right)=\ln (r_{it})$, which results in a power-law estimator for the iterative error, and $H\left(r_{it}\right)=-1/r_{it}^{p_{it}}$. If the former definition of $H\left(r_{it}\right)$ is used, three unknown constants—${\alpha }_{it}$, ${\beta }_{it}$, and ${\tilde{\varphi }}_{it,0}$—need to be determined; if the latter is used, one additional constant, $p_{it}$, needs to be determined. These constants can be solved for in the least-squares sense by minimizing the standard deviation, ${\sigma }_{it}$:

$${\sigma }_{it}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left[{\varphi }_{it}-\left({\tilde{\varphi }}_{it,0}+{\delta }_{it}\right)\right]}^2}{N-3}},$$

(10)

where $N$ is the number of simulations available for the estimation of the iterative error. In the present work, we have $N=4$. The latter form of $H\left(r_{it}\right)$ consistently produces superior fits of the available data set (lower ${\sigma }_{it}$); therefore, $H\left(r_{it}\right)=-1/r_{it}^{p_{it}}$ is used throughout for estimating ${\epsilon }_{it}$. The iterative uncertainty, $U_{it}$, can be obtained by multiplying the estimated error, ${\delta }_{it}$, to a safety factor of $F_S=1.25$ and augmented by ${\sigma }_{it}$:

$$U_{it}=F_S\left|{\delta }_{it}\right|+{\sigma }_{it}.$$

(11)

The estimation of the iterative errors and the corresponding uncertainties are illustrated in Figure 11. The resulting iterative uncertainties for all quantities in Table 6 are listed in

Table 7. Iterative uncertainties, $U_{it}$, in the magnitudes of the wave-load amplitudes.

	$\mathbf{Difference}\mathbf{Frequency},{\mathit{f}}_{\mathit{d}}$	$\mathbf{First}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_1$	$\mathbf{Second}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_2$
$F_x$	±0.5%	±2.2%	±0.9%
$F_z$	±1.4%	±0.7%	±2.1%
$M_y$	±1.2%	±2.0%	±2.0%

.

6.1.2. Discretization Uncertainty

The total discretization error consists of contributions from both temporal and spatial discretization [38]:

$$\varphi -{\varphi }_0=\epsilon ={\epsilon }_{\mathsf{\Delta }t}+{\epsilon }_{\mathsf{\Delta }x},$$

(12)

where $\varphi$ is a scalar quantity of interest obtained with a certain time step and reference grid size, and ${\varphi }_0$ is the value of the same quantity obtained at the limit of infinite temporal and spatial resolution. The combined discretization error, $\epsilon$, is the sum of the temporal discretization error, ${\epsilon }_{\mathsf{\Delta }t}$, and the spatial discretization error, ${\epsilon }_{\mathsf{\Delta }x}$.

The two sources of discretization error can be combined and estimated together by maintaining a suitable ratio between the time step and the grid size depending on the expected order of convergence during the convergence study [25]. This is usually a more efficient approach in terms of the total number of simulations required. For the present analysis, however, we opted to estimate ${\epsilon }_{\mathsf{\Delta }t}$ and ${\epsilon }_{\mathsf{\Delta }x}$ separately. This approach provides added flexibility to use different estimators for the two sources of error depending on the observed convergence behavior, while also allowing us to better understand their relative importance.

Both temporal and spatial discretization errors can be estimated using the same approach. For brevity, $\mathsf{\Delta }t$ and $\mathsf{\Delta }x$ are dropped from the subscript in the following discussion. Whenever possible, the errors were first estimated using Richardson extrapolation with the standard power-law error estimator [24]:

$${\varphi }_i-{\varphi }_0={\epsilon }_i\approx {\delta }_i={\varphi }_i-{\tilde{\varphi }}_0=\alpha h_i^p,$$

(13)

where ${\epsilon }_i$ is either the temporal or the spatial discretization error in the quantity ${\varphi }_i$ associated with $h_i$, which represents either the $i$th time step or the $i$th grid size. Here, ${\varphi }_0$ is the value of the quantity of interest obtained as $h\to 0$, and ${\tilde{\varphi }}_0$ is the estimate of ${\varphi }_0$. The estimated error is given by ${\delta }_i$. Similar to the iterative uncertainty, the constants ${\tilde{\varphi }}_0$, $\alpha$, and $p$ can be solved for in the least-squares sense by minimizing the standard deviation defined as follows [24]:

$$\sigma =\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left[{\varphi }_i-\left({\tilde{\varphi }}_0+{\delta }_i\right)\right]}^2}{N-3}},$$

(14)

where $N$ is the total number of grids or time-step sizes considered in the convergence study. In the current work, we have $N=4$ for estimating both temporal and spatial discretization errors.

Depending on the observed convergence characteristics, the estimator given in Equation (13) might not always be suitable. If the resulting apparent order of convergence, $p$, is significantly higher than the theoretical value of 2 ($p>2.05$), the resulting error estimate might be too unconservative; therefore, a modified estimator with $p=2$ was used instead [23]:

$${\delta }_i=\alpha h_i^2.$$

(15)

A direct application of Equation (13) could also result in $p<0$, which suggests a divergent solution as the time step or cell size approaches zero. This does not mean, however, that the solution is truly divergent; rather, it is simply a consequence of anomalous convergence behaviors that render the estimator of Equation (13) unsuitable. Further, we consider the least-squares fit produced with Equation (13) or Equation (15) to be poor and unreliable if the resulting standard deviation, $\sigma$, is larger than the averaged data range ${\overline{\mathsf{\Delta }}}_M$ [38]:

$$\sigma \ge {\overline{\mathsf{\Delta }}}_M={\mathsf{\Delta }}_M/\left(N-1\right),$$

(16)

where the data range, ${\mathsf{\Delta }}_M$, is given by:

$${\mathsf{\Delta }}_m=\underset{1\le i,j\le N}{\max }\left|{\varphi }_i-{\varphi }_j\right|.$$

(17)

If either $p<0$ or $\sigma \ge {\overline{\mathsf{\Delta }}}_M$, a third error estimator in the form of a quadratic function was tried [24]:

$${\delta }_i={\alpha }_1{h}_i+{\alpha }_2{h}_i^2.$$

(18)

Unlike Equations (13) and (15), which fit only monotonic convergence, the error estimator of Equation (18) also supports non-monotonic convergence. The constant coefficients, ${\alpha }_1$ and ${\alpha }_2$, are, again, obtained by minimizing $\sigma$. A fourth linear error estimator has also been suggested as an alternative to the power-law estimator of Equation (13) when $p<0.5$ [24]:

$${\delta }_i=\alpha h_i.$$

(19)

However, we did not encounter any situation where Equation (19) provided a better fit—i.e., smaller $\sigma$—than the other alternatives; therefore, Equation (19) was never used.

Finally, if the least-squares fit was still poor with Equation (18) as the error estimator, we resorted to the following range-based error estimate [39]:

$${\epsilon }_i\approx {\delta }_{{\mathsf{\Delta }}_M}=\frac{{\mathsf{\Delta }}_M}{h_N/h_1-1},$$

(20)

where $h_N$ is the coarsest time step (cell size), and $h_1$ is the finest. The introduction of the denominator, $h_N/h_1-1$, allows the range of time steps/cell sizes investigated during the convergence study to be directly incorporated into the range-based estimate. Although ${\mathsf{\Delta }}_M$ can be made artificially small by considering only very similar time steps or cell sizes, the resulting error estimate, ${\delta }_{{\mathsf{\Delta }}_M}$, will be severely penalized, with $h_N/h_1-1$ approaching zero.

The discretization uncertainty is generally defined as the estimated discretization error multiplied to a suitable safety factor, $F_s$, along with additional minor adjustments and limiters [40]. The choice of $F_s$ and the necessary adjustments can be tailored to the error estimator used and are generally based on experience.

If Equation (13) is used as the error estimator, the discretization uncertainty, $U_i$, for the $i$th cell size or time step, following the recommendation from [23], can be estimated with a safety factor of $F_S=1.25$ as:

$$U_i=F_S\left|{\delta }_i\right|+\mathsf{\sigma }.$$

(21)

In [23], an additional upper bound of $F_S{\mathsf{\Delta }}_M$ was imposed on the uncertainty, $U_i$, when the apparent order of convergence, $p$, was significantly lower than the theoretical value of 2 ($p<0.95$) to avoid being overconservative with the uncertainty estimate. This upper bound, however, is not implemented in the current study, which will result in a more conservative estimate. If Equation (15) is the error estimator, the uncertainty can be estimated with $F_S=1.25$ using the following equation [23]:

$$U_i=\max \left(F_S\left|{\delta }_i\right|+\mathsf{\sigma },F_S{\mathsf{\Delta }}_M\right).$$

(22)

Here, $F_S{\mathsf{\Delta }}_M$ is used as an additional lower bound to avoid the uncertainty estimate from being too low. If Equation (18) is the error estimator, the uncertainty is again estimated using Equation (22). Finally, if the range-based approach is adopted, the uncertainty is estimated as:

$$U_i=F_S{\delta }_{{\mathsf{\Delta }}_M}$$

(23)

with an increased safety factor of $F_S=3$ [39].

Uncertainty from Spatial Discretization

The spatial discretization uncertainty, $U_{\mathsf{\Delta }x}$, is estimated by repeating the simulation on four different grids with reference cell sizes given by $D_B/h_{ref}$=94.8, 71.1, 53.3 (baseline), and 40.0, respectively, for a constant refinement ratio of 3/4. This choice of grid refinement ratio was used in similar convergence studies [13] and is considered adequate for the present problem. Geometric similarity between the different grids was maintained as much as possible by changing only the referencing grid size while maintaining the ratios listed in Table 4; therefore, the number of cells increased rapidly from 2.7 million for the coarsest grid to 31 million for the finest grid.

The convergence trends of the wave-load amplitudes in surge, heave, and pitch at both the difference frequency and the wave frequencies are illustrated in Figure 12. For the difference-frequency amplitudes (Figure 12a,d,g), the error estimator of Equation (18) leads to a fairly good fit in all cases, with $\sigma <{\overline{\mathsf{\Delta }}}_M$. No suitable fit can be found for the wave-frequency surge force and the heave force at the first wave frequency (Figure 12b,c,e); therefore, the uncertainty estimates are based on the range-based approach of Equation (23). For all other quantities, Equation (13) provides a reasonably good fit.

Note that the present uncertainty estimates are much higher than those suggested by the traditional approach of comparing the results from the coarsest and finest grids. Meanwhile, the fact that we are not able to find a suitable fit for three out of the nine quantities of interest also points to the difficulty in applying more robust uncertainty quantification practices to practical engineering problems. The difficulty with grid convergence encountered in the present work could be partly caused by the lack of a suitable turbulence model. It is possible that improved convergence will be observed with a modified turbulence model better suited for free-surface flow (see, e.g., [41]) should they become more widely available.

The spatial discretization uncertainties associated with the baseline grid with $h_{ref}=D_B/53.3$ are presented in

Table 8. Spatial discretization uncertainties, $U_{\mathsf{\Delta }x}$, in the magnitudes of wave-load amplitudes.

	$\mathbf{Difference}\mathbf{Frequency},{\mathit{f}}_{\mathit{d}}$	$\mathbf{First}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_1$	$\mathbf{Second}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_2$
$F_x$	±18%	±0.9%	±1.1%
$F_z$	±31%	±6.0%	±8.6%
$M_y$	±9.8%	±3.8%	±2.9%

for all quantities listed in Table 6. The uncertainties in the magnitudes of difference-frequency wave-load amplitudes tend to be much higher than those at the wave frequencies. This is not surprising considering that the difference-frequency wave loads are typically two orders of magnitude smaller than the wave-frequency components. Further, the heave force at the wave frequencies tends to have higher uncertainties than the surge force and pitch moment. This is likely related to the higher viscous contribution in heave force, which tends to be more difficult to accurately capture in numerical simulations.

Uncertainty from Temporal Discretization

To quantify the uncertainty from temporal discretization, $U_{\mathsf{\Delta }t}$, the simulation was repeated with four different fixed time steps, which, from finest to coarsest, are given by $T_2/\mathsf{\Delta }t$ = 1831, 1373, 1030 (baseline), and 773. This choice of time steps is again based on a constant refinement ratio of 3/4.

The convergence of the wave-load amplitudes with time step is illustrated in Figure 13. For temporal convergence, no suitable error estimator can be found for the difference-frequency wave loads except in surge (Figure 13a); therefore, the uncertainty estimates for the heave force and pitch moment at $f_d$ are both based on Equation (23) (Figure 13d,g). We do observe regular monotonic convergence for most wave-frequency excitations, however, except for the surge-force amplitude at the second wave frequency (Figure 13c), for which the error estimator in Equation (18) provides a good fit. For all other quantities, a suitable fit can be obtained with either Equation (13) or (15).

The temporal discretization uncertainties for the baseline time step are tabulated in

Table 9. Temporal discretization uncertainties, $U_{\mathsf{\Delta }t}$, in the magnitudes of the wave-load amplitudes.

	$\mathbf{Difference}\mathbf{Frequency},{\mathit{f}}_{\mathit{d}}$	$\mathbf{First}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_1$	$\mathbf{Second}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_2$
$F_x$	±1.3%	±0.2%	±0.8%
$F_z$	±2.0%	±3.1%	±0.9%
$M_y$	±1.2%	±1.0%	±0.5%

for all quantities listed in Table 6. The temporal discretization uncertainties in the difference-frequency wave loads are one order of magnitude smaller than the corresponding spatial discretization uncertainties, indicating mesh resolution as the limiting factor to numerical accuracy.

6.2. Modeling Uncertainty

Typical sources of modeling uncertainty include turbulence models, truncated numerical domain, and mismatched boundary conditions.

In the present work, no turbulence model or wall function is used, removing the uncertainties associated with the choices of models, parameters, turbulence initial conditions, and boundary conditions. Of course, the lack of a turbulence model renders full numerical convergence more difficult, leading to higher convergence errors.

The depth and width of the numerical domain as well as the distance between the floater and the upstream boundary/wavemaker exactly match the setup of the upcoming experimental validation campaign, largely removing the modeling errors associated with a truncated numerical domain—at least for validation purposes.

One major difference from the experiment is the use of a half domain in the simulation; however, as discussed, we expect this to have minimal effect on the solution because of the low KC number. Another major difference between the experiment and the simulation is in the length of the domain and the downstream boundary. The numerical domain is slightly longer than the physical one to accommodate the use of a momentum-source wave damping zone. The experiment, on the other hand, will use a parabolic beach to minimize wave reflection. Without further information, it is difficult to gauge the effect of the mismatch in the downstream boundary condition; therefore, modeling uncertainties are neglected in the present analysis.

6.3. Statistical Uncertainty

The statistical uncertainty, $U_{stat}$, is primarily caused by the fluctuation of wave-load amplitudes, as illustrated in Figure 9. Because we simply used the solution from the last time window as the final solution (the values corresponding to $t_s=1.5T_R$ in Figure 9), the $\pm 2{\sigma }_f$ interval (the grey bands in Figure 9)—where ${\sigma }_f$ is the standard deviation of the fluctuation—provides a reasonable estimate (95% confidence level) of the statistical uncertainty. The statistical uncertainty for each quantity in Table 6 is listed in

Table 10. Statistical uncertainties, $U_{stat}$, in the magnitudes of the wave-load amplitudes.

	$\mathbf{Difference}\mathbf{Frequency},{\mathit{f}}_{\mathit{d}}$	$\mathbf{First}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_1$	$\mathbf{Second}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_2$
$F_x$	±4.0%	±0.3%	±0.3%
$F_z$	±37%	±0.6%	±0.3%
$M_y$	±2.6%	±0.3%	±0.3%

. Except for $F_z$, the statistical uncertainties are much smaller than the numerical uncertainties.

Note that had we adopted the mean values as the solution instead of those from the last time window, it would be possible to further refine the estimates of the statistical uncertainties [42]; however, this is deemed unnecessary because the numerical uncertainties dominate the statistical uncertainties for most quantities of interest in the present analysis.

6.4. Total Uncertainty

The total uncertainty for each quantity of interest is obtained by combining the numerical and statistical uncertainties. The numerical uncertainties are approximated as the sum of the iterative uncertainties, the temporal discretization uncertainties, and the spatial discretization uncertainties. The direct summation of the temporal and spatial discretization uncertainties should yield a conservative combined discretization uncertainty. Further, the iterative error is, in general, not independent from the discretization error; therefore, direct summation rather than the root sum of squares should be used when combining the iterative uncertainty and the discretization uncertainty to obtain a conservative total numerical uncertainty [25]. The numerical uncertainty, $U_{num}$, is, therefore, given by:

$$U_{num}=U_{it}+U_{\mathsf{\Delta }t}+U_{\mathsf{\Delta }x}.$$

(24)

On the other hand, it is reasonable to treat the statistical uncertainty and the numerical uncertainty as independent, allowing the two to be combined in the root-sum-of-squares fashion. Based on this discussion, the expression for the total uncertainty is given as follows:

$$U_{tot}=\sqrt{U_{num}^2+U_{stat}^2}.$$

(25)

The total uncertainties are given in

Table 11. Total uncertainties, $U_{tot}$, in the magnitudes of the wave-load amplitudes.

	$\mathbf{Difference}\mathbf{Frequency},{\mathit{f}}_{\mathit{d}}$	$\mathbf{First}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_1$	$\mathbf{Second}\mathbf{Wave}\mathbf{Frequency},{\mathit{f}}_2$
$F_x$	±20%	±3.3%	±2.8%
$F_z$	±51%	±9.7%	±12%
$M_y$	±13%	±6.7%	±5.3%

. As before, the uncertainties are expressed as percentages of the (uncorrected) baseline solution in Table 6. Overall, the uncertainties are dominated by those from spatial discretization. The only exception is the difference-frequency heave force for which the statistical uncertainty is comparable to the spatial-discretization uncertainty.

The total uncertainty in the amplitude of the difference-frequency heave force reaches 51%; fortunately, low-frequency heave force is of little engineering significance, because the natural heave frequency of FOWT semisubmersible platforms is typically near the wave-frequency range. The uncertainties in the difference-frequency surge-force and pitch-moment amplitudes—20% and 13%, respectively—though large relative to those of the wave-frequency loads, are still acceptable, especially considering the small magnitudes of these quantities.

7. Result Analysis and Discussion

To put the uncertainty estimates in context, the global wave loads on the entire floater as predicted by the present CFD simulation are shown along with the estimated uncertainties in Figure 14. The amplitudes (magnitude) of the wave-induced surge force, heave force, and pitch moment at the difference and the two wave frequencies are shown separately. The difference-frequency wave-load amplitudes have been corrected with Equation (8). The numerical uncertainties, $U_{num}$, and the total uncertainties, $U_{tot}$, are shown using two uncertainty bands. Except for the difference-frequency heave force, the two uncertainty bands are almost visually identical, indicating negligible contribution from the statistical uncertainties in the present solution.

Further, predictions from the potential-flow theory are also included in Figure 14 for reference (as represented by the hatched areas in the bar plot). The difference-frequency wave loads are based on the full QTFs from second-order potential-flow theory (see, e.g., [43]), whereas linear wave excitations are shown for the wave-frequency loads. The surge-force amplitudes at the two wave frequencies from the CFD simulation are very close to the predictions by linear potential-flow theory, indicating negligible contribution from viscous effects for the present bichromatic-wave case. This observation is consistent with the low KC number encountered. On the other hand, the CFD model predicted significantly higher heave forces at the two wave frequencies than the potential-flow theory, even considering the wider uncertainty bands. The difference suggests significant viscous excitation in heave across the large heave plates. A similar difference between the CFD and potential-flow predictions also exists, to a lesser extent, in the wave-frequency pitch moments.

More importantly, the present CFD simulation predicted much higher nonlinear wave excitation in surge and pitch at the difference frequency than the second-order potential-flow theory for the chosen bichromatic-wave case. The differences far exceed the numerical uncertainties estimated for our CFD solution and should be considered significant. The difference is most pronounced in the surge force where the CFD prediction is more than 12 times that of the second-order potential-flow theory. The large discrepancy observed here has potential consequences for the prediction of the floater surge response because the difference frequency is very close to the surge natural frequency [2]. Potential-flow prediction of the difference-frequency pitch-moment amplitude is also approximately 30% lower than the present CFD solution.

The underprediction, compared to the present CFD solution, by potential-flow theory—especially in surge—is also qualitatively consistent with our prior experience in modeling the behavior of the semisubmersible in irregular waves (see Figure 1).

8. Conclusions

As part of the OC6 project, the nonlinear difference-frequency wave excitation on a fixed semisubmersible FOWT platform was investigated using CFD simulations at the model scale with a set of unidirectional, bichromatic incident waves. The present numerical study represents a first step toward a better understanding of the underprediction of low-frequency, nonlinear wave excitation and motion responses in surge and pitch, a prevalent issue with this type of floater geometries among mid-fidelity, engineering-level tools commonly relied on by FOWT designers and researchers. The present paper focused on the initial development of a numerical model to investigate this problem and an estimate of the associated uncertainty in the CFD solutions.

The bichromatic waves generated in the CFD simulations were carefully studied. At the difference frequency, three small wave components were identified: incident (forward propagating) and reflected (backward propagating) linear free waves and incident nonlinear bound waves. The free waves were likely generated by imperfect boundary conditions and are considered a source of contamination. A wave-splitting procedure was adopted to identify and quantify the three wave components at the difference frequency, and linear potential-flow theory was used to correct the difference-frequency wave loads from the CFD simulations by subtracting the contribution from the free waves. Because the presence of difference-frequency free waves in the CFD solution cannot be easily avoided, this correction procedure is critical to obtaining uncontaminated nonlinear wave loads that can be validated against experiments in the future. In the upcoming experimental validation campaign planned for this purpose, the physically obtained waves will be similarly analyzed with the help of multiple wave probes, and the experimentally measured difference-frequency wave loads will be corrected following the same procedure to enable a consistent comparison with the CFD results.

A systematic uncertainty analysis has been carried out based on recommended practices in the literature for the numerically predicted wave loads. Sources of uncertainty considered include numerical uncertainty, which consists of iterative uncertainty and discretization uncertainty, and statistical uncertainty. The analysis shows that, with the present CFD solution, the relative total uncertainties in the nonlinear loads, though still acceptable, tend to be much higher than those of the wave-frequency loads. Further, spatial discretization error is dominant for the difference-frequency wave loads, with mesh resolution being the limiting factor to numerical accuracy in the adopted numerical setup. It should be noted that, while the adopted procedure for estimating numerical uncertainty is generally applicable to a wide range of CFD methods, the quantitative uncertainty estimates presented here are, strictly speaking, only applicable to the specific numerical setup and wave condition described in the paper.

The next phase of the OC6 CFD effort will focus on carrying out the cross-verification of the CFD results from the various OC6 participants for the same set of load cases and on validating the CFD solutions against model-scale measurements from the upcoming experimental validation campaign. For both tasks, uncertainty estimation for the CFD solutions will be important because it enables meaningful comparisons among numerical solutions and with experimental measurements. Once we establish that CFD can provide reasonably consistent and accurate predictions of the nonlinear difference-frequency wave loads, the CFD solutions will be leveraged to provide a more detailed description of the nature of the nonlinear wave excitation. Such an analysis will help guide the ongoing effort in improving the predictive capabilities of mid-fidelity, engineering-level tools regarding the nonlinear, low-frequency wave excitation and floater resonance response in realistic sea states.