Numerical analysis of wave–structure interaction of regular waves with surface-piercing inclined plates

The Mocean wave energy converter consists of two sections, hinged at a central location, allowing the device to convert energy from the relative pitching motion of the sections. In a simplified form, the scattering problem for the device can be modelled as monochromatic waves incident upon a thin, inclined, surface-piercing plate of length L′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L'$$\end{document} in a finite depth d′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d'$$\end{document} of water. In this paper, the flow past such a plate is solved using a Boundary Element Method (BEM) and Computational Fluid Dynamics (CFD). While the BEM solution is based on linear potential flow theory, CFD directly solves the Navier–Stokes equations. Problems of this type are known to exhibit near-perfect reflection (indicated by a reflection coefficient |R|≈1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|R|\approx 1$$\end{document}) of waves at specific wavenumbers k′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k'$$\end{document}. In this paper, we show that the resonant motion of the fluid induces large hydrodynamic forces on the plate. Furthermore, we argue that this low-frequency resonance resembles Helmholtz resonance, and that Mocean’s device being able to tune to these low frequencies does not act like an attenuator. For the case where the water is deep (d′>λ′/2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d'>\lambda '/2$$\end{document}, where λ′=2π/k′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda '=2\pi /k'$$\end{document}), we find excellent agreement between our simulations and previous semi-analytical studies on the value of the resonant wave periods in deep water. We also find excellent agreement between the excitation forces on the plate computed using the BEM model, analytical results, and CFD for large inclination angles (α>45∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha > 45^\circ $$\end{document}). For α≤15∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \le 15^\circ $$\end{document}, both methods show the same trend, but the CFD predicts a significantly smaller peak in the excitation force compared with BEM, which we attribute to non-linear effects such as the non-linear Froude–Krylov force


