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

: Currently, the design of ﬂoating offshore wind systems is primarily based on mid-ﬁdelity models with empirical drag forces. The tuning of the model coefﬁcients requires data from either experiments or high-ﬁdelity simulations. As part of the OC6 (Offshore Code Comparison Collab-oration, 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 veriﬁcation and validation study of computational ﬂuid 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 co-efﬁcients needed by mid-ﬁdelity models. The linear and quadratic damping coefﬁcients and the equivalent damping ratio are chosen as metrics for validation. Large numerical uncertainties are estimated for the linear and quadratic damping coefﬁcients; 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-ﬁdelity models


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 wavebasin 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 kmodel 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 followup 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 leastsquares method of [14].
Bozonnet and Emery performed CFD simulations of the forced oscillation and freedecay 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 crossverification 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.

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.
tively 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).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/m 3 , but a full-scale density of 1025 kg/m 3 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 Tables 3 and 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).

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 Sections 3.2-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 Tables A1-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 timedomain 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.

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: 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]: 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: 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.

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 freedecay 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.

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 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 largerdiameter (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.Table 6.Free-surface mesh-refinement zones with target cell sizes as fractions of h.

Refinement Zones
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.

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-

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.

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 × 10 7 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.

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.
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   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.

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.

Estimation of the Linear and Quadratic Damping Coefficients Using PQ Analysis
The linear and quadratic damping coefficients, B 1 and B 2 , can be estimated from the motion time series using the PQ method [37], which is briefly summarized here.For the weakly damped free-decay motions, we define the amplitude decrease over a halfcycle as ∆A i = A i − A i+1 , where A i and A i+1 are, respectively, the positive amplitudes of the ith 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 ith half-cycle is approximated as A i = (A i + A i+1 )/2.The energy loss over the ith half-cycle, L i , (neglecting other nonlinear damping) is given by x )dt, 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 ith half-cycle: .
where ω is the angular frequency of the motion, the following relation can be derived from Equation (7): 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: 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 Sections 5.1.1-5.1.3,example results from the PQ 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 PQ method [37], results in poor linear regression for the pitch-decay motion (see Section 5.1.3for details).Nevertheless, we believe the PQ analysis can still provide a meaningful characterization of the pitch-decay time series that enables comparisons across different results.

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 PQ analysis to minimize any start-up effect.The experimental surge free-decay motion shows poor linear regression when the standard PQ 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 PQ method can be modified by adding an additional term in the form of B 0 . .
x to the integrand on the right-hand side of Equation (7).With this addition, Equation ( 9) becomes where O = 2B 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 , 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.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 Coulombfriction-like damping force.
To verify the modified  analysis with Equation (11) further and investigate the effect of  , MARIN repeated the surge-decay CFD simulation by applying an additional force of  = 2.3 kN full scale to the floater center of gravity in the ±-direction opposing the floater surge velocity.The results from this additional simulation are labeled MARIN-B0 to distinguish from the MARIN solution without  .Note that for consistency, all numerical results in this article other than MARIN-B0 were performed with- 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.
To verify the modified PQ 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.

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.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 PQ 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.

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.

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.

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 PQ 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 PQ analysis for pitch free decay.
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  analysis for pitch free decay.

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  , the linear equivalent damping ratio  =  /(2) can also be computed from  and  using the following expression: where the motion-amplitude scaling factor,  , is given by

