Wave Power Extraction from a Dual Oscillating-Water-Column System Composed of Heave-Only and Onshore Units

: With the aim of broadening the wave-frequency bandwidth of high-efﬁciency, a small-scaled dual oscillating-water-column (OWC) system consisting of two heave-only and onshore units was numerically investigated by a well-validated computational ﬂuid dynamics (CFD) model. Based on the popular open source package OpenFOAM, the volume of ﬂuid (VOF) method was employed to track the transformation of the air–water interface under the excitation of regular waves. The six degree of freedom (6DOF) solver was applied to duplicate the heaving motion of the ﬂoating device. The effects of the two chamber widths b 1 and b 2 , the vertical restraint force (represented by the dimensionless stiffness coefﬁcient K ), the back-lip draught d 2 of the ﬂoating device, and the gap ∆ L between the two OWCs on the hydrodynamic characteristics and the wave energy conversion efﬁciencies were examined. The numerical results show that a larger width ratio b 2 / b 1 with a relatively shallow back-lip draught is more conducive to the high-performance over a broader frequency range. The ﬂoating device with a stronger vertical restraint force is more satisfactory for the high-performance of the system. Moreover, a relatively small gap is more recommended in the stage of design and construction.


Introduction
The increasing energy demand aggravates the exhaustion of traditional fossil fuels, which in turn worsens the environment issues and further drives both scientists and engineers worldwide to seek for a renewable alternative [1,2]. Wave energy, due to its huge reserves and kindly-environment, has received much attention on how to effectively exploit and utilize it [3,4]. Especially, the oscillating water column (OWC) device is considered as one of the most promising technologies due to its simple structure configuration and elegant principle of operation [5]. In general, an OWC device mainly consists of two parts: a hollow column with an immersed lower opening and a power take-off (PTO) system. Under the excitation of ocean waves, the free water surface inside the column is forced to rise and fall, which in turn compresses and decompresses the entrapped air to feed the PTO system to generate electrical energy.
Comprehensive studies have been carried out to investigate and optimize OWC devices from the perspectives of analytical analysis [6][7][8], numerical simulation [9][10][11], and physical experiment [12][13][14] in the literature. The majority of these works have more concentrated on the single OWC device and optimizing the device geometric parameters to achieve its highly-efficient hydrodynamic performance.
the reflected wave energy. The gap between the two units was expected to afford one more resonance mechanism (due to the oscillation of water column) and further broaden the high-efficiency wave-frequency range. The effects of the width ratio between the front and rear chambers, the back-wall draught of the floating OWC device, the vertical restraint force acting on the floating device, and the gap on the wave energy extraction efficiency were examined. The remainder of this paper is organized as follows. The numerical model and the method of generating and absorbing waves are described in Section 2. The verification and validation of the dynamic mesh technology model and the hydrodynamic process of a dual OWCs system are presented in Section 3. The effects of all the parameters mentioned above are discussed in Section 4. This is followed by conclusions and future works in Section 5.

