**1. Introduction**

The energy demand will increase by more than 1.0% per year up to 2040, increasing gases emission rates [1]. Moreover, the costs of commodities for energy generation from fossil fuels have increased significantly, leading to economic difficulties, risks associated with energy security, and geopolitical conflicts around the world [2]. Considering this scenario, there is a growing search for a better comprehension of the development of technologies and the use of devices and economic impacts of different renewable sources of energy such as wind, solar, geothermal, and ocean energy [2–11].

One of the important sources of renewable energy with high potential, but not frequently explored worldwide, is the conversion of ocean energy into electricity [12–15], i.e., wave energy conversion. Despite several signs of progress in technological development, there is no dominant main operational principle. Several devices have been proposed

**Citation:** dos Santos, A.L.; Fragassa, C.; Santos, A.L.G.; Vieira, R.S.; Rocha, L.A.O.; Conde, J.M.P.; Isoldi, L.A.; dos Santos, E.D. Development of a Computational Model for Investigation of and Oscillating Water Column Device with a Savonius Turbine. *J. Mar. Sci. Eng.* **2022**, *10*, 79. https://doi.org/10.3390/jmse 10010079

Academic Editors: Davide Tumino and Antonio Mancuso

Received: 29 November 2021 Accepted: 5 January 2022 Published: 7 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and investigated based on various ways to convert wave energy such as point absorbers, attenuators, oscillating surge converters, overtopping, submerged plates, and oscillating water columns (OWCs) [12–14]. The efficiency and survivability of the wave energy converters (WECs) are important issues to make them more competitive and viable [16]. In this context, the OWC device has advantages as its simplicity and maintenance, e.g., the moving parts are located outside of the water, increasing the lifetime material of the power take-off (PTO) system, and the structures of buildings are robust [12–14]. Therefore, several studies and prototypes using the OWC as the main operational principle have been developed around the world: Sakata–Japan (60 kW), Mutriku–Spain (296 kW), Pico–Portugal (400 kW), Tofteshallen–Norway (500 kW), Islay island–Scotland (500 kW), and Lewis island–Scotland (4.0 MW) [17–22].

Important experimental works have sought to improve the comprehension of the fluid dynamic behavior of water/air flow in the OWC device and investigate the influence of some parameters over its performance. For instance, experiments in the laboratory and large-scale domains analyzed the influence of the inclination of the frontal wall, entrance areas of the OWC chamber, and water depth on the device efficiency, reflection, and loading of an OWC for different wave conditions [23–25]. Recently, the experimental progress extended to understand the hydrodynamic of the fluid flow into dual chambers OWC [26,27].

The numerical simulation of OWC devices has also been worth investigating. Several works have been performed since the development of computational models to represent the main operating principle of the device and the investigation of several parameters regarding the performance. For the former purpose, the representation of fluid flow in a laboratory and large-scale devices has been done without considering the effect of a turbine, as in Maciel et al. [28]. Other studies have considered it using orifice plates, plate-baffle, obstacles, or actuator disks to simulate the head loss caused by the turbine over the airflow in the hydropneumatic chamber and air duct of the device [29–33]. Recently, an interesting approach employed the numerical simulation of water and air in the chamber and considered the effects of Wells and impulse turbines by means of analytical thermodynamic models [34,35]. All models mentioned above have been used to obtain recommendations about parameters such as the depth and inclination of the frontal wall, height and length of the chamber, the diameter of the turbine, ramp placed in the seabed below the OWC chamber, and, recently, the design of multiple coupled chambers regarding device performance for different scales and wave conditions [24,32,35–40]. It is also worth mentioning the efforts made to represent the sea state in a more trustworthy form. In this field, some studies modeled the irregular waves using a wave spectrum such as JONSWAP, and others obtained the sea state from geophysical models such as TOMAWAC and used this as an input for the modeling of a channel with the device to be investigated [41–43].

Some significant advancements regarding the PTO of OWC have also been reported in the literature. For instance, Britto-Melo et al. [44] numerically investigated the influence of the aerodynamic parameters of the Wells turbine and the influence of guide vanes and the bypass valve on the pressure drop, torque, and the overall performance. In this work, the conversion of pneumatic energy into electrical energy was estimated with a computational model based on the results extrapolated from aerodynamic tests on a scale model and empirical approximations for the generator losses. Recently, Rodríguez et al. [45] proposed a computational study to understand the behavior of OutFlow Radial (OFR) turbines in both direct and reverse modes, simulating an axisymmetric domain with the flow between a blade-to-blade arrangement of the turbine. The authors observed that the outer blade angle had poor performance in reverse mode despite the improvement of the global performance due to the rotor efficiency gain in direct mode. However, there are few studies related to the numerical simulation of the OWC device turbine considering the rotational domain with the intruded turbine rotor. Prasad et al. [46] performed similar work in this direction and developed a numerical model with a Savonius turbine immersed in a channel under regular wave flow, simulating a hydrokinetic turbine. The authors also investigated the influence of some parameters such as the submergence and the rotational speed for different blade entry angles regarding the rotor power.

