Next Article in Journal
First-Principles Study of Amorphous Al2O3 ALD Coating in Li-S Battery Electrode Design
Next Article in Special Issue
Development of Combined Load Spectra for Offshore Structures Subjected to Wind, Wave, and Ice Loading
Previous Article in Journal
A Review of Modelling of the FCC Unit—Part II: The Regenerator
Previous Article in Special Issue
Modeling the TetraSpar Floating Offshore Wind Turbine Foundation as a Flexible Structure in OrcaFlex and OpenFAST
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

OC6 Phase Ia: CFD Simulations of the Free-Decay Motion of the DeepCwind Semisubmersible

by
Lu Wang
1,*,
Amy Robertson
1,*,
Jason Jonkman
1,
Jang Kim
2,
Zhi-Rong Shen
2,
Arjen Koop
3,
Adrià Borràs Nadal
4,
Wei Shi
5,
Xinmeng Zeng
5,
Edward Ransley
6,
Scott Brown
6,
Martyn Hann
6,
Pranav Chandramouli
7,
Axelle Viré
7,
Likhitha Ramesh Reddy
7,
Xiang Li
8,
Qing Xiao
8,
Beatriz Méndez López
9,
Guillén Campaña Alonso
9,
Sho Oh
10,
Hamid Sarlak
11,
Stefan Netzband
12,
Hyunchul Jang
13 and
Kai Yu
14
add Show full author list remove Hide full author list
1
National Renewable Energy Laboratory, Golden, CO 80401, USA
2
Front Energies, Houston, TX 77084, USA
3
Maritime Research Institute Netherlands, 6708 PM Wageningen, The Netherlands
4
IFP Energies Nouvelles, 92852 Rueil-Malmaison, France
5
State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China
6
School of Engineering, Computing and Mathematics, University of Plymouth, Plymouth PL4 8AA, UK
7
Faculty of Aerospace Engineering, Delft University of Technology, 2629 HS Delft, The Netherlands
8
Department of Naval Architecture, Ocean and Marine Engineering, University of Strathclyde, Glasgow G1 1XQ, UK
9
National Renewable Energy Centre (CENER), 31621 Sarriguren, Spain
10
Renewable Energy Department, Nippon Kaiji Kyokai, Tokyo 102-8567, Japan
11
Aero- and Fluid-Dynamics Section, Department of Wind Energy, Technical University of Denmark, 2800 Lyngby, Denmark
12
Fluid Dynamics and Ship Theory, Hamburg University of Technology, 21073 Hamburg, Germany
13
Technip Energies, Houston, TX 77079, USA
14
American Bureau of Shipping, Houston, TX 77389, USA
*
Authors to whom correspondence should be addressed.
Energies 2022, 15(1), 389; https://doi.org/10.3390/en15010389
Submission received: 19 November 2021 / Revised: 14 December 2021 / Accepted: 21 December 2021 / Published: 5 January 2022

Abstract

:
Currently, the design of floating offshore wind systems is primarily based on mid-fidelity models with empirical drag forces. The tuning of the model coefficients requires data from either experiments or high-fidelity simulations. As part of the OC6 (Offshore Code Comparison Collaboration, Continued, with Correlation, and unCertainty (OC6) is a project under the International Energy Agency Wind Task 30 framework) project, the present investigation explores the latter option. A verification and validation study of computational fluid dynamics (CFD) models of the DeepCwind semisubmersible undergoing free-decay motion is performed. Several institutions provided CFD results for validation against the OC6 experimental campaign. The objective is to evaluate whether the CFD setups of the participants can provide valid estimates of the hydrodynamic damping coefficients needed by mid-fidelity models. The linear and quadratic damping coefficients and the equivalent damping ratio are chosen as metrics for validation. Large numerical uncertainties are estimated for the linear and quadratic damping coefficients; however, the equivalent damping ratios are more consistently predicted with lower uncertainty. Some difference is observed between the experimental and CFD surge-decay motion, which is caused by mechanical damping not considered in the simulations that likely originated from the mooring setup, including a Coulomb-friction-type force. Overall, the simulations and the experiment show reasonable agreement, thus demonstrating the feasibility of using CFD simulations to tune mid-fidelity models.

1. Introduction

Offshore wind energy is still a largely untapped source of renewable energy. One major barrier to the further utilization of offshore wind energy is the additional cost of the substructure supporting the wind turbine. This is especially the case in deep-water turbine farms when multiple floating substructures are required. Further cost reduction is necessary to make offshore wind energy economically competitive, and one way to achieve lower costs is through continued design optimization of the platform that supports the floating offshore wind turbine (FOWT). Presently, optimization is still mainly performed with computationally efficient mid-fidelity engineering models, such as the OpenFAST tool developed by the National Renewable Energy Laboratory (NREL) [1]; however, such tools require a priori knowledge of platform hydrodynamic characteristics to tune the model coefficients. Traditionally, the required information is obtained through wave-basin experiments with physical models. While generally reliable and efficient in terms of test scope, the wave-basin experiments are not well-suited to accommodate the rapid design changes common to optimization processes. It is therefore advantageous also to be able to obtain the necessary hydrodynamic properties through alternative means, such as high-fidelity computational fluid dynamics (CFD) simulations. This is one focus of the Offshore Code Comparison Collaboration, Continued, with Correlation, and unCertainty (OC6) project.
Previously in OC6 Phase I, we focused on using CFD to obtain estimates of the wave diffraction loads on the OC6-DeepCwind semisubmersible in a fixed condition, especially the nonlinear, low-frequency loads that potentially excite the surge and pitch resonance motion of a semisubmersible FOWT [2,3,4]. Tuning the engineering models also requires knowledge of the motion-damping characteristics of the floater, which are strongly influenced by viscous effects (see [5]); therefore, we now focus on the validation of CFD simulations of the free-decay motion of the OC6-DeepCwind semisubmersible. The primary objective is to evaluate whether the CFD setups adopted by the various OC6 Phase I participants for free-decay simulations can provide reasonable estimates of the calm-water hydrodynamic damping coefficients for tuning mid-fidelity engineering models. To this end, the numerical results from the OC6 participants are compared with each other and against experimental measurements for validation.
Because of the engineering importance of hydrodynamic damping in the design of floating offshore structures, extensive research exists on this subject, including studies that are based on the same DeepCwind semisubmersible design used in the Offshore Code Comparison Collaboration Continuation (OC4 [6]), Offshore Code Comparison Collaboration, Continued with Correlation (OC5 [7,8]), and OC6 projects [9,10] (Note that across these projects, the geometry of the semisubmersible has not changed, but the mass and inertia properties were modified to address small changes in the design of the wind turbine it supported).
Tran and Kim [11,12] carried out fully coupled aero-hydrodynamic simulations of the OC5-DeepCwind semisubmersible with a wind turbine using CFD and a catenary mooring solver. The results were compared with the OC5 experimental data [7,8]. Both free-decay and regular-wave-excited motions were investigated. It was observed that the surge-decay motions from CFD simulations with both the shear stress transport (SST) k-ω turbulence model and the Spalart–Allmaras model were in excellent agreement with those obtained without a turbulence model. On the other hand, the standard k-ϵ model resulted in excessive damping when compared to the other numerical solutions. Tran and Kim [12] also obtained good predictions of the motion response amplitude operators (RAOs) of the OC5-DeepCwind platform in regular waves; however, large errors were found in the mooring-line tension when compared to the OC5 experiment [7,8], demonstrating the difficulty with modeling the catenary mooring system.
More recent investigations paid more attention to the estimation of numerical uncertainties in the CFD results. Burmester et al. [13] investigated three representative problems in ocean engineering with increasing complexity, including the surge decay of the OC5-DeepCwind semisubmersible, with a special focus on quantifying the uncertainties in the numerical solutions. The numerical results were validated against the experimental measurements from the OC5 campaign [7,8]. Three different methods for estimating the discretization uncertainties were adopted and compared in an extensive convergence study. The least-squares approach of Eça and Hoekstra [14] was found to produce more conservative (larger) discretization errors in most cases compared to the other approaches. We therefore adopted this same approach for the present investigation. It was also found that direct uncertainty estimations for the damping coefficients led to excessively large uncertainty ranges; therefore, the uncertainty estimation was instead performed for the maximum (positive) and minimum (negative) surge displacements of the platform over the second surge period and the oscillation range (i.e., the difference between the maximum and the minimum). Burmester et al. recommended that numerical uncertainties should be estimated for quantities requiring as little postprocessing as possible to minimize uncertainties introduced by the postprocessing itself [13].
Wang et al. [15] performed CFD simulations of the OC5-DeepCwind semisubmersible undergoing pitch decay and compared the numerical solution to the OC5 experimental campaign [7,8]. Very small discretization uncertainties on the order of 1% were estimated for the maxima in floater motion and motion period based on the grid convergence index; however, substantial differences in the resulting linear and quadratic damping coefficients between the simulation and experiment were observed. Some important factors that potentially contributed to this discrepancy are the complex catenary mooring used in the experiment, the drag force on the mooring lines, the influence of a power cable attached to the model, the aerodynamic damping from the wind turbine and tower mounted on top of the model, and the three-degrees-of-freedom (3DoF) motion assumed in the CFD simulation rather than the full six-degrees-of-freedom (6DoF) motion in the experiment. In a follow-up study, Wang et al. [16] carried out a formal uncertainty analysis for the numerically predicted linear and quadratic pitch damping coefficients. An interesting approach was adopted where the numerical discretization uncertainties were first estimated for the floater motion time series and were subsequently propagated to the estimated linear and quadratic damping coefficients. It was again observed that the incorporation of a dynamic mooring model for the catenary mooring lines significantly improves the CFD predictions, especially in terms of the pitch period. Overall, successful validation was reported for the pitch period and the linear pitch damping coefficient; however, the quadratic damping was underpredicted [16].
Li and Bachynski-Polić [17] investigated the low-frequency radiation characteristics of the OC6-DeepCwind floater [9,10] by simulating free-decay motions and forced oscillation in surge, heave, and pitch. The heave and pitch decay from the CFD simulations were generally consistent with the experimental measurements, but the surge decay showed substantial differences, with the OC6 experiment having much faster decay. The increased damping was attributed to additional mechanical damping in the experimental mooring setup. (In this article, further analysis of the mechanical damping in the experimental setup is presented with an attempt to quantify and separate out at least part of the effect.) In contrast, similar issues were not found with forced oscillation, and the surge-damping coefficients obtained from the experiment and the CFD simulation showed good agreement in this case [17]. It was also concluded that at the low natural frequencies of the platform, the wave-radiation damping is generally negligible compared to the viscous damping. Because of viscous effects, the heave added mass estimated from the CFD results was also found to be substantially higher than potential-flow predictions [17].
Other examples of CFD investigations on this topic include those of Burmester et al. that investigated the surge-decay motion of the DeepCwind FOWT platform [18,19] with a focus on the effects of the various computational settings. In addition to quantifying formally the numerical discretization uncertainties with independent grid and time step convergence studies [18], the effects of various other aspects of the numerical setup, including different numerical schemes, the inclusion of a free surface (in comparison to double-body simulations), motion coupling, scaling effects, the computational domain configuration, the inclusion of a wave-absorption zone, the choice of turbulence models, and the catenary-mooring-line models, were all investigated [19]. It was discovered that the inclusion of a wave-absorption zone improves the appearance of the flow field but has very limited effects on the hydrodynamic damping of the structure. On the other hand, the parameters of the dynamic mooring-line model, more specifically the line weight and the transverse drag coefficient, have a strong impact on the motion of the structure, and proper tunning is required to achieve good agreement with the experiment.
Wang et al. simulated the motion of the DeepCwind platform in regular waves using CFD [20]. Good agreement with the experiment in terms of the surge RAO was achieved, whereas the heave and pitch RAOs were underpredicted by the CFD simulations. The errors in the heave and pitch RAOs were again attributed to the lack of a satisfactory nonlinear mooring model (In [20], the mooring loads were modeled using a linearized stiffness matrix). A follow-up study by Wang et al. [21] also formally estimated the numerical discretization uncertainties in the motion RAOs with a successful application of the least-squares method of [14].
Bozonnet and Emery performed CFD simulations of the forced oscillation and free-decay motion of a vertical cylindrical column with a thin heave plate, a common component of semisubmersible FOWT platforms [22]. The intention was to derive the relevant drag and damping coefficients that can be used with mid-fidelity potential-flow-based engineering models. For this purpose, extensive CFD simulations were performed to investigate the effects of motion amplitude, frequency, and geometric dimensions on the effective drag coefficient of the heave plate. As with [17], it was also noted that the added mass from the CFD simulations can be higher than the potential-flow predictions due to viscous effects [22], which can impact the resulting motion period.
Zhang and Kim [23] carried out fully coupled aero-hydrodynamic CFD simulations of the DeepCwind semisubmersible with the NREL 5-MW baseline wind turbine and compared their CFD results against the experimental measurements from the OC5 project [7,8].
The present investigation continues the research into the low-frequency damping characteristics of the DeepCwind offshore wind semisubmersible with the following key features aimed at addressing some of the limitations of prior studies:
(1)
The CFD simulations are validated against the measurements from a new experimental campaign from OC6 Phase Ia [9,10] specifically designed to minimize uncertainties in the experiment and to focus on the hydrodynamic problem better. The new campaign used a simplified linear taut-spring mooring system instead of the catenary mooring system in the prior OC5 project [7,8], which was frequently identified as one of the major obstacles preventing a successful validation [15,19,20]. The linear-spring mooring setup was developed to provide approximately the correct natural periods of floater motion while greatly reducing the uncertainties and difficulties associated with the numerical modeling of the mooring lines. The wind turbine and tower in the OC5 experiment were also replaced with a rigid bar and a block mass with similar inertial properties to minimize the effects of wind loading and tower flexibility. Compared to the OC5 experiment, the new design of the OC6 experimental campaign minimizes potential physical differences between the experimental and the numerical CFD setups to facilitate validation.
(2)
All CFD simulations in the present investigation have a full 6DoF floater motion. Effort was made to replicate the motion of the floater in the experiment in all directions, including the ones not directly relevant to the estimation of damping coefficients.
(3)
The uncertainty analysis for the numerical solutions is directly based on the linear and quadratic damping coefficients and the equivalent linear damping ratios estimated from the free-decay motion of the structure. This is different from prior investigations, which typically examined the uncertainty of more basic quantities, such as the maxima in the displacement of the structure. While these basic quantities are less affected by postprocessing and thus are better-suited for uncertainty analysis, it is nevertheless important to attempt to obtain uncertainty estimates for the more generalized hydrodynamic properties of the system such as damping coefficients and damping ratios, because they are of the most practical value for mid-fidelity engineering models.
(4)
The present collaborative validation study includes numerical solutions provided by many different organizations participating in the OC6 CFD investigation. Different software and CFD setups were used to produce these solutions, enabling a cross-verification study to obtain a qualitative sense of the variability that can be expected from the CFD predictions.
The rest of this article is organized as follows: the physical setup of the problem, including the geometry of the structure and the mooring configuration, is described in Section 2. The CFD setup is summarized in Section 3 with a detailed description of the baseline numerical setup adopted by NREL, which was used as a reference by the other OC6 CFD participants. A preliminary comparison of the free-decay motion from the experiment and the numerical simulations is presented in Section 4. The procedures for estimating the damping coefficients and equivalent linear damping ratios (the key metrics for validation) from the floater motion are given in Section 5. Section 6 documents the uncertainty analysis for selected CFD results, and Section 7 compares the numerical predictions and experimental measurements for verification and validation. Finally, Section 8 summarizes conclusions drawn from the study.