Governing Equations
Based on the third-party toolbox waves2Foam in OpenFOAM, the solver waveDyMFoam was applied to solve the Reynold Averaged Navier-Stokes (RANS) equations, including mass and momentum conservative equations (Equations (1) and (2)). With the assumptions of the incompressible, viscid fluid and rotational flow, the governing equations can be expressed as: ∂ρ U ∂t + ∇ · (ρ U U) − ∇ · (µ e f f ∇ U) = −∇p * − g X · ∇ρ + ∇ U · ∇µ e f f + σκ∇α. (2) Here, U is the velocity vector, ρ is the fluid (air and water) density, p * is the pressure in excess of the hydrostatic, g is the vector of gravity acceleration, X is the position vector in Cartesian coordinates, µ e f f is the efficient dynamic viscosity (comprising the molecular and turbulent viscosity components), and the last term in Equation (2) indicates the effects of surface tension: σ and κ separately denote the surface tension coefficient and interface curvature. The second to last term on the right-hand side of Equation (2) represents the viscous diffusion due to the spatial distribution of turbulent viscosity. It is worth mentioning that the turbulent model is introduced to close the RANS equations, and the Buoyancy-modified k-ω SST model was applied in the present work (for more details, the readers are referred to [24]).
The evolution and propagation of the air-water interface is crucial for the evaluation of wave energy conversion. The VOF technique [25] was employed to capture the interface. The volume fraction α, representing the fluid volume ratio in each grid cell (α = 1 for water, α = 0 for air), satisfies a transport equation, as expressed in Equation (3): To mitigate the smearing phenomenon at the interface, an artificial compression term is added into the above advection equation [26]: where U r is the compression velocity normal to the interface. In addition, to ensure the boundedness of Equation (4) (α should always be bounded between 0 and 1), a specially designed solver named MULES (multidimensional universal limiter for explicit solution) [27] was applied in the simulations.
Consequently, the mixed density ρ and viscosity µ can be weighted in terms of the phase function α: In OpenFOAM, The PIMPLE algorithm, which is a mixture of PISO (Pressure Implicit with Splitting of Operator) and SIMPLE (Semi-Implicit Method for Pressure Linked Equation) algorithms, is commonly employed to uncouple the velocity-pressure coupling. In this study, the PISO algorithm was adopted with the default of no under-relaxation. For more detailed implementation procedures, readers are referred to [28].

Wave Generation and Absorption in OpenFOAM
To simultaneously achieve the wave generation and dissipation, special treatments (or boundary conditions) are required to construct the numerical wave tank. Waves2Foam [29], one of the most widely used libraries in OpenFOAM for the studies on coastal/ocean engineering, is developed on the basis of the solver "interFoam" and provides a series of pre-processing and post-processing tools for the simulations of regular or irregular waves.
To prevent or mitigate the interference of reflected waves from both the marine structures and the outlet boundary on the incident/transmitted wave fields, two relaxation zones were introduced at the inlet and outlet boundaries, respectively. More detailedly, the velocity U in the relaxation zone is a weighted combination of the numerical and theoretical (target) solutions, as expressed in Equation (7): where Here, α R is a weighting factor dependent on the location in the relaxation zone, and χ R is a scaled distance over the relaxation zone. χ R = 1 means the end walls of the wave tank and χ R = 0 means the interface of the relaxation and non-relaxation zone.

Dynamic Mesh Technology
The dynamic mesh technology can be used to simulate the variation of grid topology or mesh deformation with time due to the motion of rigid boundaries. The form of moving boundary can either be predefined motion, i.e., the velocity or angular velocity is specified before calculation, or be non-predefined, in which case the boundary motion is determined by the calculated results of previous step. The grid updating process is automatically completed according to the boundary changes in each iteration step. When using the dynamic mesh model, it is necessary to define the initial mesh, the mode of boundary motion and the region of mesh deformation. There exists two major motion modes for the discrete meshes, which, respectively, are mesh deformation and topological changes [30].
In this study, the six degree of freedom (6DOF) solver coupled with waveDyMFoam was utilized to model the motion of floating structures. The term 6DOF refers to the three translational and three rotational displacements in three-dimensional space. To implement the 6DOF solver, the initial mass center of moving object and motion mode(s) are specified a priori. It is worth noting that the artificial restraints, to simply represent the mooring system, are imposed by defining a linear spring force, which follows the Hooke's law. The constant stiffnessK and the rest length of spring need to be prescribed. By Newton's second law, the equation of motion for the heave-only mode might be written as: where M denotes the mass of floating object, z(t) (= A 0 e −iωt , where A 0 is the response amplitude and ω is the angular frequency) is the heave displacement, F h is the hydrodynamic force from the surrounding fluid, F B is the buoyant force, F t is the shear force due to viscous effect, andKz (t) stands for the linear vertical restraint force due to the mooring system.