Despite the several above-mentioned contributions, to the authors' knowledge, the development of computational models for the simulation of OWC devices considering the rotational turbine in the air duct is an approach still little explored in the literature. Recently, however, Liu et al. [47] presented the validation of an integrated three-dimensional numerical model considering an axial-flow impulse turbine coupled with an OWC inserted in a numerical wave tank (NWT). The present work aims to perform the first step in this direction. Initially, turbulent air flow over a free turbine inserted in a long and large channel (commonly used to represent the numerical modeling of wind turbines) was simulated to verify/validate the present model (*Case 1*). The effect of the tip speed ratio (*λ*) on the time and spatial averaged power coefficient (*CP*) of a Savonius turbine was compared with numerical and experimental results available in the literature [48,49], investigating the tip speed ratios in the range 0.75 ≤ *λ* ≤ 2.00. An enclosure domain that mimics an OWC device was simulated considering a constant velocity imposed at the inlet of the domain (*Case 2*) and sinusoidal velocity that simulated the alternate flow in an OWC device (*Case 3*). For all cases, the influence of *λ* on the turbine performance and aerodynamic coefficients was investigated.

Moreover, incompressible, two-dimensional, unsteady, and turbulent flows with ReD = 867,000 were considered. URANS (Unsteady Reynolds-Averaged Navier Stokes) modeling was applied to all cases. Time-averaged equations of the conservation of mass, the balance of momentum, and transport equations of the *k*–ω SST model (used in the closure of turbulence) were solved with the finite volume method (FVM) [50–55], using the commercial code Ansys FLUENT 14.5 [56].

#### **2. Mathematical Modeling**

Incompressible, two-dimensional, unsteady, turbulent flows with constant thermophysical properties were considered. It is worth mentioning that one of the characteristics of turbulence is the three-dimensional structure of the flow [57]. Despite this fact, the present simulations represent most of the characteristics of turbulent flows, and the modeling properly predicts the transient behavior of parameters such as drag, lift, moment, and power in the Savonius turbine.

#### *2.1. Description of the Studied Cases*

Figure 1 illustrates the computational domain of *Case 1* used for verification/validation of the present model. The domain consists of a long and large channel with an inserted free Savonius turbine. The dimensions are similar to those investigated in the work of Akwa et al. [49]. In this case, the fluid flow of air is caused by the imposition of a constant velocity (*V*∞) at the inlet (left side surface). On this surface, there is also an imposed turbulent intensity of *IT* = *u*2/*u* = *σu*/*u* = 1.0%, where *u'* is the fluctuation of the velocity field, *u* is the time-averaged velocity field, and *σ<sup>u</sup>* represents the variance of *u'*. At the exit of the channel (right side surface), there is an imposed null gauge pressure (*pg* = 0 atm). At the upper and lower surfaces, there is a free slip and impermeability boundary condition (also called symmetry). In the turbine region, there is an imposed constant angular velocity (*n*) in the gray region named the "Rotational Domain" simulating the effect of wind action over the turbine. In the turbine walls, there is an imposed no-slip and impermeability boundary condition (*u* = *v* = 0 m/s) related to the rotational domain. Figure 1 illustrates the details of the turbine region with the geometric variables used to design it. Table 1 presents the parameters of the fluid flow, thermo-physical properties, dimensions of the computational domain, and turbine variables used here for the four different tip speed ratios of the rotor (*λ* = *nD*/2*V*∞) of *λ* = 0.75, 1.00, 1.25, and 2.00; *n* is the angular velocity of the turbine and *D* is the turbine diameter. For the unsteady analysis, there was a time interval of *tf* = 3.5 s, with the last 1.75 s being analyzed for computation of the drag, lift, momentum, and power coefficients.

**Figure 1.** Schematic representation of the computational domain of turbulent air flow over a free Savonius turbine used for verification/validation of the present numerical model (*Case 1*).

Figure 2 illustrates the computational domain used in *Case 2* and *Case 3*. These cases were defined based on *Case 1*. Therefore, the thermo–physical properties and dimensions of the turbine were the same, presented in Table 1 (*Case 1*). The main differences here are the domain dimensions (*H* = *h* = 6.0 m, *L* = 10.0 m, and *l* = 2.0 m) and the insertion of the Savonius turbine in an enclosure domain, leading to the imposition of a no-slip and impermeability boundary condition (*u* = *v* = 0 m/s) in the device walls. At the lower surface, there was an imposed constant velocity of *V*(*t*) = 1.4 m/s and *IT* = 1.0% for *Case 2*. The magnitude of the imposed velocity was defined in such a way to have a value of 7.0 m/s at the inlet of the air duct (in the contraction from the chamber to the air duct), leading to a Reynolds number in the turbine similar to that reached for *Case 1*. For *Case 3*, the sole difference in comparison with *Case 2* was the imposition of a sinusoidal function that mimicked the oscillating behavior in the OWC chamber:

$$V(t) = \frac{H\_W \pi}{T\_w} \cos\left(\frac{2\pi t}{T\_w}\right) \tag{1}$$

