Efficient response of an onshore Oscillating Water Column Wave Energy Converter using a one-phase SPH model coupled with a multiphysics library

In this paper the numerical modelling of an Oscillating Water Column (OWC) Wave Energy Converter (WEC) is studied using DualSPHysics, a software that applies the Smoothed Particle Hydrodynamics (SPH) method. SPH is a Lagrangian meshless method used in a growing range of applications within the field of Computational Fluid Dynamics (CFD). The power take-off (PTO) system of the OWC WEC is numerically modelled by adding a force on a plate floating on top of the free surface inside the OWC chamber. That force is implemented in the multiphysics library Project Chrono, which avoids the need of simulating the air phase that is computationally expensive in the SPH methods. Validation is carried out with experimental data received from the Korea Research Institute of Ship and Ocean Engineering (KRISO) and Ocean Energy Systems (OES) of the International Energy Agency (IEA) Task 10. The numerical and experimental water surface elevation at the centre of the OWC WEC chamber and the airflow speed through the orifice are compared for different wave conditions and different PTO systems (different orifice diameters at the top part of the chamber of the OWC WEC). Results show that DualSPHysics is a valid tool to model an OWC WEC with and without PTO system, even though no air phase is included.


Introduction
Wave energy is a potential source of clean electricity that can make a significant contribution to the de-carbonization of the world's electricity supply. The Oscillating Water Column (OWC) wave energy converter (WEC) is one type of WECs that has been studied extensively in literature (Falcão and Henriques, 2019;Sheng et al., 2013;Elhanafi et al., 2017;Zhu et al., 2020). The OWC WEC consists of a fixed or floating structure that is open to incoming sea waves, with air trapped above the free surface inside the OWC WEC chamber. The up-and downward motion of the free surface inside the OWC WEC chamber causes the air to flow through an air turbine, thereby generating electricity. Examples of full-scale OWC WECs that have been installed in the past are the LIMPET plant near the island of Islay, Scotland, the PICO plant in the Azores, Portugal, and the Mutriku plant in Spain .
Most studies concerning WECs and OWC WECs employ potential flow theory based on the linearized form of the Navier-Stokes equations. Applying linear potential flow theory allows numerical modelling of WECs in the time or frequency domain (Folley, 2016), enabling fast calculations of the WEC's motion. This method needs to assume small amplitude oscillations of the WEC and the fluid to be incompressible, inviscid and with an irrotational motion. This results in an underestimation of the wave-induced forces on WECs under highly-nonlinear sea states (Windt et al., 2018). Since previous simplifications are too restrictive when modelling WECs under non-linear conditions, it may be more appropriate to employ higher fidelity models, generally known as CFD (computational fluid dynamics) methods. The most commonly used CFD methods in hydrodynamics are mesh-based. Different OWC WECs have been studied numerically using these methods: in Elhanafi et al. (2017) the Star-CCM+ code is used, in Kamath et al. (2015) the open-source model REEF3D is applied, in Zhang et al. (2012) a two-phase model using a level-set immersed boundary method is used and in Vyzikas et al. (2017) the open-source package OpenFOAM was used. Despite being very robust mathematically and computationally, mesh-based methods face important challenges when capturing the free surface and the rapidly-evolving nonlinearities. Due to their ability to overcome these drawbacks, meshless CFD methods have gained attention during the last years, being the smoothed particle hydrodynamics (SPH) method one of the most widely used (Violeau and Rogers, 2016 The SPH method has been used in various research areas, e.g. coastal engineering (Luo et al., 2021;Gotoh and Khayyer, 2018) or waterrelated natural hazards modelling (Manenti et al., 2019). In contrast to the mesh-based methods, in SPH the fluid is discretized in a series of points, named particles, that move with the velocity calculated from the Navier-Stokes carrying all physical properties with them (Monaghan, 1992). Its meshless Lagrangian formulation makes the SPH method a very interesting alternative when simulating free-surface flows with wave-structure interaction (González-Cao et al., 2019), such as the case at hand. In this research, the open-source software Dual-SPHysics ) (available at www.dual.sphysics. org) is employed to model the OWC WEC in different sea states. The DualSPHysics code can be executed either on CPU and GPU (Graphical Processing Unit) where the parallel power of the graphics cards can be exploited to run simulations. Furthermore, DualSPHysics has been recently coupled (Canelas et al., 2018) with Project Chrono (Tasora and Anitescu, 2011), enabling the effect of the PTO system to be included in the simulations. The compression and decompression of the air in the OWC WEC chamber causes the free-surface elevation to be damped (PTO system damping). This software has proven to be a valuable tool in the modelling of wave-structure interaction in general and of WECs in particular: in Ropero-Giralda et al. (2020, Tagliafierro et al. (2019) and Quartier et al. (2021) point-absorber WECs were studied, in Brito et al. (2020) an Oscillating Wave Surge Converter (OWSC) was modelled, an OWC WEC was simulated in Crespo et al. (2017) and in Verbrugghe et al. (2018Verbrugghe et al. ( , 2019 a coupling methodology was developed in order to compute the modified wave field effects surrounding a point-absorber WEC. The application of SPH to industrial problems is one of the SPHERIC Grand Challenges  and, in particular, the simulation of wave energy converters is one of the most promising application fields (Stansby, 2018). The studies of Rafiee et al. (2013) and Edge et al. (2014) were pioneers presenting SPH simulations of oscillating wave surge devices. The works of Omidvar et al. (2013) and Yeylaghi et al. (2015) were the first ones to deal with extreme waves interacting with point-absorbers using SPH. In fact, other SPH-based codes were coupled with Project Chrono to simulate OWSCs such as Wei et al. (2019) and Liu et al. (2020). The SPH method has been applied before for the modelling of OWC WECs Wen et al., 2018;Zhu et al., 2020). However, these previous studies possess some limitations since in Crespo et al. (2017) and in Wen et al. (2018) the free surface inside the OWC WEC chamber was assumed undamped and in Zhu et al. (2020) the OWC WEC was only modelled in 2D while only assuming incompressible air in the OWC WEC chamber. In the current study the OWC WEC is modelled in 2D and 3D. Furthermore, a novel implementation of the PTO system force is added using the DualSPHysics Project Chrono coupling, including the effect of air compressibility. The experimental data from Bingham et al. (2021) has been used for validation in the current research. Therefore, the novelty of this work lies in the simulation of air compressibility of the OWC WEC but using a one-phase numerical model. The paper is organized as follows: Section 2 describes the basic theoretical principles of SPH and its implementation in DualSPHysics; Section 3 gives an overview of the experimental campaign carried out at KRISO and used for validation in the current research; Section 4 provides a description of the methodology followed to model the OWC WEC in DualSPHysics, with and without PTO system; Section 5 discusses the main results obtained using DualSPHysics; finally, Section 6 presents an overview of the conclusions.