Introduction
In recent years, there has been a renewed interest in the R&D of wave energy conversion technologies in the hope of commercialisation of the technology. This is evidenced by the various Wave Energy Converter (WEC) developers who are in The WEC is slack moored to the seabed at the bottom of the forward section, which lines it up to meet the incoming waves, similar to a hinged raft. At the end of each section, the device features so-called "wave channels" (see Fig. 1b). The resulting fluid motion inside the wave channels produces a resonant peak in the water surface motion and the corresponding hydrodynamic forces, which is beneficial to power capture. The geometry of Mocean's device resulted from an optimisation of the hydrodynamic forces on the device to give optimum energy production (McNatt and Retzler 2020a). To distinguish the operating principle of Mocean's WEC from a simple hinged raft, we consider the power absorption of the device using the Relative Capture Width (RCW, see McNatt et al. 2020), defined as where P is the power capture (response amplitude operator), r 0 = 3 3V /4π is the equivalent radius if the wetted volume V were a sphere and e f is the unit energy flux of the wave, defined as e f = 0.5ρ g c g where g is the gravitational acceleration, ρ is the density of water, and c g is the incident wave's group velocity. The RCW of a Mocean WEC and a simple hinged raft both of length, = 21 m and equal submerged volumes (corresponding to r 0 = 3.11 m) and equal widths (W = 2.9 m) are shown in Fig. 1d for wavelengths λ normalised by . There are several notable differences between Mocean's device and a simple hinged raft. The peak in the RCW for Mocean's device is four times that of a simple hinged raft, it has a broader banded response, and it peaks at wavelengths longer than the length of the device. Evidently, the device has a peak in the RCW at a wavelength in excess of six times the wave channel length L . In some designs of the device, the wave channels have sidewalls, so one might think that the resonance is connected to sloshing motion between the vertical walls. Sloshing-induced resonant behaviour has been observed to occur for a different class of WEC known as an oscillating wave surge converter (Renzi and Dias 2012). If they were present, we would expect sloshing modes to occur at incident wavelengths comparable to the characteristic width of the wave channel. However, the resonant response of the device can occur at wavelengths that are two-to-three times longer than the overall length of the machine (McNatt and Retzler 2020a) and an order of magnitude longer than the wave channel length or width. Hence, this resonance appears to be unrelated to the sloshing modes of the fluid contained within the wave channels, but a satisfying description of this resonance has yet to be established. Mocean's device does not fit neatly into an established classification of WEC (attenuator, point absorber, terminator, or flap type) due to the presence of the wave channel resonance. Instead, to accurately capture the working prin-ciple of the device, we argue that the device experiences the low-frequency resonant response similar to the Helmholtz Resonance (HR) of a sloping harbour (Mattioli 1978). The Helmholtz "cavity" of Mocean's device is the fluid wedge residing vertically above the inclined plate (see Fig. 1c) and the external forcing is the incoming waves. We will show that this description holds for small angles of a two-dimensional (2D) plate, i.e., when the vertical dimension of the cavity's opening (or mouth) is small compared with the plate's length (i.e., the characteristic size of the cavity). This situation resembles the HR of an infinitely long filled cavity that communicates with the external flow through a narrow and infinitely long slit, which has been considered previously by Bigg and Tuck (1982). The low-frequency resonance of the free surface within the wave channels of Mocean's device (at wavelengths an order of magnitude longer than the length of the wave channels) is also the most energetically favourable of all the resonant modes, which is consistent with the observations of HR within cavities such as harbours (Maas 1997).
Mathematical models that rely on the validity of linearity of the waves are often prone to fail in estimating the power capture of such resonant systems (Cummins and Dias 2017;Tan et al. 2020). This is due in part to the large-wave amplification near to the device resulting in additional forces such as the non-linear Froude-Krylov force. To obtain sensible power capture estimates, we need to develop better models that do not rely on the validity of this assumption. As a first step, here, we will consider the case of a 2D inclined plate. Although this is not a complete representation of Mocean's 3D device shown in Fig. 1, it does allow us to explore the nature of the resonance in much greater detail. Indeed, our preliminary studies have shown that the resonance observed in the full 3D problem does not differ qualitatively from the 2D problem. Rather, the resonant period T of oscillation is simply shifted towards shorter periods in the full problem. Furthermore, in 2D, this resonance has been examined analytically previously, which provides us with robust verification data for our numerical results. Once we have identified the working principle of the device in 2D, future work can identify more precisely how this translates into the full 3D problem.
In a simplified 2D form, the device involves a thin, inclined, surface-piercing plate of length L in a finite depth d of water that faces the incident waves (see Fig. 1c for an early iteration of the device). The problem of water waves travelling past a plate has been considered by many researchers over the past century. The first study on the problem considered 2D monochromatic waves of amplitude A 0 and wavelength λ travelling past a vertical plate with its top edge submerged a depth s and bottom edge submerged to a depth L beneath the free surface and under the assumption of small-amplitude waves ( A 0 λ ) where the water was considered to be of infinite depth (d > λ /2) (Dean 1945) Fig. 2). Later, the case of a surface-piercing plate (i.e., s = 0) was considered by Ursell (1947); here, the author was interested in the reflection coefficient, |R| (the ratio of wave heights of the reflected and incident waves far from the plate), and he derives an analytical expression for |R| in terms of k L , where k = 2π/λ is the wavenumber of the incident wave. Various extensions of the problem were considered in the following decades. In the 1970 s, the case of s ≤ 0 was considered afresh by Evans (1970) for deep water, who derived expressions for the forces and moments (both first and second order) on the plate. For the surface-piercing case (s = 0), the author recovered the results in the early studies (Ursell 1947;Haskind 1959). The problem of an inclined plate in deep water was first considered by John (1948) and later for the case of a nearly vertical plate by Shaw (1985). Let α denote the angle that the plate makes with respect to the water surface (see Figs. 1c and 3), so that a vertical plate corresponds to α = π/2. The case of α = 0 corresponds to a subproblem of the so-called dock problems (Friedrichs and Lewy 1948;Garrett 1971;Linton 2001); for which it is known that the reflection coefficient increases monotonically as a function of k L .
For small (but non-zero) α, the reflection coefficient for the deep-water problem exhibits a 'spiky' behaviour (Parsons and Martin 1994); indicating that at certain frequencies, the plate completely reflects the incoming waves, and with an overall trend of increasing reflection as k L increases. The authors attribute this effect to a so-called quasi-resonance between the wedge of fluid trapped above the plate and the incoming waves (Parsons and Martin 1994). Later, Parsons and McIver (1999) derived the following equation for the resonant period (with n = 0, 1, 2, . . . the mode number): The resonance originates as a balance of the energy flux into the wedge's mouth and the change in the potential energy in the wedge's free surface. Note that the n = 0 resonance is low frequency: it satisfies k 0 L 1, where k 0 is the wave number associated with the n = 0 mode so that the resonant wavelength is significantly longer than the cavity's length L , while the cavity's length is much larger than that of the mouth. In this paper, we show that this lowfrequency resonance shares many characteristics with HR. The similarities are best understood in analogy to the wellstudied topic of resonance in harbours.