Decay Test
As mentioned previously, the solver waveDyMFoam that utilizes the 6DOF technique to duplicate the dynamic response of a solid object was introduced. For brevity, the solver is termed waveDyMFoam-6DOF, which was the major computational tool for the present work to distinguish with other solvers. Undoubtedly, it is indispensable to perform a series of benchmark cases to validate the correctness and accuracy of the numerical model. Figure 1 presents a sketch of the heave decay test for a horizontal cylinder. Initially, the cylinder was immersed into the water with a draught of one-sixth cylinder diameter. During the falling process, the cylinder was only affected by gravity and buoyancy due to the hydrostatic water tank. Moreover, it is noted that the density of the cylinder is half of the water, and hence the equilibrium state is half rather than fully submerged.
As shown in Figure 1, the length of the hydrostatic tank L is 6 m, and total height is 1.44 m with a static water depth h of 1.24 m. The diameter of the floating cylinder D is 0.1523 m, and the offset displacement d is 0.0254 m. Additionally, the two ends of the tank are equipped with a one-meter long relaxation zone to absorb waves generated by the falling cylinder. To reveal the effects of boundary layer around the cylinder and capture other relevant details, the body-fitted grid was used to perform the discretization of the solution domain. As displayed in Figure 2a-c, three grids with different levels of fitness on the periphery of cylinder were designed, respectively, named coarse, medium and fine, to verify the mesh convergence. Specially, the vertical displacement, vertical velocity, and vertical force acting on the heaving cylinder during the attenuation process were focused on when comparing with the previous results, and the corresponding results are presented in Figure 3. Overall, it can be concluded that the waveDyMFoam-6DOF solver can well demonstrate the time-varying process of vertical displacement, vertical velocity and vertical force of the horizontal cylinder in the process of water entry at an initial deviation position. Moreover, the present simulated results for three grids are all in good correspondence with the experimental [31] and theoretical [32] results, and the numerical results [33]. It is worth mentioning that the instability of the dynamic grid technique usually occurred due to the velocity/pressure oscillations. In the present simulation, no obvious discrepancies were observed for the vertical velocity, and no force oscillation occurred except there existed a clear force fluctuation at the start point.

Validation of Waves Interaction with a Heave-Only Box
The simulation of the heave response of a box excited by regular waves were validated against the analytical results by Maruo [34], physical model tests by Nojiri [35], and numerical results by Koo and Kim [36] and Luo et al. [37]. Keeping the same geometric parameters as Luo et al. [37], the schematic diagram of the heave-only box interacting with incident waves is shown in Figure 4, where λ is the wavelength, A i the incident wave amplitude, A b the heave response amplitude deviating from the still water level, and B the breadth of the object.
As shown in Figure 5, the relative heave amplitude A b /A i as a function of dimensionless wave frequency ω 2 B/2g (where, ω is the angular frequency) was compared to the previous results [34][35][36][37]. Generally speaking, the present results predicted by waveDyMFoam-6DOF solver are in good accordance with analytical and numerical results, except around the system resonance frequency, which further validated the feasibility and restrictions of the dynamic mesh technology in this paper. Moreover, the interaction between waves and floating objects (such as heave-only OWCs) can be accurately simulated.