Equivalent Linear Damping Ratio
By matching the energy loss from the full linear quadratic damping model with an equivalent linear damping model with damping coefficient B e , the linear equivalent damping ratio ζ = ωB e /(2k) can also be computed from P and Q using the following expression: where the motion-amplitude scaling factor, F A , is given by with m and n being the numbers of the first and last peaks/troughs, respectively, used in the PQ 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 PQ 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 PQ method works poorly, lending more confidence to the estimated damping ratios.The linear and quadratic damping characterized by P and Q obtained using the PQ 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.

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 PQ analysis.The pitch-damping coefficients estimated using the PQ 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.

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 xmomentum 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 40 does not result in a meaningful further reduction in the residual.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.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. (a)

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 ℎ.The 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 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.

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  as the numerical solution of  from the simulation with the  th grid size and time step.The hypothetical numerical solution of  with infinite numerical resolution, i.e., with  → 0, is given by  .The dis- 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 ith 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]: where λ i is the refinement ratio associated with the ith 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]: 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]: 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 leastsquares 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 , , and  of the four repeated simulations in Table 11 are listed in Table 12.Large uncertainties are found for , which is proportional to the linear damping coefficient.This is consistent with Figure 12a, which shows a substantial decrease in  as the simulation is refined.Even with the finest simulation, the discretization uncertainty in  can be as high as ±48%.On the other hand, the value of , which is proportional to the quadratic damping coefficient, increases with the refinement (Figure 12b).The uncertainty in  is more moderate at ±19% with the finest simulation.Interestingly, the effect of the decreasing  is approximately balanced by the increase in  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 /600, leading to a temporal resolution even finer than that of the fine case in Table 11.The resulting values of , , 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 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 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].
Energies 2022, 15, x FOR PEER REVIEW 26 of 40 heave decay,  does not show monotonic convergence, and Equation ( 15) fails to provide a valid fit; therefore, the alternative range-based estimate of 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  = 3 was used [40].
(a) (b) 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  is generally very close to zero; therefore, the absolute uncertainties are shown instead of percentage uncertainties.A slightly negative linear damping, , is obtained with all but the finest CFD simulation.This is likely just a consequence of convergence error because the uncertainty ranges of  extend above zero in all cases.In contrast, the quadratic damping  is always positive.Similar to surge decay,  and  again show opposite trends as the simulation is refined.The contribution from the increasing  with finer numerical resolution is mostly canceled by that from the decreasing , 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.
Table 13.Discretization uncertainties of the heave-damping coefficients and the equivalent linear damping ratio.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.
Table 13.Discretization uncertainties of the heave-damping coefficients and the equivalent linear damping ratio.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 Tables 12 and 13 are generally an order of magnitude greater than the estimated iterative uncertainty in Tables 9 and 10 for the baseline 20 iterations per time step.We therefore neglect the contributions from iterative convergence and retain only the discretization uncertainty.

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 freedecay 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.Table 14.Sensitivity of the damping coefficients and equivalent linear damping ratio to the choice of turbulence models.

Surge Damping
Heave Damping 1 SA-DES = Spalart-Allmaras detached-eddy simulation. 2 SST = shear stress transport. 3The maximum difference in P is not given as percentages because P is very close to zero for heave decay.

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.

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.

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,  , can lead to a slight decrease in the estimated linear damping (from  = 0.066 to 0.045) and an increase in the estimated quadratic damping (from  ⋅  = 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  ) using the quadratic fit of Equation ( 12) instead of the linear fit with Equation ( 9) with the term associated with  removed, in other words, with  forced to zero.This analysis results in a reduced estimate of the linear damping of  = 0.055 and an increased quadratic damping of  ⋅  = 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  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  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

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.
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  .Furthermore, the  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  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  using Equation ( 9) is lower than what is shown, reducing the gap between the experimental results and CFD predictions.Overall,  and  ⋅  are of comparable magnitude, meaning both the linear and quadratic damping are important for the present surge free-decay load case.The equivalent linear damping ratio, , is calculated using Equation (13).The amplitude factor  used to normalize  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.
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  = 0.2 ( + ), where , , 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 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 Sections 6.1 and 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].

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.
sistent, 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 ratio, , is calculated using Equation ( 13).The amplitude factor  used to normalize  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.
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  = 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.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.

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 PQ 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.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 surgedamping coefficients with the OC6-DeepCwind floater.

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  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.Comparison of the linear and quadratic damping in pitch, characterized by the values of  and , respectively, and the equivalent linear damping ratio, , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH).The  and  values are calculated from the motion time series using the  analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients,  and  , respectively, following Equation (10).The equivalent linear damping ratio, , is calculated using Equation (13).The amplitude factor  used to normalize  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.

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  analysis and adopted as metrics for validation.The present investigation has two unique features: 1.A constant Coulomb-friction-type mechanical damping force was identified from the experimentally measured surge motion using a new modified  analysis that also

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 PQ analysis and adopted as metrics for validation.The present investigation has two unique features: 1.
A constant Coulomb-friction-type mechanical damping force was identified from the experimentally measured surge motion using a new modified PQ 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 PQ method.2.
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: 1.
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.

2.
The grid resolution is the leading contributor to the discretization error, whereas the temporal discretization error is secondary, at least for surge decay.

3.
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.

4.
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.

5.
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: 1.
With all three load cases investigated, wave radiation plays a negligible role in the overall motion damping because of the low natural frequencies.

2.
The damping in the heave and pitch directions is predominantly quadratic, most likely coming from the drag force on the heave plates.3.
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: 1.
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.

2.
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.

3.
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.4.
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/10 5 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.