2. Overview of the Physical Problem

The present investigation is based on the OC6 Phase Ia model-scale validation campaign carried out at the Concept Basin of the Maritime Research Institute Netherlands (MARIN) within the framework of the MaRINET2 project [24]. Measurements from this campaign were previously used to validate engineering-level tools in the OC6 Phase Ia project [9], and we are now comparing them to CFD simulations. The experiments were performed with a simplified version of the OC5-DeepCwind FOWT, called the OC6-DeepCwind FOWT, with the wind turbine replaced by a rigid tower and block mass, allowing the turbine/floater system to be treated as a single rigid body. The motion of the floater was recorded in all 6DoF with an optical tracking system. The model was placed in the basin about 40 m (model scale) from the wave maker (which was disabled for the free-decay cases) and in the center of the basin widthwise. For reference, the MARIN Concept Basin is 220 m long, 4 m wide, and 3.6 m deep.
Three different load cases (LC) were considered: calm-water free-decay motions in surge, heave, and pitch, labeled LC 4.2, LC 4.4, and LC 4.6, respectively, in the OC6 Phase Ia project. Both the experiment and CFD simulations were performed at 1:50 scale (the NREL simulations were equivalently performed at full scale but with increased viscosity to match the model-scale Reynolds number), but all geometric dimensions and results are presented at full scale based on Froude scaling for comparison and analysis.
The floater (at equilibrium) and the adopted coordinate system for the CFD simulations and analysis are illustrated in Figure 1. The primary components of the floater include three outer columns having a diameter of 12 m and a height of 14 m below the still water level arranged to form an equilateral triangle. The center-to-center distance between two outer columns is 50 m. Below each outer column, a heave plate having a diameter of 24 m and a height of 6 m is attached, resulting in a total draft of 20 m. A main column with a diameter of 6.5 m and a draft of 20 m is also present at the center of the triangle. The various columns and heave plates are connected to each other with pontoons and braces, collectively referred to as cross members. A complete description of the floater geometry can be found in [6]. The origin of the earth-fixed coordinate system coincides with the center of the main-column waterplane area at equilibrium. The upstream outer column is centered on the −x-axis.
Wind loading is not considered; however, the inertial properties of the wind turbine and tower are included and combined with those of the floater by treating the turbine/floater system as a single rigid body. Important full-scale dynamic properties of the combined turbine/floater system are summarized in Table 1. The model-scale water density is 998.6 kg/m3, but a full-scale density of 1025 kg/m3 is assumed when computing full-scale quantities. The dynamic viscosity of water at model scale, which determines the Reynolds number, is 8.89 × 10−4 Pa∙s.
The mooring system consists of three taut-spring mooring lines. In the experiment, each thin mooring line was redirected by a pulley at an equivalent anchor point under water to a mechanical spring above water [10]. This means the mooring lines can be effectively treated as linear springs between the fairleads and the equivalent anchor points in the CFD simulations. The positions of the fairleads with the structure at equilibrium and the equivalent anchor positions are provided in Table 2. Each mooring line has an unstretched length of 55.432 m (measured from the equivalent anchor point) and a spring constant of 48.9 kN/m full scale based on the OC6 model-scale validation campaign [10]. In the present numerical simulations, a small adjustment is consistently made to the spring constant to match the surge-decay period better (see Section 3.5 for details).
In the OC6 free-decay experiments, the floater was manually pushed to an initial offset position and subsequently released to oscillate freely. Because of practical limitations in the experimental setup, the initial floater displacement could not be precisely controlled, and the exact point in time that the floater was released could not be recorded. For the CFD simulations, we assume that the floater was released at the first peak/trough in the recorded motion time history, and we used the position and velocity of the floater at this time instant as the initial condition for the CFD simulations. In this fashion, the initial offsets and velocities of the floater in all 6DoF for the three free-decay load cases were estimated and are listed in Table 3 and Table 4. The translational offsets are defined by the displacement of the center of gravity. Note that in the experiment, nonnegligible initial displacements and/or velocities were sometimes present in directions other than the primary direction of interest. For example, a sizeable initial offset in the sway direction was recorded in the surge free-decay experiment (LC 4.2).

3. Numerical Setup

Several organizations participating in the present OC6 collaborative validation study carried out CFD simulations of some or all three free-decay load cases and provided numerical solutions for validation. The numerical setups, including mesh resolution, time step, numerical schemes, and turbulence models, differed among the CFD participants depending on the capabilities of the software used and the experience of the participants. For brevity, only the baseline numerical setup developed by NREL for the CFD software STAR-CCM+ is described in detail in Section 3.2, Section 3.3 and Section 3.4. This baseline setup was shared with the OC6 participants at the beginning of the investigation as a reference. Selected aspects of the numerical setups adopted by the OC6 participants are briefly summarized in Appendix A (see Table A1, Table A2, Table A3 and Table A4). Finally, tuning of the mooring spring constant for the CFD simulations is discussed in Section 3.5. The mooring spring constant was tuned to match the experimental surge period better and was uniformly adopted by all the participants with all load cases for consistency.
Participants in this study include the American Bureau of Shipping (ABS), National Renewable Energy Centre, ClassNK, Technical University of Denmark, Dalian University of Technology, IFP Energies nouvelles, MARIN, NREL, Delft University of Technology, the University of Plymouth, and the University of Strathclyde. All CFD results presented are based on the finite-volume method with the volume-of-fluid formulation as outlined in Section 3.1; however, several different software packages were used, including STAR-CCM+ [25], OpenFOAM [26], and ReFRESCO [27]. The complete floater geometry, including all the cross members, are included in all the CFD simulations. Additionally, the Hamburg University of Technology (TUHH) provided simulation results from a time-domain three-dimensional lower-order panel method, panMARE [28], which solves the potential flow field and free-surface elevation at each time step. In the TUHH model, the columns and the heave plates are discretized into panels, and the contributions from Morison drag are included. Empirical drag forces in the heave direction are also evaluated and applied to the faces of the heave plates. The cross members are modeled with the Morison equation only [29].
While the baseline numerical setup was shared with all CFD participants at the beginning of the project, a subgroup of participants comprising ABS, MARIN, and NREL underwent closer coordination with a frequent cross-comparison of results to ensure the physical problems were implemented in the numerical simulations as consistently as possible. Therefore, the results from this subgroup tend to be somewhat more consistent with each other compared to those of the rest of the participants, who performed the simulations more independently.

3.1. Mathematical Formulation

All CFD simulations in the present study are based on the finite-volume method and the volume-of-fluid model for multiphase flows [30]. Both water and air are treated as incompressible Newtonian fluids; therefore, the flow is described by the incompressible Reynolds-averaged Navier–Stokes equation and continuity equation:
· u = 0 ,
( ρ u ) t + · ( ρ u u ) = ( p + 2 3 ρ k ) + · [ ( μ + μ t ) ( u + ( u ) T ) ] + ρ g ,
where u is the flow velocity vector; p is pressure; and g is gravitational acceleration. The turbulence eddy viscosity is given by μ t , and the term involving the turbulent kinetic energy k might be neglected with certain turbulence models, such as the one-equation Spalart–Allmaras model. The local fluid density and dynamic viscosity, ρ and μ , are given by the volume-of-fluid model [30]:
ρ = ( 1 α w ) ρ a + α w ρ w ,
μ = ( 1 α w ) μ a + α w μ w ,
where the subscripts a and w denote air and water, respectively. The volume fraction of water, α w , is given by the following scalar transport equation:
α w t + · ( α w u ) = 0 .
Depending on the implementation and numerical setting, an additional interface compression term of the form · [ u r α w ( 1 α w ) ] might be included on the left-hand side of Equation (5) to help maintain a sharp water-air interface, with u r being an artificial velocity field aligned with the normal to the interface [31]. Note that this compression term is nonzero only when close to the interface with 0 < α w < 1 . This additional interface compression term was not used by all participants. For instance, the STAR-CCM+ simulations performed by NREL do not make use of artificial interface compression.
The baseline setup uses the Spalart–Allmaras detached-eddy simulation (SA-DES) [32] to compute the eddy viscosity with the improved delayed detached-eddy simulation (IDDES) formulation [33], all- y + wall treatment [25], and low-Re correction [34]. With the all- y + wall treatment of STAR-CCM+, the wall shear stress is computed as in laminar flows when the near wall mesh is fine enough to resolve the viscous sublayer with y + around unity or less and is estimated with boundary-layer modeling when the mesh is coarse with y + values of near wall cells in the log-law layer ( y + > 30 ). The all- y + wall treatment also produces reasonable results when y + values of near wall cells fall within the intermediate buffer layer with the help of blended wall functions [25]. With the baseline mesh described in Section 3.3, we intended to achieve a y + value of unity or less; however, higher y + values in the buffer-layer range were occasionally encountered, especially near the sharp corners of the floater. Therefore, the all- y + treatment should be used.

3.2. Numerical Domain, Initial Conditions, and Boundary Conditions

The baseline numerical domain is shown in Figure 2 with each boundary labeled. The exact locations of the boundaries are listed in Table 5 along with the corresponding boundary conditions, which are specified following the STAR-CCM+ setup. While the free-decay motion is expected to generate negligible radiated waves [17], wave-damping zones 50 m wide are nevertheless included next to the upstream and downstream boundaries. The depth and width of the numerical domain match the physical basin exactly; therefore, the side walls are treated as free-slip walls, and no wave-damping zone is employed next to the side boundaries. Compared to the cylindrical domains used in prior investigations that are intended to model open water [16], the rectangular domain used in the present study is better suited to model the narrow wave basin in which the experiment was performed. The flow field is initialized with calm water, and the floater is released from the estimated initial offset position (see Table 3) with the estimated initial velocity in the CFD simulations.

3.3. Computational Grid

The computational grid for the free-decay simulations was designed based on the expected characteristics of the flow. The Keulegan–Carpenter (KC) number for surge free decay (LC 4.2) can be estimated as
KC = 2 π A D   ,
where D is the diameter of the member, and A is the initial surge offset. For the offset columns with a diameter of 12 m, the KC number is 2.7 and half that for the larger-diameter (24 m) heave plates. The KC number based on the diameter of the central main column, 6.5 m, is higher at 4.9. The cross members should experience significant flow separation with a KC number of approximately 20; however, the contribution to the global loads is likely limited because of their small diameter of just 1.6 m. The relatively low KC numbers, especially for the larger offset columns and heave plates, suggest that viscous drag associated with the wall boundary layer can potentially have nonnegligible contributions to the total surge damping [35]; therefore, it is important to resolve the shear layer on the floater surface properly with an adequately fine prism-layer mesh. Flow separation from the corners of the heave plates also contributes to the surge damping by exerting a transverse drag force on the thick heave plates. This conjecture is supported by a preliminary round of grid sensitivity study, which shows that the mesh resolution near the corners has a strong influence on the surge-damping characteristics of the floater; therefore, fine mesh resolution near the sharp corners of the floater is also necessary. The fine mesh at the corners should also benefit the prediction of the viscous damping during heave and pitch decay, which is likely dominated by flow separation on the heave plates rather than any viscous effect on the upper columns.
Further, because of the long surge and pitch natural periods, the radiated waves are likely to be extremely small for surge and pitch free decay; we believe it is neither feasible nor necessary to resolve the radiated waves. The heave free-decay motion might generate slightly more radiated waves, but heave damping is dominated by the viscous drag on the large heave plates. We therefore anticipate wave radiation to play a very limited, if not negligible, role for the selected free-decay scenarios, and the mesh resolution near the free surface was designed simply to have a reasonably sharp water-air interface. The expectation of negligible radiation damping is also supported by [17].
Based on a preliminary grid sensitivity study, a baseline grid was developed. For convenience, a reference cell size of h = 6 m full scale is used to describe the grid resolution. In the far-field boundaries, a maximum isotropic cell size of 2   h is targeted. Near the free surface, three mesh-refinement zones spanning the full length and width of the domain are used to resolve the interface better. The refinement zones along with the targeted cell sizes in all three directions are listed in Table 6. The extents of the refinement zones are defined using the coordinate system of Figure 1. In each zone, the aspect ratio of the cells is maintained at Δ z / Δ x = 0.5 to minimize wave-dispersion error based on the recommendation of [36], even though the surface waves are expected to play a very limited role in the present investigation.
Several box-shaped mesh-refinement zones are employed to better resolve the flow near the floater. Table 7 lists the extent of each refinement zone, when the floater is at the equilibrium position (see Figure 1), and the targeted isotropic cell sizes. The target patch size of the floater surface mesh is h / 16 ; however, much smaller patch sizes are allowed where the geometry is complicated, such as at the joints between members.
Finally, to achieve even higher grid resolution near the corners of the heave plates and the bottom of the main column, three cylindrical or ring-shaped mesh refinement zones are defined. The extents of these refinement zones, expressed in terms of the ranges of the radial distance, R , measured from the column centerlines and ranges of the z coordinate, are listed in Table 8 along with the targeted isotropic cell sizes.
A prism-layer mesh with a total thickness of 0.48 m full scale divided into 15 layers with a first layer thickness of 2 mm is generated from the floater surface. This configuration resulted in a near wall y + below 1.5 in most places. Higher y + values were occasionally encountered locally at the sharp corners of the floater or at the seams where the various pontoons and braces meet.
The baseline mesh built with the trimmer meshing tool of STAR-CCM+, which generates mostly hexahedral elements, has 12.9 million cells. A y -plane cross section of this mesh is shown in Figure 3. The increased mesh resolution near the structure, especially near the corners of the heave plates and the bottom of the main column, can be seen.

3.4. Numerical Schemes and Settings