Helmholtz and sloshing modes of a rectangular basin
Consider the case of free longitudinal oscillations in a rectangular basin of length L and depth D such as that shown in Fig. 4a,b; it can be shown that the wavelength λ n of the nth mode (normalised by L ) is given by for the case of a closed rectangular basin (see Fig. 4a). Note that the basin's width does not feature in this, since this formula assumes that the motion is 2D. The corresponding wave periods T n (normalised by 2L / g D ) can be estimated using Merian's formula (Miles 1974) T n = T n g D /2L = 1 n , n = 1, 2, 3, . . . .
For the lowest frequency (n = 1) mode, the basin is said to be acting as a half-wave oscillator (Korgen 1995). Now, consider the effect of opening one end of the basin (Fig. 4b).
Notice that there is an n = 0 mode present in the openbasin case that is not present in the closed-basin case. The n = 0 resonance is also known as the Helmholtz mode, the zeroth mode, the gravest mode, and the pumping mode in the harbour literature (Rabinovich 2017). Unlike the closedbasin case, the lowest frequency (Helmholtz) mode of an open basin is said to behave as a quarter-wave oscillator. We can imagine continuously transforming from the open-to closed-basin case by reducing the size of the opening of the basin (Fig. 4c). As the channel's opening becomes smaller, the Helmholtz mode ultimately degenerates into a mode of zero frequency and uniform displacement (Murty 1977) for the closed basin. In this way, it is distinct from the sloshing modes as pointed out by the previous authors (Miles 1974).
In analogy with the rectangular basin ( Fig. 4a-c), the inclined plate at small angles of inclination corresponds to a basin of slowly varying depth D ≈ L sin α (see Fig. 4d). Notice that if we normalise the period of resonance from Parsons and McIver (1999) in the same way that we did for (5) and (7), we find that the normalised version of (2) reads We can immediately see that the period of the n = 0 mode of the open rectangular basin T 0 = 2 is three quarters that of the fundamental mode for the inclined-plate problem T 0 = 8/3. A similar shift towards longer periods happens for harbours, where sloped harbours (with their deepest point at the harbour's mouth) resonate at longer period waves than rectangular basins with the same mouth depth (Mattioli 1978). We hypothesise that the n = 0 resonance described by (8) is analogous to the n = 0 (Helmholtz) resonance observed in a gently sloping harbour (Mattioli 1978). We remark here that the corresponding 3D problem for (8) is a sufficiently wide plate, while the associated 3D problem for (6) is a sufficiently narrow basin. Mathematically, however, both (8) and (6) are derived within the same 2D framework, which allows us to sensibly draw an analogy between them. We mention for completeness that there have been numerous investigations of the case of small-amplitude waves scattering from a submerged plate. These studies include scattering by a nearly vertical plate (Mandal and Kundu 1990) and oblique scattering (Mandal and Das 1996) in deep water. In finite-depth water, the scattering by a thick vertical barrier has been considered for the case of both partially submerged and surface-piercing barriers (Kanoria et al. 1999). The case of a fully submerged inclined plate was considered by Midya et al. (2001). More recently, there have been experimental tests of inclined plates (Yagci et al. 2014), which could be useful for breakwaters used in coastal defences (Cho and Kim 2008;Rao et al. 2009). Despite the wide variety of studies investigating the inclined-plate problem, to the authors' knowledge, there has yet to be a comparison of these studies with CFD. 1 In this paper, we model the problem using linear potential flow, solved with the commercial BEM code WAMIT (WAMIT 2016), where we recover some classical results and we compare this directly to our results from CFD. We perform our CFD simulations for L = 0.1 m, corresponding to a 20th scale wave energy device (McNatt and Retzler 2020a), and in water of depth d = 4L = 0.4 m. Our CFD modelling is capable of reproducing the phenomenon of perfect reflection at specific values of k L found in the previous studies that relied on the validity of linear potential flow. By limiting the study to small-amplitude linear conditions, the results allow a verification of our CFD simulations with analytical solutions, thereby proving the robustness of the numerical model and providing a reliable computational framework to understand the non-linear resonant response of WECs of this kind in the future. We further argue that the spikes in the excitation forces on Mocean's device correspond to the 3D counterpart of the resonance described by (8). In particular, focussing on the n = 0 mode of (8), we will argue that Mocean's device exhibits HR, similar to that found in gently sloping harbours (Mattioli 1978).
The paper is organised as follows: in Sect. 2, the governing equations are presented and non-dimensionalised. The underlying theory and numerical schemes are outlined for CFD and BEM in Sects. 2.1 and 2.2, respectively. In Sect. 3, we discuss computational considerations and convergence studies, and in Sect. 4, we outline the specifics of the case studied in the paper. The results are presented in Sect. 5 for the excitation forces in Sect. 5.1 and the reflection coefficient in Sect. 5.2. Finally, we discuss our findings in Sect. 6, in particular with reference to an application to modern-day WECs.

Mathematical model
We define the reference system of coordinates O (x , y , z ) with x pointing in direction of the incoming waves, and let the y -axis lies along the width of the plate and let the z -axis points vertically upwards from the still water level (see Fig.  5). The origin O is located in the center of the plate at the still water level.
The governing equations to describe wave-structure interaction (WSI) are based on the conservation of mass and momentum, expressed in the form of the Navier-Stokes equations for incompressible flows as follows: and the Reynolds and Froude number Re = ρ U L /μ and Fr = U / g L , respectively. Here, ρ is the density of the fluid mixture (defined later in (13)), μ is the dynamic viscosity of the fluid mixture (defined later in (14)), and U is a characteristic velocity scale. n g denotes the unit vector in the direction of the gravitational acceleration. The non-dimensional force acting on the submerged plate, F = F /ρ U 2 L 2 , is calculated as the integral over the plate surface of the normal and shear stresses, as follows: where I is the unit tensor.
where F x , F y , and F z are the surge, sway, and heave forces acting on the plate, and n is a unit normal vector on the plate's surface.

Computational fluid dynamics
Due to the presence of non-linearities, no general analytical solutions are available for the governing conservation equations (10) and (11), which necessitates the development of a numerical solution. For this study, the open-source CFD toolbox OpenFOAM (Weller et al. 1998) is employed. Open-FOAM implements the finite-volume method for the spatial discretisation of the governing equations. An implicit Euler scheme is used for the temporal discretisation in this study, following Windt et al. (2019); Larsen et al. (2019).
To capture the multi-phase nature of the surface-piercing plate problem, the volume-of-fluid (VOF) method (Hirt and Nichols 1981) is implemented in OpenFOAM. In the VOF method, the volume fraction α VF is introduced, such that the governing equations of the fluid mixture can be developed for a single fluid. α VF is limited to the range 0 ≤ α VF ≤ 1, where α VF = 0 represents the air phase and α VF = 1 represents the water phase. The properties of the fluid mixture (density ρ and viscosity μ ) can be expressed as and where a subscript indicates the phase of the density and viscosity. In this study, the following values for the air and water density are used: ρ water = 1000 kgm −3 and ρ air = 1 kgm −3 . The viscosity for the air and water is set to: μ water = 1 × 10 −3 kgm −1 s −1 and μ air = 1.48 × 10 −5 kgm −1 s −1 . The transport equation for α VF is as follows: Throughout this study, the common approach is employed whereby α VF = 0.5 defines the location of the free-surface interface. To retain sharp water-air interfaces, relatively fine spatial discretisation in the interface region is required, and can be supported by the numerical interface capturing method. To that end, in OpenFOAM, a compression term (16) where the compression velocity u r is determined as representing the relative velocity between the liquid and the gaseous phase. In Eq. (17), c r controls the interface compression (Devolder et al. 2019) and is set equal to 1 throughout the simulations for this study. The MUlti-dimensional Limiter for Explicit Solution (MULES) algorithm (OpenFOAM 2020) is employed herein to discretise equation (16), ensuring boundedness and conservativeness. For more information on the performance of the interface compression method, the interested reader is referred to Deshpande et al. (2012).

Numerical wave generation
Different methodologies are available to implement numerical wave generation and absorption in CFD-based Numerical Wave Tanks (NWTs) (Miquel et al. 2018;Windt et al. 2019). Throughout this study, the relaxation zone method waves2Foam, proposed by Jacobsen et al. (2012), is employed, which exhibits good overall performance for various wave conditions (Windt et al. 2019). In addition, it is paramount to avoid re-reflections from the wave generation boundary in our study, since we are operating in a 2D domain. This can be achieved with the relaxation zone method (see Windt et al. (2019)). It should be pointed out that the relaxation zone method comes at a relatively high computational cost (e.g., compared to a static boundary method). However, for this study, the additional computational cost is a reasonable trade-off for the achieved accuracy.
In the relaxation zone method, a target solution, γ target , is blended with the computed solution, γ computed , for field for u and α V F . In the wave absorption relaxation zone (of length L a ), to the right of the simulation zone, waves are absorbed by setting the target solution for the velocity to zero and free-surface elevation to the still water level variables such as γ = u or γ = α VF , within pre-defined relaxation zones (see Fig. 6), via χ R represents a weighting function which varies smoothly along the relaxation zone with boundary values of zero at the NWT boundary and unity at the interface with the simulation zone, thereby ensuring a gradual transition in the blending of the target and computed solutions. The target solutions for u target and α VF,target at the wave generation boundary are obtained from wave theory. At the wave absorption boundary, u target is zero, and α VF,target = 0.5 defines the location of the still water level.

Linear potential flow theory
We can significantly reduce the computational overhead during the numerical solution process by representing the fluid flow as inviscid and irrotational and assuming smallamplitude waves. Under these assumptions, many WSI problems in engineering can be solved rapidly. It is well known that these assumptions break down in resonant conditions, where there can be significant viscous losses (as indicated by a Keulegan-Carpenter number K C greater than unity) (Cummins and Dias 2017) or significant non-linear free-surface effects such as the non-linear Froude-Krylov force. This can result in the theory giving overestimates of the power capture of wave energy devices in such conditions. In the present problem, the resonance relates to the resonant period of the wedge of fluid above the inclined plate and, in resonance, large free-surface elevations (A 0 ≥ λ ) can occur, which are typically overpredicted by linear inviscid potential flow. In this paper, we will see that the free-surface elevation can be four times that of the incident wave height, resulting in a K C ≈ 0.5, which indicates that viscous effects may be appreciable in such resonant conditions. Assuming an irrotational fluid flow, the velocity can be defined through the velocity potential, (x , y , z , t ), as follows: Substituting equation (19) into the continuity equation (10) yields the Laplace equation (in dimensional variables) The momentum equation (11) can be recast as Bernoulli's equation in dimensional variables, which can be readily solved once a solution for the velocity potential is found; these equations are also supplied with the appropriate boundary conditions pertaining to the free surface and to the solid-fluid interfaces in the problem. The Laplace and Bernoulli equations are commonly applied in the ocean and offshore engineering field in the form of linear potential flow theory (Faltinsen 1990;Birk 2019). Assuming small-wave amplitudes ( A 0 λ ) and small body motion (relative to the body dimensions), together with a set of (linearised) boundary conditions for the seafloor, body hull, and free surface, solutions for the velocity potential can be determined numerically, e.g., through BEM solvers, such as WAMIT. We obtain the first-order forces acting on the plate by integrating the pressure computed using the linearised version of Bernoulli's equation The force acting on the plate can then be computed using (12) without the viscous term and with pressure specified by (21).

Reduction to the 2D problem
Mocean's device is 3D, and hence to fully understand the resonance it exhibits, we would need to perform 3D simulations.
In this paper, we seek to identify if the resonance observed by Mocean's device can be understood by considering the 2D resonance observed in the flow past an inclined plate. To that end, we wish to compute the non-dimensional surge force on The modelling error relates to the value of k L corresponding to the peak in |X (1) (t)|. A best-fit power curve is fit through k L versus 1/ , which is used to predict a value of k L = 0.799 for an infinitely wide plate, and this value is used to compute the modelling error for each appearing in the final column (see  for more details) the plate per unit width X (1) (t), which we compute in different ways for WAMIT and CFD, since the models are 3D and 2D, respectively. In WAMIT, the non-dimensional surge force per unit width X (1) (t) is approximated by computing the surge force on a very wide 3D plate W = 50L where F x (t ) is the surge force on the plate of width W , = W /L , and we have chosen U = 2 A 0 g where A 0 is the incident wave's amplitude. The value of = 50 was fixed throughout the BEM studies, and its choice was based on a balance of geometrical and numerical considerations. A summary of the -dependence is presented in Table 1 alongside verification data from analytical solutions for a 2D vertical plate (Evans 1970). To calculate the modelling error associated with the choice of finite , we compute the value of k L associated with the peak force using three values of and use this information to estimate the value of k L for an infinitely wide plate using the methodology of . The modelling error for each value of is summarised in Table 1.
All subsequent WAMIT results presented in this paper are for = 50, which has a modelling error of 0.16% associated with the value of k L (see the row associated with = 50 in Table 1). The relative error between our computed value of k L = 0.8003 (for = 50) and available analytical results for k L = 0.7991 in 2D is 0.15%, which is within the modelling error (i.e., < 0.16%) of the method, and hence is a good approximation of a 2D plate.
In CFD, the non-dimensional surge force per unit width is given by where f x (t ) is the force per unit width computed using (12), where dS is the infinitesimal surface element of the 2D plate.

Computational considerations for computational fluid dynamics
To quantify the influence of the problem discretisation on the solution accuracy, convergence studies must be performed for both the temporal and spatial discretisation (Roache 1997). Following the convergence analysis described by , a figure of merit, the relative discretisation uncertainty, U, can be evaluated based on the solution from three different discretisation levels, i.e., fine, medium, and coarse. To illustrate the methodology, we compute the grid uncertainty and convergence uncertainty in Fig. 7. In Fig. 7, the relative discretisation uncertainty is determined independently for the spatial (grid size) and the temporal (time step size) discretisation. For the grid size uncertainty, the grid size in the vertical direction, z , is parametrised by the wave height, 2 resulting in three different discretisation levels, i.e., 10 cells per wave height (CPH), 20 CPH, and 40 CPH. The blue crosses in Fig. 7a show the computed value of the wave height normalised by the reference simulation value (φ g = H /H ref ) at each of these values of z (again normalised by the reference value of z ref = 20 CPH) for the case k L = 1.53, i.e., the smallest wave period considered in this study. A quadratic curve fit is constructed through the values of φ g , with standard error of the fit σ fit . This curve fit is used to calculate φ 0 , which is the value of φ g for an infinitely fine grid. The difference between the value of φ g = 1 on the reference grid and φ 0 is then used to compute the uncertainty through U g = 100 × (1.25 × |1 − φ 0 | + σ fit ) ≈ 3.16%. Note that the number of CPH refers to the discretisation size in the interface region and in the vicinity of the plate. Towards the top and bottom domain boundaries, larger cell sizes are used.
For the temporal discretisation, variable time stepping is used, where the time step is adjusted throughout the course of the simulation, based on the maximum allowable Courant number, Co max . Here, three values of Co max are considered, i.e., 0.500, 0.250, and 0.125. During the spatial convergence study, Co max is set to 0.250, while for the temporal convergence study z is held constant at 20 CPH, and a reference Courant number Co max,ref = 0.250 is used. following the same procedure as used for U g , we can compute U c ≈ 0.44%.
For both the spatial and temporal convergence studies, the measured wave height at the intended plate location during wave-only simulation is considered as the input. Results of the convergence studies for k L = 1.53 are listed in and Co max = 0.250. We use a best-fit curve to compute φ 0 , which corresponds to an infinitely fine grid resolution. b shows the same, but for three different values of Co max at fixed z = 20 CPH  Table 2. Based on the results of the convergence study, a spatial discretisation of 20 CPH (in the interface region and in the vicinity of the plate) with a temporal discretisation of Co max,ref = 0.250 is used for all simulations throughout this study. With that, a combined, relative, discretisation uncertainty of 3.60% is found, which will be used, amongst others, to define the error bars in Sect. 5. With the required discretisation sizes, our CFD model requires 24,000 s per wave period on ten cores of an Intel Xeon Gold 6132 with 2.60GHz, while our WAMIT model computes one-wave period in 0.1 s on six cores of an Intel i7-8750 H with 2.20 GHz.

Case study
Throughout this study, a half-submerged, surface-piercing, plate L = 0.1m is considered. The water depth is set to d = 4 L = 0.4 m. For the setup of the CFD-based NWT, it has to be ensured that the thickness of the plate does not significantly influence the resulting forces. 3 The plate thickness has been determined during preliminary studies and is set to L /20 = 0.005 m throughout this study: this corresponds to a thickness-to-length ratio of 5 × 10 −2 1. Throughout the study, we keep the incident wave height H the same; the associated Keulegan-Carpenter number is 1. This suggests that viscous effects are not appreciable, provided that the incident wave height H is a fair estimate of the wave height near to the plate. As we will see later, this holds for k L away from resonance.
While WAMIT operates in the frequency domain and allows sweeping over a wide parameter space in terms of wave periods and, thus, k L , CFD needs a finite set of defined waves for efficient and timely simulation. To enable   Table 3. For reference, the locations of the regular waves in the Le Méhauté (1969) diagram are shown in Fig. 8.

Excitation forces
We computed the magnitude of the surge force per unit width acting on an inclined plate at different inclinations angles 5 • ≤ α ≤ 90 • and for a range of 0 < k L < 2 using WAMIT with = 50. For α = 90 • , we see excellent agreement with the 2D analytical solution by Evans (1970), as shown in Fig.  9a, for the case of infinite depth. As α decreases, the peak in the surge force per unit width decreases in magnitude, it becomes more narrow banded, and its position shifts to lower values of k L . In Fig. 9b, we consider the effect of the finite depth on the force experienced by a vertical surface-piercing plate. Experimental tests are always conducted in a flume or tank of finite depth; hence, in this paper, we will focus on a specific depth d = 4 L . For d ≥ 4 L , there is very little difference between the force experienced by a vertical plate in finite depth and the force experienced by a vertical plate in infinite depth. In the remainder of the paper, our simulated results are all for the specific depth d = 4L .
We compute the surge force using CFD for the input waves, geometry, and scheme described in Sect. 3. We show the time trace of the non-dimensional excitation force per unit width on the plate in Fig. 10 for each of the values of α considered in Fig. 9 for k L = 0.30 (top) and k L = 1.53 (bottom), representing the upper and lower limits on k L . The non-dimensional excitation force traces are plotted over the non-dimensional time t * = t /T . A different non-dimensionalisation is chosen in this case to aid the comparison between different α and λ . As can be seen in Fig. 10, the periodic motion is established by the time that t * = 10, and so, in all of our phase-averaged results, we perform our averaging on the interval 10 < t * < 20.
In Fig. 11, we show the phase-averaged surge force per unit width of the plate for different angles of the plate as computed using WAMIT (solid curve) and CFD (data crosses), all for the case of finite depth d = 4L . The error bars on the CFD data in Fig. 11 indicate the numerical discretisation uncertainty (U = 3.6%) together with the relative standard deviation, σ r , of the peak-to-trough surge force amplitude between consecutive periods. The relative standard deviation stems from a phase-averaging procedure, using data in the temporal interrogation window between t * = 10 − 20. Based on a zero-down crossing analysis, the force time trace is cut into individual periods and the peak-to-trough surge (a) (a) Fig. 9 The magnitude of the non-dimensional excitation force (surge) per unit width of the plate |X (1) (t)| as computed using WAMIT and (22) with = 50 (dashed curves). a The case of infinite depth with varying inclination angles. The horizontal position of the peak increases as α increases, until it coincides with the analytical solution for a 2D verti-cal surface-piercing plate (Evans 1970), indicated by the black curve. force amplitude is determined for each period, allowing the computation of the mean force amplitude and its standard deviation, σ , which can be normalised by the mean value to get σ r . The surge forces as computed with CFD agree very well with those computed using BEM for α ≥ 45 • . For α ≤ 30 • , CFD computes a smaller surge force than that calculated using BEM. However, even for α = 30 • , there is still good agreement between OpenFOAM and WAMIT for the value of the peak k L . Interestingly, the predicted location of the peak is different in CFD and BEM, with CFD predicting the occurrence of the peak at a slightly lower value of k L than BEM. As shown, the difference in the prediction of CFD and BEM cannot be accounted for by the numerical error. The peak predicted by CFD is slightly offset relative to the peak predicted by linear theory. The offset is towards lower k L , which is consistent with our previous experimental measurements for a class of WECs known as oscillating wave surge converters (Crooks et al. 2015). There, we found that the peak in measured values of the hydrodynamic force peaked at lower values of k L (longer wave periods) than that predicted by numerical models that relied on linear potential flow theory. Potential reasons for the offset are non-linear effects close to the plate and domain effects such as 2D/3D effects and the fact that we have a finite horizontal extent in our CFD model.
Although K C 1 would suggest that viscous effects are not appreciable, in resonant conditions, this may not be true. Even though K C corresponding to the incident wave train is 0.12, if we consider the wave height in the definition of K C to be the wave height of the resonant response, K C can actually be as high as 0.5 for our wave conditions, which is no longer much smaller than unity. In the related harbour literature, the amplification factor as computed from linear potential flow theory significantly overestimates the resonant response of harbours near the Helmholtz mode, as measured in experiments (Cossalter et al. 1982), despite the incident wave conditions not departing significantly from linear wave theory. In this case, previous authors suggest that fluid viscosity may play a role in damping the resonant oscillations. Non-linear effects due to viscous losses are known to have a "softening" effect on fluid resonances such as these (Miles 1981), which would cause the hydrodynamic response to shift towards longer wave periods, which is consistent with what we observe here. Indeed, recent semianalytical studies on oscillating surface-piercing plates show that viscosity tends to shift their resonant peaks towards longer wave periods (Cummins and Dias 2017). A similar discrepancy arises between the linear prediction and physical experiments involving the Helmholtz response of moonpools (Miloh 1983). There have been numerous attempts to modify the potential flow theory in the region close to the moonpool to account for this discrepancy (Chen et al. 2011;Chen and Duan 2012;Tan et al. 2020). The fact that the resonant wave amplification becomes much greater when the moonpool's entrance changes from a sharp corner to a round corner (Tan et al. 2020) suggests that viscous effects are important there. In the case of a WEC, one straightforward way to accomplish better match between potential flow and CFD is to use the semi-analytical modelling used successfully for OWSCs, in which viscous dissipation is introduced near the edges of the plate (Cummins and Dias 2017).
Another form of non-linearity that is neglected in the linearised model is the changing wetted surface area of the plate. In the linear model, the wetted surface area of the plate is fixed in all wave conditions, while in CFD, it changes within a single-wave cycle and between wave conditions (see Fig.  12a,b), which results in a non-linear Froude-Krylov force. Non-linear effects of this kind can shift the resonant peak (Terra et al. 2004), and furthermore, this type of non-linearity would be amplified at smaller angles. In Fig. 11, the shift to longer wave periods is more prominent at smaller α, which is consistent with this effect.
We show the relative wave height at a non-dimensional position x = − 1 2 cos α in Fig. 13 as computed using Open-FOAM and WAMIT. For α ≥ 30 • , the wave height has a single peak in the interval 0 < k L < 2, as shown in Fig.  13a. For α > 30 • , the wave height above the plate is typically 1.5-3 times larger than the incident wave height. However, at smaller angles, the wave amplitudes become "spiky". Indeed, WAMIT predicts wave heights in excess of five times that of the incident wave height for α = 15 • . This agrees qualitatively with the predictions of OpenFOAM; however, the Fig. 11 The non-dimensional magnitude of the surge force per unit width of the plate |X (1) (t)| as computed with WAMIT using (22) with = 50, d = 4 L solid red curves: a α = 90 • , b α = 45 • , c α = 30 • , and d α = 15 • . The black crosses indicate the same computed using (23) with our 2D OpenFOAM model with d = 4L for each of these angles. The error bars indicate the numerical discretisation uncertainty (U = 3.6%) together with the relative standard deviation of the peak-to-trough surge force amplitude between consecutive periods, determined via phase averaging amplitude is much reduced here owing presumably to the same effects as for the hydrodynamic forces, described in the preceding paragraph. Indeed, the wave amplitude exhibits the same offset in the peak value towards longer wave periods. The trend in the wave amplitude is similar to that found for the wave forces. The only exception is at k L near 1.5 for small angles, as shown in Fig. 13b near k L = 1.5. This will be explored in the next section, where we describe the reflection coefficient of inclined plates. We remark here that the shift of the resonant peak towards lower k L and the increase in the resonant response is consistent with the existing numerical and analytical models of sloping harbours (Mattioli 1978;Wang et al. 2011).

Reflection coefficient
To determine the reflection coefficient, |R|, from the numerical CFD data, various methods are available (McKee et al. 2018). In this study, the method by Tuck (1971) is adopted, which allows the determination by measuring the freesurface elevation at two distinct locations in the domain, omitting the need to disentangle the incident and reflected waves via, e.g., frequency-domain analysis. First, we compute the transmission coefficient, |T |, which is defined as the ratio between the wave height of the transmitted wave, 2 A Trans , and the incident wave, 2 A 0 (at great distances from the plate) While the incident wave height 2 A 0 is determined at the intended plate location during preliminary wave-only tests, the transmitted wave is measured at a position x = max [λ ]/10 in the down-wave direction. With the transmission coefficient computed, |R| is subsequently computed as follows: To compute the error bars, shown in Fig. 14, the numerical grid uncertainty U = 3.60%, as well as the relative standard deviation of the incident and transmitted wave height for the consecutive periods in the interrogation window, stemming from the phase-averaging procedure, are considered. For a vertical plate, Fig. 14a, we find that the reflection coefficient monotonically increases with increasing k L above k L = 0.5 representing the fact that most of the inci- The relative wave height inside the wave channel |H wc (t)|/2 A 0 at the middle (y = 0) of the plate at a distance of x = 1 2 L cos α upwave of the plate as computed with WAMIT with = 50, d = 4 L red curves: a α = 30 • , b α = 15 • . The black crosses indicate the same computed using our 2D OpenFOAM model with d = 4L for each of these angles. The error bars indicate the numerical discretisation uncertainty (U = 3.6%) together with the relative standard deviation of the wave heights between consecutive periods, determined via phase averaging (a) (b) Fig. 14 The reflection coefficient |R| for a 2D inclined plate: a α = 90 • , b α = 45 • , c α = 30 • , d α = 15 • . The red curves correspond to the numerical solution for infinite depth, redrawn from Parsons and Martin (1994). The black crosses indicate the current 2D OpenFOAM model with d = 4L for each of these angles, with error bars. The error bars indicate the numerical discretisation uncertainty (U = 3.6%) together with the relative standard deviation of the incident and transmitted wave height between consecutive periods, determined via phase averaging (a) (b) (c) (d) dent wave energy at large k L is reflected back by a vertical plate (Kanoria et al. 1999). The lowest values of k L are not consistent with this trend in Fig. 14a. One might attribute this to the fact that CFD was performed with finite depth, while the red curves in Fig. 14 are computed for infinite depth. This explanation seems unlikely, since the reflection coefficient for a vertical plate in finite depth also increases monotonically with increasing k L (Mei and Black 1969;Kanoria et al. 1999). Another possible source that may account for this discrepancy is the finite thickness of the plate in CFD. It is known that the thickness of the plate strongly affects the reflection coefficient, especially at low k L (Mei and Black 1969); however, the plate is so thin with respect to its length in our CFD model that this is also unlikely to be the source of the discrepancy.
For an inclined plate, we see excellent agreement with the reflection coefficient computed using 2D linear potential flow theory of previous authors (Parsons and Martin 1994): Fig. 14b-d. Indeed, our CFD model captures the previously observed phenomenon of perfect reflection at specific values of k L for α < 90 • . To the authors' knowledge, this is the first time this phenomenon has been detected using CFD. We note that for very small α, we also detect a second peak in the reflection coefficient on the interval 0 < k L < 2, which was also previously reported for inviscid potential flow modelling. Our CFD model also captures the higher harmonics of the resonance, as shown in Fig. 14d. For small angles, the peak near k L = 1.5 corresponds to perfect reflection and a local maximum in wave height (as shown in Fig. 13b), but corresponds to almost zero in the force's amplitude, as shown in Fig. 11d.

Conclusion
In this paper, we have considered the problem of regular waves incident on an inclined surface-piercing plate using CFD and BEM numerical methods. The excitation forces computed using each method show good agreement, provided that the angle of inclination is not too small. When the angle of inclination is ≥ 30 • , the agreement is fair, and the position of the peak k L is in good agreement in both methods. For α ≤ 15 • , both methods show the same trend, but the