where *Hw* = 0.4 m, *Tw* = 0.875 s, allowing the reproduction of the piston-type movement generated by the regular waves incidence over the OWC. The range of magnitudes was limited by (−1.4 m/s ≤ *V*(*t*) ≤ 1.4 m/s). The use of real configurations of a sea wave, for example, leads to long periods of simulation, requiring a high computational effort. As the purpose of this case is to compare the imposition of the sinusoidal velocity profile with the case with constant velocity to investigate the effect of oscillating flow over the device parameters and performance, the idealized imposed velocity variation was adequate for the desired investigation. For *Case 2* and *Case 3*, the same four magnitudes of the tip speed ratio studied in the verification/validation case (*λ* = *nD*/2*V*1) of *λ* = 0.75, 1.00, 1.25, and 2.00 were investigated. It is worth mentioning that a mean velocity was measured in the air duct and before the turbine (*V*1) for the calculation of *λ*.


**Table 1.** Parameters used in the verification/validation case (*Case 1*).

In the present work, results the influence of *λ* over the drag, lift, moment, and power coefficients (*Cd*, *Cl*, *CT,* and *CP*) are based on [57]:

$$C\_d = \frac{F\_d}{1/2\rho A\_r V^2} \tag{2}$$

$$C\_l = \frac{F\_l}{1/2\rho A\_r V^2} \tag{3}$$

$$C\_T = \frac{T}{1/2\rho A\_r V^2 r} \tag{4}$$

$$\mathcal{C}\_P = \frac{P}{P\_{\text{available}}} = \frac{T}{1/2\rho A\_r V^2 r} \frac{r\mathcal{U}}{V} = \mathcal{C}\_T \lambda \tag{5}$$

where *Fd* is the drag force (N), *Fl* is the lift force (N), *T* is the rotor moment (N·m), *Ar* is the projected area of the Savonius rotor (*Ar* = *D*·*W*), *W* is the depth of the domain in the *z*-direction (m), *V* is the upstream turbine velocity (*V* = *V*<sup>∞</sup> for the *Case 1* and *V* = *V*<sup>1</sup> for *Case 2* and *Case 3*) (m/s), *r* is the radius of the rotor (m), *P* is the turbine power (W), and *Pavailable* is the available power of the wind upstream of the rotor (W).

**Figure 2.** Schematic representation of *Case 2* and *Case 3* that simulate the OWC with a Savonius turbine.

It is worth mentioning that the use of Equation (5) is valid only for the prediction of power coefficients in turbines subjected to open flow conditions (seen in tidal or wind power devices). In OWC devices, the air flow is driven by the pneumatic power, i.e., the pressure drop in the device must also be taken into account. Therefore, for power coefficients predicted for the enclosure domain of *Case 2* and *Case 3*, the *CP* is calculated as:

$$C\_P = \frac{P}{P\_{puncumtic}} = \frac{P}{\left(\Delta p + \frac{\rho V^2}{2}\right)Q} \tag{6}$$

where Δ*p* is the pressure drop between the OWC chamber and the exit of the chimney (Pa), *V* is the air velocity at the OWC turbine-duct (m/s) (*V* = *V*1), and *Q* is the volumetric flow rate of the air (m3/s).

The time-averaged magnitudes of the drag, lift, moment, and power coefficients are obtained as follows:

$$\overline{\mathcal{C}\_{d,l,T,P}} = \frac{1}{t\_f} \int\_0^{t\_f} \mathcal{C}\_{d,l,T,P} dt \tag{7}$$

#### *2.2. Governing Equations of Turbulent Flows*

For all simulations, the modeling of incompressible, two-dimensional, unsteady, and turbulent flows is given by the time-averaged conservation equation of mass and balance of momentum in the *x* and *y* directions and can be written as [50,51]:

$$\frac{\partial \overline{u}}{\partial x} + \frac{\partial \overline{v}}{\partial y} = 0 \tag{8}$$

$$
\rho \left[ \frac{\partial \overline{\boldsymbol{\Pi}}}{\partial t} + \overline{\boldsymbol{\Pi}} \frac{\partial \overline{\boldsymbol{\Pi}}}{\partial \mathbf{x}} + \overline{\boldsymbol{\sigma}} \frac{\partial \overline{\boldsymbol{\Pi}}}{\partial \boldsymbol{y}} \right] = -\frac{\partial \overline{\boldsymbol{\mu}}}{\partial \mathbf{x}} + (\mu + \mu\_l) \left( \frac{\partial^2 \overline{\boldsymbol{\Pi}}}{\partial \mathbf{x}^2} + \frac{\partial^2 \overline{\boldsymbol{\Pi}}}{\partial \mathbf{y}^2} \right) \tag{9}
$$

$$
\rho \left[ \frac{\partial \overline{\boldsymbol{\sigma}}}{\partial t} + \overline{\boldsymbol{\mu}} \frac{\partial \overline{\boldsymbol{\sigma}}}{\partial \mathbf{x}} + \overline{\overline{\boldsymbol{\sigma}}} \frac{\partial \overline{\boldsymbol{\sigma}}}{\partial \boldsymbol{y}} \right] = -\frac{\partial \overline{\boldsymbol{\sigma}}}{\partial \boldsymbol{y}} + (\mu + \mu\_l) \left( \frac{\partial^2 \overline{\boldsymbol{\sigma}}}{\partial \mathbf{x}^2} + \frac{\partial^2 \overline{\boldsymbol{\sigma}}}{\partial \boldsymbol{y}^2} \right) \tag{10}
$$