The baseline setup uses second-order discretization schemes in both time and space. More specifically, time integration is based on a second-order implicit backward-differencing scheme. The momentum advection terms are discretized in space using the hybrid upwind and bounded central-differencing scheme, and the second-order hybrid Gauss-LSQ method with Venkatakrishnan’s limiter is used for gradient reconstruction. The advection of the phase volume fraction is based on the High-Resolution Interface-Capturing scheme.
Based on a preliminary sensitivity study, a time step of T / Δ t = 400 was chosen where T is the motion period, i.e., the period of the surge, heave, or pitch free-decay motion. This baseline temporal resolution is already twice as fine as the time resolution of T / Δ t   200 used in [12], which was found to produce converged unsteady solutions for a similar problem. Pressure-velocity coupling is achieved with the semi-implicit method for pressure-linked equations (SIMPLE) with 20 iterations per time step. The under-relaxation factors for velocity and pressure are 0.8 and 0.4, respectively. The maximum number of iterations of the 6DoF rigid-body motion solver is also 20 per time step. The sensitivity of the results to the number of iterations is investigated in Section 6.1. To accommodate the motion of the floater, b-spline mesh morphing is used [25]. With an implicit algorithm for fluid-structure coupling, the mesh morphing is updated at each iteration [25]. Each morphing operation starts from the same initial mesh with the floater at the equilibrium position to avoid gradual deterioration of the mesh quality over time.

3.5. Model Parameter Tuning

Due to the uncertainties in the experimentally measured model parameters [10], it is sometimes necessary to tune the floater/tower and mooring dynamic properties used in the CFD simulations slightly. In the present investigation, the mass of the floater and the mooring spring constant were adjusted simultaneously to achieve the target design draft and the surge free-decay period measured in the experiment. In the end, a final combined floater/tower mass of 1.4046 × 107 kg (a 0.2% decrease) and a mooring spring constant of 52.32 kN/m were consistently used in all numerical simulations. The adjusted mooring stiffness represents a 7% increase from the experimentally measured value of 48.9 kN/m. This adjustment is partially justified by the fact that there is approximately a ±10% uncertainty in the experimental mooring spring constant [10]. All other dynamic properties of the floater/tower system and the mooring system used in the CFD simulations are consistent with those from the experiment [10] with no further adjustments.

4. Comparison of the Floater-Motion Time Series

The instantaneous position of the center of gravity in the earth-fixed coordinate system is used to define the translational motion of the structure in all three directions. The positive direction for each rotation is given by the right-hand convention.
As a preliminary step towards formal validation, the free-decay motion time series of the floater from the numerical simulations are compared to the experimental measurements (EXP) in Figure 4. The NREL CFD solutions shown are obtained with the baseline numerical configuration described in Section 3.
In Figure 4a, most CFD solutions (colored curves) show a qualitatively consistent surge decay motion; however, the experimentally measured surge motion decays visibly faster especially towards the end of the time series when the motion amplitude is small. The potential-flow solution of TUHH, on the other hand, significantly underpredicts the damping compared to the CFD solutions. Both issues are further explored in subsequent analyses.
With heave and pitch decay, the CFD solutions and the TUHH potential-flow solution agree with the experiment visually for the most part. Some CFD solutions slightly underpredict the heave and pitch periods, resulting in a small phase shift relative to the experiment after several periods toward the end of the time series. Despite the phase shift however, the periods of the decaying motion, as predicted by the CFD simulations, are all very close to the experiment. A detailed comparison of the periods is presented in Section 7.1, which shows that all motion periods from the numerical simulations are within ± 2 % of the experiment.
Overall, it appears that all numerical solutions capture the decaying motion reasonably well. In the rest of this article, further validation is carried out based on the damping coefficients and equivalent linear damping ratio estimated from the motion time series. These parameters are of practical engineering interest and are more sensitive to minor changes in the time history, leading to a more stringent validation.

5. Method of Analysis

The methods for computing the linear and quadratic damping coefficients and the equivalent linear damping ratio—key metrics for validation—from the motion time series are presented in this section.

5.1. Estimation of the Linear and Quadratic Damping Coefficients Using P Q Analysis

The linear and quadratic damping coefficients, B 1 and B 2 , can be estimated from the motion time series using the P Q method [37], which is briefly summarized here. For the weakly damped free-decay motions, we define the amplitude decrease over a half-cycle as Δ A i = A i A i + 1 , where A i and A i + 1 are, respectively, the positive amplitudes of the i th peak (trough) at time t = t i and the immediate next trough (peak) at t i + 1 in the surge-, heave-, or pitch-motion time series. The mean motion amplitude over the i th half-cycle is approximated as A ¯ i = ( A i + A i + 1 ) / 2 . The energy loss over the i th half-cycle, L i , (neglecting other nonlinear damping) is given by
1 2 k ( A i 2 A i + 1 2 ) = L i = t i t i + 1 ( B 1 x ˙ ( t ) 2 + B 2 x ˙ 2 | x ˙ | ) d t ,
where k is the total stiffness of the system for the mode of motion of interest, and x ˙ is the instantaneous translational or angular velocity of the floater. With the following approximation over the i th half-cycle:
x ˙ ( t ) ± A ¯ i ω sin ( ω ( t t i ) ) ,
where ω is the angular frequency of the motion, the following relation can be derived from Equation (7):
Δ A i A ¯ i = Δ A ¯ i = P + Q A ¯ i ,
which states that the normalized amplitude decrease Δ A ¯ is a linear function of the mean amplitude A ¯ . The y -intercept, P , and slope, Q , which can be determined from linear regression, are related to B 1 and B 2 by the following relations:
B 1 = 2 k π ω P   and   B 2 = 3 k 4 ω 2 Q .
Because P and Q are directly proportional to B 1 and B 2 , we use them as surrogates of the actual linear and quadratic damping coefficients in the subsequent analysis and sometimes refer to them as such for brevity. When performing the above analysis, the motion of the floater center of mass was used; however, very similar results were also obtained if the motion of a body-fixed origin at the center of the calm-water plane was used. In Section 5.1.1, Section 5.1.2 and Section 5.1.3, example results from the P Q analysis are shown. The relation in Equation (9) generally describes the surge-decay and heave-decay motions from the simulations well; however, the presence of surge-pitch coupling, which is not considered in the P Q method [37], results in poor linear regression for the pitch-decay motion (see Section 5.1.3 for details). Nevertheless, we believe the P Q analysis can still provide a meaningful characterization of the pitch-decay time series that enables comparisons across different results.

5.1.1. Surge Free Decay

Examples of the application of Equation (9) to the surge free-decay motion of the platform are given in Figure 5a. The CFD simulations performed by ABS, MARIN, and NREL are shown. Note that the NREL simulation is based on the baseline numerical setup of Section 3. The first half-period of surge oscillation was consistently excluded from the P Q analysis to minimize any start-up effect. The experimental surge free-decay motion shows poor linear regression when the standard P Q method of Equation (9) is applied. As illustrated in Figure 5b, the data points from the measurement show increased Δ A ¯ when the mean amplitude A ¯ is small, which cannot be described by the linear-plus-quadratic damping model. This behavior suggests the presence of a Coulomb-friction-type damping force of constant magnitude B 0 opposite the surge velocity in the experimental setup. To obtain an estimate of this force, the standard P Q method can be modified by adding an additional term in the form of B 0 x ˙ 2 / | x ˙ | to the integrand on the right-hand side of Equation (7). With this addition, Equation (9) becomes
Δ A ¯ i = O A ¯ i + P + Q A ¯ i ,
where O = 2 B 0 / k . The additional term O / A ¯ leads to increasing Δ A ¯ as A ¯ 0 as is observed with the experimental surge-decay motion. By multiplying Equation (11) on both sides by A ¯ i , a quadratic fit with Δ A i as a function of A ¯ i ,
Δ A i = O + P A ¯ i + Q A ¯ i 2 ,
can be carried out to determine the constants O , P , and Q . The resulting best fit in the form of Equation (11) shows good agreement with the experimental results as shown in Figure 5b. The magnitude of the friction-like force based on this analysis is 2.2 kN at full scale or only 0.017 N at model scale. This additional force, which likely originated from the pulley system placed under water, acts as a small external mechanical damping force when the floater is moving in the surge direction.
To verify the modified P Q analysis with Equation (11) further and investigate the effect of B 0 , MARIN repeated the surge-decay CFD simulation by applying an additional force of B 0 = 2.3 kN full scale to the floater center of gravity in the ± x -direction opposing the floater surge velocity. The results from this additional simulation are labeled MARIN-B0 to distinguish from the MARIN solution without B 0 . Note that for consistency, all numerical results in this article other than MARIN-B0 were performed without B 0 to focus just on the hydrodynamic damping. A comparison of the MARIN CFD solutions with and without B 0 and the experimental measurements are shown in Figure 6.
In Figure 6a, the MARIN CFD solution with added B 0 decays much faster, especially toward the end of the simulation when the motion amplitude is small, showing improved agreement with the experiment. The regression analyses for estimating the surge-damping coefficients are shown in Figure 6b. With the added B 0 , the MARIN-B0 solution is also showing a higher normalized decrease in surge amplitude as the mean surge amplitude approaches zero; the best fits for the MARIN-B0 solution and the experiment show good agreement in this region. More importantly, the modified regression analysis with Equation (12), when applied to the MARIN-B0 solution, recovers the value of B 0 = 2.3 kN prescribed in the CFD simulation, thus demonstrating the validity of the modified P Q analysis. The MARIN solution without B 0 is well-described by the linear regression with Equation (9). Interestingly, the estimated linear and quadratic damping coefficients of the MARIN-B0 solution are slightly different from those without B 0 ; however, this difference might not be physical and is likely just a consequence of the slightly different behaviors exhibited by the linear fit with Equation (9) and the quadratic fit of Equation (12). See Section 7.2 for further comparison and discussion. Finally, note that the experimental setup might also have additional mechanical linear damping; however, unlike the Coulomb-friction-type damping of B 0 , it is not possible to separate out the linear mechanical damping from the hydrodynamic contribution based on the motion time history alone.

5.1.2. Heave Free Decay

With heave free decay, no clear indication of the presence of a Coulomb-friction-like damping force is observed with the experimental measurements; therefore, Equation (9) is consistently used for the experiment and all numerical simulations. Heave damping is strongly influenced by the vortex shedding from the corners of the heave plates, which leads to strong flow-memory effects. As a result, the first few periods of heave oscillation generally do not follow the linear trend given by Equation (9). To avoid this issue and to minimize the effects of the mismatch in the initial conditions between the experiment and the numerical simulations, the first 3.5 periods of the heave oscillation were consistently discarded from the computation of the damping coefficients. In general, reasonable linear fits are observed with the numerical results of heave free-decay motion. Examples are shown in Figure 7a. The experiment, on the other hand, shows more scatter (see Figure 7b). It is not possible to determine exactly what is causing the heave-decay motion from the experiment to show slightly more scatter about the fitted linear relation compared to the CFD results. The normalized decrease in motion amplitude is highly sensitive to subtle changes in the motion time series, and any minor perturbations, such as minor nonideal behaviors displayed by the pulley-spring mooring system in the experiment, can cause the increased scatter. Nevertheless, the linear fit in Figure 7b still describes the trend of the experimental data reasonably well, and the resulting estimates of the linear and quadratic damping coefficients should be physically meaningful.

5.1.3. Pitch Free Decay

With pitch free decay, we typically observe poor linear regression with Equation (9). This is primarily caused by the coupling from the low-frequency surge motion during the pitch-decay test. As shown in Table 3, the pitch-decay load case also has a relatively large initial offset in surge, and the presence of surge motion, through motion coupling, leads to an irregular pitch-decay time history that is not well described by Equation (9), which is derived based on a single-degree-of-freedom motion as explained in Section 5.1. The poor linear regression is simply a consequence of the limitation of the postprocessing technique, the P Q analysis, rather than any issue with the experimental measurements or the numerical solutions. To avoid this problem in future investigations, care must be taken to minimize unwanted surge motion during pitch-decay experiments.
Nevertheless, because effort was made to replicate, in the CFD simulations, the surge motion present in the pitch free-decay experiment, the scatter of the data points in Figure 8 caused by the coupling from surge is qualitatively consistent between the experiment and the CFD solutions. In this article, we continue to use Equation (9) to characterize the pitch decay despite the poor linear regression in the hope that the effect of surge coupling cancels out over multiple periods. As with surge free decay, the first half-period of pitch oscillation was discarded when performing the P Q analysis for pitch free decay.

5.2. Equivalent Linear Damping Ratio

By matching the energy loss from the full linear and quadratic damping model with an equivalent linear damping model with damping coefficient B e , the linear equivalent damping ratio ζ = ω B e / ( 2 k ) can also be computed from P and Q using the following expression:
π ζ = P + F A · Q ,
where the motion-amplitude scaling factor, F A , is given by
F A = i = m n A ¯ i 3 i = m n A ¯ i 2   ,
with m and n being the numbers of the first and last peaks/troughs, respectively, used in the P Q damping analysis. Note that we are only interested in the hydrodynamic damping; therefore, the Coulomb-friction-like forces estimated from the experimental surge decay and that prescribed in the MARIN-B0 simulation are not included when computing ζ .
Alternatively, ζ can be estimated from the logarithmic-decrement method. To be consistent with the P Q method, a separate damping ratio should be computed for each half-cycle using the logarithmic-decrement method followed by an A ¯ 2 -weighted averaging. In the absence of the Coulomb-friction-like force, this approach was found to produce very similar estimates of ζ compared to those from Equation (13), even with pitch decay where the linear regression of the P Q method works poorly, lending more confidence to the estimated damping ratios.
The linear and quadratic damping characterized by P and Q obtained using the P Q method as well as the equivalent linear damping ratio ζ are used to evaluate numerical convergence in Section 6 and for validation against the experiment in Section 7.

6. Estimation of Numerical Uncertainty

The numerical error in the CFD solutions primarily consists of the discretization error associated with finite temporal and spatial resolutions and the iterative convergence error coming from a segregated and iterative solution of the relevant equations at each time step. With double-precision computation, the truncation error can be neglected [38]. To estimate the resulting numerical uncertainty, convergence studies were performed, but were limited to only the surge-decay and heave-decay motions that show good linear regression with the P Q analysis. The pitch-damping coefficients estimated using the P Q method likely contain substantial uncertainties from the linear regression itself; therefore, the estimation of numerical uncertainty is omitted for pitch decay. The modeling uncertainty in the solution associated with the choice of turbulence models was also briefly investigated. The uncertainty analysis was performed for the NREL CFD solutions only.

6.1. Iterative Uncertainty