Validation of Wave Energy Conversion Efficiency ξ
The emphasis of this study was to predict the energy-conversion performance of the coupled OWC system, thus it was crucial to validate the calculations of wave energy conversion efficiency against the previous results. Fortunately, by using the commercial software Star-CCM+, Elhanafi et al. [18] investigated the hydrodynamic performance of a dual-chamber offshore-stationary OWC device, as depicted in Figure 6. Some constant parameters are presented in the graph, i.e., incident wave amplitude (A i = 25 mm), water depth (h = 1500 mm), and freeboard (h c = 150 mm). The draughts of the front and rear walls are 25 and 300 mm, respectively. Two draughts of the mid-wall are chosen: 200 and 300 mm. The slot opening ratios for the front (e 1 ) and rear chambers (e 2 ) are both 1.67%, and the tested wave periods (T) range from 0.9 to 2.0 s. In general, the amount of wave energy available for conversion by the OWC device essentially lies on both the heaving motion of the water column and the air pressure drop inside the air chamber. The time-averaged hydrodynamic energy extracted by an OWC device, E OWC , can be calculated by: where P a is the pneumatic pressure drop between the inside and outside of the chamber, η is the inside vertical motion displacement of the water column, the overdot indicates the differentiation with respect to the corresponding argument, b is the air chamber width, w is the width of wave tank (for 2D problem, w = 1), and T is the incident wave period. The wave energy absorption efficiency, ξ, which is defined by the ratio of the absorbed energy to the incident wave power (P inc = ρ w gA 2 i ω 4k 1 + 2kh sinh 2kh , where ρ w is the water density, ω the angular frequency, and k the wavenumber), can be calculated by It is worth noting that the wavenumber k is related to the angular frequency ω through the dispersion relation ω 2 = gk tanh(kh).
For the OWC device with dual chambers, the total conversion efficiency is the sum of the efficiencies from the front and rear devices. In Figure 7, the comparison of wave energy conversion efficiency ξ between the present results and those in [18] is displayed. It can be easily concluded that a good agreement was achieved for both relatively shallow (d 2 = 200 mm) and deep (d 2 = 300 mm) draughts, which further demonstrated the feasibility and correctness of the method to calculate the efficiency.

Convergence Tests
To obtain the converged results with a relatively high computational efficiency, the grid convergence for this work was necessary. Figure 8 depicts three grids with coarse, medium and fine spatial resolutions around OWC devices, where the minimum horizontal grid sizes were, respectively, 0.005 m, 0.004 m and 0.002 m. The smallest horizontal grid size was deployed on the inner sides of the OWC device (i.e., inside the column), and then became coarser gradually until equaling to horizontal unrefined grid size. The unrefined grid sizes in the horizontal and vertical directions were H/16 (H the wave height) and λ/90 (λ the wavelength), respectively. As for the temporal discretization, a fixed time step of ∆t = 0.001 s was chosen in this study, which could ensure both the numerical accuracy and stability, as recommended by Simonetti et al. [38].
and opening ratios e 1 = e 2 = 1%. It can be seen that, except for the peak and trough, there was almost no discrepancy for the three sets of grids, which indicates that the mesh convergence for the dual-OWC system was achieved with a maximum error less than 5%. In the following simulations, the medium spatial resolution was used. It is worth mentioning that, to mitigate the negative impact of the non-physical pressure fluctuations occurring in the process of coupling the CFD and 6DOF solvers, the treatment of low-pass filtering strategy to the air pressure was first performed before calculating the mean conversion efficiency over several wave periods. 22

Results and Discussion
The hydrodynamic performance of the dual OWC system, affected by the geometrical configuration and vertical restraint force, was examined numerically. Unless otherwise specified, the numerical wave tank was 20 m long, the constant water depth was h = 0.637 m, the incident wave height was H = 0.02 m, and the incident wave periods T ranged from 0.9 s to 1.8 s, as listed in Table 1. As for the locations of the OWC devices, the rear wall of the fixed one coincided with the right boundary of the wave tank, and the initial interval ∆L between OWCs was 0.075 m. As shown in Figure 10, nine numerical wave gauges (G 1 -G 9 ) were used to measure the instantaneous surface elevations at different positions. G 1 closer to wave-maker was used to monitor the generated incident wave. The time series sampled by G 2 and G 3 were used to separate the incident and reflected waves. The two sets of wave gauges (i.e., G 4 -G 6 and G 7 -G 9 ) were used to measure the surface elevations separately inside the front and rear chambers, and the averaged surface displacement was selected to calculate the vertical free-surface velocity. Moreover, four pressure probes (S 1 -S 4 ) were allocated to monitor the air pressure variation, and the average value was regarded as the air pressure drop between the inside and outside of the chambers. It is worth mentioning that only one relaxation zone at the inlet boundary was employed to undertake the functions of wave generation and dissipation, and the relevant physical quantities were notated by a subscript 1 or 2, representing the quantities in the front or rear chamber, respectively. From the perspective of the practical engineering applications, the floating OWC device was pile-supported, and hence only the heaving motion was allowed. In addition, to limit the vertical response and improve the survivability of the structure, the mooring system was applied by adding a linear spring force (with a stiffness coefficientK) in the following simulations. Here, let us introduce a non-dimensional stiffness coefficient K as: where S b is the lower plane area of the OWC device.