DualSPHysics: a solver based on smoothed particle hydrodynamics
DualSPHysics applies the SPH method, which is a meshless Lagrangian method used in a growing range of applications within the field of CFD. In the current study DualSPHysics is run on a Graphics Processing Unit (GPU). In SPH the fluid is discretized in a set of particles for which the position, velocity, density and pressure are computed by solving the Navier-Stokes equations and by interpolation of the values of neighbouring particles. The contribution of each of the neighbouring particle depends on the Kernel function , which has an area of influence defined by the smoothing length ℎ - (Monaghan, 1992). In the DualSPHysics code the smoothing length is defined starting from the initial inter-particle distance ( ), Domínguez et al. (2021).

Governing equations
The Navier-Stokes equations written in their SPH notation are solved at each timestep and for each of the particles. The momentum and continuity equation are given in Eqs.
(1) and (2), respectively. In the following equations the physical properties of particle are calculated, with each of its neighbouring particles: with the mass, the velocity, the pressure, the density and the gravity acceleration. represents the Kernel function and depends on the distance between particles and . In this work a Quintic Kernel (Wendland, 1995) is applied. The diffusion term introduced in the right hand side of the continuity equation (Eq. (2)) acts as a numerical noise filter thereby improving the numerical stability and smoothing the density and pressure field (Antuono et al., 2012). In this paper the recently proposed density diffusion term of Fourtakas et al. (2019) is applied, since it has proven to produce more accurate results for the pressure field near boundaries while keeping the computational cost limited. The term takes the following form: where = − with the position of particle and the dynamic density, equal to the difference of the total ( ) and hydrostatic ( ) density: = − .
represents the artificial viscosity term, proposed in Monaghan (1992) with: with = − , the velocity of particle k, = ℎ ⋅ 2 + 2 , = 0.5( + ) is the mean speed of sound, 2 = 0.1ℎ 2 and is a coefficient tuned for proper dissipation -see Domínguez et al. (2021). Throughout this paper will be set equal to 0.01 . Since the fluid is weakly-compressible in DualSPHysics, an equation of state, Eq. (5), is used to calculate the fluid pressure as function of the density instead of solving a Poisson-like equation. The speed of sound is artificially lowered such that a reasonable time step can be employed while ensuring that density variations are kept lower than 1% during the simulation.