where *x* and *y* are the spatial coordinates (m), *u* and *v* are the velocity components in the *x* and *y* directions, respectively (m/s), *p* is the pressure (N/m2), *μ* is the dynamic viscosity (kg/m·s), *μ<sup>t</sup>* is the turbulent viscosity (kg/m·s), and the overbar represents the time-averaged operator.

For the *k*–*ω* SST closure model, the turbulent viscosity (*μt*) is [52,53]:

$$
\mu\_t = \frac{\overline{\rho} \alpha\_1 k}{\max(\alpha\_1 \omega, SF\_2)} \tag{11}
$$

The transport equations of turbulent kinetic energy (*k*) and its specific dissipation rate (*ω*) are as follows:

$$\frac{\partial k}{\partial t} + \frac{\partial (\overline{u}\_i k)}{\partial x\_i} = \widetilde{P}\_{\mathbf{k}} - \frac{k^{\frac{3}{2}}}{L\_T} + \frac{\partial}{\partial x\_i} \left[ (\mu + \sigma\_k \mu\_t) \frac{\partial k}{\partial x\_i} \right] \tag{12}$$

$$\frac{\partial \boldsymbol{\omega}}{\partial t} + \frac{\partial (\overline{\boldsymbol{u}}\_{i} \boldsymbol{\omega})}{\partial \mathbf{x}\_{i}} = \frac{\boldsymbol{\alpha}}{\mu\_{t}} \widetilde{\boldsymbol{P}}\_{k} - \beta \boldsymbol{\omega}^{2} + \frac{\partial}{\partial \mathbf{x}\_{i}} \left[ (\mu + \sigma\_{\omega^{i}} \mu\_{t}) \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_{i}} \right] + 2(1 - F\_{1}) \frac{\sigma\_{\omega 2}}{\omega} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{i}} \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_{i}} \tag{13}$$

where *P <sup>k</sup>* is a function that prevents the turbulence generation in stagnation regions, *i* represents the direction of fluid flow (*i* = 1 represents the *x* direction and *i* = 2 represents the *y* direction), *β* = 0.09, *α*<sup>1</sup> = 5/9, *β*<sup>1</sup> = 3/40, *σ<sup>k</sup>* = 0.85, *σ<sup>w</sup>* = 0.5, *σ<sup>2</sup>* = 0.44, *β<sup>2</sup>* = 0.0828, *σk2* = 1, *σw2* = 0.856 are ad hoc constants used in [52]. *F1* and *F2* are blending functions defined as follows:

$$F\_1 = \tanh\left\{ \left\{ \min \left[ \max \left( \frac{k^{1/2}}{\beta^\* \omega y}, \frac{500v}{y^2 \omega} \right), \frac{4\rho v\_{\omega 2} k}{\overline{\mathbb{C}D\_{k\omega} y^2}} \right] \right\}^4 \right\} \tag{14}$$

$$F\_2 = \tanh\left\{ \left[ \max\left( \frac{2k^{1/2}}{\beta^\* \omega y}, \frac{500\upsilon}{y^2 \omega} \right) \right]^2 \right\} \tag{15}$$

In Equation (13), the term *CDk<sup>ω</sup>* is calculated as follows:

$$\mathbb{C}D\_{k\omega} = \max\left(2\rho\sigma\_{\omega 2}\frac{1}{\omega}\frac{\partial k}{\partial \mathbf{x}\_i}\frac{\partial \omega}{\partial \mathbf{x}\_i}, 10^{-10}\right) \tag{16}$$

#### **3. Numerical Modeling**

The solution of the governing equations was performed with the FVM using the commercial package Ansys FluentTM [54–56]. To tackle the advective terms of the balance of momentum and transport of *k* and *ω*, the second-order upwind interpolation function was employed. The algorithm SIMPLE (Semi-Implicit Method for Pressure Linked Equations) was used for pressure–velocity coupling. The simulations were considered to have converged when the residuals for continuity, the balance of momentum, *k* and *ω* transport equations were less than 10−5. Moreover, the maximum number of iterations per time step was 200. Concerning the time advancement, there was an implicit time advancement scheme and a fixed time step of <sup>Δ</sup>*<sup>t</sup>* = 1.75 × <sup>10</sup>−<sup>3</sup> s. All simulations were performed using desktops with six core Intel® Core ™ i7 5820K @ 3.30 GHz processors and 16 Gb of RAM memory. The processing time for the simulation of *t* = 3.5 s of physical time was nearly <sup>20</sup> × <sup>10</sup><sup>3</sup> s.