First, the iterative convergence of the surge free-decay simulation was checked by repeating the simulation with 5, 10, 20 (baseline), and 40 SIMPLE iterations per time step using the baseline grid and time step described in Section 3. For consistency, the maximum number of iterations of the 6DoF motion solver was also set to the same number. To illustrate the influence, the final unnormalized root-mean-square residual of the x -momentum equation at the end of each time step, as computed by NREL using STAR-CCM+, is shown in Figure 9. A significant reduction in the residual is achieved when increasing the number of iterations from 5 to 10, whereas doubling the number of iterations from 20 to 40 does not result in a meaningful further reduction in the residual.
The time series of the surge and heave free-decay motions obtained with different numbers of iterations shown in Figure 10 corroborate the behavior of the residual shown in Figure 9. In Figure 10a, the surge-decay motion obtained with only 5 iterations shows large differences with the other results, whereas the solutions obtained with 10, 20 (baseline), and 40 iterations per time step are almost identical to each other. The heave-decay simulation is only repeated twice: once with the baseline setup of 20 iterations per time step and once with 40 iterations per step. The two solutions are visually identical to each other.
In terms of the damping coefficients and damping ratio, summarized in Table 9, the simulation with only 5 iterations per step shows poor iterative convergence, while 10, 20, and 40 iterations all provide similar results. The linear damping coefficient is the most sensitive to the number of iterations. As it is no longer feasible to reduce the residuals meaningfully beyond those of 40 iterations, we simply estimate the iterative uncertainties in the damping coefficients and the equivalent linear damping ratio as the difference relative to the results obtained with 40 iterations. This is the same procedure adopted in [39]. The iterative uncertainties in the heave-damping coefficients and damping ratio are similarly estimated and provided in Table 10. Based on the observed changes, the baseline of 20 iterations per time step appears to be a cost-effective choice, leading to an uncertainty of 0.4% in ζ due to iterative convergence. Assuming the pitch-decay simulation shows a similar iterative convergence as surge and heave, we expect 20 iterations per step to provide reasonable predictions for pitch damping as well.

6.2. Discretization Errors and Uncertainties

To estimate the discretization error, the surge- and heave-decay simulations were repeated four times with different mesh resolutions and time steps. As both temporal and spatial discretization schemes are formally second order, the time step and cell size were simultaneously adjusted following the same refinement ratio during the convergence study. Further, to maintain the geometric similarity of the mesh as much as possible, the mesh was refined and coarsened globally by simply adjusting the reference size h . The details are listed in Table 11, in which the baseline case is considered to have a refinement ratio of 1. Based on the iterative convergence analysis, 20 SIMPLE iterations per time step were used throughout the grid and time-step convergence study.
In Figure 11, an informal comparison of the motion time series obtained from the CFD simulations with the four different levels of numerical refinement described in Table 11 is shown. The time series are visually very consistent with each other despite the wide range of refinement levels considered, indicating a degree of numerical convergence.
In the rest of this section, a more formal estimation of the discretization uncertainties is carried out for the damping coefficients and damping ratios using a least-squares approach, which involves fitting a suitable convergence trend to the available data from the convergence study and extrapolating the solution to an infinite numerical resolution. The difference between the extrapolated and the actual numerical solutions provides a measure of the discretization error, which can be used to construct a discretization uncertainty interval.
For any given scaler quantity of interest, ϕ , we define ϕ i as the numerical solution of ϕ from the simulation with the i th grid size and time step. The hypothetical numerical solution of ϕ with infinite numerical resolution, i.e., with λ 0 , is given by ϕ 0 . The discretization error associated with ϕ i is defined as ϵ i = ϕ i ϕ 0 . With the simultaneous and consistent refinement of cell size and time step, the convergence error ϵ i can be assumed to follow the standard power-law error estimator (Richardson extrapolation) [38]:
ϵ i = α λ i p ,
where λ i is the refinement ratio associated with the i th time step and cell size (see Table 11). The constant coefficient α , the apparent order of convergence p , and the value of ϕ 0 can be determined if ϕ i from, at a minimum, N = 3 simulations with different cell size/time step are available; however, it is more reliable to have N > 3 simulations and solve for the three unknowns in the least-squares sense. More specifically, the values of α , p , and ϕ 0 are estimated by minimizing the standard deviation [38]:
σ = i = 1 N [ ϕ i ( ϕ 0 + α λ i p ) ] 2 N 3     .
The resulting minimized value of σ also provides a gauge of how well the convergence follows the assumed trend of Equation (15). If σ is too large, the error estimate should be considered invalid, and an alternative approach should be adopted. On the other hand, if σ is small and the apparent order of convergence is in the range 0.95 < p < 2.05 (assuming formally second-order schemes), a corresponding discretization uncertainty can be developed based on the estimated discretization error, ϵ i , with the help of a suitable safety factor, F S [40]:
U i = F S · | ϵ i | + σ .
Finally, if the apparent order of convergence from the least-squares fit is greater than the theoretical order of 2, say, p > 2.05 , the error estimates might be under-conservative; therefore, we simply set the value of p to 2 in this case when performing the least-squares fit, and an additional lower bound of F S · Δ M can be imposed on the discretization uncertainty where Δ M is the range spanned by ϕ i across all cell sizes and time steps investigated [40]. If σ is small, the widely adopted value of F S = 1.25 from the literature can be used (see [38,40]).
While, in principle, the above technique for estimating discretization uncertainty can be applied to any scaler quantity, it is generally better to apply it to quantities that require as little postprocessing as possible, such as the force/moment on the structure, or the instantaneous displacement of the floater. In this way, the errors are mostly coming from the numerical discretization rather than postprocessing itself [13]. Quantities such as P , Q , and ζ are therefore not preferred candidates for the estimation of discretization uncertainty. From the point of view of practical application however, the primary goal of performing the CFD simulations is to obtain estimates of the generalized damping characteristics that can be used in the mid-fidelity engineering models; therefore, it is more important to obtain uncertainty estimates for P , Q , and ζ rather than quantities that are specific to a particular realization of free-decay motion. We therefore attempt to estimate the discretization uncertainties directly for the damping coefficients and damping ratio in the present article, while fully acknowledging that they are not best suited for this type of analysis. In addition to the use of heavily postprocessed quantities, several other factors also complicate the estimation of the discretization errors. The unstructured mesh generated with the trimmer meshing tool does not guarantee full geometric similarity between grids with refinement, violating the assumption of the Richardson extrapolation. The mesh could also influence the behavior of the SA-DES turbulence model; therefore, the change in the solution is not due to the discretization error alone. Despite these problems, we still believe the least-squares method is the most robust way to obtain estimates of the discretization uncertainty, and we applied it in the present study whenever possible.
The values of P , Q , and ζ for surge free decay obtained from the CFD simulations with the four different refinement ratios in Table 11 are shown in Figure 12. The fitted convergence trends based on Equation (15) are also included, along with the estimated discretization uncertainties for each CFD simulation. For all three damping parameters, the best-fitting apparent order of convergence, p , is greater than 2.05; therefore, the least-squares fits shown are based on p = 2 instead, and the lower bound of F S · Δ M is imposed on the estimated discretization uncertainties as in [40]. In general, the fitted convergence trends describe the actual CFD solutions reasonably well, even with the value of p constrained to 2, lending confidence to the estimated discretization uncertainties.
The estimated percentage discretization uncertainties in P , Q , and ζ of the four repeated simulations in Table 11 are listed in Table 12. Large uncertainties are found for P , which is proportional to the linear damping coefficient. This is consistent with Figure 12a, which shows a substantial decrease in P as the simulation is refined. Even with the finest simulation, the discretization uncertainty in P can be as high as ± 48 % . On the other hand, the value of Q , which is proportional to the quadratic damping coefficient, increases with the refinement (Figure 12b). The uncertainty in Q is more moderate at ± 19 % with the finest simulation. Interestingly, the effect of the decreasing P is approximately balanced by the increase in Q with refinement, resulting in a more consistent equivalent linear damping ratio, ζ , which only has a discretization uncertainty of ± 10 % for both the baseline and fine simulations.
The uncertainties in Table 12 are the combined spatial and temporal discretization uncertainties. To obtain a measure of the relative importance of the temporal and spatial discretization errors, the surge-decay simulation was repeated one more time with the baseline mesh but a smaller time step of T / 600 , leading to a temporal resolution even finer than that of the fine case in Table 11. The resulting values of P , Q , and ζ are 0.058, 0.0270 m−1, and 4.2%, respectively. Compared to the baseline solution in Table 12, the changes are much smaller than those between the fine and the baseline cases, which means that the limited spatial resolution is the primary contributor to the discretization error, and the temporal discretization error is secondary.
The discretization uncertainties in the heave-damping coefficients were estimated following a similar procedure. The uncertainties along with the best-fitting error estimators are shown in Figure 13. Note that only the convergence of P and Q are shown, for which the error estimator of Equation (15) describes the convergence trends well. For heave decay, ζ does not show monotonic convergence, and Equation (15) fails to provide a valid fit; therefore, the alternative range-based estimate of
U = F S · Δ M ,
suggested by [40], was used to provide a discretization uncertainty for ζ . Because of the reduced confidence associated with the range-based approach, an increased factor of safety of F S = 3 was used [40].
The damping coefficients and equivalent linear damping ratio along with the estimated discretization uncertainties for heave decay are given in Table 13. The value of P is generally very close to zero; therefore, the absolute uncertainties are shown instead of percentage uncertainties. A slightly negative linear damping, P , is obtained with all but the finest CFD simulation. This is likely just a consequence of convergence error because the uncertainty ranges of P extend above zero in all cases. In contrast, the quadratic damping Q is always positive. Similar to surge decay, P and Q again show opposite trends as the simulation is refined. The contribution from the increasing P with finer numerical resolution is mostly canceled by that from the decreasing Q , resulting in an almost constant equivalent linear damping ratio, which only has a 7% discretization uncertainty. It appears that it is quite challenging to obtain fully converged linear and quadratic damping coefficients, but the resulting equivalent linear damping ratio is generally very consistent and not as sensitive to the numerical setup. It can therefore be argued that the equivalent linear damping ratios predicted by CFD simulations can be relied on with a higher degree of confidence compared to the separate linear and quadratic damping; however, the drawback is that ζ does not provide information on the exact composition of linear and quadratic damping.
The total numerical uncertainty can be obtained by combining the discretization uncertainty and iterative uncertainty. In the present analysis, however, the estimated discretization uncertainties in Table 12 and Table 13 are generally an order of magnitude greater than the estimated iterative uncertainty in Table 9 and Table 10 for the baseline 20 iterations per time step. We therefore neglect the contributions from iterative convergence and retain only the discretization uncertainty.

6.3. Uncertainties from the Choice of Turbulence Models

In addition to the numerical iterative and discretization errors, the CFD solutions also contain modeling errors. While there are many potential sources of modeling errors, we only consider here those associated with the choice of, or a lack of, turbulence models. The impact of the turbulence model was gauged by repeating the surge and heave free-decay simulations with the baseline numerical setup but with different turbulence models or no turbulence model activated. In addition to the SA-DES model [32] of the baseline configuration, the commonly used Menter SST k - ω model [41] was also tried. It is important to test the different turbulence models with both surge and heave free-decay motions because the nature of the damping in surge and heave can be very different. Surge damping contains a significant linear contribution coming from the viscous boundary layer, while heave damping is dominated by the quadratic damping from the drag on the heave plates (see Section 7). As a result, surge and heave decay could show different levels of sensitivity to the turbulence model.
The results from this sensitivity study are presented in Table 14. The quadratic damping Q in surge shows the most sensitivity to the choice of turbulence model, showing an 8% increase with no turbulence model compared to the baseline solution with the SA-DES model. In other cases, the variations in P and Q caused by the turbulence model are at least an order of magnitude smaller than the discretization uncertainty. With ζ , the impact of the choice of turbulence model is of the same order of magnitude as the corresponding discretization uncertainty but still less than half in value. The relatively weak influence from the turbulence models on the damping coefficients observed here is consistent with [12].
Based on the above observations, we consider the uncertainty associated with the turbulence model to be secondary to the discretization uncertainty; therefore, in the ensuing comparison and discussion of the numerical solutions and the experimental results, we only retain the discretization uncertainties for the NREL solutions, while acknowledging that the actual uncertainties in the CFD solutions could be slightly higher because of other contributions, such as those from the iterative error and the turbulence model investigated here.

7. Cross-Verification and Validation of the Numerical Results

In this section, a comparison of the damping coefficients P and Q , equivalent linear damping ratio ζ , and the periods of the free-decay motions T from all numerical simulations supplied by the OC6 participants and the experiment is carried out. Apart from the TUHH solutions, which are based on a potential-flow model with empirical drag forces, all other numerical solutions are based on finite-volume CFD simulations.
To compare the relative importance of the linear and quadratic damping directly, characterized by the values of P and Q , respectively, Q is nondimensionalized by the motion-amplitude scaling factor F A from Equation (14). For example, if P F A · Q , the linear damping and quadratic damping contribute roughly equal energy dissipation. The exact value of F A , of course, depends on the motion and varies slightly among the numerical simulations and the experiment. For the sake of a consistent comparison of Q however, we simply prescribe an approximated value of F A , denoted as F ˜ A , for each free-decay load case that will be consistently applied across all solutions to normalize Q . More specifically, we have F ˜ A = 2.8 m for surge decay (LC 4.2), 0.48 m for heave decay (LC 4.4), and finally 0.053 radian for pitch decay (LC 4.6). These prescribed values are representative of the actual scaling factor F A observed with each CFD simulation or experimental measurement, which falls within ± 15 % of F ˜ A . Note that the equivalent linear damping ratio ζ shown in this section was computed using Equation (13) with the exact F A of each simulation or experiment.
For validation against the experiment, the discretization uncertainties of the NREL solutions estimated in Section 6.2 are included whenever available; however, we do not have robust uncertainty estimates for the damping coefficients determined from the experiment. The validation performed here should therefore be considered an informal one.

7.1. Motion Periods

The motion periods are first compared in Figure 14. All numerical solutions show mostly consistent periods in the surge, heave, and pitch directions. With the increased mooring spring constant, the surge periods from the numerical simulations also match that of the experiment well. The system stiffness in heave and pitch, on the other hand, primarily comes from the hydrostatic restoring force/moment with the mooring playing a more limited role. The heave and pitch periods are also consistent between the numerical solutions and the experiment, suggesting that the hydrodynamic added mass/moment of inertia was successfully captured by the numerical simulations. Finally, the four repetitions of the surge- and heave-decay simulations performed by NREL for the discretization convergence study all show effectively identical motion periods. This means the periods are not sensitive to the numerical setup and therefore only provide a weak check on numerical convergence. In other words, other physical quantities, such as the damping coefficients, could be far from convergence even if the motion period has already converged. The two MARIN CFD results also show that the inclusion of B 0 in the simulation has little influence on the surge period. Overall, the agreement in the motion periods is deemed satisfactory with most numerical predictions falling within ±2% of the experiment.

7.2. Surge Damping

