Ocean Energy Systems Wave Energy Modeling Task 10.4: Numerical Modeling of a Fixed Oscillating Water Column

: This paper reports on an ongoing international effort to establish guidelines for numerical modeling of wave energy converters, initiated by the International Energy Agency Technology Collaboration Program for Ocean Energy Systems. Initial results for point absorbers were presented in previous work, and here we present results for a breakwater-mounted Oscillating Water Column (OWC) device. The experimental model is at scale 1:4 relative to a full-scale installation in a water depth of 12.8 m. The power-extracting air turbine is modeled by an oriﬁce plate of 1–2% of the internal chamber surface area. Measurements of chamber surface elevation, air ﬂow through the oriﬁce, and pressure difference across the oriﬁce are compared with numerical calculations using both weakly-nonlinear potential ﬂow theory and computational ﬂuid dynamics. Both compressible-and incompressible-ﬂow models are considered, and the effects of air compressibility are found to have a signiﬁcant inﬂuence on the motion of the internal chamber surface. Recommendations are made for reducing uncertainties in future experimental campaigns, which are critical to enable ﬁrm conclusions to be drawn about the relative accuracy of the numerical models. It is well-known that boundary element method solutions of the linear potential ﬂow problem (e.g., WAMIT) are singular at inﬁnite frequency when panels are placed directly on the free surface. This is problematic for time-domain solutions where the value of the added mass matrix at inﬁnite frequency is critical, especially for OWC chambers, which are modeled by zero-mass elements on the free surface. A straightforward rational procedure is described to replace ad-hoc solutions to this problem that have been proposed in the literature.


Introduction
Wave energy converters (WECs) represent a small, but potentially significant, segment of the global renewable energy mix. In order to be competitive with offshore wind or solar nonlinearity of about 15% of the maximum stable wave steepness. In reality, such a device would be exposed to irregular wave forcing, which complicates the analysis. In particular, for a weakly-nonlinear frequency-domain approach, a fixed wave-steepness value must be associated with each frequency component, which adds some uncertainty to the approach; whereas for CFD methods, properly resolving a large range of frequency components over a domain suited to the longer waves substantially increases the computational effort. This topic should be addressed in future work. Contributions to the study include several weakly-nonlinear potential flow solvers applying both incompressible-and compressibleflow models for the air phase, as well as several CFD solutions applying both compressibleand incompressible-flow assumptions. Unfortunately, we were not able to attract participants using fully-nonlinear potential flow. The calculations generally compare well with the experimental measurements for the amplitude of the primary-frequency, meanchamber motion response and the pressure difference across the orifice plate, as well as the mean absorbed power. Air compressibility is found to be relatively unimportant for the flow through the orifice plate, but important for the interaction between the air volume in the chamber and the response of the internal chamber surface. This introduces a phase shift between the pressure and the elevation, which is captured by the weakly-nonlinear time-domain model; however, only the fully-nonlinear, compressible-flow CFD model is able to fully capture the detailed nonlinear chamber response. This highlights the need for a multi-level analysis strategy, where fast potential flow solutions are used to map out the response and target specific conditions for detailed analysis by much more expensive CFD calculations. Whether a fully-nonlinear potential flow model, with a compressible air-flow treatment, could also capture these details is an interesting question for future investigation.