Effects of Chamber Widths
Generally, the chamber width is crucial for the high-performance of the OWC device. The effects of the two chamber widths, where b 1 denotes the front and b 2 the rear, were investigated first for a proper size arrangement. Seven chamber width ratios (i.e., b 2 /b 1 = 1/2, 12/18, 13/17, 1, 17/13, 2 and 5) were considered with the same total chamber width (i.e., b 1 + b 2 = 0.3 m), as listed in Table 2. It is noteworthy that, as demonstrated by He et al. [13] and Elhanafi et al. [18], a relatively large back-chamber width is profitable for high-performance of a dual-OWC device, thus the minimum value of b 2 /b 1 = 1/2 (namely, the front-chamber width is two times of the back) was chosen. The same opening ratios (e 1 = e 2 = 1%) were chosen as well, which is considered to be an optimal opening for high-performance [39]. Elhanafi et al. [40] recommended that the draught of the front-lip equivalent to the half of wave height is more beneficial for the transmission of wave energy. However, to avoid air leakage underneath the front-lip and obstruction of wave energy, all draughts (i.e., d 1 , d 2 and d 3 ) were equal to 0.05 m. The non-dimensional spring coefficient K was equal to 2 for the modeling of the modest vertical restraint force. Table 2. Cases for different front (b 1 ) and rear (b 2 ) chamber widths. The variations of the relative surface elevations (η j /A i , j = 1, 2, where η j is the mean amplitude of surface elevation), relative pressure drops (P j /ρgA i , j = 1, 2, where P j is the mean amplitude of pressure drop), and extraction efficiencies (ξ j , j = 1, 2) versus ω 2 h/g are plotted in Figure 11a-f. In Figure 11a, the curves display the similar variation tendency except for the relatively small front-chamber width (i.e., b 2 /b 1 = 1/2), where the relative amplitude has a relatively large value over a wide range of wave frequencies. According to Boyle's law, the air pressure inside a closed chamber is inversely proportional to the air volume. In other words, the more violent is the internal water-free surface, the more dramatic is the inside pressure fluctuation. However, due to the heave response of the floating OWC device as well as the phase lag between the air volume flux and the air pressure fluctuation, the relatively drastic surface fluctuation did not correspondingly induce a violent pressure oscillation, as shown in Figure 11c. The variation trends of the relative surface elevation and pressure drop for the rear device are depicted in Figure 11b,d. Similarly, the relatively violent pressure fluctuation did not precisely correspond to the strong surface elevation due to the phase lag. Meanwhile, it was found that a relatively small chamber width, referring to the curves of b 2 /b 1 = 1/2, is helpful to mitigate the phase lag effects. Figure 11e shows the variation of the front extraction efficiency ξ 1 against the dimensionless wave frequency ω 2 h/g. The curves have multiple peaks with increasing ω 2 h/g, which was due to the existence of multiple resonance mechanisms (such as the two water columns and heaving-free structure) [17,37]. However, for the stationary device, the efficiency curves shown in Figure 11f have only one peak for all the chamber widths, which means that, except for the resonance of the rear OWC, the effects of the other resonance mechanisms were negligible. By comparing the peak values in Figure 11e-f, it can be concluded that the majority of the total extracted energy was contributed to the rear device, while the energy extracted from the front was only a good complement. It is worth mentioning that, for b 2 /b 1 = 5, the energy extracted from the front device reduced by about 37% compared to the case of b 2 /b 1 = 2, whereas the energy from the back was almost equivalent. Moreover, it was demonstrated that the inside free surface displacement significantly determined the extraction efficiency (a similar conclusion was presented by Ning et al. [19]).  The total extraction efficiency, which is the sum of the efficiencies of the front and rear devices (i.e., ξ tot = ξ 1 + ξ 2 ), is depicted in Figure 12. A similar variation tendency to the efficiency of the rear device was observed. Under the condition of a constant total chamber width (i.e., b 1 + b 2 = 0.3 m), a relatively large width ratio (i.e., b 2 /b 1 = 2) was significantly beneficial for the performance improvement of the present system, approximately increasing by 30% via comparing the peak values of b 2 /b 1 = 2 and 1/2. A similar conclusion was found by He et al. [17]. However, continuously increasing the width ratio (e.g., from b 2 /b 1 = 2 to 5) could not considerably lift the total efficiency over the whole wave-frequency range. Especially, a reduced efficiency could be noticed in the intermediateand short-wave regimes. Therefore, the chamber width ratio of b 2 /b 1 = 2 is recommended for high-performance. More specifically, as the chamber width ratio was b 2 /b 1 = 2, the maximum extraction efficiency even approached 0.87, increasing by about 30% compared to the ratio of b 2 /b 1 = 1/2. It should be noted that the present finding is somehow contradictory with those of Ning et al. [19], who concluded that the total efficiency is not sensitive to the width ratio. The reason for this contradictory is probably due to the existence of the gap.