Boundary conditions and floating bodies
In DualSPHysics the boundaries are described by a set of particles for which the same equations (Eqs. (1) and (2)) as for fluid particles are solved. However, the particles belonging to the boundary do not move according to the forces acting on them: boundary particles remain either fixed or move according to a predefined movement. These boundary conditions are called Dynamic Boundary Conditions (DBC) and have the advantage of being able to deal with complex geometries (Zhang et al., 2018) and being computationally efficient (Crespo et al., 2007). However, due to excessive repulsive forces near the boundary between a structure and the fluid, a gap appears of the order of magnitude of the smoothing length ℎ -see English et al. (2020). Furthermore at this same boundary unphysical values of the pressure are observed. These issues are dealt with in the modified DBC (mDBC) implementation recently added to DualSPHysics, as described in English et al. (2020). For each boundary particle of mDBC a ℎ is mirrored into the fluid across the boundary interface in a similar procedure to Marrone et al. (2011). Boundary particles then receive the properties of the fluid at the position of the ghost node. The density is calculated using a first order consistent SPH interpolation, as proposed in Liu and Liu (2006). The mDBC method has proven to show results with more realistic physical values for the pressure at the boundaries and an elimination of the gap between fluid and boundary without a significant extra computational cost (Quartier et al., 2021).
DualSPHysics DualSPHyics has the capability to accurately simulate fluid-driven objects which will be used extensively in this paper. Validation of sinking and buoyant objects can be found in Canelas et al. (2015) while floating box under regular waves are compared with experimental data in Domínguez et al. (2019) (with experimental data from Ren et al., 2015).
The force per unit mass acting on one boundary particle of the floating body is calculated by summing up the contributing forces exerted on this boundary particle by fluid particles within the compact support of the Kernel: with the force per unit mass acting on boundary particle of the floating body and the force per unit mass acting on boundary particle exerted by fluid particle , calculated with Eq. (1), Canelas et al. (2015). Once the force acting on the floating body is computed, the body's motion can be determined assuming it is rigid: with the body's total mass, the moment of inertia, the velocity, the rotational velocity, the centre of mass, the mass of floating boundary particle and the position of floating particle . Eqs. (7) and (8) are integrated over time in order to compute the velocity of the floating particle:

Coupling DualSPHysics-Project Chrono
Project Chrono is an open-source software package that enables the numerical modelling of mechanical constraints and collision between objects (Tasora and Anitescu, 2011). It has recently been successfully coupled to DualSPHysics (Brito et al., 2020;Canelas et al., 2018) and can be used in this case for the modelling of the PTO system of a WEC.
Already implemented in the coupling was the linear damping, which is useful for the simulation of a WEC with a linear damping PTO system -see Eq. (10) with , the linear PTO system damping coefficient and ( ) the WEC's heave velocity. The DualSPHysics-Project Chrono is here extended with the PTO system force of the OWC, acting on the free surface, see Section 4 and Eqs. (17), (20).

