An impressive number of wave energy technologies have been developed over the last 25 years [
27]. Among the various types of WEC devices that convert sea wave energy into electrical energy, whether in the design phase, test phase, or pre-commercial phase, the OWC type device is possibly the most extensively studied, due to its simplicity and conversion potential [
10]. An OWC WEC device is basically composed of a partially submerged hollow structure, namely a hydropneumatic chamber, which can be fixed or floating. In addition, the device contains a turbine duct, a turbine, and a generator. The turbines used in OWC WEC devices are usually self-rectifying axial flow air turbines [
28,
29]. There are three types of air turbines, namely Wells turbines, impulse turbines, and Denniss–Auld turbines, of which the most commonly used are Wells turbines and impulse turbines [
10,
16]. The hydropneumatic chamber of the device possesses an open bottom below the free water surface, retaining the trapped air inside the chamber. The movement of the waves causes the air to be continuously compressed and decompressed, which causes it to flow through a turbine coupled to a generator [
10].
Figure 1 depicts this principle.
Several theories have been developed in an attempt to describe ocean waves, where several factors are taken into account. However, not all theories are applicable to every type of situation, since some assumptions and physical simplifications are made. While the use of regular waves is a widely adopted approach to represent ocean waves, irregular waves represent the sea state in a less idealized and more realistic way [
30]. Therefore, in this study, the focus was on the incidence of realistic irregular waves over the OWC WEC device. However, the characteristics of representative regular waves were used as parameters for the spatial discretization of the computational model as well as for a final comparison regarding the OWC WEC behavior. The representative regular waves were achieved through a statistical analysis of spectral data, which is explained in detail in
Section 2.2.1. Thus, the linear wave characteristics are shown in
Figure 2. Another important characteristic is the wave period (
T, in s), which represents the time required for two successive crests (or two successive troughs) to reach a specific point [
31].
2.1. Mathematical and Numerical Modeling
In order to reproduce an OWC WEC device subjected to realistic sea states, the commercial computational code ANSYS Fluent version 2019 R2 was employed. This is a type of computational fluid dynamics (CFD) software, based on the finite volume method (FVM), which allows broad modeling of faculties for incompressible and compressible fluids as well as laminar and turbulent flow in complex geometries [
32].
For this study, a non-linear multiphase model was employed to tackle the water–air interactions. For the fluids in the model, the governing equations were the continuity and momentum conservation equations, given by [
33,
34]:
in which
is the fluid density (kg/m
),
t is the time (s),
is the velocity vector (m/s),
p is the pressure (N/m
),
is the absolute viscosity coefficient (kg/m
s), and
S is the damping sink term.
The multiphase model used in this study employed the volume of fluid (VOF) technique [
35], which enables the reproduction of two or more immiscible fluids and the monitoring of the volume fraction (
) associated with these fluids within individual computational cells across the domain. Consequently, within a domain comprising water and air, there can be three distinct scenarios based on the volume fraction values:
= 0 denotes a cell devoid of water and filled entirely with air;
= 1 indicates a cell completely filled with water; and any other
value between 0 and 1 means the presence of the interface between the two phases within the cell. The volume fraction for these dual phases across the two-dimensional domain was governed by a transport equation, defined as [
34,
36]:
The values of the density and absolute viscosity coefficient, both essential fluid properties within the momentum conservation equation (Equation (
2)), were established using the following equations [
34]:
Regarding the flow regime, Gomes et al. [
20] performed numerical simulations with OWC WEC devices, comparing them with experimental results. In their study, the mathematical and numerical models employed considered laminar flow, and the numerical results showed good agreement with the experimental results. The authors stated that, in the flow of water and air through the hydropneumatic chamber, the predominant factor is the pressure difference, in which case, for the purpose of analyzing the airflow, turbulence can be disregarded. Therefore, only laminar flow was considered in this study.
Due to the potential interference of wave reflections with the surface elevation within the channel, a numerical beach was introduced in the domain. In this scenario, a damping sink term (
S) was incorporated into the momentum equation for the cell zone near to the pressure outlet boundary located at the edge of the wave channel [
32]:
in which
and
are, respectively, the linear (s
) and quadratic (m
) damping coefficients,
V represents the velocity in the vertical direction (m/s), and
z stands for the vertical distance from the free surface level (m). Additionally,
denotes the free surface coordinate along the vertical axis (m) and
is the bottom coordinate also along the vertical axis (m). Furthermore,
x denotes the horizontal coordinate (m),
indicates the beginning of the numerical beach (m), and
indicates its ending (m), all in the horizontal direction.
As previously mentioned, regular waves representative of the sea state were also generated in order to evaluate the hydrodynamic performance of the device in comparison to realistic irregular waves. Therefore, the horizontal (
u) and vertical (
w) wave velocity components and the water surface elevation were calculated by [
31]:
in which
g is the gravitational acceleration (m/s
),
k represents the wave number (m
), and
is the wave angular frequency (rad/s).
In addition, some numerical parameters were considered for the simulations conducted in this study. The numerical simulations employed a pressure-based solver. To address the challenge of the velocity–pressure linear interdependence, the pressure-implicit with splitting of operators (PISO) coupling scheme was used. Spatial derivatives within the momentum equations were discretized using a first-order upwind scheme, and time discretization was achieved through a first-order implicit formulation. The geo-reconstruct method was employed for the volume fraction calculation, while for the pressure interpolation at the volume faces, we employed the pressure staggering option (PRESTO!) scheme. This scheme is similar to the staggered-grid schemes used with structured meshes [
37,
38].
Finally, to analyze all the different geometries of the OWC WEC device, the hydropneumatic power (
), which depends on the static pressure, air velocity, and cross-sectional area of the duct, was calculated using the equation [
17,
39]:
in which
is the air static pressure at the turbine duct (N/m
),
is the air mass flow rate inside the turbine duct (kg/s), and
is the air velocity through the turbine duct (m/s), which, in turn, is obtained as
in which
is the two-dimensional area of the turbine duct (m
).
2.2. WaveMIMO Methodology
Machado et al. [
22] introduced the WaveMIMO methodology, which transforms sea state spectra into free surface elevation series. These data can be subsequently utilized as boundary conditions within a hydrodynamic model (such as Fluent software), allowing us to numerically generate irregular waves.
It is important to highlight that the WaveMIMO methodology allows us to always numerically generate the same sequence of irregular waves, making it useful for studies focused on the geometric configuration analysis of WECs (or other naval, coastal, and oceanic structures). Therefore, the application of the WaveMIMO methodology makes it possible to investigate different geometric configurations of a WEC under the incidence of the same sequence of irregular waves. This aspect promotes a geometric optimization study with a fairer performance comparison. Other approaches, for instance, the JONSWAP spectrum or the Pierson–Moskowitz spectrum, cannot generate the same sequence of irregular waves in different simulations. In addition, the WaveMIMO only needs the water free surface elevation of the sea state, i.e., if the elevation of the free surface is known (obtained in situ by an acoustic Doppler current profiler (ADCP) or a pressure sensor or by means of a numerical simulation (for example, in
Tomawac), this sea state can be accurately reproduced with a CFD code. On the other hand, a limitation of the WaveMIMO methodology is that the high-frequency waves are not fully propagated along the wave channel. However, concerning the fluid dynamics behavior of WECs (the focus of the present work) this limitation in reproducing high-frequency waves does not significantly affect the analysis, since low- and medium-frequency waves are the most relevant for wave energy conversion. More details about the advantages, limitations, and applications of the WaveMIMO methodology can be found in Machado et al. [
22].
Therefore, in the present work, the water free surface elevation of the sea state of interest was obtained from the
Tomawac spectral model version 7.1 (part of the Open
Telemac-Mascaret modeling suite) [
40].
Tomawac is a scientific software package that models the changes in the power spectrum of waves and the wave agitation in the oceanic domain in both intracontinental seas and coastal zones [
41]. To obtain the sea state data used in this study, the WaveMIMO methodology was employed, which corresponds to a sequence of procedures performed to obtain velocity profile data equivalent to realistic irregular waves, and then subsequently apply these in numerical simulations. An illustrative flowchart depicting all the steps necessary to employ the WaveMIMO methodology can be observed in
Figure 3.
The first step in applying the methodology is a sea state simulation, which was performed with
Tomawac by Oleinik et al. [
43]. The computational domain of the
Tomawac simulation was the continental shelf extending from the southern region of the state of Rio Grande do Sul at the southern border of Brazil to the northern part of the state of Santa Catarina.
In this simulation, three types of boundary conditions were used: surface, bottom, and oceanic. For the surface wind boundary condition, data from the Reanalysis 1 project from the National Oceanic and Atmospheric Administration (NOAA) were used [
44,
45]. For the bottom boundary condition, two types of bathymetric data were used: those for the seabed on the continental shelf and those for the deep ocean. The bathymetry of the continental shelf was obtained from nautical charts published by the Directorate of Hydrography and Navigation of the Brazilian Navy [
46] and digitized by Cardoso [
47], and the bathymetric data for the deep ocean were acquired from the General Bathymetric Chart of the Oceans (GEBCO) [
48,
49]. The ocean boundary condition was achieved using historical data sourced from the
Wave Watch III wave model [
45].
For the simulation with
Tomawac, a mesh that encompassed the Brazilian coast was employed, totaling approximately 735 km of coastline. The mesh also extended approximately 300 km into the ocean. The shape of the coastline used was obtained from the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG), developed by Wessel and Smith [
50].
In the oceanic region, the edges of the elements were approximately 10 km in length, gradually decreasing to 1.3 km on the continental shelf, and down to 250 m on the coastline. The mesh had a total of 165,548 nodes and 323,929 triangular elements.
Lastly, the
Tomawac simulation was performed from 1 January to 31 December 2014. Calibration and validation of the
Tomawac wave model for this region were performed by Oleinik et al. [
43,
45].
For the present work, the point corresponding to node 1545 of the mesh of triangular elements was chosen for the analysis (shown in
Figure 4), because it is a point located on Cassino beach in the city of Rio Grande, in the state of Rio Grande do Sul, Brazil. This point is located at coordinates 52°17
47.25
W, 32°22
30.95
S, at a distance of 2 km from the coast and with a depth of 9.52 m.
2.2.1. Spectral Data Conversion and Temporal Location
The software used for the numerical simulation of the OWC WEC device does not use spectral data as an input; therefore, conversion of these data was necessary before they were imposed as boundary conditions. For this, wave spectra were converted into a time series of the free surface elevation, which was processed by employing the methodology outlined in previous studies [
22,
42] and subsequently converted into orbital velocity components. The resulting time series encompassed significant wave heights and mean wave periods for the year 2014, as shown in
Figure 5 and
Figure 6.
Utilizing the spectral data extracted from Tomawac, it became possible to obtain wave parameters for Cassino Beach for the year 2014. Nonetheless, conducting numerical simulations of an OWC WEC device demands substantial computational resources, requiring the simulation of a relatively brief time span that was incompatible with the volume of available data spanning an entire year. Therefore, it was necessary to carry out a statistical analysis to determine the wave parameters and the time interval that best represented the studied year.
In this context, a bivariate histogram representing the time series of the significant wave height () and mean wave period () for the year 2014 was outlined. This allowed the identification of the most frequent values for these variables, ensuring an accurate portrayal of the sea state.
The histogram presented in
Figure 7 shows, for every pair of variables (
,
), the number of occurrences of that sea state in the year 2014, where it was possible to identify the consistency of the sea state attributes over the studied time. It can be observed that the predominant significant wave height was 0.66 m and the wave period was 6.30 s, with an occurrence of over 2000 times, which represents 6% of the total for the year 2014.
Although the values for and , which represent the wave climate for the coastal region of Rio Grande in the state of Rio Grande do Sul, Brazil, were calculated, it was not yet known when they occurred in the year 2014.
It is worth noting that each spectrum generated by
Tomawac encapsulates a 15 min data span. Consequently, numerous variance spectra for the year 2014 were examined in a previous study [
49]. The aim was to determine the time interval that most accurately portrayed the wave climate for the year 2014. Subsequently, the spectrum that closely matched
≈ 0.66 m and
≈ 6.30 s was selected. Remarkably, this spectrum corresponded to 23 August at 00:00 (
Figure 8).
Therefore, this is the time interval that is representative of the sea state for the city of Rio Grande which was adopted in the numerical simulations of the OWC WEC device. As previously mentioned, the time series of the free surface elevation was appropriately processed using a methodology outlined in previous studies [
22,
42]. This processing resulted in the transformation of the data into the orbital velocity components of water particles. These components were subsequently employed as boundary conditions within the hydrodynamic model.
2.3. Numerical Model Verification
Before the hydrodynamic model could be used to simulate an OWC WEC device with a realistic sea state, it needed to be verified and validated. It should be noted that Maciel et al. [
51] accomplished the verification and validation of the WaveMIMO methodology, while Maciel et al. [
52] performed a further investigation regarding the aspects involved in applying this methodology to generate numerical waves. However, an additional verification of the generation and propagation of waves for the present study was essential. For this, only the wave channel was numerically simulated and evaluated for 5 min of the simulation. In this case, wave velocity data were imposed as boundary conditions in Fluent, and then the results were compared with the realistic irregular waves obtained from the
Tomawac wave model.
The length of the wave channel was determined by calculating the characteristic wavelength of the realistic sea data for the 15 min of data used in the numerical simulations. Machado et al. [
22] computed the waves by counting and measuring them using the zero-up crossing method [
31] and then obtained their periods and wavelengths by solving the dispersion relation for
k, which is given by [
30]:
As recommended in a previous study [
20], the length of the wave channel was five times the characteristic wavelength of the realistic sea data for the 15 min period, which is
= 51.6 m. Hence, the computational domain was 258 m long, 12 m high, and 9.52 m deep, which is the same as the water depth of the point at Cassino beach, corresponding to node 1545 (see
Figure 4) of the mesh used in the spectral model simulation [
22,
43,
45].
As for the numerical beach, Machado et al. [
22] applied the WaveMIMO methodology and compared two lengths for this region of the domain: 2
and 2.5
. They obtained better results with the latter length; therefore, 2.5
on the end of the wave channel was assigned to the numerical beach, as shown in
Figure 9 as a means to dampen waves and prevent reflection. As for the calibration of the numerical beach, the values for the linear (
) and quadratic (
) damping coefficients were adopted as recommended in previous studies [
51,
53], in which
s
and
m
.
It should be noted that, in the WaveMIMO methodology, the prescribed velocity boundary condition is assigned to the line segments (or velocity profile divisions) at the entrance of the channel, which requires the discretization of the number of subdivisions that provide the proper generation of numerical waves. Previous studies have adopted different numbers of line segments [
22,
51]; however, Maciel et al. [
52] performed an investigation of the number of subdivisions that more adequately represent the generation of realistic irregular waves for the wave characteristics adopted in this study, in which the best results were obtained with 10 subdivisions at the entrance of the channel. Therefore, 10 velocity profile divisions were adopted for the boundary condition at the entrance of the wave channel in this study.
An additional aspect recently investigated regarding the generation of realistic irregular waves pertained to the process of spectral data conversion using the WaveMIMO methodology. During this process, the orbital wave velocity was calculated at the center of each line segment and then converted into horizontal and vertical wave velocity components (
u and
w). Nevertheless, some uncertainty remained regarding whether this location was suitable for accurate wave generation, or if alternative positions (such as the top or bottom of each line segment) could yield improved results [
52]. Addressing this, an investigation was conducted to determine the most appropriate location for calculating the wave velocity components alongside each line segment. This analysis revealed that better results for the wave characteristics of the present investigation were achieved by computing the wave velocity components at the bottom of each subdivision, a method that was then adopted here.
Concerning the spatial discretization of the domain, four different meshes were evaluated to ensure that proper wave generation and propagation occurred. In each mesh, mesh elements with different horizontal and vertical sizes were analyzed for the region containing only water, the free surface elevation, and the region containing only air. In the free surface region, a refinement zone was applied along the entire length of the channel, based on the stretched mesh technique [
54]. This refinement zone presented a height twice the value of the greatest free surface elevation during the 15 min of realistic sea data,
H =
= 0.66 m. The dimensions adopted in each mesh are shown in
Table 1. As for the temporal discretization, Liu et al. [
55] guaranteed numerical accuracy when
/50. Therefore, a time step of
= 0.05 s was adopted, as also performed in previous works [
22,
56].
Figure 10 illustrates the free surface elevation at
x = 0 m over a simulation period of 5 min and the free surface elevation values derived from spectral data obtained via
Tomawac conversion. In
Figure 11, the same results are presented, but with a shorter time interval (160 s
220 s), highlighting the more pronounced differences in the free surface elevation.
As can be observed in
Figure 10 and
Figure 11, the free surface elevation at the crests and troughs varies considerably, a behavior that is characteristic of realistic irregular waves. It is also possible to notice that the numerical wave model adopted for this study adequately reproduced the free surface elevation values acquired from
Tomawac by employing the WaveMIMO methodology. To assist with the analysis of the numerical model capability to represent the wave behavior, the RMSE values were calculated for each simulation and are presented in
Table 2. All analyzed meshes presented similar RMSE values, ranging from 0.0689 to 0.0686. Since meshes 3 and 4 presented the lowest values, mesh 3 was adopted for this study. This verification test ensured proper wave generation and propagation in the computational domain.
2.4. Numerical Simulations and Geometric Optimization
Once all parameters had been configured, numerical simulations were conducted to generate realistic irregular waves, which accurately reproduced the sea state conditions within the channel. Subsequently, the simulation results were compared with the time series data obtained from Tomawac. Following the tests to verify the accuracy and reliability of wave generation and propagation, the OWC WEC device was inserted into the wave channel.
The OWC WEC device used in this study was based on the onshore device located on the island of Pico in the Azores archipelago, Portugal. The report of the Non-Nuclear Energy Programme [
57] presents a technical drawing of the OWC device at Pico island. The inner cross-sectional area of the hydropneumatic chamber is 12 m × 12 m at the bottom, the turbine axis is located 10 m from the mean water level, the turbine duct has a diameter of 2.8 m and a length of 11.3 m, and the draught of the front wall is 3.4 m from the mean water level, as shown in
Figure 12.
In order to insert the OWC WEC device into the numerical wave channel shown in
Figure 9, the geometry of the device was simplified, as a means to make the mesh generation easier, as performed by [
49]. In the simplified geometry, the hydropneumatic chamber presented a rectangular shape, with the turbine duct placed vertically and the walls constructed with the minimal thickness. In the turbine duct, the dimensions of 11.3 m × 2.8 m were kept, but the duct was positioned vertically, above the hydropneumatic chamber. The width of the chamber was kept at 12 m, although in the real device, there is a narrowing towards the top, and, lastly, the distance from the mean water level to the entrance of the turbine duct was kept at 10 m. It is important to mention that, alongside this simplification, an adaptation regarding the device was made. The Pico island device is a fixed onshore device, where the rear wall of the hydropneumatic chamber is a solid barrier that completely blocks the passage of waves, and the one adopted for this study was a floating offshore device. The reason for this adaptation was to avoid wave reflection from the device and, thus, to allow a comparison with the input data.
In addition, it should be noted that energy losses due to the power take-off (PTO) system were not considered; hence, the investigation concerning the OWC performance was developed considering its available power. This simplification has been adopted in several studies regarding the OWC WEC device [
21,
22,
23,
58,
59]. Therefore, considering the turbine was not the original goal of this study.
Figure 13 shows the OWC WEC device inserted into the domain, and
Figure 14 shows the domain, indicating the location of the device.
As can be seen in
Figure 14, the OWC device was positioned at a distance of 1.2 wavelengths from the entrance of the wave channel, which corresponds to one wavelength, plus a mesh transition zone and a refined region before the device, with a length of
/5. Regarding the spatial discretization inside the device, a refined discretization was adopted, which was the same as that used for the free surface region, i.e., 20 volumes per wave height. It is also important to mention that there is a region of transition from the wave channel before the device with rectangular volumes to the OWC region with square volumes, which was made using a hybrid mesh of triangular and quadrangular volumes (see
Figure 13).
In order to measure physical quantities from the flow, numerical probes (NP) were inserted into the numerical wave channel located within the domain, as indicated in
Figure 15. Probes NP1, NP2, and NP3 measured the free surface elevation, with NP1 located at
x = 0 m, NP2 located at
x = 50 m, and NP3 located at the center of the OWC device, as calculated by the relation
x = 1.2
+ (
2). Positioned in the middle of the turbine duct, NP4 served as a measurement point for the static pressure and mass flow rate, which are essential variables for calculating the available hydropneumatic power, as indicated in Equation (
10).
Upon defining the device geometry, the constructal design method could be applied. This draws its principles from the constructal theory established by Adrian Bejan [
60]. This theory posits that the emergence of flow structures commonly observed in nature (such as river basins, lungs, atmospheric circulation, and vascularized tissues, among others) should be regarded as a physical phenomenon. Therefore, the principles of flow configuration evolution can be employed to obtain an optimization of the arrangement in a system, and this method has been widely used in various scientific fields.
In this study, the constructal design method was applied to explore alternative device geometries to analyze how varying the geometry of an OWC WEC device affects its performance under realistic sea state conditions and to estimate the hydrodynamic power available within the study area. The constructal design method requires some parameters to be applied, such as geometric constraints, degrees of freedom, and performance indicators. The procedural steps for implementing this method are outlined in
Figure 16.
The degree of freedom
, which is the ratio of the height and length of the hydropneumatic chamber, was varied in this case. Meanwhile, the degrees of freedom
(ratio of the height and length of the turbine duct) and
(submergence depth of hydropneumatic chamber) were kept constant.
Figure 17 illustrates these dimensions.
The hydropneumatic chamber volume (
) and the total volume of the OWC WEC device (
) are the geometric constraints of the problem and, as mentioned, they were kept constant. The volumes
and
are described by:
in which the third dimension,
, was kept constant and possessed a value equal to one, so the problem was regarded as two-dimensional. The equations for the length and height of the hydropneumatic chamber can be derived from Equation (
13), which changed accordingly to the degree of freedom
:
Likewise, using Equation (14), the equations for both the length and height of the turbine duct could be obtained:
Nevertheless, the ratio
was kept constant for all simulations, with
being equal to 11.3 m and
l being equal to 2.8 m, as previously mentioned. Next, the values for the ratio
were defined, which allowed the definition of values for both the lengths and heights of the hydropneumatic chambers to be simulated. Therefore, to define the volume restriction for the device chamber, dimensions akin to those of the OWC WEC device on Pico island, Azores, Portugal, were adopted, where the length of the hydropneumatic chamber is 12 m and its height is 13.4 m, as shown in
Figure 12. This represents an
value of 1.1167. For the minimum value of
investigated, the adopted height of the hydropneumatic chamber was half the height of the turbine duct, which led to
L being equal to 28.46 m,
being equal to 5.65 m, and
being equal to 0.1985. For the maximum value of
considered, the adopted hydropneumatic chamber length was three times the turbine duct length, which led to an
L value equal to 8.40 m, an
value equal to 19.14 m, and an
value equal to 2.2789. All values of
considered in this study are presented in
Table 3.
Throughout this study, the performance metric analyzed was the root mean square (RMS) of the available hydropneumatic power generated by the OWC WEC device, and the objective was to maximize this indicator.
Upon obtaining the optimal geometry for the ratios analyzed, this device configuration was subjected to regular waves with wave parameters representative of that sea state, i.e.,
H of 0.66 m and
T of 6.30 s (see
Figure 7 and
Figure 8). This analysis was performed in order to evaluate how the available hydropneumatic power changed in comparison to the use of realistic irregular waves.