Effects of the Back-Lip Draught d 2 of the Heave-Only Device
The immersion depth of the OWC device is a significant factor affecting the hydrodynamic characteristics. Many studies on the effects of the wall draughts have been performed (e.g., [17,41,42]), and the consistent conclusion that a relatively smaller front-lip draught is more helpful for the highly-efficient extraction is drawn. In this study, a constant front-lip draught (i.e., d 1 = d 3 = 0.05 m) was chosen, and five different back-wall draughts (i.e., d 2 = 0.05, 0.1, 0.2, 0.3, and 0.4 m) of the heave-only device were introduced to consider the effects on the hydrodynamic performance. It is worth mentioning that, referring to the above findings in Section 4.1, the chamber width ratio of b 2 /b 1 = 2 was used in experiments presented here and in forthcoming subsections. Figure 13 shows the effects of the back-wall draught d 2 of the front floating device on the front ξ 1 and rear ξ 2 conversion efficiencies and the total efficiency ξ tot . Figure 13a shows that gradually increasing the draught from d 2 = 0.05 to 0.2 m improved the performance of the front device in the intermediate wave-frequency range tested, while further increasing, namely from d 2 = 0.2 to 0.4 m, reduced the performance over the entire tested range. However, different from the performance of the front device, the peak of efficiency curve of the rear device shifted to the lower frequency range, and the highly-effective wave-frequency bandwidth became narrow with the increase of the back-wall draught, as shown in Figure 13b. Although the back wall was not directly connected to the rear fixed device, it seemed to play the same role as the front-lip of the rear device. As a result, a larger d 2 seriously weakened the performance of the rear OWC as most of the wave energy was reflected from the back-wall, especially for short waves. In addition, the increase of the back-wall draught prolonged the water particle trajectory within one wave period such that the peak frequency shifted to the low frequency zone [15].
The total extraction efficiency as a function of ω 2 h/g for five different back-wall draughts is depicted in Figure 13c. Since the efficiency Was mainly contributed to the rear device, the overall variation tendency of the total efficiency was the same as that shown in Figure 13b. To improve the competitiveness of the dual-OWC system, the minimum draught d 2 = 0.05 m is recommended. In this way, the entire system can operate with a maximum efficiency of ξ tot ∼ 0.87 and over a wider highly-effective bandwidth. In the practical engineering applications, the shallower immersion depth is beneficial for the high-performance of the system, but the air leakage must be taken into account in the design stage. Figure 14 illustrates the effects of the back-wall draught on the reflection coefficient (which is defined by the ratio of the reflected wave height to the incident one). It further demonstrates that increasing the draught could prevent more of the wave energy from transmitting into the rear OWC, and cause the total-reflection frequency to shift to the low wave-frequency area.