Experimental campaign
As part of the International Energy Agency (IEA) Technology Collaboration Program for Ocean Energy Systems (OES) a WEC modelling task was included (Wendt et al., 2019;Nielsen et al., 2018). This task aims at assessing the accuracy and reliability of numerical tools for the modelling of WECs (Kramer et al., 2021), with one of those WECs being the OWC WEC. The experimental data used for validation in the current research is a 1/4 scaled onshore OWC WEC model tested at the ocean basin located in the Korea Research Institute of Ships and Ocean Engineering (KRISO) (Sewan Park et al., 2019). This experimental data is used also for comparison to different numerical models in Bingham et al. (2021), however in Bingham et al. (2021) no comparison with results obtained from an SPH model was carried out. The OWC WEC chamber was provided with an air duct at the top part that has an orifice inside to account for the damping effect of the air turbine. Experiments were performed with different orifice diameters: 0.4D and 0.5D, with D the duct diameter of 0.2 m. The geometry and dimensions of the physical model are given in Fig. 1.
The variables measured during the experimental campaign were (i) the water surface level inside the chamber, measured at the centre of the OWC WEC chamber, (ii) the airflow speed through the orifice and (iii) the differential pressure between just below and above the orifice. A more detailed description of how the experiments were conducted and how the measurements were carried out can be found in Bingham et al. (2021) and Sewan Park et al. (2019).