With the periods of the free-decay motion validated, the damping coefficients and the equivalent linear damping ratios are compared next. As already demonstrated in Section 6, the damping coefficients are very sensitive to the adopted numerical setup; thus, a comparison of the damping coefficients should provide a stronger check on the quality of the numerical solutions. The damping coefficients are also the primary quantities of interest that we would like to extract from the numerical simulations.
The damping coefficients and equivalent linear damping ratio in surge are compared in Figure 15. We first compare the MARIN and MARIN-B0 solutions, which show that the inclusion of the Coulomb-friction-type damping, B 0 , can lead to a slight decrease in the estimated linear damping (from P = 0.066 to 0.045 ) and an increase in the estimated quadratic damping (from F ˜ A · Q = 0.069 to 0.088 ). As discussed earlier, however, these relatively small changes might not be physical and are likely just a consequence of the slightly different behaviors exhibited by the linear fit with Equation (9) used for the MARIN solution and the quadratic fit with Equation (12) used for the MARIN-B0 solution. To demonstrate this, we repeated the analysis of the MARIN solution (without B 0 ) using the quadratic fit of Equation (12) instead of the linear fit with Equation (9) with the term associated with B 0 removed, in other words, with O forced to zero. This analysis results in a reduced estimate of the linear damping of P = 0.055 and an increased quadratic damping of F ˜ A · Q = 0.081 , thus at least partially explaining the differences observed between the MARIN and MARIN-B0 solutions shown in Figure 15. Interestingly, the two MARIN simulations, despite having different linear and quadratic damping coefficients, yield effectively the same equivalent linear damping ratio, thus demonstrating a level of consistency between the analysis with Equation (9) and that with Equation (12). Again, note that we did not include the contribution from B 0 when computing ζ as shown in Equation (13) to only focus on the hydrodynamic damping.
Next, we compare all the numerical solutions in Figure 15 to the experiment for validation. The NREL solutions show that the linear damping characterized by P tends to decrease as the simulation is refined in space and time, with sizeable discretization uncertainty even for the finest CFD solution. Substantial variation is also observed among the other CFD solutions for the linear surge damping, suggesting a high degree of sensitivity to the numerical setup. The linear damping from the experiment is mostly consistent with the NREL CFD solutions, falling within the NREL discretization uncertainty ranges; however, this agreement should be interpreted with care because the experimental linear damping could potentially be influenced by two factors. First, the use of the quadratic fit with Equation (12) to accommodate the Coulomb-friction-type damping in the experiment could potentially reduce the estimated linear damping as demonstrated by the two MARIN CFD simulations. In other words, the linear damping obtained with an idealized experimental setup without B 0 using Equation (9) could be higher. Second, the experimental setup might also exert additional mechanical linear damping on the structure, which cannot be distinguished from the linear hydrodynamic damping based on the available information. This means the actual hydrodynamic damping can also be lower than what the experiment suggests.
The CFD simulations tend to provide more consistent estimates of the quadratic damping characterized by Q . Furthermore, the Q value from the NREL simulations gradually increases as the mesh and time step is refined; however, it is unclear whether a fully converged CFD solution would agree with the experiment because the Q value from the experiment falls just outside the discretization uncertainty ranges. Most of the CFD solutions, except those from ABS and the University of Plymouth, also underpredict the quadratic damping when compared to the experiment. One possible explanation of this discrepancy is again provided by the two MARIN simulations, which show that the use of the quadratic fit of Equation (12) can lead to a slightly higher estimate of the quadratic damping. It is therefore possible that the actual hydrodynamic quadratic damping estimated with an idealized experimental setup without B 0 using Equation (9) is lower than what is shown, reducing the gap between the experimental results and CFD predictions. Overall, P and F ˜ A · Q are of comparable magnitude, meaning both the linear and quadratic damping are important for the present surge free-decay load case.
Finally, the majority of the CFD solutions provide consistent estimates of the equivalent linear damping ratio ζ between 4% and 5%, with the NREL, ABS, MARIN, and the University of Plymouth solutions showing very good agreement with each other. This is approximately half of the 10% surge damping ratio derived from B 1 = 0.2 k ( m + μ ) , where k , m , and μ are the surge mooring stiffness, floater mass, and surge added mass, respectively, suggested by Bureau Veritas for semisubmersibles [42], which is not surprising because surge damping can be very sensitive to the floater geometry and the motion amplitude. The low damping ratio encountered here is likely a consequence of the circular cross section of the columns and heave plates and the relatively low KC number. The NREL solution converges away from the experiment, with the experimental ζ falling outside the NREL discretization uncertainty ranges, suggesting the presence of some mismatch between the actual physical setup and the numerical setup, such as additional linear mechanical damping in the experiment.
Whereas it would certainly be of interest to link the differences in the CFD solutions of the various groups to the differences in the adopted numerical setups documented in the appendix directly, it is generally not possible to do so with any degree of confidence because the CFD configuration adopted by each group differs from the baseline setup and from each other in multiple respects, and the effects of the different settings on the results can either reinforce or cancel each other, further obscuring the influence of each individual aspect of the numerical setup. Out of all the participants other than NREL, the MARIN setup is the closest to the baseline setup described in Section 3, and some limited comparisons can be made with the NREL baseline solution. Apart from the different software packages and numerical discretization schemes (Table A1), the MARIN and NREL CFD setups are mostly consistent with each other with the same formal order of convergence (Table A1), the same SIMPLE pressure-velocity coupling scheme (Table A2), and a mostly similar computational mesh and identical time step (Table A3). The primary differences are in the mesh resolution near the corners of the heave plates and the bottom of the main column ( h / 64 in the NREL baseline simulation and h / 16 in the MARIN simulation as shown in Table A3), the level of iterative convergence (the MARIN simulation likely has better iterative convergence with up to 50 iterations per time step compared to the fixed 20 iterations per step used in the NREL baseline setup as shown in Table A2), and the different turbulence models (Table A4). As a result, it is unsurprising that the MARIN and NREL baseline solutions are mostly consistent with each other as shown in Figure 15 with almost perfect agreement in ζ . Based on the iterative convergence study and the turbulence-model sensitivity study performed in Section 6.1 and Section 6.3, the effects of iterative convergence and the choice of turbulence model are likely not enough to explain the remaining differences between the NREL and MARIN solutions in the linear and quadratic damping coefficients; therefore, the differences are most likely caused by the different mesh resolutions near the corners of the heave plates and at the bottom of the main column with possible contributions from the different numerical schemes used by the different software packages as well.
The potential-flow result of TUHH mostly matches the quadratic damping from the experiment and the CFD simulations by tuning the empirical transverse drag coefficient, C D , on the structure ( C D = 0.48 on the offset columns, 0.44 on the main column, 1.25 on the heave plates, and 0.5 on the cross members); however, it shows negligible linear wave damping, which is consistent with our expectation that wave radiation plays a negligible role in the present surge free-decay load case. It also suggests that the linear damping observed with the experiment and the CFD solutions is also caused by fluid viscosity. The linear viscous damping associated with the boundary layer on the surface of the structure typically becomes significant and comparable to the quadratic drag force for circular cylinders when the KC number is small (estimated to be 2.7 based on the diameter of the offset upper columns) and flow separation is weak; therefore, the original form of the Morison equation [29] without a linear drag term is fundamentally incorrect [35]. In this flow regime, potential-flow models with empirical drag forces need to either include an additional global linear damping matrix or use an alternative generalized form of the Morison equation with a linear drag term [35].

7.3. Heave Damping

The heave-damping coefficients and damping ratio are compared in Figure 16. Unlike surge damping, heave damping is dominated by the quadratic damping from the drag force on the heave plates. This behavior is captured by most of the numerical solutions. The linear damping from wave radiation is again secondary, confirming our expectation outlined at the beginning of Section 3.3.
The P and Q values from the experiment both fall within the discretization uncertainty ranges of the NREL CFD solutions, which gradually approach the experimental results with more refinement; in fact, the extrapolated values of P = 0.013 and F ˜ A · Q = 0.055 as the refinement ratio λ 0 are both in good agreement with the experiment. This observation suggests that the experiment and the CFD simulations are mostly consistent, and unlike surge decay, no clear signs of a significant mismatch between the experimental and CFD setups can be identified. In fact, we have several CFD solutions showing very good agreement with the experiment in terms of the quadratic damping. Out of all the CFD solutions, ABS and IFP Energies nouvelles are the closest to the experiment, matching both the linear and quadratic damping coefficients well. The potential-flow solution from TUHH, which was obtained with an empirical axial/normal drag coefficient of 4.8 for the faces of the heave plates, is also mostly consistent with the experiment and the CFD predictions. The same drag coefficient is also used by TUHH in the pitch-decay simulation as well.
The equivalent linear damping ratios from the numerical simulations are mostly consistent between 2% and 2.5%; quite a few predictions fall within the discretization uncertainty ranges of the NREL solutions, whereas the experimental ζ , unfortunately, falls just outside of the uncertainty ranges. It is possible that the less reliable range-based estimate of Equation (18) underpredicted the discretization uncertainty, even with the increased factor of safety of F S = 3 . Other sources of uncertainty in the numerical solution and the experimental uncertainties, which are not readily available, could also help explain the disagreement. Nevertheless, the agreement among the numerical simulations and with the experiment is deemed acceptable, suggesting that CFD simulations can, indeed, be used to obtain meaningful estimates of the damping characteristics of the floating platform, especially in terms of the equivalent linear damping ratio.
Interestingly, the dominating quadratic damping coefficients in heave from the MARIN simulation and the NREL baseline simulation are highly consistent, more so than for surge decay, despite the vastly different mesh resolutions near the corners of the heave plates and the bottom of the main column. The equivalent linear damping ratios are similarly in excellent agreement. This observation suggests that the heave-damping coefficients might be less sensitive to the mesh resolution in this region compared to the surge-damping coefficients with the OC6-DeepCwind floater.

7.4. Pitch Damping

Finally, the damping coefficients and damping ratio in pitch are shown in Figure 17. Because of the coupling from the surge motion, pitch decay generally shows poor linear regression when the P Q analysis is applied (see Figure 8). We therefore relegate the comparison to a qualitative one, without formal convergence and uncertainty analysis. Interestingly, despite the poor linear regression, the numerical solutions generally show good agreement with the experiment with several CFD simulations matching all three damping parameters compared in Figure 17. As with heave damping, the pitch damping is also predominantly quadratic, coming from the drag force on the heave plates, confirming our initial expectation. The potential-flow model of TUHH also shows good agreement with the experiment with the help of the empirical quadratic drag forces on the heave plates.

8. Conclusions

As part of the OC6 Phase I project, a collaborative verification and validation study for CFD simulations of the free-decay motion of the DeepCwind offshore wind semisubmersible was carried out. Three load cases focusing on the surge, heave, and pitch free decay were investigated. Several organizations provided CFD results for some or all three load cases to be validated against the corresponding experimental data obtained as part of the OC6 Phase Ia experimental campaign. The motion periods, linear and quadratic damping coefficients, and equivalent linear damping ratios were estimated from the experimental and numerical floater motion time series using the P Q analysis and adopted as metrics for validation. The present investigation has two unique features:
  • A constant Coulomb-friction-type mechanical damping force was identified from the experimentally measured surge motion using a new modified P Q analysis that also includes a constant damping force. The validity of the procedure was confirmed numerically by repeating the surge-decay CFD simulation with an added constant damping force, which was then successfully recovered from the numerical motion time series using the modified P Q method.
  • For selected CFD solutions, the discretization uncertainties were directly evaluated for the linear and quadratic damping coefficients and the equivalent linear damping ratio, which are the primary metrics that we would like to obtain from the CFD simulations to support the engineering modeling effort. This contrasts with previous efforts, which have typically focused on uncertainty estimates for more fundamental physical quantities, rather than the heavily postprocessed quantities such as damping coefficients.
Overall, the present validation study indicates that the CFD simulations performed by the OC6 participants can, in fact, produce meaningful estimates of the hydrodynamic damping characteristics of a floating offshore wind semisubmersible, especially in terms of the motion periods and equivalent linear damping ratios, that are mostly consistent with each other and in reasonable agreement with the experiment; however, it is important to pay attention to numerical convergence. The key conclusions derived from the convergence study and uncertainty analysis for the selected CFD solutions are as follows:
  • In the present analysis, even with the finest mesh of 25.6 million cells and the smallest time step of 533 steps per period, the estimated discretization uncertainty remains substantial, especially for the separate linear and quadratic damping coefficients.
  • The grid resolution is the leading contributor to the discretization error, whereas the temporal discretization error is secondary, at least for surge decay.
  • The numerical iterative errors and modeling uncertainties associated with the choice of turbulence model were also estimated but were found to be secondary to the discretization uncertainty for the load cases investigated and the numerical setup adopted.
  • Interestingly, the linear and quadratic damping tend to show opposite trends as the simulation is refined, resulting in a very consistent equivalent linear damping ratio. This trend can also be observed when comparing across the CFD solutions of the OC6 participants, with the equivalent linear damping ratio being the most consistent across simulations.
  • As a result, the equivalent linear damping ratio is the preferred metric against which mid-fidelity engineering models based on potential-flow theory and/or the Morison equation should be tuned; however, it does not provide information on the exact composition of the linear and quadratic damping.
Based on the numerical and experimental results presented in this study, we made the following key observations regarding the behavior and characteristics of the hydrodynamic damping of the DeepCwind semisubmersible:
  • With all three load cases investigated, wave radiation plays a negligible role in the overall motion damping because of the low natural frequencies.
  • The damping in the heave and pitch directions is predominantly quadratic, most likely coming from the drag force on the heave plates.
  • On the other hand, the surge damping can have comparable contributions from linear and quadratic damping when the KC number is small, with the former primarily coming from the viscous boundary layer.
The above observations on the characteristics of the hydrodynamic damping have several important implications for the development and tuning of mid-fidelity engineering models based on potential-flow theory and empirical drag forces, such as the panMARE model of TUHH discussed in this paper:
  • Because wave-radiation damping is negligible compared to the viscous damping, the accuracy of the mid-fidelity models in predicting free-decay motion is almost exclusively determined by the tunning of the drag coefficients.
  • The predominantly quadratic heave and pitch damping can be modeled satisfactorily with a quadratic empirical drag force on the heave plates, as demonstrated by the TUHH solution in the present investigation.
  • For surge decay, however, the conventional form of the Morison equation without a linear drag term is fundamentally incorrect especially for low KC numbers. Either additional linear damping or a generalized Morison equation with a linear drag term should be used with mid-fidelity engineering models to model surge decay properly.
  • Apart from not being able to capture the linear damping in surge, the mid-fidelity potential-flow model of TUHH generates mostly consistent predictions compared to the CFD simulations and shows a similar level of agreement with the experiment. The computing time, in terms of core hours, of a free-decay simulation using a typical mid-fidelity model is approximately 1/105 that of the corresponding CFD simulation. However, mid-fidelity potential-flow models are not fully predictive, and their accuracy strongly depends on the empirical drag coefficients used, which need to be tuned against some reference data. Therefore, to leverage the high efficiency of the mid-fidelity models fully, it is crucial to be able to obtain valid reference data from high-fidelity CFD simulations to help tune the model drag coefficients, especially during the initial design stage when the floater design is subject to frequent changes. This is, in fact, the primary goal of the present investigation.