Effects of Non-Dimensional Spring Coefficient K
In practice, the offshore structure is almost impossible to be fixed on the water surface, especially for the heave motion. Next, the effects of the non-dimensional spring coefficient K (representing the vertical restraint force) defined by Equation (12) on the hydrodynamic efficiencies, such as the front, rear and total efficiencies, were examined. Invoking the previous findings, the present geometrical configuration was updated as follows: Additionally, four different non-dimensional spring coefficients K, namely K = 2, 4, 8 and ∞, were considered. It is noted that the front OWC device was fixed as K → ∞. Figure 15 illustrates the variations of the heave-response amplitude of the front OWC device against the dimensionless wave frequency ω 2 h/g. As expected, the smaller was the spring stiffness, the stronger was the heave response under the excitation of regular waves. To further improve the performance of the floating OWC device, a phase control method could be adopted by tuning the spring coefficient to achieve that the internal water-free surface moves up, while the base structure falls down. As shown in Figure 16a, the variation of spring coefficients has noticeable influences on the front floating device, and generally a larger vertical restraint force was more helpful for the high-performance of the front device over a wider tested frequency range. Physically, the heaving response of the front device Was more sensitive for a smaller coefficient K. Hence, by properly tuning the phase difference between the front water column and the heaving motion, the highly-effective bandwidth can be broadened. From the efficiency curves of K = 2, 4 and 8, multi peaks were observed, which conformed to the feature of the floating OWC device.
The effects of the non-dimensional spring coefficients K on the efficiency ξ 2 of the rear fixed OWC device are described in Figure 16b. Although the spring force was not directly imposed on the rear device, the influence could not be ignored. For the intermediate-and long-wave regimes tested, the efficiency reduced with increasing the spring stiffness. As for the total efficiency shown in Figure 16c, the completely inverse variation was observed compared to that in Figure 16b. From the perspective of highly-effective operations, a larger non-dimensional spring coefficient is more beneficial. However, although a completely stationary marine structure is more satisfactory, the corresponding construction cost is high and the feasibility is quite low. Therefore, the non-dimensional spring coefficient K = 8 was chosen for the following simulations.