Concerning spatial discretization, hybrid triangular and rectangular finite volumes were used with the domain in simulations of *Case 1*, *Case 2*, and *Case 3*. Figure 3a–c illustrates the mesh generated with software GMSH [58] for the free Savonius turbine (*Case 1*), the configuration similar to the OWC device (*Case 2* and *Case 3*), and a detail of the mesh in the blades of Savonius turbine, respectively. In detail, it is possible to observe a region around the blades with the refined rectangular mesh. The dimensions of the rectangular volumes were defined as a function of a grid independence study and the parameter for representation of the boundary layer (*y+*), which must be *y+* ≤ 1.0 in the walls; *y+* was defined as follows [51–53]:

$$y^{+} = \frac{y\sqrt{\pi\_w/\rho}}{\nu} = \frac{yu\_\pi}{\nu} \tag{17}$$

where *y* is the normal distance to the wall (m), *τ<sup>w</sup>* is the surface tension in the wall (N/m2), *ν* is the kinematic viscosity (m2/s), and *u<sup>τ</sup>* is the friction velocity (m/s).

For the grid independence study, four different meshes were simulated, and the results for the time-averaged power coefficient for the free Savonius rotor (*Case 1*) with *ReD* = 867,000, *λ* = 1.25, and bucket overlap ratio of *RS* = *s*/*c* = 0.15 are presented in Table 2. The mesh was considered independent when the relative difference between the results of *CP* for two successive grids met the criterion given by:

$$dev\left(\%\right) = 100 \times \left| \frac{\overline{\mathbb{C}\_p}^j - \overline{\mathbb{C}\_p}^{j+1}}{\overline{\mathbb{C}\_p}^j} \right| \le 5.0 \times 10^{-1} \tag{18}$$

where *j* represents the result obtained with the coarser mesh, and *j* + 1 represents the result obtained with the next successive refined mesh. Based on the results of Table 2, the mesh with 369,653 volumes was used in the independent grid. The same parameters used in this mesh were applied to the spatial discretization of the domain of *Case 2* and *Case 3*.

**Figure 3.** Employed mesh for the studied cases: (**a**) turbine in the open channel (*Case 1*), (**b**) turbine in an OWC chamber domain (*Case 2* and *Case 3*), (**c**) detail of grid refinement the surfaces of the blades.

**Table 2.** Grid independence test for a free turbine with *ReD* = 867,000, *λ* = 1.25, and *RS* = 0.15.


#### **4. Results and Discussion**

This section is divided into two parts, the verification/validation of the developed computational model (*Case 1*) and the investigation of the influence of the enclosure model and imposition of the sinusoidal velocity inlet over the aerodynamic and performance coefficients (*Case 2* and *Case 3*).

#### *4.1. Verification/Validation of the Computational Model (Case 1)*

Figure 4 illustrates the instantaneous drag, lift, and moment coefficients as a function of time for *Case 1* with *ReD* = 867,000, *λ* = 2.0, and *RS* = 0.15. For the first instants of time, mainly for *t* ≤ 0.5 s, the coefficients had a strong variation due to the incidence of the fluid flow and the imposition of angular velocity in the rotational domain region. Therefore, the present model must be used to predict power in the turbine when the flow is stabilized, and the first cycles of rotation of the turbine were disregarded for the analysis of coefficients

and power. Here, for the computation of time-averaged parameters, only the results in the range of time 1.75 s ≤ *t* ≤ 3.5 s were used. Despite the complexity of the fluid flow, the results also demonstrated a regular oscillation in the magnitudes of *Cd*, *Cl*, and *CT*, which had similar behavior to that previously obtained in Akwa et al. [49]. Moreover, the crest and cave magnitudes were also similar, showing the generation of regular wakes of vortices behind the turbine.

**Figure 4.** Transient coefficients of the lift (*Cl*), drag (*Cd*), and torque (*CT*) for the free turbine case with *ReD* = 867,000, *Rs* = 0.15, and *λ* = 2.0.

Figure 5 shows the comparison of power coefficients (*CP*) as a function of the tip speed ratio (*λ*) obtained with the present computational model and the numerical predictions of Akwa et al. [49] obtained with other commercial code (Star-CCM+) also based on the FVM and the experimental results of Blackwell et al. [48]. The results predicted with the present method are in close agreement with those previously obtained in the literature, verifying and validating the method used here. Even for *λ* = 1.00 and 2.00, where the highest differences between the present results and the experimental ones were obtained, the deviations were lower than the uncertainty of the experiment. The results indicated that the highest magnitudes of *CP* were reached at the lowest values of *λ*. With the increase in *λ*, the magnitude of *CP* had a slight decrease in the range 1.00 ≤ *λ* ≤ 1.25 and a step decrease for *λ* ≥ 1.25. For *λ* = 2.00, negative magnitudes of *CP* were obtained, indicating that the device supplies energy for the fluid flow and not the contrary.

**Figure 5.** Comparison between the power coefficient (*CP*) as a function of *λ* obtained in the present work and that from the literature [48,49].

Figures 6 and 7 illustrate the velocity and pressure fields, respectively, for *Case 1* with *ReD* = 867,000 and *λ* = 1.25 for five different instants of time: (a) *t* = 0.0 s, (b) *t* = 0.53 s, (c) *t* = 1.05 s, (d) *t* = 1.58 s, and (e) *t* = 2.10 s, which represent different angular positions of the Savonius rotor. For the initial time step, in Figures 6a and 7a, the velocity and pressure fields noticed were generated by the imposition of angular velocity in the rotational domain by the computational model. This behavior does not represent the real condition of fluid flow over the turbine since the turbine should not be in motion before the incidence of the fluid flow. Therefore, the present computational method must be used when the flow is stabilized, which happens a few cycles later at the beginning of the fluid flow (as shown in Figure 4 for monitoring the coefficients). As the time advances, mainly for *t* > 1.0 s, Figures 6c–e and 7c–e show an increase in the pressure field magnitude on the concave side of the advancement blade and a pressure drop as a consequence, being the main reason for the pressure drag force that drives the turbine. The results also show the fluid flowing in the region between the two blades. The results also show the generation of wakes behind the rotor and vortices generated in the tip region of the blade. It is also worth mentioning that the behavior found here is similar to that described in previous literature, e.g., in Prasad et al. [45] and Blackwell et al. [49].

**Figure 6.** Velocity fields around the free turbine (*Case 1*) with ReD = 867,000 and *λ* = 1.25: (**a**) *t* = 0.0 s, (**b**) *t* = 0.53 s, (**c**) *t* = 1.05 s, (**d**) *t* = 1.58 s, and (**e**) *t* = 2.10 s.

**Figure 7.** Pressure fields around the free turbine (*Case 1*) with ReD = 867,000 and *λ* = 1.25: (**a**) *t* = 0.0 s, (**b**) *t* = 0.53 s, (**c**) *t* = 1.05 s, (**d**) *t* = 1.58 s, and (**e**) *t* = 2.10 s.

#### *4.2. The Results of the Savonius Turbine Inserted in an OWC Domain (Case 2 and Case 3)*

For the turbulent air flow in the enclosure domain, the instantaneous aerodynamic coefficients (*Cd*, *Cl*, and *CT*) were obtained as a function of time to better understand the behavior of the turbine in an OWC domain. Figure 8 depicts the coefficients for the flow with ReD = 867,000, *λ* = 2.0, and *RS* = 0.15 for a constant imposed velocity at the inlet (*Case 2*). In general, the results show an average increase in *Cd* (black line) for all instants of time compared to the free turbine simulations (*Case 1*) due to the insertion of the turbine in the enclosure domain. The increase in the drag coefficient (*Cd*) also led to an augmentation of the moment coefficient (*CT*) represented in red color in a similar form reached for the *Cd*, while the lift coefficient (*Cl*) did not suffer important modifications in its mean magnitude. The results also indicate that the transient behavior of the coefficients was strongly modified compared to *Case 1*, showing *Cd* and *Cl* with sharp peaks and smooth troughs, i.e., each cycle did not behave in a sinusoidal form as noticed in Figure 4 for the free turbine configuration. It is worth mentioning that, for other magnitudes of *λ*, similar behavior for the instantaneous coefficients was obtained. Therefore, for the sake of brevity, the instantaneous coefficients for other magnitudes of *λ* are not presented.

**Figure 8.** Instantaneous coefficients of the drag (*Cd*), lift (*Cl*), and moment (*CT*) for the case of OWC with constant imposed velocity at the inlet (*Case 2*) with ReD = 867,000, *λ* = 2.0, and *RS* = 0.15.

As previously mentioned, the air flow in the OWC chamber is subjected to the pistontype oscillatory motion of the water column (hydropneumatic chamber). In order to simulate this effect, the results of instantaneous power coefficient obtained for *Case 3* and *Case 2*, considering the same conditions (ReD = 867,000, *λ* = 2.0, and *RS* = 0.15), are presented in Figure 9. It is important to reinforce that, for prediction of *CP* for *Case 2* and *Case 3*, Equation (6) was used instead of Equation (5), and the time-averaged magnitudes were predicted with Equation (7). The results reveal a strong similarity between the instantaneous *CP* for both cases, indicating that the imposition of sinusoidal velocity at the inlet of the domain did not have a significant influence over the pattern of the *CP* investigated here, which is not intuitively expected since the mean imposed momentum decreased for *Case 3* in comparison with *Case 2*. The results also demonstrate a slight increase in the power coefficients when the sinusoidal velocity was imposed, with a difference of nearly 9.0% on average. Possible explanations for the behavior found here include the synchronization of velocity augmentation with the rotation of the turbine and the imposition of the rotational domain, as well as the acceleration of the fluid caused by the variation of the imposed velocity at the inlet for *Case 3*, which did not happen for *Case 2*. The magnitude of the imposed mean velocity at the inlet of *Case 2* and *Case 3* and acceleration at the inlet for *Case 3* are illustrated in Figure 10 to make this visualization easy. It is worth mentioning that, despite different imposed inlet velocities, the transient fields of velocity and pressure had only slight differences. Future investigations should be performed with other magnitudes of amplitude and periods of imposed velocity to corroborate this hypothesis and the development of other models where the inertia moment of the turbine is taken into account.

**Figure 9.** Effect of imposed velocity at the inlet of the OWC domain over the instantaneous power coefficient in the Savonius turbine, considering ReD = 867,000, *λ* = 2.0, and *RS* = 0.15.

**Figure 10.** Magnitude of the imposed velocity at the inlet of the OWC domain for *Case 2* and *Case 3* and acceleration of the fluid at the inlet of the domain for *Case 3*.

Figures 11 and 12 illustrate the velocity and pressure fields, respectively, for the simulation of *Case 3* with ReD = 867,000, *λ* = 2.0, and *RS* = 0.15 for six different time steps: (a) *t* = 0.0 s, (b) *t* = 0.26 s, (c) *t* = 0.53 s, (d) *t* = 0.79 s, (e) *t* = 1.05 s, and (f) *t* = 1.31 s. As the behaviors of *Case 3* and *Case 2* were almost the same, only the fields for one of the cases are illustrated. As previously mentioned, for the free turbine case (see Figures 6 and 7), the initial time steps should be disregarded due to the artificial imposition of the angular velocity in the turbine region. In spite of that, when the flow was stabilized, the model properly represented the physical problem. In turn, one can note in Figures 11 and 12 an increase in the pressure difference between the concave and convex sides of the advancement blade for *Case 2* and *Case 3* compared to *Case 1*, which explains the augmentation of the pressure drag in comparison with *Case 1*.

Moreover, the pressure difference between the upstream and downstream regions of the turbine was also augmented. This influence over the water oscillation was not investigated, but the results pointed out that this aspect is worth investigation when the turbine is taken into account in the problem in order to avoid the restriction of the water column, mainly when coupling with a wave channel is performed. The results also indicated the variation of the pressure magnitude in the chamber for different instants of time, explaining the non-symmetric differences between the peaks and troughs with respect to the mean magnitudes of the coefficients. The momentum of the fluid flow intensified due to the insertion of the turbine in the enclosure domain. Moreover, the wakes generated in the turbine and secondary vortices in the tip of the blades could not be spread in the spanwise direction of the main flow due to the limitation imposed by the air duct walls. Therefore, as expected, the insertion of the turbine in the enclosure domain affected the fluid dynamic behavior of the flow considerably.

**Figure 11.** Velocity fields for the OWC case with sinusoidal velocity at the inlet for different instants of time for the case with *ReD* = 867,000, *λ* = 2.0, and *RS* = 0.15: (**a**) *t* = 0.0 s, (**b**) *t* = 0.26 s, (**c**) *t* = 0.53 s, (**d**) *t* = 0.79 s, (**e**) *t* = 1.05 s, and (**f**) *t* = 1.31 s.

**Figure 12.** Pressure fields for the OWC case with sinusoidal velocity at the inlet for different instants of time for the case with *ReD* = 867,000, *λ* = 2.0, and *RS* = 0.15: (**a**) *t* = 0.0 s, (**b**) *t* = 0.26 s, (**c**) *t* = 0.53 s, (**d**) *t* = 0.79 s, (**e**) *t* = 1.05 s, and (**f**) *t* = 1.31 s.

To summarize, Figure 13 shows the effect of *λ* on the time-averaged magnitudes of *Cd*, *Cl*, *CT* and *CP* for *Case 2* and *Case 3*, and Figure 14 illustrates the same effect considering only *CP* for *Case 1*, *Case 2*, and *Case 3*. The results indicate that the differences between the values of the aerodynamic and power coefficients obtained for *Case 2* and *Case 3* were not significant, with the highest difference being lower than 10.0%. The effect of *λ* on the coefficients was also similar for both cases. Concerning the magnitude of *Cd*, it decreased in the range 0.75 ≤ *λ* ≤ 1.00, and after this point, the magnitude was almost constant. The *Cl* and *CT* magnitudes decreased with the increase in *λ*. However, for *CT* the decrease was more significant only in the range 1.00 ≤ *λ* ≤ 1.25, being almost constant in other intervals of *λ*. For the *CP*, there was an increase with the augmentation of *λ*. For *Case 3*, for example, the magnitude increased from *CP* = 0.2922 when *λ* = 0.75 to *CP* = 0.5054 when *λ* = 2.0. Figure 14 indicates that in the region 0.75 ≤ *λ* ≤ 1.25, a similar trend of *CP* was noticed when *Case 1* and the enclosure cases (*Case 2* and *Case 3*) were compared. However, for *λ* = 2.0, contrary to the behavior noticed for *Case 1*, *Case 2* and *Case 3* did not have a reduction of *CP* for the highest magnitude of *λ* investigated. Therefore, the results indicated that changes in the domain where the turbine is placed are important to define its application range and the effect of *λ* over the problem performance. Further studies should be performed to define the range of application of λ for the enclosure domain.

**Figure 13.** Averaged coefficients as a function of *λ* for the cases of OWC with constant- and sinusoidalimposed velocity.

**Figure 14.** Power coefficients as a function of the tip speed ratio for the three studied cases (*Case 1*, *Case 2*, and *Case 3*).

#### **5. Conclusions**

The present work developed a computational model to investigate turbulent flows in enclosure domains with an inserted Savonius turbine simulating the air flow in an OWC/WEC device. Initially, a free turbine configuration (*Case 1*) was investigated in order to perform the verification/validation of the present method. Then, two different study cases were investigated (*Case 2* and *Case 3*) with constant and sinusoidal velocity imposed at the domain inlet. For all cases, *ReD* = 867,000, *RS* = 0.15, and four different magnitudes of *λ* were studied (*λ* = 0.75, 1.00, 1.25, and 2.00). For the simulation of incompressible, two-dimensional, transient, and turbulent flows, the time-averaged equations of mass conservation, the balance of momentum, and transport equations of the *k–ω* SST model were solved with the FVM.

The proposed computational model was verified and validated using a comparison of *CP* with the numerical and experimental results of the literature [48,49], reproducing the effect of *λ* on *CP* for the turbulent flow over a free Savonius turbine. Moreover, the velocity and pressure fields demonstrated that the driven force of the Savonius turbine was dominated by the pressure difference between the concave and convex sides of the advancement blade, which the present model adequately predicted.

After the verification/validation of the computational model, new recommendations were reached for the turbulent flows over the Savonius turbine inserted in the enclosure domain representing an OWC. The results demonstrated that the new configuration had a strong influence over the behavior and magnitudes of the instantaneous coefficients such as *Cd* and *CT*, also affecting the instantaneous values of *CP*. For the investigated values of *λ*, the results indicated that the insertion of a turbine in the enclosure domain led to an overall augmentation of *CP* for all values of *λ*. In addition, the effect of *λ* on *Cd*, *Cl*, *CT*, and *CP* was strongly affected by the domain change, indicating that the geometric configuration can be important to define the range of applicability of the turbine and the design of the PTO in the OWC device. The comparison between the imposition of sinusoidal velocity (*Case 3*) and constant velocity (*Case 2*) at the domain inlet led to similar performance and aerodynamic coefficients. It is important to highlight that due to the two-dimensional approach, some characteristics of the air flow phenomenology through the Savonious turbine cannot be completely revealed from the obtained results, being a limitation of the proposed computational model.

Future investigations should be performed considering a model where the inertia moment is taken into account to verify if this behavior continues to be observed. It is recommended to investigate other parameters, such as the effect of the overlap and spacing between the blades (*s* and *a*) and other tip speed ratios (*λ*) on the aerodynamic and performance of the OWC and the imposition of irregular velocity variation to represent the real sea state movement. It is also recommended to investigate the coupling between the present model and the wave channel for adequate simulation of the interaction between the oncoming waves, the structure of the device, and air flow over the turbine inserted in the OWC air duct.

**Author Contributions:** Conceptualization, A.L.d.S., R.S.V., J.M.P.C. and E.D.d.S.; methodology, A.L.d.S., A.L.G.S. and R.S.V.; software, A.L.d.S., A.L.G.S. and R.S.V.; validation, A.L.d.S. and A.L.G.S.; formal analysis, A.L.d.S., J.M.P.C. and E.D.d.S.; investigation, A.L.d.S. and A.L.G.S.; resources, C.F., L.A.O.R., L.A.I. and E.D.d.S.; data curation, J.M.P.C., R.S.V., L.A.I. and E.D.d.S.; writing—original draft preparation, A.L.d.S., A.L.G.S. and E.D.d.S.; writing—review and editing, C.F., L.A.O.R. and L.A.I.; visualization, C.F. and J.M.P.C.; supervision, R.S.V., J.M.P.C. and E.D.d.S.; project administration, L.A.O.R., L.A.I. and E.D.d.S.; funding acquisition, C.F., L.A.O.R., L.A.I. and E.D.d.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Brazilian National Council for Scientific and Technological Development—CNPq (Processes: 306012/2017-0, 307791/2019-0, 306024/2017-9, 131487/2020 and 440010/2019-5) and the Research Support Foundation of the State of Rio Grande do Sul—FAPERGS (Public Call FAPERGS 07/2021—Programa Pesquisador Gaúcho—PqG).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy reasons.

**Acknowledgments:** The author A.L.d.S. thanks CNPq for the Master of Science scholarship (Process: 131487/2020). The author A.L.G.S. thanks CNPq for the doctorate scholarship (Process: 440010/2019-5). The authors L.A.I., L.A.O.R., and E.D.d.S. thank CNPq for the research grant (Processes: 306012/2017-0, 307791/2019-0 and 306024/2017-9). The authors A.L.G.S. and E.D.d.S. thank CNPq for the financial support in the CNPq/Equinor Energia Ltd., Call Nº 38/2018 (Process: 440010/2019-5). L.A.I. thanks FAPERGS (Public Call FAPERGS 07/2021—PqG). J.M.P.C. thanks the Portuguese Foundation for Science and Technology funding through UNIDEMI (Process: UID/EMS/00667/2020) and the Sabbatical Leave Fellowship (Process: SFRH/BSAB/150449/2019).

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

### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Journal of Marine Science and Engineering* Editorial Office E-mail: jmse@mdpi.com www.mdpi.com/journal/jmse

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-3115-1