In the future, we propose to leverage the observations made in the present investigation to determine how best to tune the engineering models of offshore wind turbines using CFD simulations, especially for complex novel floater geometries during the preliminary design phase. Another important question requiring further investigation is the applicability of the damping coefficients obtained from calm-water free-decay motion in a wave-like environment, because the ultimate goal is to predict the motion of the structure in realistic sea states accurately.

Author Contributions

Conceptualization, A.R., J.J. and L.W.; methodology, A.K., A.R., H.J., J.K., L.W., S.N. and Z.-R.S.; software, A.K., H.J., J.K., S.N. and Z.-R.S.; validation, L.W.; visualization, L.W.; formal analysis, L.W.; investigation, A.B.N., A.K., A.R., A.V., B.M.L., E.R., G.C.A., H.J., H.S., J.J., J.K., K.Y., L.R.R., L.W., M.H., P.C., Q.X., S.B., S.N., S.O., W.S., X.L., X.Z. and Z.-R.S.; data curation, L.W.; writing—original draft preparation, L.W.; writing—review and editing, A.B.N., A.K., A.R., A.V., B.M.L., E.R., G.C.A., H.J., H.S., J.J., J.K., K.Y., L.R.R., L.W., M.H., P.C., Q.X., S.B., S.N., S.O., W.S., X.L., X.Z. and Z.-R.S.; supervision, A.R. and J.J.; project administration, A.R. and L.W.; funding acquisition, A.R. and M.H. All authors have read and agreed to the published version of the manuscript.

Funding

Funding for this work was provided in part by the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Wind Energy Technologies Office under Contract No. DE-AC36-08GO28308. The University of Plymouth team would like to acknowledge the ongoing support of the Engineering and Physical Sciences Research Council (EPSRC) via project EP/T004177/1. MARIN would like to acknowledge that their contribution is partly funded by the Dutch Ministry of Economic Affairs and Climate Policy. The simulations by the Delft University of Technology team made use of the Dutch national e-infrastructure with the support of the SURF Cooperative with grant no. EINF-1649. The simulations by the Dalian University of Technology team were supported by the National Natural Science Foundation of China with grant no. 52071058.

Data Availability Statement

The numerical solutions analyzed in the present study will be made publicly accessible from the Data Archive and Portal of the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy at https://a2e.energy.gov/data/oc6/oc6.phase1a (accessed on: 24 September 2021).

Acknowledgments

This work was authored in part by the National Renewable Energy Laboratory, managed and operated by Alliance for Sustainable Energy, LLC for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes. A portion of this research was performed using computational resources sponsored by the Department of Energy’s Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A

Key aspects of the CFD numerical setups adopted by the participants are summarized in Table A1, Table A2, Table A3 and Table A4.
Table A1. CFD software and numerical discretization scheme.
Table A1. CFD software and numerical discretization scheme.
GroupCFD SoftwareTemporal SchemeSpatial SchemesInterface Capturing Scheme
ABSOpenFOAM Ver. 2006VOF: Crank-Nicolson; Others: 2nd order implicitGradient: Gauss linear; Divergence (momentum/turbulence): Gauss limitedLinear; Divergence (other): Gauss linear; Laplacian: Gauss linear correctedMULES
CENEROpenFOAM Ver. 18121st order implicitGradient: cellLimited Gauss linear 1; Divergence (momentum/turbulence): Gauss linearUpwind; Divergence (other): Gauss linear; Laplacian: Gauss linear correctedMULES; Gauss vanLeer
CLNKOpenFOAM Ver. 20062nd order implicitGradient: cellLimited Gauss linear 1; Divergence (momentum/turbulence): Gauss linear; Divergence (other): Gauss linear; Laplacian: Gauss linear correctedMULES; Gauss vanLeer
DTUOpenFOAM Ver. 1912Equal blending of first and second order implicitGradient: cellLimited Gauss linear 1; Divergence (momentum): Gauss vanLeer; Divergence (turbulence): Gauss upwind; Laplacian: Gauss linear correctedMULES; Gauss vanLeer
DUTSTAR-CCM+ Ver. 14.06.0132nd order implicit2nd order hybrid-BCD for convection; 2nd order hybrid Gauss-LSQ for gradient with Venkatakrishnan’s limiter; 2nd order upwind for turbulence quantitiesHRIC
IFPENOpenFOAM Ver. 18121st order implicitGradient: Gauss linear; Divergence (momentum/turbulence): Gauss upwind; Divergence (other): Gauss linear; Laplacian: Gauss linear correctedMULES; Gauss vanLeer
MARINReFRESCO Ver. 2.72nd order implicit2nd order harmonic TVD scheme of Van Leer; Gaussian least-squaresReFRICS [43]
NREL 1STAR-CCM+ Ver. 13.06.0122nd order implicit2nd order hybrid-BCD for convection; 2nd order hybrid Gauss-LSQ for gradient with Venkatakrishnan’s limiter; 2nd order upwind for turbulence quantitiesHRIC
TUDOpenFOAM Ver. 20121st order implicitGradient: cellLimited Gauss linear 1; Divergence: Gauss linear; Laplacian (LC4.2): Gauss linear corrected; Laplacian (LC4.4/4.6): Gauss linear limited 0.33MULES; Gauss vanLeer
UOPOpenFOAM Ver. 19061st order implicitTurbulence quantities: Gauss linearUpwind; Other quantities: Gauss linearMULES; Gauss MUSCL
UOSOpenFOAM Ver. 4.01st order implicitGradient Gauss linear; Divergence: Gauss limitedLinearV 1; Laplacian: Gauss linear correctedMULES; Gauss vanLeer01
1 NREL setup follows the baseline setup described in Section 3.
Table A2. Pressure-velocity coupling schemes.
Table A2. Pressure-velocity coupling schemes.
GroupAlgorithmNumber of PISO IterationsResidual ToleranceMaximum Outer Iterations
ABSPIMPLE2N/A6
CENERPIMPLE1N/A3
CLNKPISO3N/AN/A
DTUPISO3N/AN/A
DUTSIMPLEN/AN/A20
IFPENPISO3N/AN/A
MARINSIMPLEN/A10−450
NREL 1SIMPLEN/AN/A20
TUDPISO3N/AN/A
UOPPISO3N/AN/A
UOSPIMPLE2U: tol 10−9, relTol 10−3; p_rgh: tol 10−5, relTol 10−310
1 NREL setup follows the baseline setup described in Section 3.
Table A3. Time step and computational mesh (all dimensions at full scale).
Table A3. Time step and computational mesh (all dimensions at full scale).
GroupTime Step 1Cell Sizes at Free Surface 2 as Fractions of h = 6 mIsotropic Cell/Patch Sizes Near the Floater as Fractions of h = 6 mPrism Layers
Δ x = Δ y Δ z Near Floater 3Columns and Heave Plates 4Heave-Plate Corners and Btm. of Main Column 5Floater SurfaceNo. of LayersFirst LayerTotal Thickness
ABST/2501/61/151/41/81/321/1660.04 m0.45 m
CENERAdaptive: Co ≤ 15 at interface, Co ≤ 25 elsewhere1/41/81/41/161/641/16100.01415 m0.48 m
CLNKAdaptive: Co ≤ 1 everywhere1/0.61/101/41/401/401/4050.145 m1.22 m
DTUAdaptive: Co<51/81/81/61/81/321/1670.005 m0.48 m
DUTT/4001/41/81/41/161/641/16150.002 m0.48 m
IFPENAdaptive: Co ≤ 0.5 at interface, Co ≤ 2.5 elsewhere1/81/81/161/161/321/3270.00566 m0.185 m
MARINT/4001/41/81/41/161/161/16150.002 m0.48 m
NREL 6T/4001/41/81/41/161/641/16150.002 m0.48 m
TUDAdaptive: Co ≤ 1.51/41/43/41/81/321/1690.05 m0.48 m
UOPT/600 to T/7501/81/81/41/641/641/64No prism layers
UOSTmin/6001/41/81/81/81/161/1680.02 m0.4 m
1 Fixed time steps are expressed as fractions of the period, T , of the motion or as fractions of the heave period T m i n . Adaptive time step is based on the Courant number, Co. 2 Approximately corresponds to the first mesh refinement zone of Table 5. 3 Approximately corresponds to the first mesh refinement zone of Table 6. 4 Approximately corresponds to the rest of the mesh refinement zones of Table 6. 5 Approximately corresponds to the mesh refinement zones of Table 7. 6 NREL setup follows the baseline setup described in Section 3.
Table A4. Numerical domain, floater motion, and other settings.
Table A4. Numerical domain, floater motion, and other settings.
GroupDomain Size (Full Scale)Wave-Damping Zone
(Upstream/Downstream)
6Dof Motion Solver ParametersDynamic MeshTurb. Model
x (m) y (m) z (m)Half Domain
( y -Sym.)
Acceleration Relaxation
(OpenFOAM)
Acceleration Damping
(OpenFOAM)
DOF
ABS[−150, 190][−145, 145][−145, 145]No72 m/100 mN/A (ABS in-house)N/A (ABS in-house)6Overset (ABS in-house)SST k-ω
CENER[−200, 200][−100, 100][−180, 180]NoNone1.01.06MorphingSST k-ω
CLNK[−200, 200][−100, 100][−180, 100]No50 m/50 m0.61.06MorphingSST k-ω
DTU[−200, 200][−100, 100][−180, 180]NoNone0.95N/A6MorphingSST k-ω
DUT[−200, 200][−100, 100][−180, 180]No50 m/50 mN/AN/A6MorphingSpalart–Allmaras DES
IFPEN[−210, 210][−120, 120][−180, 150]NoNone1.01.06MorphingRNG k-ε
MARIN[−200, 200][−100, 100][−180, 180]No50 m/50 mN/AN/A6MorphingSST k-ω
NREL 1[−200, 200][−100, 100][−180, 180]No50 m/50 mN/AN/A6MorphingSpalart–Allmaras DES
TUD[−200, 200][−100, 100][−180, 180]NoNone1.01.06MorphingSST k-ω
UOP[−200, 200][−100, 100][−180, 180]No50 m/50 m0.71.06MorphingSpalart–Allmaras DES
UOS[−200, 200][−100, 100][−180, 180]No50 m/50 m0.71.06MorphingNone
1 NREL setup follows the baseline setup described in Section 3.