Figure 1 .
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 .
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 2 .
Figure 2. Numerical domain for CFD simulations.The three taut-spring mooring lines attached t the heave plates are also shown.

Figure 2 .
Figure 2. Numerical domain for CFD simulations.The three taut-spring mooring lines attached to the heave plates are also shown.

Figure 3 .
Figure 3. Vertical section of the baseline mesh.

Figure 4 .
The NREL CFD solutions shown are obtained with the baseline numerical configuration described in Section 3.

Figure 5 .
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 Coulombfriction-like damping force.

Figure 5 .
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 6 .
Figure 6.Comparison of the MARIN CFD solutions for surge decay (LC 4.2) with and without  = 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 .
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 7 .
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 8 .
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 .
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 2022 , 40 Figure 9 .
Figure 9. Unnormalized root-mean-square (RMS) residual of the -moment equation at the end of each time step.The time series of the surge and heave free-decay motions obtained with different numbers of iterations shown in Figure10corroborate the behavior of the residual shown in Figure9.In Figure10a, 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 Table9, 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 Table10.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.

Figure 9 . 40 Figure 9 .
Figure 9. Unnormalized root-mean-square (RMS) residual of the x-moment equation at the end of each time step.The time series of the surge and heave free-decay motions obtained with different numbers of iterations shown in Figure10corroborate the behavior of the residual shown in Figure9.In Figure10a, 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.

Figure 10 .
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 2022 ,Figure 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.

Figure 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.

Figure 12 .
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 .
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 13 .
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 .
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 14 .
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 Table11.

Figure 14 .
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 Table11.

Figure 15 .
Figure 15.Comparison of the linear and quadratic damping in surge, characterized by the values of  and , respectively, and the equivalent linear damping ratio, , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH).The  and  values are calculated from the motion time series using the  analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients,  and  , respectively, following Equation (10).The equivalent linear damping ratio, , is calculated using Equation(13).The amplitude factor  used to normalize  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 Figure4.The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table11.The uncertainty bands attached to the NREL solutions are the numerical discretization uncertainties estimated in Section 6.2.

Figure 15 .
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 PQ 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 Table11.The uncertainty bands attached to the NREL solutions are the numerical discretization uncertainties estimated in Section 6.2.

Figure 16 .
Figure16.Comparison of the linear and quadratic damping in heave, characterized by the values of  and , respectively, and the equivalent linear damping ratio, , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH).The  and  values are calculated from the motion time series using the  analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients,  and  , respectively, following Equation(10).The equivalent linear damping ratio, , is calculated using Equation (13).The amplitude factor  used to normalize  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 Figure4.The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table11.The uncertainty bands attached to the NREL solutions are the numerical discretization uncertainties estimated in Section 6.2.

Figure 16 .
Figure16.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 PQ 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 Figure4.The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table11.The uncertainty bands attached to the NREL solutions are the numerical discretization uncertainties estimated in Section 6.2.
Figure16.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 PQ 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 Figure4.The four NREL solutions correspond to the four different numerical setups used in the convergence study listed in Table11.The uncertainty bands attached to the NREL solutions are the numerical discretization uncertainties estimated in Section 6.2.

Figure 17 .
Figure 17.Comparison of the linear and quadratic damping in pitch, characterized by the values of  and , respectively, and the equivalent linear damping ratio, , from the experiment (EXP), the CFD simulations, and the potential-flow simulation (TUHH).The  and  values are calculated from the motion time series using the  analysis discussed in Section 5 and are proportional to the linear and quadratic damping coefficients,  and  , respectively, following Equation(10).The equivalent linear damping ratio, , is calculated using Equation(13).The amplitude factor  used to normalize  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 Figure4.

Figure 17 .
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 PQ 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 Figure4.

Table 2 .
Geometry of the mooring lines at equilibrium.

Table 3 .
Initial displacements in all six degrees of freedom (6DoF) from the selected free-decay experiments.

Table 4 .
Initial velocities in all 6DoF from the selected free-decay experiments.

Table 5 .
Baseline domain boundary locations and conditions for the computational fluid dynamics (CFD) simulations.

Table 7 .
Box-shaped 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.

Table 9 .
Variation in surge-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.

Table 11 .
Numerical configurations for convergence study.

Table 9 .
Variation in surge-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.

Table 11 .
Numerical configurations for convergence study.

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.
Discretization uncertainties of P are expressed in terms of actual values instead of percentages because P is very close to zero. 1

Table A3 .
Time step and computational mesh (all dimensions at full scale).