Experimental Measurements
The device was tested at a scale of 1:4 relative to a proposed full-scale deployment as part of a breakwater in a water depth of 12.8 m. The experimental measurements were carried out at the Korea Research Institute of Ships and Ocean Engineering (KRISO) ocean basin in 2019. Preliminary results along with a comparison to CFD calculations were presented in Reference [17]. The geometry and dimensions of the basin are shown in Figure 1, and the wave probe locations are indicated in Figure 2.  The wave conditions in terms of period T and height H are given in Table 1. As shown, three wave heights are used at each wave period, which are denoted by a WaveID of 'Low', 'Med' and 'High'. The wave steepness values for each condition are shown in Figure 3 (left), which plots H/λ, with λ the wavelength according to linear theory, along with the theoretical limiting value, H max , given by the parameterized data of Reference [18]. Figure 3 (right) shows the steepness as a percentage of the maximum wave height, plotted against the relative water depth, kh, where k = 2π/λ, showing that all of the wave conditions are nearly linear.
The geometry and dimensions of the chamber are shown in Figure 4. Also indicated in the figure are the dimensions of the air duct, which was fitted with orifice plates of three sizes to model the energy extraction from an air turbine. The chamber was fitted with nine wave gauges to measure the motion of the water column surface. Their locations are indicated in Figure 5, which also shows the locations of the six pressure sensors and the two hot wire anemometers that were installed in the duct. The wave gauges are capacitance-type gauges, manufactured in-house, with an accuracy of ±1 mm. The pressure gauges were connected to an Aplisens (Warsaw, Poland) model APR-2000ALW differential pressure transmitter that has a range of ±2000 Pa, with errors up to ±0.075% of the maximum differential pressure. The airflow speed through the duct was measured by two omnidirectional TSI (Shoreview, MN, USA) model number 8465 hot-wire anemometers. They have a range of 0 to 50 m/s with an uncertainty of up to ±2.5% of the maximum flow speed. All data were logged at 50 Hz.  . The test wave conditions in terms of nonlinearity, H/λ, relative water depth, kh, and the maximum possible wave steepness, H max /λ, as given by Reference [18].   Tests were performed with no orifice plate installed in the duct, and with three orifice plate sizes corresponding to 0.3D, 0.4D, and 0.5D, where D = 0.2 m is the duct internal diameter. Table 2 shows the conditions tested along with their TestID numbers. These conditions are meant to cover a wide range of possible operating conditions for the fullscale device. All of the tests shown in Table 2 were first performed without the chamber in place to establish the undisturbed incident wave conditions. For these tests, an additional probe was placed at the location of the center of the internal chamber, as indicated in Figure 5. Three examples of the measured incident wave signal in the open basin for Tests 106, 107, and 108 are shown in Figure 6, along with the chosen analysis section, and the target wave amplitude. The wavelengths here run from λ = 7.8 m for WaveID 02 to λ = 17.8 m for WaveID 08. Given that the length of the wave absorber is 9 m, we should expect wave reflections to become increasingly significant for the longer-wave cases. Also shown in this figure is an estimated time, t r , for each wave to travel the 74 m from the chamber to the end of the beach and back again, based on the linear group velocity, where c is the phase velocity, and ω = 2π/T is the wave frequency. Considering the combination of λ and t r for each case, WaveID 02 to 06 all provide a reliable, three to four wave period steady-state section for analysis without significant reflection effects.
To extract the wave amplitude components, we make a linear least-squares fit of the signal over the analysis section to a sum of sinusoidal components at the primary wave frequency and the first four higher harmonics. In contrast to using a discrete Fourier transform, this technique is very robust and requires very few wave periods to accurately extract the harmonic components. This is verified by checking the residual of the fit and ensuring that it is at least several orders of magnitude smaller than the signal. For WaveID 07 and 08, reflections are clearly apparent during the chosen analysis section, but we have chosen to include them in the comparisons anyway, keeping in mind that they include additional uncertainty. This is justified by the observation that the departure from the target amplitude is relatively small during the analysis section, and that these cases are the most energetic of the series, thus being of interest in terms of air compressibility and nonlinear effects.    Figure 7 shows the internal-chamber, surface-elevation measurements and the measured pressure difference across the orifice plate for Test 106. Here, the chamber elevation, η, represents the vertical position of the water surface, i.e., the projection of the raw wave gauge distance to the vertical by multiplication with sin θ, where θ = 33.69 • is the chamber slope angle. In Figure 7a, we see that the three probe signals along each horizontal line are nearly identical, indicating that there is no significant transverse motion of the chamber water surface. On the other hand, we can identify a significant front-to-back standing wave oscillation from the relatively large differences in amplitude of the signals along each row of probes. In Figure 7b, we plot the pressure difference between each pair of pressure probes, so that Gauge 1 represents DP1-DP2, etc. A positive pressure difference, p, represents a greater pressure on the chamber side of the orifice plate than on the outside. Ideally, we would like to know the pressure difference between the chamber and the atmosphere, but the pressure differences between the chamber and the duct (and the duct and the atmosphere) are relatively small compared to the drop across the orifice plate. This conclusion is motivated by the following analysis.
(a) Measured vertical chamber motions at the nine wave gauges shown in Figure 5 for Test 106.
(b) Measured pressure difference across the orifice plate at the three gauges shown in Figure 5 for Test 106. First, we consider the effects of air compressibility. As long as the pressure differences are small compared to the atmospheric pressure, p a , then a good approximation for the air density is given by where γ = 1.4 is the heat capacity ratio, ρ a = 1.225 kg/m 3 , and p a = 101, 325 Pa are the density and atmospheric pressure of air near sea level at 20 • C. Figure 8 shows the variation in relative air density for Test 118, as estimated from Equation (2) and the measured pressure drop across the orifice plate. This case represents the largest measured pressure from the series and shows a maximum density variation of about 1.5%. From Figure 8, we might be tempted to conclude that compressibility effects should not be very significant for these experiments. However, as discussed further in Section 5.2, while air compressibility is not very important for the flow through the orifice, it leads to significant effects from the large-volume air chamber on the response of the chamber water surface. As air flows from the chamber into the duct, it encounters a sudden pipe contraction, as discussed by Bullen et al. [19]. This is shown in Figure 9a, where A c is the area of the chamber water surface and U c , p 1 are the mean air flow velocity and pressure in the chamber, whereas A d is the area of the duct, with U d and p 2 the corresponding mean velocity and pressure in the duct. Both of these control points are located away from any local disturbances caused by the contraction, so that the flows at these points are essentially homogeneous. Ignoring compressibility effects, and following Reference [19], we can write the Bernoulli equation for the pressure as where K c is a pressure loss coefficient representing the head loss through the contraction.
Using this in Equation (3), we can express the pressure drop across the contraction as where σ = A d /A c is the area ratio, and we have introduced the coefficient Considering a typical case from these experiments, the peak chamber velocity is U c = 0.05 m/s, which gives a chamber Reynolds number of R = U c L 1 /ν ≈ 8500. Using the model derived in Reference [19] gives C d = 0.8356 for this case. This gives a maximum pressure drop of about 20 Pa for the most energetic tests, which is less than 1% of the maximum pressure drop across the orifice plate with a diameter of 0.4D. We can expect a similar effect as the flow passes from the duct to the atmosphere. Thus, for the cases that we focus on here, we can conclude that the measured pressure drop is very close to the pressure drop between the chamber and the atmosphere.  For the flow across the orifice plate, the situation is sketched in Figure 9b. Again, considering two control points located away from local disturbances due to the orifice plate, we can write the Bernoulli equation as where, in this case, σ = A o /A d is the ratio of the orifice area to the duct area and the total pressure loss across the orifice is expressed by the pressure loss coefficient, K c . From the experimental measurements, the estimated value corresponds to C d = 0.82. Using this value, and mass conservation, we can compare estimates of the volume flux in the chamber from each of the three measurements: 1. Q = A cū , whereū = ∂ tη , withη the mean of the nine internal wave gauge measurements and ∂ t the time derivative (here taken by a fourth-order finite-difference scheme). 2.
where u a is the measured air flow velocity through the orifice. Here, it was determined that the most accurate result was found by taking the value from the upper gauge for the outflow and from the lower gauge for the inflow.

3.
From Equation (6), we can write Q = 2(C d A o ) 2 /ρ |p| sign(p), wherep is the mean of the three pressure difference measurements.
Two examples are shown in Figure 10, where we can see reasonably good agreement between the predictions from the three measurements. The absorbed power is given by P = Q p, so there are nine possible ways to estimate absorbed power from the three experimental measurement quantities. In order to minimize the use of the empirical pressure-flux relation, Equation (6), we chose to focus on the three of these given by the previously described estimates of Q and the average of the measured pressures,p. Figure 11 shows two examples for Tests 103 and 106. The mean absorbed power, per wave cycle, is indicated by the horizontal lines.