OWC WEC modelling
In experimental tests the damping effect caused by the air turbine present in the OWC WEC is modelled by providing a small orifice at the top of the OWC WEC chamber. This causes air above the free surface to compress and decompress when the free surface moves upwards or downwards, respectively. In the DualSPHysics code, this can be modelled by adding an air phase (Mokos et al., 2015(Mokos et al., , 2017. However, modelling the air with SPH particles implies a significant increase in computational time due to three main reasons: (i) the air volume must be discretized into particles which increases drastically the number of particles involved in the simulation, (ii) many more time steps are needed since the speed of sound required to model the air is several times lower than the one of the fluid, (iii) the problem to solve is more complex and involves more computational operations. Therefore the execution runtime increases significantly as shown in the simulation of a dam break impacting an obstacle studied in Mokos et al. (2015), where the execution runtime using the multiphase code is approximately 278 times higher in comparison with the single-phase DualSPHysics code. On the other hand, the computation runtime required by Chrono is not affected by the number of particles and it becomes negligible in high-resolution 3D simulations. In this paper the damping effect of the orifice is modelled by adding a plate floating on top of the free surface on which a PTO force is acting, representing the compressing/decompressing force of the air phasesee Fig. 2. When the free surface inside the OWC WEC chamber rises from 1 at time 1, to 2 at time 2, an air pressure develops. In the N. Quartier et al. current study this air pressure is taken into account by implementing acting on the floating plate on top of the free surface - Fig. 2. Note that this model will be effective when the free surface is uniform inside the air chamber as is shown in Fig. 2.
The PTO force is numerically modelled using the DualSPHysics-Project Chrono coupling, modified to represent a correct . The floating plate itself is a thin plate with a low mass, designed so that it almost perfectly follows the water surface elevation at the centre of the OWC WEC chamber when no is acting on it. This will be true when the heave natural period of the plate is lower than the period of the incoming waves. It is ensured that the floating plate can only move along the direction parallel to the front or back side of the OWC WEC chamber. The air pressure acting on the floating plate can be calculated by assuming incompressible air -see Section 4.1 or compressible air inside the OWC WEC chamber -see Section 4.2. Several studies have shown the importance of air compressibility in the air chamber and how this is related to the model scale, so that the air compressibility should be considered in large-scale models but can be ignored in small-scale ones (Simonetti et al., 2018;Chen et al., 2021).

Incompressible air inside OWC WEC chamber
Initially, the air inside the OWC WEC chamber is assumed incompressible, which is acceptable considering the small scale (Sheng et al., 2013;Dimakopoulos et al., 2017). This means that the air flow rate driven by the water surface (see Eq. (11)) equals the air flow rate through the PTO system (or in this case through the orifice).

= −
with a constant depending on the air density = 1.2 kg∕m 3 , a calibration constant equal to 0.82 for the currently studied experiments (see Bingham et al., 2021), and the orifice area : The calibration constant is calibrated with experimental data. If no experimental data were available could be estimated analytically using Eq. (14) as suggested in : with the pressure loss coefficient representing the head loss through the contraction and 0 the area of the OWC WEC chamber. from Eq. (11) can now be written as follows: witḣ( ) the heave velocity of the floating plate. It is important to note thaṫ( ) represents the heave velocity of the floating plate, which is equal to its total velocity multiplied with ( ), = ( 1 1.5

)
. Filling in Eqs. (13) and (15) in Eq. (12) leads to the following equation: It is assumed that the pressure is the same in 2D as in 3D, but the will be different since the pressure in a 2D simulation only acts along the width of the plate. The given in Eq. (17) is implemented in the DualSPHysics-Project Chrono coupling in 2D as follows: N. Quartier et al. Considering air pressure from Eq. (16) we will get: with 0,2 the 2D area of the floating plate being equal to its width. Since the air is considered incompressible the airflow speed through the orifice can be calculated as in Eq. (19) When simulating the OWC WEC in 3D the same dimensions as in the 2D simulation will be used. It is assumed that the pressure is the same in 2D as in 3D, but the PTO force is different since the pressure now acts on a plate with another area 0,3 -see Eq. (20).

Compressible air inside the OWC WEC chamber
It has been concluded in previous research that assuming incompressible air inside the OWC WEC chamber may lead to an overestimation of the average absorbed power (Sheng et al., 2013). Therefore it was decided to add the effect of compressibility in the DualSPHysics simulations by slightly altering the formula. When compressibility is taken into account the air flow rate through the orifice does no longer equal , but changes due to the variable air density and pressure inside the OWC WEC chamber (as mentioned in Sheng et al., 2013)-see Eq. (21): where is the specific heat ratio of air (=1.4) and 0 the air volume in the OWC WEC chamber in calm water. This is a simplified formula assuming the air volume changes due to the oscillating water surface inside the OWC WEC chamber are small compared to the initial air volume in calm water 0 . No distinction is made between the air flow rate for inhalation and exhalation. The pressure difference is then again calculated similarly as for the incompressible case: The difficulty in implementing this modified equation for the pressure difference lies in the derivative of , appearing in the formula for (Eq. (21)). The pressure depends on , but in order to calculate , needs to be known at the current time step. Therefore an approximation for is used by deriving Eq. (22): In Eq. (24), is simplified as , with from Eq. (15). In Eq. (24) the z derivatives are computed with the experimental water surface elevation and a good agreement is obtained when compared with the experimental value of -see Fig. 3. It is expected that this approach is less accurate when considering large-scale models since in Eq. (26) the assumption is made that the airflow caused by the free surface inside the OWC chamber is the same as the airflow through the orifice . This is no longer the case when air compressibility becomes more significant.
Eq. (24) is inserted in Eq. (21), resulting in the following expression for : Eq. (25) is now used to calculate the air pressure with Eq. (22): and finally to calculate the PTO system force -see scheme in Fig. 4.

Numerical validation and discussion
The OWC WEC with the same geometrical characteristics as given by OES TASK 10 is simulated in DualSPHysics. The testcases are described first and the results are discussed in the following paragraphs. Two different GPU cards have been used for the simulations: GeForce GTX Titan Black for the 2D cases and GeForce RTX 2080 Ti for the 3D ones.

Testcases and numerical setup
The simulations that are carried out in DualSPHysics are summarized in Table 1. Table 1 includes two different wave conditions: Case A with wave height, , of 0.0881 m and wave period, , of 3.25 s and Case B with =0.089 m and =3.5 s. Different PTO systems will be also considered: an orifice with a diameter of 50% of the duct diameter (0.5 ), an orifice with 40% of duct diameter (0.4 ) and no PTO system with the completely open duct (1.0 ). Fig. 5 shows a schematic representation of the numerical wave flume in DualSPHysics. The OWC WEC structure is located at one wavelength ( ) from the wave generation piston and an Active Wave Absorption System (AWAS) is applied in order to reduce reflection from the wave generation piston . A dissipative beach is included from the start of the beach to the end of the tank to minimize wave reflection. In fact, a numerical damping layer  is provided in combination with the beach in order to guarantee that, for the simulated wave conditions (listed in Table 1) N. Quartier et al.   a reflection coefficient of 3% was obtained. The numerical tank will be 3.2 m wide while the chamber width is 2.5 m, therefore periodic boundaries (Gomez-Gesteira et al., 2012) and damping layers  of 0.2 m wide are applied at the lateral walls, in order to reduce sideways reflection . The interparticle distance is defined in DualSPHysics as . Following the outcome of a convergence study (see Section 5.2), a value of = 0.022 m (= ∕4) will be used in the simulations. Note that the value of smoothing length ℎ will be then equal to ℎ=0.0458 m, being 2ℎ the interaction distance using the kernel function. The OWC WEC's front wall has a thickness of 0.25 m (see Fig. 1) resulting in a thickness of more than 10 particles, thereby ensuring that fluid particles on both sides of the walls do not influence each other. Since the OWC WEC chamber has a rectangular shape it is acceptable to start with numerical simulations in 2D.

Convergence study
Initially simulations are carried out in 2D and without PTO system. The free-surface elevation inside the OWC WEC chamber is measured and compared to the experimental one for CaseB-1.0D of Table 1 (undamped case without PTO system). A convergence study is carried out in order to choose a suitable interparticle distance for the simulations, using three different values of equal to ∕2, ∕4 and ∕8. Fig. 6 displays the results, together with the experimental data for the freesurface elevation at the centre of the OWC WEC chamber. In order to quantify the accuracy of the numerical results compared to the experimental values, an index of agreement is calculated using Eq. (27),

Table 2
Interparticle distance, number of particles, GPU runtime (GeForce GTX Titan Black) and index of agreement 1 for different resolutions, for CaseB-1.0D simulated in 2D.

Particles
Runtime with the value of the DualSPHysics prediction and the value of the experimental observation. In the current research the index of agreement is calculated using the data of the free-surface elevation. It is observed that an interparticle distance of = ∕4 provides accurate results, while still keeping the number of particles reasonably low, as also shown in Table 2. Therefore = ∕4 is the interparticle distance used in the simulations.

Compressible versus incompressible
Next, simulations are carried out both with incompressible and compressible air (in 2D) using a case with PTO system (CaseB-0.4D), where is computed using Eqs. (18) and (20) respectively. The results of the free-surface elevation are displayed in Fig. 7, comparing experimental results with the numerical ones.
The index of agreement 1 equals 0.85 for the simulations with compressible air and 0.71 for the simulations with incompressible air. Since 1 lies closer to 1 (which refers to 100% agreement) for the simulation with compressible air the following simulations in this manuscript will be carried out assuming that approach. This result was expected since modelling compressible air is the physically more correct assumption. In addition, the computational cost of simulations with compressible and incompressible air is the same. The difference in free-surface elevation and airflow speed through the orifice between simulations with compressible and incompressible air is expected to increase with scale (Dimakopoulos et al., 2017). It was shown in Simonetti et al. (2018) and Chen et al. (2021) that the air compressibility effect should be taken into account in large-scale models, whereas it can be neglected in small-scale models.

Validation
First, the cases without PTO system are analysed in terms of freesurface elevation only. Fig. 8 shows the free-surface elevation at the centre of the OWC WEC chamber computed with DualSPHysics in 2D and in 3D, and obtained from experiments for the undamped cases where there is no PTO system present (CaseA-1.0D and CaseB-1.0D). The index of agreement for each of the cases is given in Table 3.

Table 3
Index of agreement 1 for CaseA-1.0D and CaseB-1.0D in 2D and 3D.  Fig. 8 and Table 3, it can be concluded that more accurate results are obtained for CaseA-1.0D than for CaseB-1.0D. This can be related to the fact that the CaseB included additional uncertainties compared to CaseA, as mentioned in Bingham et al. (2021) where it is stated that the experimental free-surface elevation in an empty wave basin showed more reflection in CaseB than in CaseA. Furthermore, the 3D simulations show a better correspondence with experiments than the 2D simulations, indicating non-negligible 3D effects. When modelling an oscillating water column in 2D it has been shown that standing waves can form in front of the chamber, leading to overestimated performance of the OWC WEC (Sewan Park et al., 2019).
Next, the free-surface elevation and airflow speed through the orifice is calculated with DualSPHysics for the cases with PTO system damping: CaseA-0.4D, CaseB-0.4D, CaseA-0.5D and CaseB-0.5D (according to Table 1). These simulations are carried out both in 2D and in 3D for a more complete comparison. All indices of agreement (Eq. (27)) are calculated using the free-surface elevation at the centre of the OWC WEC chamber. Since now there is PTO system damping, the heave motion of the floating plate is used as the numerical value of the freesurface elevation at the centre of the OWC WEC chamber. The floating plate is restricted to translational motion, so, the heave of the plate is the same along the width of the OWC WEC chamber. Several instants of the same period simulating the CaseB-0.5D in 2D are depicted in Fig. 9. Horizontal velocity field (left panels) and pressure field (right panels) are displayed. Fig. 10 displays the 2D validation for the cases with PTO system damping. For all cases it is clear that the airflow speed through the orifice corresponds better with experiments for the positive values (that corresponds to air flowing out the chamber) than for the negative ones (air flowing inside the chamber). This is because the calibration constant appearing in Eq. (13) is calibrated only experimentally with the positive values of the airflow -see Bingham et al. (2021). Fig. 11 displays the 3D validation for the cases with PTO system damping. Results of 3D simulations show much better agreement than results from 2D simulations. In addition the airflow speed fits better with experimental data when its value is positive.
In Table 4 the average number of particles and the average runtime is given for the 2D and 3D cases, where the simulation of PTO force solved by Project Chrono takes less than 0.2% of the total runtime in the 3D cases. Since the number of particles and the runtime differs only slightly among the different cases, but depends mostly on the number of dimensions (2D or 3D), only average values are mentioned. In Table 5 the indices of agreement for the 2D and 3D cases are summarized. An    overall good agreement is achieved, but it can be concluded that the 3D simulations show better correspondence with the experiments than the 2D simulations. Therefore, the 3D effects are not negligible here and future numerical analysis conducted using this one-phase approach of DualSPHysics should be conducted in 3D. A possible explanation for the higher accuracy in 3D simulations is the fact that standing waves can occur in front of the OWC WEC chamber when simulating in 2D, as mentioned in Sewan Park et al. (2019).

Conclusions
In this study the numerical modelling of an Oscillating Water Column WEC with DualSPHysics has been investigated, considering both incompressible and compressible air inside the OWC WEC chamber. The numerical results of the free-surface elevation inside the OWC WEC chamber as well as the airflow speed through the orifice have been N. Quartier et al. Fig. 9. Instants of the same period of the horizontal velocity (left panels) and pressure (right panels) for CaseB-0.5D simulated in 2D.  validated with experimental data, provided by the tests performed at the KRISO.
The PTO system has been modelled in DualSPHysics by adding a plate floating on the free water surface inside the OWC WEC chamber, on which a PTO force acts implemented with the coupling between DualSPHysics and Project Chrono. Results show that DualSPHysics is a valuable tool for OWC WEC modelling, even when only a onephase solver is used. When considering the air as compressible, the Table 5 Index of agreement 1 for the 2D and 3D cases. motion of the free water surface inside the OWC WEC chamber agrees slightly better with experimental results than using the incompressible approach.
In order to reduce the computational cost, a multi-resolution technique (Vacondio et al., 2013;Barcarolo et al., 2014;Zhang et al., 2021) should be applied to this type of problems where the higher resolution N. Quartier et al. is employed to the area around the OWC WEC while less resolution will be used in the rest of the domain. This will be part of the future work.
Since the method developed in the current research allows accurate modelling of an onshore OWC WEC, the method can be extended to an application with floating OWC WECs. DualSPHysics is suited to model floating OWC WECs since it can handle extreme motions and the mooring line dynamics. However, disadvantages of the applied method in the current study is that the air pressure inside the OWC WEC chamber can only be calculated with theoretical formulas, and that the free-surface elevation is the same inside the entire OWC WEC chamber since the floating plate does not pitch.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.