References

  1. Jonkman, J.M. The new modularization framework for the FAST Wind Turbine CAE Tool. In Proceedings of the 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Grapevine, TX, USA, 7–10 January 2013. [Google Scholar] [CrossRef] [Green Version]
  2. Wang, L.; Robertson, A.; Jonkman, J.; Yu, Y.-H. Uncertainty assessment of CFD investigation of the nonlinear difference-frequency wave loads on a semisubmersible FOWT platform. Sustainability 2021, 13, 64. [Google Scholar] [CrossRef]
  3. Wang, L.; Robertson, A.; Jonkman, J.; Yu, Y.-H.; Koop, A.; Borràs Nadal, A.; Li, H.; Bachynski, E.; Pinguet, R.; Shi, W.; et al. OC6 Phase Ib: Validation of the CFD predictions of difference-frequency wave excitation on an FOWT semisubmersible. Ocean Eng. 2021, 241, 110026. [Google Scholar] [CrossRef]
  4. Wang, L.; Robertson, A.; Jonkman, J.; Yu, Y.-H.; Koop, A.; Borràs Nadal, A.; Li, H.; Shi, W.; Pinguet, R.; Zhou, Y.; et al. Investigation of nonlinear difference-frequency wave excitation on a semisubmersible offshore-wind platform with bichromatic-wave CFD simulations. In Proceedings of the ASME 2021 3rd International Offshore Wind Technical Conference, Boston, MA, USA, 16–17 February 2021. [Google Scholar] [CrossRef]
  5. Böhm, M.; Robertson, A.; Hübler, C.; Rolfes, R.; Schaumann, P. Optimization-based calibration of hydrodynamic drag coefficients for a semisubmersible platform using experimental data of an irregular sea state. J. Phys. Conf. Ser. 2020, 1669, 012023. [Google Scholar] [CrossRef]
  6. Robertson, A.; Jonkman, J.; Masciola, M.; Song, H.; Goupee, A.; Coulling, A.; Luan, C. Definition of the Semisubmersible Floating System for Phase II of OC4; National Renewable Energy Laboratory: Golden, CO, USA, 2014. Available online: https://www.nrel.gov/docs/fy14osti/60601.pdf (accessed on 24 September 2021).
  7. Robertson, A.; Wendt, F.; Jonkman, J.; Popko, W.; Dagher, H.; Gueydon, S.; Qvist, J.; Vittori, F.; Azcona, J.; Uzunoglu, E. OC5 Project Phase II: Validation of global loads of the DeepCwind floating semisubmersible wind turbine. Energy Procedia 2017, 137, 38–57. [Google Scholar] [CrossRef]
  8. Coulling, A.; Goupee, A.; Robertson, A.; Jonkman, J.; Dagher, H. Validation of a FAST semi-submersible floating wind turbine numerical model with DeepCwind test data. J. Renew. Sustain. Energy 2013, 5, 023116. [Google Scholar] [CrossRef]
  9. Robertson, A.; Gueydon, S.; Bachynski, E.; Wang, L.; Jonkman, J.; Alarcón, D.; Amet, E.; Beardsell, A.; Bonnet, P.; Boudet, B.; et al. OC6 Phase I: Investigating the underprediction of low-frequency hydrodynamic loads and responses of a floating wind turbine. J. Phys. Conf. Ser. 2020, 1618, 032033. [Google Scholar] [CrossRef]
  10. Robertson, A.; Bachynski, E.; Gueydon, S.; Wendt, F.; Schünemann, P. Total experimental uncertainty in hydrodynamic testing of a semisubmersible wind turbine, considering numerical propagation of systematic uncertainty. Ocean Eng. 2020, 195, 106605. [Google Scholar] [CrossRef]
  11. Tran, T.T.; Kim, D.H. The coupled dynamic response computation for a semi-submersible platform of floating offshore wind turbine. J. Wind Eng. Ind. Aerodyn. 2015, 147, 104–119. [Google Scholar] [CrossRef]
  12. Tran, T.T.; Kim, D.H. Fully coupled aero-hydrodynamic analysis of a semi-submersible FOWT using a dynamic fluid body interaction approach. Renew. Energy 2016, 92, 244–261. [Google Scholar] [CrossRef]
  13. Burmester, S.; Vaz, G.; el Moctar, O. Towards credible CFD simulations for floating offshore wind turbines. Ocean Eng. 2020, 209, 107237. [Google Scholar] [CrossRef]
  14. Eça, L.; Hoekstra, M. A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies. J. Comput. Phys. 2014, 262, 104–130. [Google Scholar] [CrossRef]
  15. Wang, Y.; Chen, H.-C.; Vaz, G.; Burmester, S. CFD simulation of semi-submersible floating offshore wind turbine under pitch decay motion. In Proceedings of the ASME 2019 2nd International Offshore Wind Technical Conference, St. Julian’s, Malta, 3–6 November 2019. [Google Scholar] [CrossRef]
  16. Wang, Y.; Chen, H.-C.; Koop, A.; Vaz, G. Verification and validation of CFD simulations for semi-submersible floating offshore wind turbine under pitch free-decay motion. Ocean Eng. 2021, 242, 109993. [Google Scholar] [CrossRef]
  17. Li, H.; Bachynski-Polić, E. Experimental and numerically obtained low-frequency radiation characteristics of the OC5-DeepCwind semisubmersible. Ocean Eng. 2021, 232, 109130. [Google Scholar] [CrossRef]
  18. Burmester, S.; Vaz, G.; Gueydon, S.; el Moctar, O. Investigation of a semi-submersible floating wind turbine in surge decay using CFD. Ship Technol. Res. 2020, 67, 2–14. [Google Scholar] [CrossRef]
  19. Burmester, S.; Vaz, G.; el Moctar, O.; Gueydon, S.; Koop, A.; Wang, Y.; Chen, H. High-fidelity modelling of floating offshore wind turbine platforms. In Proceedings of the ASME 2020 39th International Conference on Ocean, Offshore and Arctic Engineering, Fort Lauderdale, FL, USA, 3–7 August 2020. [Google Scholar] [CrossRef]
  20. Wang, Y.; Chen, H.-C.; Vaz, G.; Burmester, S. CFD simulation of semi-submersible floating offshore wind turbine under regular waves. In Proceedings of the 30th International Ocean and Polar Engineering Conference, Shanghai, China, 11–16 October 2020. [Google Scholar]
  21. Wang, Y.; Chen, H.-C.; Vaz, G.; Mewes, S. Verification study of CFD simulation of semi-submersible floating offshore wind turbine under regular waves. In Proceedings of the ASME 2021 3rd International Conference on Ocean, Offshore and Arctic Engineering, Boston, MA, USA, 16–17 February 2021. [Google Scholar] [CrossRef]
  22. Bozonnet, P.; Emery, A. CFD simulations for the design of offshore floating wind platforms encompassing heave plates. In Proceedings of the 25th International Ocean and Polar Engineering Conference, Kona, HI, USA, 21–26 June 2015. [Google Scholar]
  23. Zhang, Y.; Kim, B. A fully coupled computational fluid dynamics method for analysis of semi-submersible floating offshore wind turbines under wind-wave excitation conditions based on OC5 data. Appl. Sci. 2018, 8, 2314. [Google Scholar] [CrossRef] [Green Version]
  24. Gueydon, S.; van Alfen, R. MARINET2 OC6 Model-Tests: A Series of Tests Focusing on the Hydrodynamics of the OC5 Semisubmersible; Maritime Research Institute Netherlands: Wageningen, The Netherlands, 2018. [Google Scholar]
  25. Siemens PLM Software. Simcenter STAR-CCM+ User Guide; Version 13.06.012; Siemens PLM Software: Plano, TX, USA, 2018. [Google Scholar]
  26. Weller, H.G.; Tabor, G.; Jasak, H.; Fureby, C. A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput. Phys. 1998, 12, 620–631. [Google Scholar] [CrossRef]
  27. ReFRESCO. Available online: www.refresco.org (accessed on 21 May 2021).
  28. Ferreira González, D.; Göttsche, U.; Netzband, S.; Abdel-Maksoud, M. Advances on simulations of wave-body interactions under consideration of the nonlinear free water surface. Ship Technol. Res. 2021, 68, 27–40. [Google Scholar] [CrossRef]
  29. Morison, J.; O’Brien, M.; Johnson, J.; Schaaf, S. The force exerted by surface waves on piles. J. Pet. Technol. 1950, 2, 149–154. [Google Scholar] [CrossRef]
  30. Hirt, C.W.; Nichols, B.D. Volume of Fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 1981, 39, 201–225. [Google Scholar] [CrossRef]
  31. Cifani, P.; Michalek, W.; Priems, G.; Kuerten, J.; van der Geld, C.; Geurts, B. A comparison between the surface compression method and an interface reconstruction method for the VOF approach. Comput. Fluids 2016, 136, 421–435. [Google Scholar] [CrossRef] [Green Version]
  32. Spalart, P.; Jou, W.-H.; Strelets, M.; Allmaras, S. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES. In Advances in DNS/LES; Greyden Press: Columbus, OH, USA, 1997; pp. 137–147. [Google Scholar]
  33. Shur, M.; Spalart, P.; Strelets, M.; Travin, A. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities. Int. J. Heat Fluid Flow 2008, 29, 1638–1649. [Google Scholar] [CrossRef]
  34. Spalart, P.; Deck, S.; Shur, M.; Squires, K.; Strelets, M.; Travin, A. A new version of detached eddy simulation, resistant to ambiguous grid densities. Theor. Comput. Fluid Dyn. 2006, 20, 181–195. [Google Scholar] [CrossRef]
  35. Ren, C.; Lu, L.; Cheng, L.; Chen, T. Hydrodynamic damping of an oscillating cylinder at small Keulegan-Carpenter numbers. J. Fluid Mech. 2021, 913, A36. [Google Scholar] [CrossRef]
  36. Kim, J.; Baquet, A.; Jang, H. Wave propagation in CFD-based numerical wave tank. In Proceedings of the ASME 2019 38th International Conference on Ocean, Offshore and Arctic Engineering, Glasgow, UK, 9–14 June 2019. [Google Scholar] [CrossRef]
  37. Helder, J.A.; Pietersma, M. UMaine-DeepCwind/OC4 Semi Floating Wind Turbine Repeat Tests; Maritime Research Institute Netherlands: Wageningen, The Netherlands, 2013. [Google Scholar]
  38. Eça, L.; Vaz, G.; Toxopeus, S.L.; Hoekstra, M. Numerical errors in unsteady flow simulations. J. Verif. Valid. Uncert. 2019, 4, 021001. [Google Scholar] [CrossRef]
  39. Wang, W.; Wu, M.; Palm, J.; Eskilsson, C. Estimation of numerical uncertainty in computational fluid dynamics simulations of a passively controlled wave energy converter. Proc. Inst. Mech. Eng. Part M J. Eng. Marit. Environ. 2018, 232, 71–84. [Google Scholar] [CrossRef] [Green Version]
  40. Eça, L.; Hoekstra, M. Evaluation of numerical error estimation based on grid refinement studies with the method of the manufactured solutions. Comput. Fluids 2009, 38, 1580–1591. [Google Scholar] [CrossRef]
  41. Menter, F. Two-equation eddy-viscosity turbulence modeling for engineering applications. AIAA J. 1994, 32, 1598–1605. [Google Scholar] [CrossRef] [Green Version]
  42. Bureau Veritas. Classification of Mooring Systems for Permanent and Mobile Offshore Units Rule Note NR 493; Bureau Veritas: Paris, France, 2021; Available online: https://erules.veristar.com/dy/data/bv/pdf/493-NR_2021-07.pdf (accessed on 21 September 2021).
  43. Klaij, C.M.; Hoekstra, M.; Vaz, G. Design, analysis and verification of a volume-of-fluid model with interface-capturing scheme. Comput. Fluids 2018, 170, 324–340. [Google Scholar] [CrossRef]