Numerical Modeling of OWC Chambers
This section describes the numerical models that have been applied in this study to predict the KRISO OWC chamber response and power absorption. Three levels of modeling have been applied, corresponding to: weakly-nonlinear potential flow theory in the frequency domain, weakly-nonlinear potential flow theory in the time domain, and CFD. Both compressible-and incompressible-flow models are applied to the air phase. These models are outlined in the following subsections.

Weakly-Nonlinear Potential Flow Theory in the Frequency Domain
Detailed presentations of linearized, potential flow theory for the analysis of wavestructure interaction can be found in the classic textbooks by Newman, Falnes, or Faltinsen [20][21][22]. Here, we will summarize the theory and define the notation used. In order to apply radiation/diffraction theory to OWC chambers, there are two standard approaches. The first approach follows that outlined in References [21,[23][24][25], in which new degrees of freedom are introduced to represent the pressure distribution applied to the interior free surface of each chamber by the air turbine. This approach is implemented in the software WAMIT [26], under the label Free-Surface Pressure (FSP) modes, and it was applied to the analysis of a fixed OWC chamber by Reference [11]. A second approach, which is more straightforward to apply in the time domain, is to treat the interior free surface of the chamber as a massless element of the structural boundary, which is free to move in a set of predefined generalized modes. This approach is described in Reference [27] and is also implemented in WAMIT. When generalized modes are included, the equations of motion for a floating structure can be written as Here, ξ j (ω) is the generalized body-response phasor in 6 + M g degrees of freedom, where j = 1, 2, 3 correspond to translational motions in surge, sway, and heave, j = 4, 5, 6 represent rotational motions in roll, pitch, and yaw, and j = 7, 8, . . . , 6 + M g describe the motions of the interior chamber surface. Each generalized mode is defined by the boundary conditions where φ j is the radiation potential in mode j, w j (x) defines the vertical displacement of the interior free surface due to unit amplitude motion in mode j, S i is the internal free surface, S b is the submerged body surface, and ∂/∂n represents the derivative in the direction normal to S b . M jk is the linearized body inertia matrix, A jk (ω) and B jk (ω) are the radiation added mass and damping coefficient matrices, c jk is the hydrostatic restoring coefficient matrix, and X j (ω) is the diffraction exciting force coefficient vector. The assumption here is that the incident wave can be described by a superposition of linear Stokes waves at frequency ω of the form Here, {} indicates the real part of a complex quantity, x = [x, y] is a horizontal position vector, and the z-axis is oriented vertically upward, with z = 0 at the still-water level and z = −h at the fluid bottom. The free-surface elevation is η 0 and φ 0 is the velocity potential, with g the gravitational acceleration. The incident wave has a period T = 2π/ω, a length λ = 2π/k, an amplitude A = H/2, and it propagates in the direction defined by the angle β, measured from the positive x-axis. The dispersion relation, Equation (1), relates the wave period and length, and defines the wave phase and group velocities c and c g . All the coefficients appearing in Equation (7) can be computed for any desired floating structure using, for example, the frequency-domain radiation/diffraction code WAMIT [26], except for the matrix B 0 jk , which represents the damping applied to the chamber motion by the power-extracting air turbine, or in this case the orifice plate. The evaluation of B 0 jk is discussed in Section 3.3, and, because the applied damping is a nonlinear function of the chamber flux, its inclusion in the otherwise linear equations makes them weakly-nonlinear.

Weakly-Nonlinear Potential Flow Theory in the Time Domain
As originally noted by Cummins [28], the frequency-domain formulation described earlier has an equivalent statement in the time domain. A weakly-nonlinear version of this formulation can be written as where x j (t) is the time history of the motion response in mode j, A ∞ jk = A jk (∞) is the infinite frequency limit of the added mass, and the over-dots represent time differentiation.
The radiation impulse response functions, K jk , and the diffraction exciting forces, F jD , are related to the frequency-response functions by the Fourier transforms where η is a particular incident wave elevation measured at the origin of the coordinate system. Here, K jD is the diffraction impulse response function, and the two equivalent forms of F jD indicate how the wave excitation force can be computed in either the time or the frequency domains, when the incident wave elevation signal is known for all time.
Other external forces applied to the body, for example by the air turbine or a mooring system, are represented by F j0 (t) and can be nonlinear.

Potential Flow Modeling of the Orifice Plate Damping-Incompressible Flow
Under the assumption that air compressibility effects are negligible, Equation (6) defines the relation between the flux in the chamber, Q, and the pressure difference across the orifice plate, p, Defining generalized mode j = 7 to be the piston mode (uniform vertical motion of the chamber surface), the flux is Q = A cẋ7 , with A c the chamber free-surface area. The force applied by the orifice plate is then As long as all other mode shapes are defined to have zero mean value, they will not contribute to the flux, and thus not directly experience any forcing from the orifice plate. Equation (14) can be applied directly in the time-domain model, but in the frequency domain an equivalent linearized damping coefficient must be developed. For example, as shown in Reference [11], assuming a sinusoidal flux at frequency ω, and a linear relationship between the pressure and the flux gives as the equivalent linear damping coefficient that ensures the same power extraction per cycle as that defined by Equation (14). Because this damping coefficient depends on the chamber motion response, ξ 7 , the solution to Equation (7) is also weakly nonlinear and must be found iteratively for each combination of frequency and incident wave steepness, = 2A/λ.

Potential Flow Modeling of the Orifice Plate Damping-Compressible Air Volume with Incompressible Orifice Flow
As discussed in Section 2.1, the effects of air compressibility in these experiments are small but significant, suggesting that, at full scale, such compressive effects cannot be neglected. This is also discussed, for example, by Kelly et al. [29] and Sheng et al. [30]. In the current work, the influence of the air compressibility on the potential flow models of an OWC is considered from the perspective of the volume of air contained within the OWC above the water column and separately from the perspective of the air flow between the air volume and atmosphere through the orifice. A relationship between the pressure in the air chamber with respect to both the mass flow of air into or out of the chamber and the variation in the water column displacement, and in which the compressibilty of air within the chamber is considered, may be obtained from the First Law of Thermodynamics, and is given [31] where V is the volume of air in the chamber above the water column, and c s is the speed of sound in air. Assuming air is an ideal gas, c 2 s = γp/ρ, and the change in air density due to the compression of air may be obtained from Equation (2). The air flow direction, either into or out of the chamber, is dependent on the value of p. If the flow through the orifice is considered incompressible, then orifice theory may be used to model the mass flow of air into or out of the OWC chamber, as given in Equation (17) The force acting on the water column due to the pressure differential between the air trapped in the OWC chamber above the water column, and atmospheric pressure is F 70 = pA c . The pressure differential p is assumed to act uniformly over the surface of the water column and is treated as a state variable in the time-domain state-space model discussed in Section 4.3.

Potential Flow Modeling of the Orifice Plate Damping-Compressible Air Volume with Compressible Orifice Flow
If the flow of air through the orifice is also considered compressible, Equation (16) may still be used to model the pressure variation within the OWC chamber, while Equation (17) may be modified to allow for the expansion of gases flowing from a region of higher pressure to a region of lower pressure. Under unchoked conditions, the mass flow rate through an orifice may be given [32] bẏ where σ = A o /A d , p 1 , and ρ 1 are the pressure and density of the fluid in the region of higher air pressure, and p 2 is the fluid pressure in the region of lower air pressure.

Numerical Potential Flow Solutions
All of the participants in this study who used a potential flow assumption started by computing the coefficients appearing in Equation (7) using the high-order Boundary Element Method (BEM) solver, WAMIT [26]. Complete details on this software can be found in the WAMIT user manual and associated references on the web page. Analyzing OWC chambers with this technique presents some special challenges, which are discussed in the following subsections.
For the generalized modes representing the motion of the internal chamber surface (j > 6), the inertia terms, M jk , are all zero; therefore, it is of critical importance to have an accurate estimate of A ∞ jk in order to solve the problem in the time domain. Unfortunately, when the body has panels lying directly on the free surface, the integral equation is, in fact, singular [33]. Thus, the infinite-frequency problem for OWC-type modes cannot be solved in the usual way using a standard radiation/diffraction BEM, such as WAMIT. This failure has been noted by several authors and was studied in some detail by Sheng et al. [13], who suggested replacing the zero-mass (zero-thickness) internal free surface with a light piston of small finite displacement, which makes the calculation of A ∞ jk non-singular. A number of other authors have also used this approach; recent examples include References [14][15][16]. However, there are two problems with this approach. First, it is not obvious that the finite-mass floating piston is a good representation of the true free surface, nor how heavy the piston should be made to find a good balance between modeling the chamber motion and computing A ∞ jk accurately. Second, even if this gives a good representation of the piston mode, what about the other possible modes of motion for the internal free surface? In other words, it is surely possible, with some effort, to tune this approach to a particular chamber and get reliable results, but it is an ad-hoc and theoretically questionable method.
Instead of this ad-hoc approach, we suggest a more reliable and straightforward method. For a given OWC geometry, the finite-frequency values of the added mass and damping coefficients are computed as usual. Then, a new geometry is constructed by reflecting the original submerged geometry about the z = 0 plane. The patches/panels that describe the interior free surface are then redefined as dipole patches [26], and the body is submerged to the original water depth, h, while the depth is redefined to be 2h. Now the infinite (and zero) frequency problems can be solved accurately and reliably, and the results (divided by two) combined with the originally computed finite-frequency values. This strategy also ensures that A ∞ jk is correctly computed for all generalized modes defined on the chamber free surface, as well as the coupling between all rigid-body and generalized modes. Figure 12 shows the two KRISO chamber geometries used here. These geometries have been created using the analytic geometry feature of WAMIT by including a dedicated subroutine that is linked to the WAMIT executable at run-time. The body-coordinate system has its origin at the center of the interior chamber surface, with the x-axis pointing toward the incoming waves. The first patch of the geometry (shown in green in Figure 12) represents the internal chamber free surface. To model the motions of the chamber water column, a corresponding set of generalized modes has been added to the newmodes.f file to describe the vertical motions of the chamber surface. These modes are of the form where L c = 1.25 m is the x-length of the interior free surface, and x 0 is the position of the front wall of the chamber in the body-coordinate system. After comparing solutions using four modes and then two modes, it was determined that, over the range of wave frequencies of interest here, the higher-mode responses are very small and do not influence the first two modes of motion. Thus, we include only m = 0 (the piston mode) and m = 1, the first sloshing mode. These modes correspond to j = 7, 8, respectively. Figure 13 illustrates how the two-mode approximation of the internal surface motion relates to the geometry of the chamber. The solid black lines indicate the internal walls of the chamber, and the dashed red line shows the still water surface. Although, in reality, the water surface will follow the sloped chamber walls, linear theory assumes that the motion is vertical, following the dashed black lines. The dashed green lines indicate the front and back set of wave probes, which are offset by a distance of s l = 0.46875 m from the chamber centerline. The dotted blue line indicates the piston mode, and the solid blue line the sum of the piston and the first sloshing-wave modes. Three grids have been run to ensure converged results. These are labeled Grids 1-3 and correspond to setting PANEL_SIZE = [1, 0.5, 0.25] m. The resulting added mass and damping coefficients are shown in Figure 14, where we note that the ω = ∞ limit is plotted at ω L/g = 4, with L = 1 m. Here, we can see that the damping coefficient is rapidly convergent for all frequencies, whereas the high-frequency range of the added mass is actually divergent. This is due to the singularity discussed earlier in the BEM solution as ω → ∞. This singularity also leads to increasingly ill-conditioned linear systems at high, but finite frequency, resulting in the observed slow convergence there.
(b) Reflected geometry for ω = 0, ∞.   Convergence of the wave-exciting force coefficients and the piston-mode response is shown in Figure 15, where we can see that these coefficients are also rapidly convergent. Figure 16 compares the undamped (no orifice plate) WAMIT response amplitude operator with the response computed in the time domain via the inverse cosine transform of the damping coefficients and A ∞ jk . This confirms that the previously described strategy for computing the infinite-frequency limits of the added mass works as desired.

CFD Solutions
Recently, more computationally intensive modeling approaches (models using CFD) have become more widely applied due to the advantages of treating physical effects that are not easily captured by the aforementioned potential flow theory-based models [34]. In this paper, Unsteady Reynolds-Averaged Navier-Stokes (URANS) models are implemented from the open-source framework, OpenFOAM (versions 5.0, 7, and v1912), and the commercial code, StarCCM+ 13.06. The numerical models utilize the finite-volume method to discretize the RANS equations. The interface between the two fluid phases, including air and water flows, is tracked by a volume of fluid advection scheme (see, e.g., Reference [35]). Both the air phase and the water phase are treated as viscous, isothermal, and immiscible flows. Modeling has been done using both an incompressible flow assumption and allowing for the air phase to be compressible.

DTU_FD:
The Technical University of Denmark (DTU) weakly-nonlinear frequencydomain contributions solve Equation (7) using coefficients computed by the high-order BEM solver, WAMIT [26]. The computed coefficients are read into a post-processing program, and combined with an initial guess for B 0 77 at each value of ω and , to solve Equation (7) iteratively until the maximum change in ξ j is below a tolerance of 10 −6 . This gives the weakly-nonlinear, frequency-domain solution.
DTU_TD: The DTU weakly-nonlinear time-domain contributions solve Equation (11) using the open-source package, DTUMotionSimulator [36], which is described in the associated documentation. This code is based on the solution originally described by Bingham [37]. In brief, the Fast Fourier Transform (FFT) is used to compute the coefficients in Equation (12), and the explicit fourth-order Runge-Kutta integration scheme is applied to evolve the equations in time, with the fourth-order Simpson rule used to evaluate the convolution integral over the nonzero range of K jk . Equation (14) is used to compute the incompressible-flow model of the nonlinear power-take-off damping force at each stage of the time integration.

The National Renewable Energy Laboratory Team
NREL_WECSim: The National Renewable Energy Laboratory Team (NREL) WEC-Sim model was developed in MATLAB/SIMULINK and solves Equation (7) in the time domain [38]. In WEC-Sim, the rigid-body dynamics are solved using the multibody dynamics solver, Simscape Multibody, but the additional degrees of freedom from the generalized modes are formulated in a state-space form [39]. Similar to DTU_TD, the FFT is applied to compute the coefficients in Equation (12), and Equation (7) is solved in the time domain using the fourth-order Runge-Kutta integration scheme.
NREL_CFD: The high-fidelity CFD code, STAR-CCM+,is utilized to verify and validate the wave tank testing data, with the medium wave height case (WaveID = Med06) with TestID = 106 used as an example. An implicit, three-dimensional, incompressible, and unsteady RANS model is applied for all simulations. For the turbulence closure model, the Shear Stress Transport (SST) k-ω model, with "all y+ wall" treatment, is utilized. Herein, the "all y+ wall" treatment is a hybrid wall solution that attempts to combine high y+ wall treatment (y+ > 30) and low y+ wall treatment (y+ < 1). The free surface is modeled using the Eulerian multiphase volume of fluid method, utilizing typical fluid properties (ρ = 1000 kg/m 3 , ρ a = 1.225 kg/m 3 ).
The computational domain in the wave propagation direction (as shown in Figure 17) is dimensionalized similar to the experimental wave tank size with a total length of 58 m. Due to geometrical symmetry, only half of the computational domain width, set at 15 m, is used and its height is specified to be 9.7 m. In this study, a wave forcing model with a length of 14 m is applied to damp out wave reflections from the tank inlet and outlet boundaries. A velocity inlet boundary condition is specified at the tank inlet. A pressure outlet boundary condition, with a wave damping length of 6 m, is specified at the tank outlet. A wall boundary condition, without wave forcing, is specified at the bottom of the tank, whereas a pressure outlet boundary condition, without wave damping, is specified at the top of the tank. Additionally, symmetry boundary conditions are specified at both of the x-z planes.
The trimmed cell mesh technique is used to generate a high-quality grid for the fixed OWC WEC geometry with the orifice structure. In the wave-propagation region, the number of cells per wave length (TestID = 106) is specified to be 180, whereas the grid ratio at the water surface (∆x/∆z) is specified to be 8. Approximately, this grid ratio results in a number of cells per wave height of 9. Refined-mesh regions are specified at both the interior wave region and the external wave around the OWC body, to better capture the wave run-up behavior. In addition, a refined mesh is also generated around the orifice region. Near the wall surfaces of the entire WEC system, we generate seven layers of boundary layer mesh with a first layer thickness of 0.5 mm and a progression factor of 1.25, which leads to an average y+ value of 5. The maximum and minimum grid sizes for the OWC WEC system are 0.04 m and 0.00625 m, respectively. After running two-dimensional convergence studies, including grid size (e.g., base sizes = 0.014, 0.010, 0.007, 005 m), time-step size (e.g., ∆t = T/300, T/400, T/600, T/800), grid ratio (e.g., ∆x/∆z = 2, 4, 8), and inner iteration (e.g., 10, 15, 20, and 25), several three-dimensional computational domain convergence studies were conducted, including an evaluation of the horizontal grid ratios (∆x/∆z) of the interior-wave region and the time-step size. The final grid resolution used in this study is selected based on the aforementioned convergence studies, in terms of the trade-offs between numerical accuracy and computational cost. Particularly, the total number of cells used in the NREL study is 8.3 million. The simulations were run at NREL's high-performance computing center. Each simulation took an average of 16.3 h to complete and used 4700 cpu·hr.

1-DOF state-space model:
For the single-degree-of-freedom (1-DOF) case, two models were produced by the Maynooth University/Dundalk IT (MU/DkIT) Team. The first assumes incompressible flow through the orifice and uses Equation (17) to determine the mass flow rate of air between the chamber and atmosphere, while the second model assumes compressible flow through the orifice and uses Equation (18) to the same end. Both models use Equation (16) to model the pressure within the chamber, allowing a comparison to be made between the relative importance of the different compressible elements of the two models. Both 1-DOF models operate under the assumption that the water column acts in pumping mode only, referred to here as mode 7. The first sloshing mode, referred to as mode 8, results in an antisymmetric displacement of the free surface about a line parallel to the wave front running through the centroid of the free surface. Hence, there is no overall change in the air volume above the water column due to the sloshing mode, and it does not contribute to the power absorbed by the OWC; see Figure 13. The pumping mode is, therefore, the main 'power' mode of the OWC. The equation of motion of the water column in the KRISO device, for a single degree of freedom in the time domain, can be represented by a variation of Equation (11), also known as the modified Cummins equation [28]. The coupled equations Equations (11), (16), and (17) (or, in the case of the incompressible flow model, Equation (18)) are recast in state-space form. Such a form facilitates the replacement of the convolution term with a state-space approximation. The state variables are: with the state-space model for this 1-DOF case documented in Reference [31], containing a parametric approximation for the convolution term, conv 77 = t 0 K 77 (t − τ)ẋ 7 (τ)dτ, calculated through the FOAMM toolbox [40]. 2-DOF state-space model: Following a similar procedure as for the 1-DOF models, two 2-DOF models were produced. The first 2-DOF model assumes the air flow across the orifice to be incompressible and uses Equation (17) to model the mass flow rate of air between the OWC chamber and atmosphere, while the second model assumes the air flow across the orifice to be compressible, but uses Equation (18) in place of Equation (17). Both 2-DOF models represent the motion of the water column in the KRISO using a modified version of Equation (11) for both mode 7 and mode 8, Equation (16) to model the pressure variation within the chamber, and either Equation (17) or Equation (18) to model the mass flow rate through the orifice, as appropriate. The systems of equations are then cast in state-space form. Based on the behavior of the WAMIT coefficients for the chamber, the following approximations are found to be justified:

•
The sloshing mode (mode 8) and the coupling modes (mode 78 and mode 87) are waveless (i.e., the convolution integrals for these modes can be set to zero). • The coupling added masses, a 78 and a 87 , are independent of frequency and equal to their infinite frequency values, a ∞ 78 and a ∞ 87 , respectively. The chosen state variables for the 2-DOF state-space model are: The state-space model can be represented as: Using the parametric approximation of conv 77 given by the matrices: A r , B r and C r , as described for example in Reference [40], the final state-space matrices for the 2-DOF model are The 1-DOF and 2-DOF models are simulated in MATLAB as sets of ODEs in matrix form (state space), which are solved using the 'ODE45' solver based on an explicit Runge-Kutta integration formulation. The parametric approximations are used for both the radiation and excitation forces (Equation (12)) by using the methods described in References [40,41]. By generating both 1-DOF and 2-DOF models, the importance of the first sloshing mode on the performance of the model can be assessed.

Incompressible CFD model:
The open-source CFD code, OpenFOAM (version 4.1), is utilized in preliminary simulations of four of the physical wave tank test cases (Tests 107, 207, and 407). The code solves the unsteady, incompressible, RANS equations for two isothermal, immiscible fluids using a volume of fluid interface capturing scheme based on the multi-dimensional limiter for explicit solution (MULES) algorithm [42] and typical fluid properties (ρ water = 1000 kg/m 3 , ρ air = 1.225 kg/m 3 ). A variable time-stepping approach is used based on a maximum Courant number, Co, of 1.0. Pressure-velocity coupling is achieved via the Pressure Implicit with Splitting of Operators (PISO) algorithm [43] using one outer and three inner correctors. The k-ω SST (Shear Stress Transport) turbulence closure model is used with standard wall functions, for k, ω, and ν t , on all walls (including the surface of the OWC and the tank floor).
The computational domain used in this study is similar in dimension to the physical wave tank, but due to the geometrical symmetry of the problem, only half the width of the physical problem is simulated. The computational domain is, therefore, 15 m wide, 6.2 m tall, and 59.5 m long. Wave generation and absorption is achieved, using the waves2Foam toolbox [44], via the second-order Stokes theory, expression-based boundary conditions for the phase fraction and velocity on the upstream and downstream boundaries and 5-m-long relaxation zones adjacent to both (Figure 18). A symmetry boundary condition is applied to the boundary coincident with the symmetry plane of the problem and a pressure inlet/outlet boundary condition is applied to the top boundary. All other boundaries are modeled using a no-slip condition.
A convergence study was undertaken in which four different meshes were considered using the same variable time-stepping approach and maximum Courant number, Co, of 1.0. These meshes had background mesh sizes of 0.38, 0.25, 0.16, and 0.11 m. The selected computational mesh, shown in Figure 18, has an average background cell side length of ≈0.25 m, with grading of the cells in the along-crest direction for y > 3 m. The mesh has a number of regions of refinement, achieved using the quadtree refinement strategy and OpenFOAM's snappyHexMesh utility; the mesh in the free-surface region and on the OWC surface is refined to level 2 (≈0.06 m), whereas the duct and orifice regions are refined to level 4 (≈0.015 m). The surfaces of the duct and orifice are also refined with four prismatic cell layers, with a first layer thickness of 1.6 mm and an expansion ratio of 2. The mesh has 1.03 million cells in total. The simulations were run on the University of Plymouth high-performance computing cluster. Each simulation took approximately 44 h to complete and used 700 cpu.hr. Figure 18. OpenFOAM computational domain and grid refinement regions for the wave regimes, the duct, and the orifice.

Comparison of Results
In this section, we compare all of the simulation results from the participants, including the pressure difference across the orifice plate, the flow rate, estimated power, and internal chamber surface elevation, to the KRISO experimental data. The undisturbed incident wave elevation signals obtained from the experimental measurements were directly applied in the time-domain solutions, but the frequency domain and CFD simulations used the specified wave height and wave period given in Table 1. The analysis section was selected to minimize the influence of wave reflections from the beach on the measured data, as shown in Figure 6.
The measured and simulated first-harmonic, time-averaged amplitudes of the internal chamber surface elevation, orifice flow rate, and orifice pressure difference are equal to half of the peak-to-trough values, which are calculated using the time history data over the three wave periods in each analysis section. The power is estimated based on the flow rate and the pressure difference across the orifice plate, and the average absorbed power over the analysis section can be calculated from where T s and T r are the start and end times of the analysis section. For a fair comparison, the first-harmonic amplitude values and the averaged power output are scaled by the wave amplitude, which is equal to half of the wave height that was computed directly based on the peak and trough of the wave elevation time history from the analysis section.

Potential Flow Solutions
In this section, we focus mainly on the potential flow model results, but also include one set of CFD solutions provided by the University of Plymouth team. The pressure difference is created by the presence of the orifice plate. In the potential flow model, this is the pressure difference between the internal chamber and the atmosphere, while in the experiments and the CFD simulations, it is the pressure difference across the orifice plate, as shown in Figure 5. The measured and computed pressure difference across the orifice plate are plotted in Figure 19. We note that the incompressible-flow, weakly-nonlinear, timedomain solutions computed by the three teams (DTU, NREL, MU/DkIT) all gave results that are nearly graphically identical, so only one curve/point is visible for this method. This confirms that the different numerical implementations of this formulation are all reasonably converged and consistent. In general, there are only small differences between the incompressible-flow, time-and frequency-domain solutions, indicating a consistency between the two evaluations of the weakly-nonlinear pressure term. The solutions from the compressible models give slightly higher peaks than those predicted by the incompressible methods. For the compressible-flow, air-chamber models of the MU/DkIT team, very small differences are found between the compressible-and incompressible-flow models for the flow through the orifice plate; therefore, only the 1-DOF and 2-DOF fully compressible-flow values are shown, and only the 2-DOF fully compressible-flow solutions are presented in the time history plots. It is notable that the quadratic damping term in the time-domain potential flow models results in a nonlinear behavior of the signal near the zero-crossing, in the simulated time history. However, this nonlinear behavior was not observed in the KRISO experimental measurements. Figure 20 plots the measured and simulated flow rate from the potential flow solutions and the KRISO experimental data. In the experiments, by comparing to the time derivative of the internal chamber surface elevation, it was determined that it is more accurate to estimate the flow rate by taking the value from the upper gauge (HWA1) for the outflow and from the lower gauge (HWA2) for the inflow, as mentioned in Section 2.1. For consistency, the flow rate in the CFD simulations was computed using the same upper/lower gauge approach that was used in the experiments. In the potential flow models, the simulated flow rate depends on the coefficient C d , appearing in the pressure-flux relation of Equation (6). Based on the pressure difference and the 'outflow' flow rate obtained from the experimental data, the value of C d was estimated to be 0.8271. Since the C d value was tuned based on the experimental measurements, all codes with different modeling approaches generally agree with the estimated flow rate from the KRISO data. Because the C d value was estimated using only the outflow, the potential flow solutions have a better agreement with the KRISO data for the outflow cycle, when the flow rate is positive. Figure 21 plots the measured and simulated internal surface elevations at the center of the chamber. The predicted internal chamber surface elevation amplitude is similar to the measurements, except for the high wave-height cases, where all the potential flow models underestimate the surface elevation. However, the nonlinear behavior of chamber surface elevation at the peaks, which is apparent in the experiments, is not captured by any of the potential flow models. In addition, a phase shift is observed between the compressible and incompressible solutions, in particular for Test 106. This is one expected effect of the compressibility of the air phase in the chamber, since it acts like a mechanical spring which is in phase with the chamber motion rather than the chamber velocity. This introduces a phase shift compared to the incompressible model where the pressure and the chamber surface velocity are exactly in phase. Figure 22 shows the average absorbed power, which is calculated from the measured and simulated data. In the time history plots, the high and low peaks observed in the measured data are mainly caused by the difference of flow rate between the outflow and inflow, which is not reproduced by the potential flow models because the same C d value is used for both the outflow and inflow in the simulations. Thus, it also follows that the estimated power outputs from the potential flow models agree better with the experiments for the outflow cycle.     To summarize these results, we calculated the maximum percentage difference between the KRISO measurements and the solutions from the potential flow models. The percentage difference is determined based on the maximum value of the pressure difference, chamber surface elevation, flow rate, and averaged power from the test cases, and the results are listed in Table 3. Overall, the differences between the potential flow model solutions and experimental data are within 21%, except for the 1-DOF model, which only considers the piston mode for the chamber surface motion. Although all of the results show a similar trend, with relatively small differences in the predicted first-harmonic amplitudes, only the compressible-flow model is able to capture the small phase shift in the chamber surface elevation relative to the pressure. Table 3. A summary of the difference between the measured and potential flow solutions.

CFD Simulations
As discussed in Sections 2.1 and 5.1, the expected variation of air density due to air compressibility inside the chamber is small for these experiments. However, while the influence of air compressibility on the flow through the orifice plate was found to be very small, a small but significant phase shift between the chamber surface elevation and the pressure was accurately captured by the compressible model for the air volume in the chamber. To further investigate the significance of nonlinearity in the system, verify the numerical solutions and better understand the effect of air compressibility, both incompressible and compressible CFD simulations were performed for Test 106. This case was chosen because it is the largest-amplitude case which is still free of reflection effects from the beach.
The measured and simulated pressure differences across the orifice plate for Test 106 are plotted in Figure 23. In addition to the CFD solutions, the result from the incompressible and compressible potential flow time-domain models (NREL_WECSim and MU's 2-DOF), with both C d = 0.65 and C d = 0.8271, are also included. The quadratic nonlinearity of the signal near the zero-crossing observed in the potential flow model is perfectly captured in the incompressible CFD simulation. However, this nonlinear behavior is not observed in the compressible CFD simulation, and the effect of air compressibility then seems to explain why this nonlinear behavior near the zero crossing is not observed in the experimental data. The compressible CFD simulation also captures the small phase shift of the peaks observed in the KRISO experiments. In general, the effect of air compressibility has a small influence on the predicted pressure difference. The amplitude of the pressure difference across the orifice obtained from the compressible CFD simulation is slightly higher than that predicted by the incompressible CFD simulation, which is consistent with the trend shown by the potential flow calculations.
On the other hand, estimation of the flow rate is the primary source of uncertainty in the experiments. To better understand this, we compare the measured flow rate to that obtained from the CFD simulations and the potential flow models. The results are plotted in Figure 24. Here, two different approaches are used to estimate the experimental flow rate. One is based on the previously described upper/lower air velocity gauge approach, and the other is obtained by taking the average of the two measurements. In both the incompressible and compressible CFD simulations, we directly compute the flow rate by integrating the flow velocity across a surface plane at the orifice. The predicted flow rate from the CFD simulations is closer to the experimental data using the average value from the two velocity gauges and the potential flow solution when C d = 0.65 is used. C d represents a time-averaged value, which is calculated using the correlation between the pressure difference and flow rate, as described in Section 2.1, and the time history data of the last three wave periods from the compressible CFD simulation. For the effect of air compressibility, which is similar to what we have found for the pressure difference, the peaks are slightly shifted, and the flow rate magnitude obtained using the compressible CFD simulation is slightly higher than the prediction from the incompressible CFD simulation.  To better quantify the flow rate calculation uncertainty, we compare the computed flow rate from the compressible CFD simulations using different methods including: the upper/lower gauge approach; the average value from the upper and lower gauges; the velocity at the center of the orifice; and, finally, integration of the flow velocity over a surface plane inside the duct or at the orifice, which should give the true flow rate. This comparison is presented in Figure 25, and the solution from the incompressible CFD simulation is also included for reference. Based on the results from the CFD simulations, the flow rate is overestimated when using the upper/lower gauge approach. To further investigate the differences, we plot the velocity magnitude at the upper (HWA1) and lower (HWA2) gauges in Figure 26. The CFD solutions are close to the measurements. However, the measured downstream velocity magnitude is smaller than the CFD solutions, and the measured upstream velocity magnitude is larger than the CFD results, particularly for HWA2. Moreover, the measured velocity magnitude has a higher peak value for the outflow than the inflow, whereas the calculated values are nearly the same. In addition, as shown in Figure 26, the calculations predict that the velocity near the orifice is very sensitive to the vertical location, and varies almost linearly. This probably explains why the estimated flow rate is smaller if calculated from the average value of the two velocity gauges and why there is a significant difference between the outflow flow rate and the inflow flow rate.  (c) Velocity contour on the orifice surface plane.
(d) Point measurements between HWA1 and HWA2. Typically, compressible flow is considered as being relevant for a flow velocity with a Mach number greater than 0.3. For Test 106, the Mach number is about 0.1, using the maximum flow velocity near the orifice from the CFD simulations. Based on the estimated Mach number, and the calculated air density variation, as described in Section 2.1, we initially assumed that air compressibility effects would be insignificant for these experiments. This is the case for the pressure difference and the flow rate. However, the effect of air compressibility has a significant influence on the chamber surface elevation. A comparison of the internal chamber surface elevation at the center of the chamber (wave gauge #5), using both the compressible and incompressible CFD simulations, is shown in Figure 27. The predicted chamber elevation from compressible CFD is slightly smaller but captures the nonlinear time history profile fairly well. There is a significant difference between the solutions from compressible and incompressible CFD simulations. The air compressibility also introduces a phase shift in the chamber surface elevation. In addition, like the measured data, the compressible flow CFD results show a larger time derivative near the zero-crossing compared to the incompressible CFD solution, which is the main reason why the estimated value for the peak flow rate is larger when using the chamber wave elevation for the calculation. Figure 28 shows a comparison of the measured and simulated power outputs, which are directly calculated based on the pressure difference ( Figure 23) and the flow rate ( Figure 25). Note that the experimental pressure difference is obtained from direct measurement, regardless of how the flow rate is estimated. Different ways of estimating the flow rate result in different power estimates. For Test 106, the estimated power from the compressible CFD simulation is about 13% higher than the prediction from the incompressible CFD simulation, and the estimated power from the experiment using the upper/lower gauge flow rate calculation approach is about 15% higher than the estimated power from the compressible CFD simulation. The latter is essential because the time derivative of the chamber wave elevation is often used to estimate the flow rate, when the flow velocity through the orifice is not available. Focusing on the compressible models and the measurements here, we can see the sensitivity of the power estimate to the assumed flow-rate through the orifice and the chosen value of C d , illustrating the importance of getting these values accurately from the experimental conditions. The results from this study indicate that this approach can be misleading when the effect of air compressibility is non-negligible. It is anticipated that, for smaller orifice cases, the difference is most likely to be larger because of the increase of flow velocity and a larger Mach number, where the effect of air compressibility becomes even more important.

Conclusions
A set of relatively large-scale measurements for a breakwater-mounted OWC chamber has been carefully analyzed and compared with numerical predictions based on weaklynonlinear time-and frequency-domain potential flow methods, and CFD. The air phase has been modeled using both incompressible-and compressible-flow models. The effects of air compressibility are found to be small for the pressure difference across the orifice plate, but they produce significant nonlinear effects and a significant phase shift of the internal chamber surface motion. All models capture the time-averaged quantities well, but only the compressible CFD model is able to accurately reproduce the detailed responses.
This study highlights a number of important guidelines for the production of benchmark experimental data suitable for comparing the relative accuracy of different levels of numerical modeling. Ideally, benchmark data should include repetition of each case, so that the uncertainty of all measured quantities can be quantified to produce reliable confidence bounds. This begins with the undisturbed wave conditions in the basin without the structure in place, and then all measured quantities with the structure in place. Once the realization uncertainties are established, they can be combined with measurement errors to establish 95% confidence intervals, which then allow the data to be used to precisely evaluate the accuracy of the numerical models. For OWC-type devices, special care needs to be taken to determine the pressure-flux relation for the flow through the orifice plate, or the model turbine, which is used to extract power from the system. Ideally, the pressure difference should be measured between the chamber and the atmosphere to avoid local effects near the orifice/turbine. Finally, we have illustrated a simple and accurate method for using state-of-theart BEM solutions to accurately determine the infinite frequency added mass for OWC chambers. This method should replace ad-hoc approaches which are commonly applied to address this problem.