Investigation of Turbulence Modeling for Point-Absorber-Type Wave Energy Converters

Reviewing the literature of CFD-based numerical wave tanks for wave energy applications, it can be observed that different flow conditions and different turbulence models are applied during numerical wave energy converter (WEC) experiments. No single turbulence model can be identified as an `industry standard’ for WEC modeling. The complexity of the flow field around a WEC, together with the strong dependency of turbulence effects on the shape, operational conditions, and external forces, hampers the formulation of such an `industry standard’. Furthermore, the conceptually different flow characteristics (i.e., oscillating, free surface flows), compared to the design cases of most turbulence models (i.e., continuous single-phase flow), can be identified as a source for the potential lack of accuracy of turbulence models for WEC applications. This communication performs a first step towards analyzing the accuracy and necessity of modeling turbulence effects, by means of turbulence models, within CFD-based NWTs for WEC applications. To that end, the influence of turbulence models and, in addition, the influence of the initial turbulence intensity is investigated based on different wave–structure interaction cases considering two separately validated WEC models. The results highlight the complexity of such a `turbulence analysis’ and the study suggests specific future work to get a better understanding of the model requirements for the flow field around WECs.


Introduction
The undeniable increase of the atmospheric temperature and the widespread recognition of human-induced climate change have, in recent years, fueled research and development (R&D) into novel technologies to harness renewable energy sources. Among the available renewable energy sources, ocean wave energy, once economically viable, can make a valuable contribution towards a sustainable, global, energy mix [1]. The economic viability is commonly expressed by a single, monetary value: the levelized cost of energy (LCoE). The LCoE is defined as the ratio between the overall cost of an energy conversion device/array over its lifetime and the total generated electric energy over the lifetime. From a hydrodynamic point of view, to drive down the LCoE, optimization of the hydrodynamic response of the device is required, which can be achieved through optimization of the shape, the operational degrees of freedom (DoFs), or implementation of energy maximizing control systems (EMCSs). For the hydrodynamic optimization and the synthesis and evaluation of EMCSs at low to mid technology readiness levels (TRLs) [2], researchers and developers rely on complementary experimental and numerical analysis, allowing experimentation in relatively controlled environments, compared to open ocean trials.
where u is the fluid velocity vector, ρ is the fluid density, p denotes the pressure, ν is the kinematic viscosity, and f b are body forces, such as gravity.
With Equations (1) and (2), a CFD-based NWT can capture all the relevant hydrodynamic non-linearities. However, numerical testing in a CFD-based NWT also features specific challenges with respect to the model requirements for WEC analysis. Figure 1 schematically depicts the main requirements of a CFD-based NWT for WEC applications. Specifically, these requirements are: • The capturing (and monitoring) of the multi-phase problem [7], • The determination of a correct discretization of the domain in space and time, with a focus on the spatial discretization in the free surface interface region [8], • The generation and absorption of waves [9], • The accommodation of body motion, induced by hydrodynamic and external forces, through dynamic mesh motion methods [10], • The modeling of turbulence.
This communication will focus on the investigation of the modeling of turbulence with a CFD-based NWT for WEC applications.

Wave absorption Wave Generation
Body motion Turbulence Multi-phase Discretisation Figure 1. Visualization of the requirements for a CFD-based NWT for WEC analysis. (adapted from [9]).

Turbulence Modeling
Considering viscous, as opposed to inviscid flows, the characteristics of the observed flow may show significant differences based on the prevailing conditions: laminar or turbulent.
In laminar flow, infinitesimal fluid layers move smoothly past each other and no mixing occurs between the layers. In such orderly flows, small perturbations in the flow are damped out and viscous forces are predominant over inertial forces. On the contrary, turbulent flows can be characterized as somewhat random and, following Kundu et al. [11], feature: fluctuations of the field quantities; rapid diffusion due to mixing; three-dimensional (3D) turbulent structures due to vorticity, i.e., eddies; dissipation of energy and vorticity; non-linearities. If these turbulence characteristics must be considered within the numerical model, but simulations should be performed at reasonable computational cost, turbulence models, based on the Reynolds averaged Navier-Stokes (RANS) equations, can be employed (Please note that more detailed but more costly approaches to consider turbulence in numerical models, such as direct numerical simulation (DNS) or large eddy simulation (LES) are omitted here due to their relative irrelevance for WEC modeling [12].). The RANS equations can be deduced by applying Reynolds decomposition, i.e., splitting quantities into a mean and a temporally fluctuating component, to all relevant flow quantities in the momentum Equation (2).
In the RANS momentum equation, special attention must be paid to additional terms resulting from the Reynold decomposition of the fluid velocity. In the light of the momentum exchange within fluid layers due to turbulence, these terms can be interpreted as Reynolds stresses. The occurrence of additional terms in the RANS momentum equations requires additional equations to achieve closure of the system of equations. These additional equations can be provided in the form of turbulence models.
The two-equation (standard) k − ε turbulence model, developed by Launder et al. [13], introduces two new transport equations for the kinetic turbulent energy, k t , and the dissipation rate of turbulent kinetic energy, ε t : where u , v , and w denote the variance of the velocity fluctuations, stemming from Reynolds decomposition. S ij denotes the fluctuation of the rate of strain in suffix notation. i or j = 1 denotes the x-direction, i or j = 2 denotes the y-direction, and i or j = 3 denotes the z-direction. The transport equations for k t and ε t , omitted here for brevity, feature terms for the diffusive transport of k t and ε t . The diffusive transport is connected to the eddy viscosity, which follows: where C µ is a model constant. The definition of µ t can be derived from dimensional analysis and the definition of two turbulent scales for the velocity, Υ , and the length, , where With the eddy viscosity, the Reynolds stresses are computed based on the Boussinesq hypothesis [14], following: Please note that the term 2 3 ρ k t δ ij is included to ensure physically correct results. A derivation of this term can be found in ( [15], Chapter 3). The transport equations for k t and ε t include additional model constants which are found by empirical data fitting.
Special care must be taken with the initial and boundary conditions for the standard k − ε model and, in fact, any turbulence model. To set the initial conditions for k t and ε t , the following definitions are presented in [15]: where U ref denotes a reference velocity and T i is the turbulence intensity, defined as the ratio between the fluctuation and the mean of the flow velocity. L is the characteristic length scale of the flow.
The turbulence length scale in the standard k − ε turbulence model is a function of the rate of dissipation of turbulent kinetic energy (see Equation (7)). Another definition for can be found as: with The turbulence eddy viscosity then follows as: which is the definition used in the k − ω turbulence model [19][20][21]. Transport equations for k t and ω t can be formulated and a set of tunable constants are required. Compared to the standard k − ε turbulence model, benefits for the k − ω model are found in the treatment of near wall boundary layer flows, but accuracy issues may occur in the free stream. Thus, Menter [22][23][24][25] proposed the k − ω shear stress transport (SST) model, blending the standard k − ε and k − ω models, using a blending function F C (F C = 0 at the wall; F C = 1 in the far field) for a smooth transition between the wall and the edge of the boundary layer. The blending results in accurate modeling for both far wall and near wall regions.

Wall Treatment
When employing turbulence models, care must be taken with the boundary conditions in the vicinity of solid walls. Due to the steep gradient of the velocity when moving towards the wall, the flow characteristics in the near wall region are significantly different from the characteristics in the far wall region. Four regions of the flow, in the vicinity of the wall, can be defined:

•
Outer layer: Inertial forces largely overweight viscous forces • Log-law layer: Viscous force are smaller than internal forces, but are non-negligible • Buffer layer: Viscous and inertial forces are approximately equal • Viscous sublayer: Viscous forces outweigh inertial forces Although the flow in the outer layer does not require specific treatment, the log-law, buffer, and viscous (sub)layers deserve special attention. The location of the different layers can be defined by the dimensionless distance to the wall, y + , following: where u τ denotes the friction velocity, following u τ = √ τ w/ρ, with the wall shear stress, τ w . In the viscous sublayer (y + < 5), viscous forces are dominant, implying that the mean velocity,ū, is independent of the free stream properties and can be expressed as a function of the wall shear stress τ w , the distance from the wall y, and the dynamic viscosity µ.
In the log-law layer (30 < y + < 300), viscous and inertial forces affect the flow and, following [26], the non-dimensional velocity can be expressed as: where κ is the von Karman constant (κ = 0.41) and E is a model constant (E = 9.8).
In the buffer layer (5 < y + < 30), a transition from the linear to the logarithmic relation between y + and u + takes place; however, no explicit formulations are available.
In the context of turbulence modeling, the analytical expressions of the non-dimensional velocity in the viscous sub-and log-law layer are important, enabling the reduction of the required spatial discretization at solid walls. For high Reynolds numbers, the viscous sublayer is very thin and can be neglected. Therefore, discretization of this layer can also be avoided. As a result, and providing that the location of the first cell at the wall falls within the range 30 < y + < 300, validity of the log-law (see Equation (15)) can be assumed and wall functions can be defined for the turbulence quantities. At low Reynolds numbers, the viscous sublayer cannot be neglected and the log-law cannot be applied all the way 'down to the wall'. Thus, wall functions for the turbulence quantities also loose validity, ultimately requiring discretization of the viscous sublayer. For that, y + ≈ 1 at the wall must be ensured, requiring relatively small grid sizes.

Related Studies
The k − ε model is often employed when simulating oscillating water columns (OWCs), with applications ranging from performance analysis [27][28][29][30] and (structural) optimization [31,32] to system identification [33]. Justification for the choice of the turbulence model is, however, omitted in all studies.
Stansby et al. [34,35] study the drag effects on the point-absorber type M4 WEC and compare the standard k − ε, k − ω SST, and the adopted k − ε V2F [36] models with experimental data and laminar flow simulation. The authors concluded that the k − ω SST and the k − ε V2F models give similar results for the considered free decay tests. Ultimately, in [34], the the k − ε V2F model is employed for subsequent studies, while in [35] the k − ω SST is used for the remainder of the study.
Wei et al. [37] present comparative results considering the standard k − ε, the RNG k − ε, the realizable k − , and the k − ω SST turbulence model for the analysis of viscous effects for an oscillating wave surge converter (OWSC). The negligible differences between the results, together with the reduced computational costs for the standard k − ε model, lead to its implementation in the remainder of [37].
Scarpetta et al. [38] provide a rather qualitative comparison between turbulent (k − ω turbulence model) and laminar simulations, concluding that the necessity of modeling turbulence effects must be assessed on a case-by-case basis. Further references for the k − ω turbulence model can be found in [38][39][40][41][42][43], justification for the choice of the specific turbulence model is generally omitted.
Nguyen et al. [44] perform an assessment of different RANS turbulence models (standard k − ε, RNG k − ε, standard k − ω) based on a bottom hinged OWSC, using experimental data are as a reference. Comparing, for example the rotational angle of the flap, it is found that the standard k − ω turbulence model delivers the best agreement with the experimental data. However, given the close agreement of all numerical results, and the good agreement between laminar simulation and the experimental data, the authors' conclusion is somewhat diluted.
The k − ω SST turbulence model can be identified as one of the most applied turbulence models in the WEC literature. Elhanafi et al. [45] compare the k − ε, the realizable k − ε, and the k − ω SST turbulence models when studying the 2D flow field around a fixed OWC. The two k − ε models produce virtually identical results. In contrast, the comparison between the k − ε model and the k − ω SST model reveal a closer agreement with experimental reference data, for the k − ω SST model.
Van Rij et al. [46] justify the use of the k − ω SST turbulence model in their study on the design loads of a WEC device by identifying it as a "good compromise of computational stability, cost, and accuracy between the simpler RANS turbulence models and the more complex LES".
For the simulation of point-absorber devices, Devolder et al. [47][48][49] employ their proposed buoyancy-modified k − ω SST model [50,51]. By including density in the transport equations, and adding a buoyancy term in the turbulent kinetic energy equation, the modified turbulence model tackles the problem of excessive wave damping.
A summary of the complete literature review in [12] with respect to the modeled flow conditions is shown in Figure 2.

Objectives
Reviewing the existing literature of CFD-based NWTs for WEC experiments, a general uncertainty regarding turbulence modeling (and the required wall treatment) during the simulation of WECs can be noticed. The complexity of the flow field around a WEC, together with the strong dependency of turbulence effects on the shape, operational conditions, and external forces fuels this uncertainty. Furthermore, the conceptually different flow characteristics (i.e., oscillating, free surface, flows), compared to the design cases of most turbulence models (i.e., continuous single-phase flow), can be identified as a source for a potential lack of accuracy of turbulence models for WEC applications.
To get a better understanding of the modeling requirement for, and capabilities of, industry standard turbulence models, a thorough analysis of the available turbulence models for WEC modeling in realistic conditions is desired. To that end, this study analyses the influence of turbulence models on the free surface elevation, the wave excitation forces, and the WEC dynamics, based on inviscid, viscous laminar, and viscous turbulent simulations. Results are compared for four different test cases: wave-only tests, wave excitation force test, and wave-induced motion tests with varying motion constraints. Furthermore, for the turbulent simulations, the effect of the initial turbulence intensity T i is investi-gated by considering three different values, i.e., T i = 1%, 5%, 10%, reflecting cases of low, mid, and high turbulence [52]. The study is based on two representative WEC structures, for which the numerical model (with a laminar flow assumption) is validated in [53,54].

Outline of the Paper
The remainder of the paper is organized as follows: The case study, including a description of the WEC structures, the input waves, and the test cases is presented in Section 2. Section 3 presents the numerical wave tank with a focus on the requirements of the turbulence models. Following on, Section 4 presents the results of the test cases, based on which conclusions are drawn in Section 5.

Case Study
This section presents the case study used to investigate the influence of turbulence modeling on the WEC dynamics. First, the WEC structures are introduced, followed by a description of the input waves.

WEC Structures
Two different WEC structures, W1 and W2, are considered in the following. The structures are based on the systems employed in the blind test series 2 and 3 of the collaborative code project wave structure interaction (CCP-WSI) [55,56] and resemble axisymmetric, cylindrical point-absorber-type devices. The relevant physical properties are shown in Figure 3. The vertical location of the CoMs, as well as the values of the inertia Ixx, Iyy, and Izz, are listed in Table 1. The taut-mooring of W1 and W2 is implemented with a linear spring (k mooring = 67 N m −1 ), connecting the structures to the tank floor. For this study, device motion is constrained in the numerical model to two separate scenarios: 3-DoF (i.e., heave, surge, pitch) and single DoF (i.e., heave) motion. Furthermore, a PTO system, as in [10], with a PTO force f u is implemented for the single DoF case.

Input Waves
An irregular JONSWAP wave train, with a peak period of T p = 1.94 s and a significant wave height of H s = 0.10 m, is considered in this study. The particular wave train represents realistic conditions at the AMETS test site in Bellmullet, Co. Mayo, off the West Coast of Ireland [57] at a scale of 1/30th. Figure 4 shows a time trace of the recorded free surface elevation, recorded at the intended WEC location during a preliminary wave-only simulation.

Test Cases
In this study, four different test cases are considered, following an incremental approach: (1) waves-only; (2) wave excitation (static device); (3) wave-induced motion (single DoF); (4) wave-induced motion . Throughout the test cases, simulations are performed assuming inviscid (Euler) and viscous (laminar and turbulent) flow conditions. For the turbulent simulations, the k − ω SST and the k − ε turbulence models are employed with varying turbulence intensities of T i = 1%, 5%, 10%.

Waves-Only
To investigate the effect of the inclusion of viscosity (i.e., Euler versus laminar simulation) and the use of turbulence models (i.e., wave damping, as shown in [50]), wave-only test cases are performed and results are compared by means of the normalized root mean square, nRMS, following Here, RMS(Π) is the RMS of the free surface elevation from the inviscid or viscous turbulent simulation and RMS Π is the RMS of the free surface elevation from the viscous laminar simulation. Throughout this study, the results from the viscous laminar simulations are the reference, since the numerical model has been successfully validated in [53,54], with the assumption of laminar flow.

Wave Excitation
During the wave excitation force test, the device is held fixed in its equilibrium position, while exposed to the irregular wave train depicted in Figure 4. For the comparison between viscous laminar, viscous turbulent, and inviscid simulations, the wave excitation force is post-processed and compared.

Wave-Induced Motion: Heave-Only
Wave-induced motion tests are considered by introducing non-fixed WEC structures, constrained to move in the heave DoF only. For this particular case, both controlled and uncontrolled devices are considered, investigating the influence of WEC control. In the controlled case, the energy maximizing optimal controller is synthesized in an outputfeedback form [58]. Here, both heave displacement and velocity of the device are used as measurable feedback variables. The optimal control law f u , written in a parametric form, follows with H = [k u b u ], and z andż denoting the heave displacement and velocity of the device, respectively. In particular, the optimal gain H is computed as H W1 = [−1318 N m −1 64 N s m −1 ] and H W2 = [−1385 N m −1 33 N s m −1 ] (see also [59] for more details on the WEC control). For the comparison between viscous laminar, viscous turbulent, and inviscid simulations, the WEC's heave motion is post-processed and compared.

Wave-Induced Motion: 3-DoF
Finally, wave-induced motion of structures, constrained to move in the heave, surge, and pitch DoF, are considered. For this test case, only the uncontrolled case is considered (Please note that no EMCS is available for the 3-DoF model due to a lacking control design model for this case). For the comparison between viscous laminar, viscous turbulent, and inviscid simulations, the WEC's heave, surge, and pitch motion is post-processed and compared.

Numerical Wave Tank
For this study, the numerical wave tank in [10] is adapted (see Figure 5). All simulations are performed with the open source CFD toolbox OpenFOAM, in particular version 1812 of the ESI branch [60]. Waves are generated and absorbed with the IHFOAM [61] toolbox, which is readily implemented in OpenFOAM 1812. The body motion is accommodated using the mesh morphing method (see [10] for more details).

Computational Domain
The overall domain and spatial discretization layout is adapted from the setup in [10] (see Figure 5). However, when using turbulence modeling, the discretization around the WEC is important to capture the boundary layer by means of wall functions. As stated in, e.g., [62], the condition on the y + value (i.e., 30 < y + < 300) poses some challenges for oscillating flows.
During the setup of the numerical domain, preliminary (wave excitation force) tests have been conducted to monitor the achieved y + value. During these preliminary tests, problems were encountered indicating that larger first-layer thicknesses would be required to obey the 30 < y + < 300 criterion; however, the background grid size (leading to converged solutions, see [10]) does not allow arbitrarily large first-layer thicknesses. Thus, a trade-off must be accepted with a first-layer thickness of 5 × 10 −3 m (see Figure 6), leading to y + values of 30 or lower. These preliminary studies highlight the complexity of achieving an acceptable grid layout, complying with the condition on y + .
Since the achieved spatial discretization does not strictly obey the aforementioned y + condition, a preliminary sensitivity analysis on the use of wall functions is conducted based on the single DoF, heave-only, test case. Simulations with the k − ω SST and k − turbulence model and T i = 5% are conducted using no wall functions (i.e., under-resolving the boundary layer), as well as standard and low-Re wall functions. Table 2 lists the nRMS between the different viscous and inviscid flow simulations for a controlled, heave-only, test case.
The results in Table 2 clearly indicate that the nRMS values are increased (by approx. an order of magnitude) when omitting the wall function on an under-resolved grid. Given the successful validation of the numerical model under the assumption of laminar flow [53,54], it can be assumed that these large nRMS values indicate unrealistic results. It should be pointed out, however, that no direct validation of the turbulence model can be conducted here, due to the lack of appropriate, high-fidelity, reference data. Finally, for all subsequent simulations, standard wall functions are used.

Results and Discussion
This section presents the results of the four different test cases outlined in Section 2.3.

Wave-Only
For a qualitative comparison, Figure 7 shows the time traces of the free surface elevation, measured at the intended device position, from the inviscid and viscous (laminar and turbulent) simulations. Visually, a relatively small deviation between the results can be observed. It can be seen in the zoom box in Figure 7 that by considering turbulence modeling and with increasing turbulence intensity, the numerical wave damping increases, which is consistent with findings in, e.g., [50,51].
The quantitative analysis based on the nRMS is shown in Table 3, listing the results including the two turbulence models and the different turbulence intensities. For the inviscid simulation, negligible differences of << 0.1% are calculated, indicating marginal larger free surface elevation amplitudes, compared to the laminar case. Similarly, for the k − ω SST turbulence model with T i = 1%, the deviation between the viscous turbulent and viscous laminar simulation is marginal; however, the negative sign indicates wave damping due to the inclusion of turbulence modeling. With increasing turbulence intensity, the deviation increases. For T i = 5% and 10% the two turbulence models show relatively similar deviations, with maximum values of the order of −6%. Given that any physical dissipation due to the mere presence of fluid viscosity would also be captured in the viscous laminar case, together with the findings in [50,51], the behavior of the turbulent test cases indicates that turbulence modeling in the considered wave-only case induces artificial numerical dissipation, which should be avoided.  Table 3. nRMS between the different viscous and inviscid flow simulations for the wave-only test case.

Euler
Turbulent In Figure 8a-d, the nRMS for the inviscid simulations indicates that the heave and surge excitation forces are slightly larger compared with the viscous laminar simulation. For all viscous turbulent simulations, the nRMS suggests lower excitation forces, compared with the laminar case, which is expected in the light of the result in Section 4.1.

Wave Excitation
Furthermore, the results indicate similar trends in the nRMS for W1 and W2. However, for both structures, the nRMS in the heave force shows maximum values > −4%, the nRMS in the surge force shows maximum values > −10%, which are similar values to the nRMS in the free surface elevation. Increased viscous effects in the surge DoF, compared to the other DoFs, has also been highlighted in [63]. Comparing the results for W1 and W2, it can be seen that the (significant) difference in the shape of the two structures is not reflected in the calculated nRMS values, which are of a similar magnitude for W1 and W2.
For both heave and surge forces, the nRMS increases with increasing turbulence intensity T i . Comparing the results for the k − ω SST (square markers in Figure 8) and the k − turbulence model (diamond markers in Figure 8), larger nRMS values can be observed for the k − ω SST model with T i = 5% and 10%.
It is interesting to observe that particularly for the heave forces, the nRMS is smaller than the nRMS of the surface elevation. In Figure 8 this is highlighted by the fact that the data points for the viscous turbulent simulation fall above the gray line, indicating a linear correlation. These results show that a reduction of the excitation force due to wave damping (shown by the nRMS η) is, to some extend, counteracted by forces on the structure induced by the inclusion of turbulence modeling. This effect is enhanced when considering higher turbulence intensities.

Wave-Induced Motion-Heave-Only
For the case of wave-induced motion and structures constrained to the heave DoF, Figure 9a,b show the nRMS of the heave displacement for W1 and W2 over the nRMS of the free surface elevation, respectively. At first sight, it can be seen that the heave motion nRMS follows the trend of the heave excitation force nRMS. Furthermore, both WEC structures show similar results. Relatively small maximum nRMS values of > −3% and > −2% can be found for W1 and W2 with the k − ω SST turbulence model and T i = 10%, respectively. To investigate the effect of WEC control and, thus, exaggerated device motion on the influence of turbulence modeling, simulations of an uncontrolled device, constrained to the heave DoF, are performed. Please note that for brevity, only cases with T i = 5% are considered.  Based on the time traces in Figure 10 a,b, Figure 11a,b show the nRMS of the heave displacement over the nRMS of the free surface elevation for W1 and W2, respectively, for both the controlled and uncontrolled cases. The results indicate an increased nRMS for the uncontrolled WECs, compared to the controlled cases, most significantly for W1. These results are somewhat unexpected, since the enhanced device motion in the controlled case is assumed to 'trigger' turbulence effects. However, due to the coupling between the effect of turbulence modeling on the device motion and the free surface elevation, an exact determination of the differences between the viscous laminar and turbulent simulations is not straightforward and requires additional analysis (see Section 5).

Wave-Induced Motion-3-DoF
Finally, wave-induced motion of the uncontrolled structures W1 and W2, constrained to move in the heave, surge, and pitch DoFs only, is considered. Figure 12a-f show the nRMS of the heave, surge, and pitch displacements over the nRMS of the free surface elevation, for both W1 and W2. The corresponding time traces of the device motion are plotted in Figure 13a-f. Please note that again for brevity, only results of the k − ω SST and the k − turbulence model with T i = 5% are considered in the following. For the heave DoF (Figure 12a,b), similar values of the nRMS are found for W1 and W2 as for the uncontrolled heave-only test case (see Figure 11), where overall larger nRMS values are found for structure W1. Similar results are found for the pitch DoF (Figure 12e,f). Again, slightly larger nRMS values are visible for W1, compared to W2, and the nRMS for the k − ω SST turbulence model is larger compared to the k − turbulence model. Overall, for both the heave and pitch displacement, relatively small maximum nRMS values of −4% are found. More significant differences can be seen for the surge DoF (Figure 12c,d). For both W1 and W2, the nRMS of the device motion is significantly larger, compared to the displacement in the heave or pitch DoF, with maximum values of up to −16%. The results are consistent with the results for the heave and pitch DoF, such that the k − ω SST turbulence model shows larger deviations compared to the k − turbulence model. These observations are underlined by the time traces shown in Figure 13c,d.
To some extent, the results for the surge displacement of the non-fixed WEC structures W1 and W2 are consistent with the findings for the wave excitation force test (see Figure 8b,d), where the wave excitation forces in the surge DoF are also notably larger than for the heave DoF.

Conclusions
From the presented results, it can be stated that the use of turbulence models has a significant influence on the modeled free surface elevation, the wave excitation forces, and the device dynamics. The presented results suggest that the degree of influence depends on the specific turbulence model, as well as the turbulence intensity, employed in the NWT. However, for all tested WSI cases, the excitation forces or device motion in the heave DoF show relatively small nRMS values of > −4%, the excitation forces and device motion in the surge DoF shows larger nRMS values of up to −16%. These results are, to some extend, consistent with the findings in [63].
However, more in-depth analysis is required to get a better understanding of the effect and accuracy of turbulence modeling for WEC applications. In future work cases with regular wave excitation should be considered and the recently proposed, modified, turbulence models by Devolder et al. [50,51] should be considered to separate the influence of artificial wave damping and the behavior of the WSI, when turbulence modeling is included in the CFD-based NWT. Furthermore, dedicated, high quality, validation data (including e.g., pressure measurements or particle image velocimetry data) should be acquired to validate the numerical models. Until then, there remains an uncertainty regarding the inclusion and capabilities of 'industry standard' turbulence models for WEC application. In the light of the successful validation of several different WEC devices for a range of different test cases assuming laminar flow conditions, such assumptions for the flow conditions may be suitable, avoiding uncertainty in the selection of a 'correct' turbulence model.