Effects of Devices Interval ∆L
As aforementioned, in addition to the common tasks of absorbing and extracting wave energy, the two devices will inevitably interact with each other. With the aim of broadening the highly-effective bandwidth, another water column, denoted by a gap ∆L in between the front and rear devices, was introduced. The effects of the gap on the hydrodynamic efficiencies were examined. More specifically, ten intervals Were classified into two regimes: one was the relatively short interval regime (i.e., ∆L/b 1 = 0.75, 2, 3, 4, and 5) such that the spatial variation of the free surface could be ignored; and the other was the relatively long interval regime (i.e., ∆L/b 1 = 8, 10, 15, 18, and 20) such that the spatial variation of the free surface was not negligible. The detailed parameters are listed in Table 3. For the short interval regime, the effects of the interval ∆L on the front, rear and total conversion efficiencies were first investigated under the conditions of b 2 /b 1 = 2, d 1 = d 2 = d 3 = 0.05 m, and K = 8, as shown in Figure 17. The efficiencies for both the floating and fixed single OWC devices with the same dimensions, in comparison with the dual-OWC system, are presented in Figure 17a,b, respectively. Figure 17a shows that the increase of the interval reduced the front efficiency in the long-wave range but increased the efficiency in the intermediate-and short-wave range. In addition, the advantage of the floating single device only reflected within the long-wave regime, which is not desirable in the engineering practice. However, as shown in Figure 17b, reducing the interval (from ∆L/b 1 = 5 to 0.75) could significantly improve the entire performance with a maximum increment of about 100%. The fixed single device was more satisfactory only for the short-wave regime. Figure 17c shows the variations of the total efficiency ξ tot . A similar tendency to that in Figure 17b was found. Moreover, as ∆L/b 1 = 0.75, the maximum efficiency value of 0.86 was obtained, and the efficiency larger than 0.6 was achieved over the range of 1 ≤ ω 2 h/g ≤ 3.1. In Figure 18a-c, the variations of the hydrodynamic efficiencies against the dimensionless frequency ω 2 h/g for five different intervals in the long-interval regime are illustrated. For comparison, the results of the single OWC devices are presented as well. Physically, when the gap was larger than half of wavelength, the standing wave was formed in between the two devices. As a result, more energy was dissipated due to the multi-transmission processes, and some of wave energy was entrapped in the water column bounded by these two devices. It was expected that multiple peaks of the efficiency curves occurred in the tested frequency range. Moreover, with increasing the gap, the peak frequency shifted to the low frequency area. For the total efficiency, the widest bandwidth of the efficiency value larger than ξ tot = 0.6 was in the range of 1.6 ≤ ω 2 h/g ≤ 2.6 as ∆L/b 1 = 8. Finally, the contour of the total conversion efficiency as a function of ω 2 h/g and ∆L is plotted in Figure 19. It was found that a relatively small interval was more beneficial for the highly-effective extraction over a wider wave-frequency range. The total extraction efficiency tot Figure 19. The contour of the total efficiency as a function of the gap ∆L and the dimensionless frequency ω 2 h/g. In all cases, b 2 /b 1 = 2, d 1 = d 2 = d 3 = 0.05 m, and K = 8.

Conclusions
In this study, the hydrodynamic performance of a dual-OWC system comprised of two heave-only and onshore stationary units was investigated numerically, using the open source package OpenFOAM. The influences of the front (b 1 ) and rear (b 2 ) chamber widths, the back-lip draught (d 2 ) of the floating device, the non-dimensional spring coefficient K, and the interval ∆L between two devices on the hydrodynamic efficiencies were examined. After analyzing and discussing, the following conclusions can be drawn from this study: (1) For a dual-OWC device, the configuration that the rear chamber width is larger than the front one is more helpful for high-performance over a wide wave-frequency range, and the total efficiency is mainly dominated by the rear unit.
The back-lip draught d 2 of the front heave-only unit can significantly affect the hydrodynamic efficiencies and the resonance of the rear OWC. A relatively small draught can improve the system performance to a great extent, especially for intermediate and low frequency waves.
The vertical spring restraint force significantly influences the front heave-only device, and a larger spring stiffness is more advisable for the high-performance of the whole system. (4) The interval between these two OWC devices is a significant factor that influences both the front and rear OWC devices. A small interval is recommended for practical engineering applications.
The 2D dual-OWC system was investigated in this study. To make the topic more realistic, a 3D numerical model should be considered, as some pivotal factors and characteristics are indispensable, such as the wave diffraction phenomenon. With a scaling factor of 1:20, the total width of the present prototype is 6 m. Under the conditions of water depth h = 12 m, wave period T = 6 s and wave height H i = 0.4 m, the incident wave power per unit crest length is P i = 1080 W/m, and the extracted wave power for a total extraction efficiency of ξ = 0.5 is approximately P out = 540 W/m. It is noteworthy that the extracted wave energy based on a 2D numerical model might be overestimated within some wave-frequency ranges, as reported by Elhanafi et al. [40]. Nevertheless, the hydrodynamic performance and extraction efficiency presented by the current investigations might still provide a useful guidance to the practical engineering design and applications, such as the optimal configuration of the combination of dual OWC devices.
Author Contributions: Z.D. proposed the idea and gave the essential suggestions on paper writing; C.W. performed the numerical simulations, collected and organized the results, and wrote the initial version of the manuscript; P.W. and Y.Y. contributed to the discussion of this study. All of the above authors revised the paper.

Conflicts of Interest:
The authors declare no conflict of interest.