Figure 1. The OC6-DeepCwind floating offshore wind turbine (FOWT) semisubmersible. (a) The geometry of the floater and the adopted coordinate system. The surge, sway, and heave motions are along the x-, y-, and z-directions, respectively. The roll, pitch, and yaw motions are also about the x-, y-, and z-axes, respectively. (b) Setup of the free-decay experiment in the MARIN Concept Basin (photo by Amy Robertson, NREL).
Figure 1. The OC6-DeepCwind floating offshore wind turbine (FOWT) semisubmersible. (a) The geometry of the floater and the adopted coordinate system. The surge, sway, and heave motions are along the x-, y-, and z-directions, respectively. The roll, pitch, and yaw motions are also about the x-, y-, and z-axes, respectively. (b) Setup of the free-decay experiment in the MARIN Concept Basin (photo by Amy Robertson, NREL).
Energies 15 00389 g001
Figure 2. Numerical domain for CFD simulations. The three taut-spring mooring lines attached to the heave plates are also shown.
Figure 2. Numerical domain for CFD simulations. The three taut-spring mooring lines attached to the heave plates are also shown.
Energies 15 00389 g002
Figure 3. Vertical section of the baseline mesh.
Figure 3. Vertical section of the baseline mesh.
Energies 15 00389 g003
Figure 4. Free-decay motion of the floater in (a) surge (Load Case [LC] 4.2), (b) heave (LC 4.4), and (c) pitch (LC 4.6). Legend: CENER = National Renewable Energy Centre; CLNK = ClassNK; DTU = Technical University of Denmark; DUT = Dalian University of Technology; IFPEN = IFP Energies nouvelles; TUD = Delft University of Technology; UOP = University of Plymouth; UOS = University of Strathclyde; TUHH = Hamburg University of Technology; ABS = American Bureau of Shipping; MARIN = Maritime Research Institute Netherlands; NREL = National Renewable Energy Laboratory; EXP = experiment.
Figure 4. Free-decay motion of the floater in (a) surge (Load Case [LC] 4.2), (b) heave (LC 4.4), and (c) pitch (LC 4.6). Legend: CENER = National Renewable Energy Centre; CLNK = ClassNK; DTU = Technical University of Denmark; DUT = Dalian University of Technology; IFPEN = IFP Energies nouvelles; TUD = Delft University of Technology; UOP = University of Plymouth; UOS = University of Strathclyde; TUHH = Hamburg University of Technology; ABS = American Bureau of Shipping; MARIN = Maritime Research Institute Netherlands; NREL = National Renewable Energy Laboratory; EXP = experiment.
Energies 15 00389 g004
Figure 5. Regression analyses extracting the surge-damping coefficients from (a) the CFD simulations performed by ABS, MARIN, and NREL (note that the NREL simulation is based on the baseline numerical setup of Section 3) and (b) the experimental measurement. The symbols are the data points from the surge free-decay motion of the platform, and the lines or curve are the best fits of the data points. In (a), the fits are of the form given by Equation (9), which only considers linear and quadratic damping. In (b), the fit is given by Equation (11), which also includes a constant Coulomb-friction-like damping force.
Figure 5. Regression analyses extracting the surge-damping coefficients from (a) the CFD simulations performed by ABS, MARIN, and NREL (note that the NREL simulation is based on the baseline numerical setup of Section 3) and (b) the experimental measurement. The symbols are the data points from the surge free-decay motion of the platform, and the lines or curve are the best fits of the data points. In (a), the fits are of the form given by Equation (9), which only considers linear and quadratic damping. In (b), the fit is given by Equation (11), which also includes a constant Coulomb-friction-like damping force.
Energies 15 00389 g005
Figure 6. Comparison of the MARIN CFD solutions for surge decay (LC 4.2) with and without B 0 = 2.3 kN full scale. The experimental measurement (EXP) is included for reference. (a) Surge time series and (b) regression analysis for extracting the surge-damping coefficients.
Figure 6. Comparison of the MARIN CFD solutions for surge decay (LC 4.2) with and without B 0 = 2.3 kN full scale. The experimental measurement (EXP) is included for reference. (a) Surge time series and (b) regression analysis for extracting the surge-damping coefficients.
Energies 15 00389 g006
Figure 7. Regression analyses extracting the heave-damping coefficients from (a) the CFD simulations performed by ABS, MARIN, and NREL (note that the NREL simulation is based on the baseline numerical setup of Section 3) and (b) the experimental measurement. The symbols are the data points from the heave free-decay motion of the platform, and the lines are the best fits given by Equation (9).
Figure 7. Regression analyses extracting the heave-damping coefficients from (a) the CFD simulations performed by ABS, MARIN, and NREL (note that the NREL simulation is based on the baseline numerical setup of Section 3) and (b) the experimental measurement. The symbols are the data points from the heave free-decay motion of the platform, and the lines are the best fits given by Equation (9).
Energies 15 00389 g007
Figure 8. Regression analyses extracting the pitch-damping coefficients from (a) the CFD simulations performed by ABS, MARIN, and NREL (note that the NREL simulation is based on the baseline numerical setup of Section 3) and (b) the experimental measurement. The symbols are the data points from the pitch free-decay motion of the platform, and the lines are the best fits given by Equation (9).
Figure 8. Regression analyses extracting the pitch-damping coefficients from (a) the CFD simulations performed by ABS, MARIN, and NREL (note that the NREL simulation is based on the baseline numerical setup of Section 3) and (b) the experimental measurement. The symbols are the data points from the pitch free-decay motion of the platform, and the lines are the best fits given by Equation (9).
Energies 15 00389 g008
Figure 9. Unnormalized root-mean-square (RMS) residual of the x -moment equation at the end of each time step.
Figure 9. Unnormalized root-mean-square (RMS) residual of the x -moment equation at the end of each time step.
Energies 15 00389 g009
Figure 10. Time series of (a) the surge-decay motion and (b) the heave-decay motion obtained from the CFD simulations with different numbers of SIMPLE iterations per time step.
Figure 10. Time series of (a) the surge-decay motion and (b) the heave-decay motion obtained from the CFD simulations with different numbers of SIMPLE iterations per time step.
Energies 15 00389 g010aEnergies 15 00389 g010b
Figure 11. Time series of (a) the surge-decay motion and (b) the heave-decay motion obtained from the CFD simulations with the four different levels of numerical refinement described in Table 11.
Figure 11. Time series of (a) the surge-decay motion and (b) the heave-decay motion obtained from the CFD simulations with the four different levels of numerical refinement described in Table 11.
Energies 15 00389 g011
Figure 12. Convergence of (a) the linear damping coefficient, (b) quadratic damping coefficient, and (c) equivalent linear damping ratio in surge with simultaneous time-step and grid refinement. The crosses are the CFD solutions computed with the four different refinement ratios in Table 11, and the estimated discretization uncertainties are shown as symmetric uncertainty bands attached to the crosses.
Figure 12. Convergence of (a) the linear damping coefficient, (b) quadratic damping coefficient, and (c) equivalent linear damping ratio in surge with simultaneous time-step and grid refinement. The crosses are the CFD solutions computed with the four different refinement ratios in Table 11, and the estimated discretization uncertainties are shown as symmetric uncertainty bands attached to the crosses.
Energies 15 00389 g012
Figure 13. Convergence of (a) the linear damping coefficient and (b) the quadratic damping coefficient in heave with simultaneous time-step and grid refinement. The crosses are the CFD solutions computed with the four different refinement ratios in Table 11, and the estimated discretization uncertainties are shown as symmetric uncertainty bands attached to the crosses.
Figure 13. Convergence of (a) the linear damping coefficient and (b) the quadratic damping coefficient in heave with simultaneous time-step and grid refinement. The crosses are the CFD solutions computed with the four different refinement ratios in Table 11, and the estimated discretization uncertainties are shown as symmetric uncertainty bands attached to the crosses.
Energies 15 00389 g013
Figure 14. Periods of free-decay motions in surge, heave, and pitch from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH). For interpretation of the legends, which denote different participants of the OC6 Phase Ia project, please refer to the caption of Figure 4. The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table 11.
Figure 14. Periods of free-decay motions in surge, heave, and pitch from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH). For interpretation of the legends, which denote different participants of the OC6 Phase Ia project, please refer to the caption of Figure 4. The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table 11.
Energies 15 00389 g014
Figure 15. Comparison of the linear and quadratic damping in surge, characterized by the values of P and Q , respectively, and the equivalent linear damping ratio, ζ , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH). The P and Q values are calculated from the motion time series using the P Q analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients, B 1 and B 2 , respectively, following Equation (10). The equivalent linear damping ratio, ζ , is calculated using Equation (13). The amplitude factor F ˜ A used to normalize Q is 2.8 m. For interpretation of the legends, which denote different participants of the OC6 Phase Ia project, please refer to the caption of Figure 4. The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table 11. The uncertainty bands attached to the NREL solutions are the numerical discretization uncertainties estimated in Section 6.2.
Figure 15. Comparison of the linear and quadratic damping in surge, characterized by the values of P and Q , respectively, and the equivalent linear damping ratio, ζ , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH). The P and Q values are calculated from the motion time series using the P Q analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients, B 1 and B 2 , respectively, following Equation (10). The equivalent linear damping ratio, ζ , is calculated using Equation (13). The amplitude factor F ˜ A used to normalize Q is 2.8 m. For interpretation of the legends, which denote different participants of the OC6 Phase Ia project, please refer to the caption of Figure 4. The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table 11. The uncertainty bands attached to the NREL solutions are the numerical discretization uncertainties estimated in Section 6.2.
Energies 15 00389 g015
Figure 16. Comparison of the linear and quadratic damping in heave, characterized by the values of P and Q , respectively, and the equivalent linear damping ratio, ζ , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH). The P and Q values are calculated from the motion time series using the P Q analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients, B 1 and B 2 , respectively, following Equation (10). The equivalent linear damping ratio, ζ , is calculated using Equation (13). The amplitude factor F ˜ A used to normalize Q is 0.48 m. For interpretation of the legends, which denote different participants of the OC6 Phase Ia project, please refer to the caption of Figure 4. The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table 11. The uncertainty bands attached to the NREL solutions are the numerical discretization uncertainties estimated in Section 6.2.
Figure 16. Comparison of the linear and quadratic damping in heave, characterized by the values of P and Q , respectively, and the equivalent linear damping ratio, ζ , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH). The P and Q values are calculated from the motion time series using the P Q analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients, B 1 and B 2 , respectively, following Equation (10). The equivalent linear damping ratio, ζ , is calculated using Equation (13). The amplitude factor F ˜ A used to normalize Q is 0.48 m. For interpretation of the legends, which denote different participants of the OC6 Phase Ia project, please refer to the caption of Figure 4. The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table 11. The uncertainty bands attached to the NREL solutions are the numerical discretization uncertainties estimated in Section 6.2.
Energies 15 00389 g016
Figure 17. Comparison of the linear and quadratic damping in pitch, characterized by the values of P and Q , respectively, and the equivalent linear damping ratio, ζ , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH). The P and Q values are calculated from the motion time series using the P Q analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients, B 1 and B 2 , respectively, following Equation (10). The equivalent linear damping ratio, ζ , is calculated using Equation (13). The amplitude factor F ˜ A used to normalize Q is 0.053 radian. For interpretation of the legends, which denote different participants of the OC6 Phase Ia project, please refer to the caption of Figure 4.
Figure 17. Comparison of the linear and quadratic damping in pitch, characterized by the values of P and Q , respectively, and the equivalent linear damping ratio, ζ , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH). The P and Q values are calculated from the motion time series using the P Q analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients, B 1 and B 2 , respectively, following Equation (10). The equivalent linear damping ratio, ζ , is calculated using Equation (13). The amplitude factor F ˜ A used to normalize Q is 0.053 radian. For interpretation of the legends, which denote different participants of the OC6 Phase Ia project, please refer to the caption of Figure 4.
Energies 15 00389 g017
Table 1. Principal full-scale dynamic properties of the OC6-DeepCwind FOWT.
Table 1. Principal full-scale dynamic properties of the OC6-DeepCwind FOWT.
ParametersValueUnit
Displaced   Volume ,   1.392 × 104m3
Vertical Center of Buoyancy, VCB (from SWL 1)−13.17m
Vertical Center of Gravity, VCG (from SWL)−7.53m
Mass ,   m 1.407 × 107kg
Roll   Moment   of   Inertia ,   I x x about CG 21.2898 × 1010kg-m2
Pitch   Moment   of   Inertia ,   I y y about CG1.2851 × 1010kg-m2
Yaw   Moment   of   Inertia ,   I z z about CG1.4189 × 1010kg-m2
1 SWL = still water level. 2 CG = center of gravity.
Table 2. Geometry of the mooring lines at equilibrium.
Table 2. Geometry of the mooring lines at equilibrium.
LinesEnd Points x (m) y (m) z (m)
1Fairlead (FL1)−40.870.00−14.0
Anchor (AC1)−105.470.00−58.4
2Fairlead (FL2)20.43−35.39−14.0
Anchor (AC2)52.73−91.34−58.4
3Fairlead (FL3)20.4335.39−14.0
Anchor (AC3)52.7391.34−58.4
Table 3. Initial displacements in all six degrees of freedom (6DoF) from the selected free-decay experiments.
Table 3. Initial displacements in all six degrees of freedom (6DoF) from the selected free-decay experiments.
Load CaseSurge (m)Sway (m)Heave (m)Roll (deg)Pitch (deg)Yaw (deg)
4.2—Surge Decay−5.10.90.00.10.70.3
4.4—Heave Decay0.1−0.3−2.20.20.30.1
4.6—Pitch Decay−2.10.1−0.1−0.1−5.70.0
Table 4. Initial velocities in all 6DoF from the selected free-decay experiments.
Table 4. Initial velocities in all 6DoF from the selected free-decay experiments.
Load CaseSurge (m/s)Sway (m/s)Heave (m/s)Roll (deg/s)Pitch (deg/s)Yaw (deg/s)
4.2—Surge Decay0.00.00.00.00.00.0
4.4—Heave Decay0.00.00.0−0.10.00.0
4.6—Pitch Decay−0.20.00.00.00.00.0
Table 5. Baseline domain boundary locations and conditions for the computational fluid dynamics (CFD) simulations.
Table 5. Baseline domain boundary locations and conditions for the computational fluid dynamics (CFD) simulations.
BoundaryLocationTypeVelocityPressurePhase Fraction
FloaterN/AWallNo slipZero normal gradientZero normal gradient
Upstream x = 200 mVelocity InletZeroZero normal gradient z > 0 :   air ;   z 0 : water
Downstream x = + 200 mVelocity InletZeroZero normal gradient z > 0 :   air ;   z 0 : water
Sides y = ± 100 mWallFree slipZero normal gradientZero normal gradient
Bottom z = 180 mVelocity InletZeroZero normal gradientWater only
Top z = + 180 mPressure OutletExtrapolated backflow dir.AtmosphericAir only
Table 6. Free-surface mesh-refinement zones with target cell sizes as fractions of h .
Table 6. Free-surface mesh-refinement zones with target cell sizes as fractions of h .
Refinement Zones Δ x / h Δ y / h Δ z / h
z [−5 m, 5 m]1/41/41/8
z [−40 m, 20 m]1/21/21/4
z [−80 m, 30 m]111/2
Table 7. Box-shaped mesh-refinement zones with target cell sizes as fractions of h .
Table 7. Box-shaped mesh-refinement zones with target cell sizes as fractions of h .
Refinement Zones x Range (m)y Range (m) z Range (m) ( Δ x , Δ y , Δ z ) / h
Near Floater[−60, 46][−50, 50][−40, 20]1/4
Main Column[−7, 7][−4, 4][−22, 12]1/16
Upstream Col.—Heave Plate[−45, −13][−14, 14][−23, −11]1/16
Starboard Col.—Heave Plate[−2, 30][11, 39][−23, −11]1/16
Port Col.—Heave Plate[−2, 30][−39, −11][−23, −11]1/16
Upstream Col.—Upper Part[−38.5, −19.5][−7, 7][−11, 12]1/16
Starboard Col.—Upper Part[5, 24][18, 32][−11, 12]1/16
Port Col.—Upper Part[5, 24][−32, −18][−11, 12]1/16
Floater Surface MeshN/AN/AN/A1/16
Table 8. Corner mesh-refinement zones with target cell sizes as fractions of h .
Table 8. Corner mesh-refinement zones with target cell sizes as fractions of h .
Refinement Zones R Range (m) z Range (m) ( Δ x , Δ y , Δ z ) / h
Bottom of the Main Column[0, 4.25][−21, −13]1/64
Top Edge of Heave Plates[10.75, 13.25][−15.25, −12.75]1/64
Bottom Edge of Heave Plates[10.75, 13.25][−21.25, −18.75]1/64
Table 9. Variation in surge-damping coefficients and equivalent linear damping ratio with the number of iterations.
Table 9. Variation in surge-damping coefficients and equivalent linear damping ratio with the number of iterations.
No. Iter. P = π ω B 1 2 k Q = 4 ω 2 B 2 3 k ζ
ValueDifferenceValue (m−1)DifferenceValueDifference
50.11297%0.032720%6.06%45%
100.06412%0.02653.5%4.32%2.9%
200.0593.4%0.02692.0%4.22%0.4%
400.057N/A0.0274N/A4.20%N/A
Table 10. Variation in heave-damping coefficients and equivalent linear damping ratio with the number of iterations.
Table 10. Variation in heave-damping coefficients and equivalent linear damping ratio with the number of iterations.
No. Iter. P = π ω B 1 2 k Q = 4 ω 2 B 2 3 k ζ
ValueDifferenceValue (m−1)DifferenceValueDifference
20−0.00396%0.1700.8%2.306%0.4%
40−0.0041N/A0.171N/A2.313%N/A
Table 11. Numerical configurations for convergence study.
Table 11. Numerical configurations for convergence study.
Case Refinement   Ratio ,   λ Reference   Cell   Size ,   h Number of Cells Time   Step ,   Δ t No. Iter.
Very Coarse4/38 m6.7 million T / 300 20
Coarse7/67 m8.6 million T / 343 20
Baseline16 m12.9 million T / 400 20
Fine3/44.5 m25.6 million T / 533 20
Table 12. Discretization uncertainties of the surge-damping coefficients and the equivalent linear damping ratio.
Table 12. Discretization uncertainties of the surge-damping coefficients and the equivalent linear damping ratio.
Cases P = π ω B 1 2 k Q = 4 ω 2 B 2 3 k ζ
ValueDiscretization UncertaintyValue (m−1)Discretization UncertaintyValueDiscretization Uncertainty
Very Coarse0.075±60%0.0235±38%4.5%±16%
Coarse0.071±50%0.0244±29%4.4%±13%
Baseline0.059±47%0.0269±21%4.2%±10%
Fine0.055±48%0.0276±19%4.1%±10%
Table 13. Discretization uncertainties of the heave-damping coefficients and the equivalent linear damping ratio.
Table 13. Discretization uncertainties of the heave-damping coefficients and the equivalent linear damping ratio.
Cases P = π ω B 1 2 k Q = 4 ω 2 B 2 3 k ζ
ValueDiscretization Uncertainty 1Value (m−1)Discretization UncertaintyValueDiscretization Uncertainty
Very Coarse−0.018±0.040.20±55%2.28%±7%
Coarse−0.007±0.040.18±51%2.31%±7%
Baseline−0.004±0.030.17±42%2.31%±7%
Fine+0.003±0.030.15±33%2.26%±7%
1 Discretization uncertainties of P are expressed in terms of actual values instead of percentages because P is very close to zero.
Table 14. Sensitivity of the damping coefficients and equivalent linear damping ratio to the choice of turbulence models.
Table 14. Sensitivity of the damping coefficients and equivalent linear damping ratio to the choice of turbulence models.
Surge DampingHeave Damping
P Q (m−1) ζ (%) P Q (m−1) ζ (%)
SA-DES 1 (Baseline)0.0590.02694.22−0.0040.1702.31
SST   2   k ω 0.0610.02724.29−0.0010.1692.38
No Turbulence Model0.0570.02904.34−0.0040.1712.31
Max. Difference Relative to Baseline3.3%8.0%2.9%0.003 30.9%3.1%
1 SA-DES = Spalart–Allmaras detached-eddy simulation. 2 SST = shear stress transport. 3 The maximum difference in P is not given as percentages because P is very close to zero for heave decay.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wang, L.; Robertson, A.; Jonkman, J.; Kim, J.; Shen, Z.-R.; Koop, A.; Borràs Nadal, A.; Shi, W.; Zeng, X.; Ransley, E.; et al. OC6 Phase Ia: CFD Simulations of the Free-Decay Motion of the DeepCwind Semisubmersible. Energies 2022, 15, 389. https://doi.org/10.3390/en15010389

AMA Style

Wang L, Robertson A, Jonkman J, Kim J, Shen Z-R, Koop A, Borràs Nadal A, Shi W, Zeng X, Ransley E, et al. OC6 Phase Ia: CFD Simulations of the Free-Decay Motion of the DeepCwind Semisubmersible. Energies. 2022; 15(1):389. https://doi.org/10.3390/en15010389

Chicago/Turabian Style

Wang, Lu, Amy Robertson, Jason Jonkman, Jang Kim, Zhi-Rong Shen, Arjen Koop, Adrià Borràs Nadal, Wei Shi, Xinmeng Zeng, Edward Ransley, and et al. 2022. "OC6 Phase Ia: CFD Simulations of the Free-Decay Motion of the DeepCwind Semisubmersible" Energies 15, no. 1: 389. https://doi.org/10.3390/en15010389

APA Style

Wang, L., Robertson, A., Jonkman, J., Kim, J., Shen, Z. -R., Koop, A., Borràs Nadal, A., Shi, W., Zeng, X., Ransley, E., Brown, S., Hann, M., Chandramouli, P., Viré, A., Ramesh Reddy, L., Li, X., Xiao, Q., Méndez López, B., Campaña Alonso, G., ... Yu, K. (2022). OC6 Phase Ia: CFD Simulations of the Free-Decay Motion of the DeepCwind Semisubmersible. Energies, 15(1), 389. https://doi.org/10.3390/en15010389

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop