Comparison of analytical and numerical solutions for wave interaction with a vertical porous barrier

Analytical solutions for wave interaction with a vertical porous barrier are presented. The analytical solutions are derived using two different methods for taking the depth-average of the pressure drop across the porous barrier. Both solutions assume that the evanescent modes in the wave field can be neglected. The results from the analytical models are compared to results from an iterative boundary element method (BEM) model. The BEM model shows that neglecting evanescent modes is a reasonable assumption for long waves, but that for short waves the velocity through the porous wall from the evanescent modes can be up to 25% of the velocity from the progressive modes at the free surface. However, the effect of neglecting the evanescent modes has only a small effect on the depth-averaged velocity through porous wall and the analytical models derived using depth-averaged assumptions are shown to give good agreement with the BEM model for the reflection coefficient, horizontal force and overturning moment on the porous barrier. The analytical models are used to investigate the effects of the drag and inertia coefficients of the porous barrier on the behaviour of the solution. It is shown that for fixed values of the drag coefficient, wave frequency and amplitude, the solutions for the reflection coefficient lie on approximately semicircular arcs on the complex plane, with the position on the arc determined by the inertial coefficient. This places bounds on the size of the phase change in the reflected and transmitted wave that are possible. The analytical models are also used to derive the asymptotic behaviour of the solution in long and short waves. The implications of the results for more general cases of wave interaction with porous structures are discussed.


Introduction
Porous or perforated structures are of interest in a wide range of applications in coastal and ocean engineering. Wave interaction with porous structures has been studied in contexts such as fixed and floating breakwaters (e.g. Huang et al., 2011;Dai et al., 2018), motion damping of floating structures (e.g. Williams et al., 2000;Lee and Ker, 2002;Molin, 2011;Vijay and Sahoo, 2018), wave absorbers in narrow flumes (Twu and Lin, 1991;Molin and Fourest, 1992), cages used for aquaculture (e.g. Zhao et al., 2010;Dokken et al., 2017), foundations for offshore wind turbines (Park et al., 2014) or tuned liquid dampers (Faltinsen et al., 2011;Crowley and Porter, 2012;Molin andRemy, 2013, 2015) One of the simplest cases for wave-porous structure interaction is the case of a vertical porous barrier occupying the full water column on a flat seabed. This case arises as an idealisation of a wave absorber in a flume, a slotted barrier in a harbour or closely spaced pile breakwater. However, the case is also of interest as a simple situation to investigate behaviour of the solutions. The purpose of this work is to compare analytical and numerical solutions and investigate the effect of the simplifying assumptions.
To model wave interaction with porous or perforated structures, the effect of the openings on the flow can either be explicitly resolved or parametrised in some way. Explicit modelling of the openings has been undertaken by some authors. For example Crowley and Porter (2012) consider the diffraction of waves through an infinitely thin slotted barrier in a potential flow context. Chen et al. (2019) used a Reynolds-Averaged Navier-Stokes (RANS) model to explicitly model the flow through a perforated wall. Mentzoni and Kristiansen (2019) also used a Navier-Stokes solver to explicitly model the flow through oscillating perforated plates. A smooth particle hydrodynamics (SPH) approach was investigated by Meringolo et al. (2015) and Ren et al. (2018).
Explicit modelling of the flow through the openings is computationally expensive and for many practical applications it is sufficient to treat the porous structure as a homogeneous area with a pressure drop across it, where the pressure drop is a function of the physical properties of the structure and the flow velocity. For large volume porous structures such as rubble mound breakwaters, the energy dissipated within the structure can be approximated by Darcy's law, where the pressure drop across the structure is assumed to be linearly proportional to the flow through the porous structure (Chwang and Chan, 1998). For thin porous structures, such as slotted screens or perforated sheets, flow separation at the edges of the openings causes turbulence, resulting in a force proportional to the flow velocity squared (Mei et al., 1974;Molin, 2011). The pressure drop across a thin porous barrier in oscillatory flow can be written in terms of a quadratic drag term and an inertial term due to acceleration of the flow through the openings. The solutions for wave interaction with a vertical porous wall that have been developed, can be classified firstly by whether the quadratic drag is linearised or treated explicitly and secondly by whether inertial effects are considered or not. Yu (1995) presented a solution in which the energy loss is assumed to be linearly dependent on the flow velocity. Yu presented a simple closed form solution for the velocity potential, forces and reflection coefficient in terms of the so-called 'porous effect parameter' defined by Chwang (1983). The porous effect parameter involves a linearised dissipation coefficient, which must be determined empirically. As the dissipation is a function of the flow velocity through the porous barrier, the linearised coefficient is dependent on both wave frequency and amplitude as well as the physical characteristics of the structure. Chakrabarti and Sahoo (1996) considered a nearly-vertical porous wall and used a Green function technique to solve the boundary value problem, using linearised dissipation. Sivanesan (2016, 2017) considered a porous barrier occupying a portion of the water column, with linearised dissipation across the porous barrier. Sivanesan and Manam (2019) extended this work to consider the case of a vertical porous barrier with two gaps. Kaligatla et al. (2017Kaligatla et al. ( , 2018 used the mild slope equations to investigate the effects of bottom topography on linear wave dissipation by vertical porous barriers. Geng et al. (2018) considered the use of multiple porous sheets used as a wave absorber in a wave tank, assuming linear dissipation across each sheet.
Various authors have developed solutions for wave interaction with a vertical porous barrier which treat the quadratic loss explicitly. Mei et al. (1974) used shallow water theory, where the flow velocity is assumed to be constant with depth, to solve the quadratic loss problem, including inertial effects. Hagiwara (1984) used an integral equation method to develop a solution for intermediate depth waves with quadratic loss. Bennett et al. (1992) also solved the quadratic loss problem for a thin vertical barrier, using an eigenfunction expansion method. Fugazza and Natale (1992) considered the case of one or more vertical porous walls in front of a solid back wall, representing a Jarlantype breakwater, and solved for the quadratic energy loss. However, their solution neglected the evanescent components of the wave field. Molin and Fourest (1992) considered the problem of modelling a series of vertical porous plates used as wave absorbers in a flume. They modelled the quadratic loss across the plates using an iterative eigenfunction expansion method. Kakuno and Liu (1993) used the method of matched asymptotic expansion, considering both the quadratic loss and inertia coefficient. Isaacson et al. (1998) extended the work of Bennett et al. (1992) to the case of a vertical slotted barrier occupying a portion of the water column, using an eigenfunction expansion method. However, Isaacson et al. linearised the energy loss term, using the same porous effect parameter adopted by Yu (1995). Liu and Li (2017) used an iterative boundary element method to provide a general solution for wave interaction with a thin porous structure in two dimensions.
The works mentioned above all require numerical solution. Various analytical solutions for wave interaction with a vertical porous barrier with quadratic loss have also been proposed. Hayashi et al. (1966) presented a solution for shallow water waves, but did not consider inertial effects. Kriebel (1992) extended the solution of Hayashi et al. to intermediate depth waves and also neglected inertial effects. Kim (1998) presented a solution for intermediate depth waves which accounts for inertial effects (the same solution is also presented in Suh et al. 2011a andKoraim 2011). Huang (2007) gave a solution for two adjacent slotted walls, which reduces to Kim's formula for a single wall when the distance between the two walls is zero and inertial effects are neglected.
The analytical solutions mentioned above all make the simplifying assumption that evanescent components of the wave field are zero. This is a valid assumption in the long wave limit, where the velocity profile is constant with depth. However, for waves in intermediate and deep water, the evanescent component of the wave field is non-zero (this point is discussed further in Section 4).
In this work we compare analytical and numerical solutions to assess the impact of neglecting the evanescent terms in the velocity potential. We start by deriving analytical solutions in two slightly different ways. The first approach extends that of Kriebel (1992) to account for inertial effects. In the second we derive a model that gives the same results as Kim's solution (Kim, 1998), but is solved in terms of a quadratic equation rather than the quartic polynomial used by Kim. The analytical solutions are used to explain the influence of the wave conditions and the porosity coefficients (drag and inertial terms) on the behaviour of the solution. The analytical solutions are also used to provide bounds on the range of possible solutions in terms of the reflection and transmission coefficients and to examine the limiting behaviour of the solution in long and short waves.
Secondly, we assess the impact of neglecting the evanescent modes by comparing the analytical results to results from an iterative boundary element method (BEM) model which does not require simplifying assumptions about the form of the velocity potential. A BEM model was chosen over an eigenfunction expansion method as it was found that the solution was quickly convergent for all parameter values investigated in this study. The BEM approach adopted in this study is similar to that proposed by Liu and Li (2017). An extensive validation of the iterative BEM model against experimental data was presented in Liu and Li (2017), so the focus here is on the differences between the analytical and numerical solutions.
The paper is organised as follows. A brief review of models for the pressure drop over a thin porous surface is presented in Section 2. The mathematical formulation of the problem is presented in Section 3. The analytical solutions are derived in Section 4, together with a comparison between the two solutions and an investigation of the behaviour of the solution. The BEM model is presented in Section 5 and the results of the analytical and BEM models are compared in 6. Finally, conclusions are presented in Section 7.

Pressure drop over a thin porous surface
The pressure drop, , across a thin porous surface can be modelled as (Sollitt and Cross, 1972) where is the fluid density, is the kinematic viscosity, is the component of the velocity normal to the porous wall (assumed to be the average velocity close to the wall, rather than the flow speed through the gaps), is a length scale (related to wall thickness, hole size, etc.), is a dimensionless friction or drag coefficient and is an inertial coefficient with dimension of length. The first term on the RHS of (1) is a viscous friction loss, the second term is a turbulent dissipation loss, the third term is due to phase difference in the flow across the boundary and does not dissipate energy. Sollitt and Cross (1972) noted that the linear drag term is dominant at low Reynolds number flow and the quadratic term is dominant at high Reynolds number. The Reynolds numbers for wave interaction with thin porous structures are usually high enough that linear viscous forces can be neglected.
Under the assumption that the openings in the porous structure are small relative to the wavelength, the wakes will be quickly regularised and homogenised, so that the flow away from the porous wall can be modelled using potential flow theory (Molin, 2011). The potential flow in the outer region is assumed to be time-harmonic with angular frequency and represented by a complex potential: where is the spatial component of the potential. The flow velocity normal to the porous wall is given by The quadratic pressure loss term introduces higher harmonics into the flow. However, Mei et al. (1974) showed that the higher harmonics are effectively negligible compared to the fundamental harmonic, and the time-dependence can be linearised using Lorenz's principle of equivalent work: The pressure is given by the linearised Bernoulli equation Substituting (2)-(5) into (1) and neglecting the viscous loss term gives the following boundary condition on the porous wall where 1 and 2 are the potentials on either side of the porous boundary and is a coefficient defined as Since is dependent the fluid velocity through the porous surface, ∕ , it is dependent on wave amplitude and frequency and varies over the porous surface of the structure due to the spatial variation in the fluid velocity.
The drag and inertial terms are normally considered separately, with the drag term quantified assuming steady viscous flow and the inertial term quantified assuming inviscid oscillatory flow. Models for the quadratic friction coefficient have been reviewed in Molin (2011) and Huang et al. (2011). Two commonly used models are those of Mei et al. (1974) and Molin (2011). Mei et al. defined the friction coefficient as where is the porosity of the surface, defined as the open-area ratio (i.e. the ratio of the area of the openings to the total area of the surface) and is the coefficient of contraction, approximated by = 0.6 + 0.4 2 . Molin and Fourest (1992) defined the friction coefficient as where is an empirically determined discharge coefficient, usually taking values in the range 0.4-0.5. The models (8) and (9)  , where is the drag coefficient for steady flow. For larger KC numbers, corresponding to small perforation sizes or large flow velocities, the influence of the KC number decreases and the drag coefficient tends to the value for steady flow. In the present study, the value of , is assumed to be independent of KC number. Mei et al. (1974) noted that for slotted screens, the inertial coefficient can be calculated using a long-wave approximation used in acoustics (Morse and Ingard, 1968): where is the distance between the centre of the slots. McIver (1998) developed a model for the blockage coefficient of a circular hole in a rectangular duct, also using a long-wavelength approximation. McIver compared the results to the expression of Tuck (1975), developed under the assumption that the hole size is small compared to the duct width, and developed a quadratic correction to Tuck's formula which gives 'graphically indistinguishable' results to the full solution. Reformulating equations (32)-(35) in McIver (1998) in terms of porosity and evaluating all the coefficients in the expressions gives: where is the length of the side of the square duct (equivalent to the separation between adjacent hole centres). Molin and Remy (2015) developed a similar model for the blockage coefficient of a circular hole in a circular duct, which gives results that are almost identical to those from (11). Eqs. (10) and (11) assume that the porous barrier is infinitely thin. Kakuno and Liu (1993) used a blockage coefficient derived in Flagg and Newman (1971) and Taylor (1973) to model the inertial effect of waves passing through an array of cylinders of finite thickness , given by Expressions (10)-(12) are compared in Fig. 2 for 0 < < 0.5. The inertial coefficients from (10) and (11) are similar for > 0.2, but the difference increases for lower porosities. Expressions (10) and (12) agree well for = 0, but (12) predicts a significant influence of the sheet thickness for > 0.
Expressions (10) and (11) were developed under a long-wavelength approximation. Crowley and Porter (2012) solved the full diffraction problem for an infinitely thin slotted barrier and showed that the reflection and transmission coefficients are very similar to those calculated from (10) when < 0.1, where is the wavenumber and 2 is the gap width. In our notation 2 = . Molin and Remy (2015) found that (10) and (11) gave good agreement with measurements of added mass and damping coefficients for experiments of sloshing in a tank with a slotted or perforated screen.

Problem formulation
The water is assumed to be of constant depth ℎ. A flat, rigid, vertical porous wall that extends throughout the water column is located at = 0. The thickness of the wall is assumed to be much less than the wavelength. Away from the porous wall, the fluid is assumed to be inviscid and incompressible, and its motion is irrotational, so that a velocity potential can be used to describe the fluid motion. Regular linear waves of amplitude and angular frequency , propagating from the negative direction are incident on the porous wall. The fluid domain is divided into two regions, with the potential in the region < 0 denoted 1 and the potential in the region > 0 denoted 2 . A sketch of the problem definition is shown in Fig. 3.
The fluid motion is assumed to be time harmonic so that the potentials in each domain can be written as where is the spatial component of the potential in domain . The spatial potentials in each domain satisfy the Laplace equation, the linearised free-surface condition and the no-flow condition on the seabed: where = 2 ∕ is the infinite-depth wavenumber and is the acceleration due to gravity. The porous wall is represented as a homogeneous region with infinitesimally small holes, so that the flow through the porous wall in the normal direction is assumed to be continuous on either side of the boundary, so that The pressure drop across the porous wall, results in the additional boundary condition given in (6).
Once the problem has been solved, the horizontal force and moment about the seabed per unit width, , on the porous wall are given by

Analytic solutions
The potentials in each domain can be expressed using the standard eigenfunction expansions as and are the unknown complex reflection and transmission coefficients, corresponding to the progressive components of the wave field. and are the unknown complex amplitudes of the evanescent components of the wave field. is the positive real solution of = tanh( ℎ) and are the positive real solutions of − = tan( ℎ), ordered in increasing value. The vertical eigenfunctions are given by Substituting (20) and (21) into (17), multiplying by 0 , integrating over depth and using the orthogonality of the vertical eigenfunctions gives 1 − = .
( 23) Similarly, multiplying by and integrating over depth, gives It is instructive to rewrite the dynamic boundary condition (6) as Suppose that evanescent terms are zero. Under this assumption, the terms on the LHS of (25) have a vertical dependence on cosh( ( + ℎ)) and terms on the RHS are proportional to cosh 2 ( ( +ℎ)). This condition is met in the long wave limit when = 0 or on the seabed when + ℎ = 0. However, in general, for positive frequencies and positive elevations above the seabed we have cosh( ( + ℎ)) ≠ cosh 2 ( ( + ℎ)). Therefore, there must be some evanescent terms which are non-zero.
So far, no approximations other than those of linear wave theory have been made. To derive an analytical solution, we assume that the evanescent component of the wave field is negligible and solve the system using a depth averaged approximation. Solutions can be derived in two ways, depending on the stage at which the depth-averaging is applied. The two methods of solution are outlined in Sections 4.1 and 4.2. It will be shown in Section 6, through comparison with results from the BEM model, that the depth-average velocity from the evanescent components is close to zero, and neglecting the evanescent components therefore results in a reasonable approximation. For both solutions we will make the assumption that and are constant with depth, so that the physical characteristics do not change with depth and we assume that there is no dependence of and on the amplitude of the oscillatory velocity (i.e. KC number).
Under the assumptions that the evanescent terms are zero, substituting (20) and (21) into (18) and (19) gives the non-dimensional force and moment on the porous wall as Note that when = 1 these expressions give the linear force and moment on a solid wall. Therefore, is the ratio of both the force or moment on the porous wall to the force or moment on the solid wall.

Method 1: Depth-average of dynamic boundary condition
Substituting (20) and (21) into (6), neglecting the evanescent components and integrating both sides with respect to gives a quadratic equation for the transmission coefficient and A similar approach was taken by Kriebel (1992) to solve the case where = 0. In this case is real and (28) can be solved explicitly to give In the case that > 0, an explicit solution of (28) can be obtained, as explained in Section 4.3. Alternatively, (28) can be solved numerically using standard algorithms. In this work the MATLAB function fsolve with a first guess = 0 has been used.

Method 2: Depth-average linearised drag coefficient
The quadratic drag term in the dynamic boundary condition (6) can be linearised using Lorenz's principle of equivalent work to give a linear boundary condition Substituting (20) and (21) into (34), neglecting the evanescent components and evaluating the integrals gives Substituting this back into the linearised dynamic boundary condition (33) together with the expressions for the potentials gives a quadratic equation for the transmission coefficient and . Kim (1998) derived an analytical solution in a similar way, by linearising the quadratic drag coefficient and neglecting the evanescent components. In Kim's method the equations were solved in a different way, requiring the solution of quartic polynomial. The derivation outlined above results in a quadratic rather than quartic equation to solve. As noted in the preceding section, the quadratic (37) can be either be solved numerically or explicitly as the root of a quartic polynomial, as explained below.

Comparison of methods and explicit solution
The two methods outlined above result in the same form of quadratic equation to be solved (cf. (28) and (37)), with the coefficients 1 and 2 differing only in the form of the depth dependence terms ( ℎ) and ( ℎ). The functions ( ℎ) and ( ℎ) are illustrated in Fig. 4. The two functions tend to the same value for low ℎ (i.e. long waves) and tend to constant but different values for high ℎ. The asymptotic values of the two functions for small and large values of ℎ are The two methods will therefore give the same results for long waves but will give a constant offset for short waves. The agreement in long wave conditions is to be expected, since the velocity potential is less variable with depth in long waves, the method used to take the depth average will have less influence on the results. The explicit solution to (28) and (37) for the case = 0 is given in (31) and (32). To obtain an explicit solution in the case > 0 we write the transmission coefficient using Euler's formula as = E. Mackay and L. Johanning | |(cos + sin ), where is the phase of . Substituting this into (28) and separating the real and imaginary parts gives the following where is either 1 or 2 . Substituting (45) into (44) and using the identity sec( ) ≡ ± √ 1 + tan 2 gives a quartic equation for | |: When = 0 this reduces to the expression given in (32). The general solution for positive and can be written down explicitly. However, the explicit solution involves a large number of terms and is not particularly instructive to look at. The explicit solution is therefore not copied here, but can easily be obtained by standard computer packages.

Behaviour of analytical solution
From (28) and (37) it is apparent that the reflection coefficient (and hence also the horizontal force and moment) is a function of (either 1 or 2 depending on the method used) and only. The relationship between , , and is shown in Fig. 5. The magnitude of the transmission coefficient is monotonically decreasing with increasing , corresponding to increasing (decreasing porosity), increasing wave steepness , or decreasing ℎ (longer waves). An increase in the inertia coefficient leads to a reduction in the transmission coefficient and a negative shift in the phase of the transmitted wave. The inertia coefficient has a proportionally larger effect at lower values of , as would be expected from inspection of (28) and (37).
The reflection coefficient is related to the transmission coefficient by = 1 − (this is true for the general solution including evanescent terms, as explained at the start of this section). When = 0 the phase of is zero and is monotonically increasing with . When > 0 the phase change in leads to a non-monotonic variation in | | with .
The variation of the magnitude and phase of the and with and can be better understood, by plotting the real and imaginary parts for fixed values of , as shown in Fig. 6. Each coloured line on the plot corresponds to a solution for a fixed value of and values of between 0 and ∞. When = 0 the solution to (28) and (37) is .
In fact, this is the exact solution, since the evanescent terms are zero in the case = 0. Making the substitution = 2 tan( ∕2), we obtain This is the equation of a circle in the complex plane, centred at = 1∕2 and with radius 1∕2. Since 0 ≤ < ∞, we have 0 ≤ < and Eq. (48) describes a semicircle in the fourth quadrant of the complex plane. The solutions for = 0 correspond to a purely inviscid interaction at the porous boundary (the fluid is already assumed to be inviscid away from the porous wall), where there is no energy loss due to turbulence and  the pressure on the porous wall results from acceleration of the fluid through the gaps. The range of values that and can take is bounded by semicircles in the first and fourth quadrants of the complex plane, respectively, representing the solutions for = 0. For > 0 the arcs representing the solutions for constant are not quite a semicircle, with the solutions lying slightly outside the semicircle defined by the two solutions on the real axis (for = 0 and = ∞). For short waves, the interaction with any structure with a vertical porous wall at the free surface will tend toward the present case, where the porous wall occupies the full water column. The bounds on the range of the reflection and transmission coefficients derived here, would therefore be expected to apply for short waves in the more general case as well. The results imply that for given values of ℎ, and , the maximum possible phase change in is governed by solution for = 0 which is given by (31) and (32). If we denote the solution for = 0 as 0 , then the solutions for > 0 lie approximately on a semicircle centred at (1 + 0 )∕2 with radius (1 − 0 )∕2. In contrast, the transmission coefficient can have a maximum phase lag of ∕2 for any value of , as the solution arcs always intersect the imaginary axis as → ∞. However, large phase lags are associated with small transmission coefficients. Fig. 7 shows the magnitude and phase of and for fixed values of and 0 ≤ < ∞. It is evident that for a given value of | | or | |, the maximum phase angle is bounded by the phase angle of the solution at = 0. From (47)  It is interesting to consider the behaviour of the reflection coefficient in long and short waves. The solution is dependent on the variable , which is a function of wave steepness. If wave steepness is held constant then the amplitude tends to infinity as ℎ → 0. Conversely, if amplitude is held constant then the steepness tends to infinity as ℎ → ∞. So both cases lead to non-physical results in one of the limits. Example solutions are shown in Fig. 8 for the cases of both constant steepness and constant amplitude, assuming = 100 and = 0. Note that since = 0, is real. The solutions shown have been calculated using method 1, but the behaviour is the same for solutions from method 2.
In the case of constant amplitude, we have In this case, the reflection coefficient will tend to a constant value for long waves and tend to one for short waves. In the case of long waves, the horizontal velocity tends to a constant value, √ ∕ℎ, and hence the drag resistance of the porous wall will also tend to a constant value, leading to a constant reflection coefficient. The physical interpretation of the short wave limit is that since the velocity through the holes is proportional to , if is constant then the velocity will tend to infinity in short waves, causing the flow through the porous barrier to experience increasing drag. In the limit, the flow through the barrier will decrease to zero, causing all energy to be reflected. In reality, the wave will break when the steepness exceeds the Miche limit and other effects will become relevant before this limit is reached. In particular, E. Mackay and L. Johanning the assumption that the porous wall can be modelled as a homogeneous surface will become less valid and diffraction effects will need to be accounted for (see e.g. Crowley and Porter, 2012).
In the case of constant steepness we have In this case, the reflection coefficient tends to a constant value for short waves and tends to one for long waves. For the short wave case, maintaining a constant steepness makes the problem scale-invariant and hence a constant reflection coefficient is expected under the assumption that the porous wall can be treated as a homogeneous surface. The physical interpretation for the long wave limit is the same as that for the short wave limit in the constant amplitude case. In the case of constant steepness waves, the horizontal velocity tends to infinity as ℎ → 0, so the drag resistance of the porous wall will also increase to infinity. In Fig. 8, there appears to be a distinct change in the behaviour of the solution around ℎ = 1. The convergence of the solutions toward their asymptotic values is dependent on the behaviour of the depth functions ( ℎ) and ( ℎ). From Fig. 4 it can be seen that for ℎ < 1, both ( ℎ) and ( ℎ) are well-approximated by 1∕ ℎ and for ℎ > 3, both ( ℎ) and ( ℎ) are approximately constant.
In the case that > 0, the asymptotic behaviour of the reflection coefficient remains the same for long waves (i.e. tending to a constant value for constant amplitude or tending to one for constant steepness), but for short waves the reflection coefficient will tend to one for both cases of constant amplitude and constant steepness. This can be seen from (28) and (37), since → ∞ as → ∞ the only solution can be = 0. As before, the assumption that the porous wall can be treated as a homogeneous barrier becomes less valid as the ratio of the wavelength to the size of the openings decreases and diffraction effects become more important.

BEM solution
The multi-domain BEM approach used here is similar to that presented in Liu and Li (2017). To simplify the analysis we define dimensionless potentials as where the superscript indicates a dimensionless quantity. Two domains are defined to solve the BEM problem. The first domain, located on the up-wave side of the porous screen (LHS in Fig. 3), is bounded by the porous wall, the free surface, the sea floor and a control surface at = − . The second domain, located on the down-wave side of the wall (RHS in Fig. 4), is the mirror of the first domain in the porous wall, with the control surface at = . The normal vectors to the boundaries, , are defined to point outward from each domain. Application of Green's third identity to the two domains yields integral equations for 1 and 2 : where is the boundary of domain and is a Green function, satisfying the Laplace equation, defined as and ( , ) are the coordinates of an arbitrary point on . The boundary conditions for the integral equations on the mean surface and sea bed are given by (15) and (16) The dynamic boundary condition on the porous wall is given by (6) and (54) as The boundary conditions on the control surfaces are derived as follows.
If is sufficiently large then the evanescent components tend to zero, so Assuming that is sufficiently large for (60) and (61) to hold, expressions for in terms of 1 and 2 at = ± are obtained by integrating both sides of (60) and (61) over to give E. Mackay and L. Johanning Substituting (62) and (63) into (60) and (61) and taking partial derivatives gives the boundary conditions on the control surfaces as The system is solved by discretising the boundary into discrete segments and assuming the potential is constant on each panel. The boundary conditions are substituted into the discretised form of (55) giving a linear system of equations for the potential (see Appendix), where the strength of the dissipation coefficient ( ) on each panel of the porous wall is an additional unknown (note that ( ) varies with water depth). The system is solved iteratively, with a first guess of ( ) = 10 5 , corresponding to a high porosity. After the system is solved, the normalised velocity through the porous wall, ( ) = 1 (0, )∕ , is recalculated from (58) and the new velocity is defined as +1 ( ) = ( ( ) + −1 ( ))∕2, where the superscripts denote the iteration number and 0 ( ) = 0. The new velocity is then used to define a new value of ( ) on each panel using (59) and the system is resolved. The iterations are terminated when max Once the system has been solved, the reflection coefficient can be calculated from (62) and the force and moment per unit width on the wall can be calculated from (18) and (19).
A convergence study was conducted to determine an adequate discretisation length and position for the control surfaces, for various values of , , ℎ and . It was found that using a constant discretisation length of ℎ∕200 and locating the control surfaces at = 5ℎ was sufficient to obtain a stable solution.

Velocities from evanescent components
Before comparing the BEM and analytical solutions directly, it is instructive to consider the size of the evanescent components of the potential relative to the progressive components. The total dimensionless flow velocity through the porous screen is given by The dimensionless velocity through the porous screen associated with the progressive components of the potential is given by The dimensionless velocity through the porous screen associated with the evanescent components is therefore given by Figs. 9 and 11 show the ratio | ( )∕ (0)|, where and have both been calculated from the BEM model. The denominator in the ratio | ( )∕ (0)| is the maximum value of the velocity from the progressive components of the wave field, which occurs at the free surface. The cases shown are for = 0.1, = 10, 50, 100, 200 and = 0, 0.5. Results are shown for various values of ℎ, with the colour of the line indicating the value of ℎ. As mentioned in Section 4, in the long wave limit, when ℎ = 0, the evanescent component of the wave field tends to zero. Figs. 9 and 11 show that the relative size of the evanescent waves increases with ℎ. The magnitude of the evanescent components also increases with the friction coefficient , with a peak value of | ( )∕ (0)| of around 25% at the free surface for the range of variables shown here. The magnitude of the evanescent components decreases slightly with increasing inertial coefficient .
The phase of the velocity from the evanescent components relative to the phase of the progressive component is shown in Fig. 10 for E. Mackay and L. Johanning  the cases corresponding to those in Fig. 9 for = 0. At the top of the water column, the evanescent components have a phase lead over the progressive component, whilst further down the water column the trend is reversed. The pattern is similar for the cases with = 0.5 and is not shown here.

Reflection, transmission and forces
Figs. 13 and 14 present a comparison of the magnitude of the reflection coefficient and normalised surge force and moment from the BEM and analytical models. The force and moment on the porous wall have been normalised by the force and moment on a solid wall for the same wave conditions, denoted 0 and 0 respectively. For the two E. Mackay and L. Johanning  analytical models, we have = ∕ 0 = ∕ 0 . However, for the BEM model the presence of the evanescent components of the wave field results in differences between these quantities, so it is interesting to compare all three values. In all cases, the agreement is good in longer waves, with ℎ < 1. This is consistent with the observation that the relative size of the evanescent waves decreases with ℎ, implying that neglecting the evanescent components is a reasonable approximation for the analytical models. At higher values of ℎ the magnitude | | from the BEM model agrees better with the analytical model using method 2, with method 1 giving slightly lower values in the case = 0. In the case = 0.5, the agreement is closer between all three models, which is consistent with the observation above that the evanescent waves are smaller in this case.
For the surge force, the analytical solution using method 1 agrees better with the BEM model for = 10, but the force from the BEM model is slightly lower than both analytical solutions for higher and = 0, with method 2 giving a larger discrepancy. For the moment about the seabed, the BEM model gives a value in between the two analytical solutions for 1 < ℎ < 3, with method 1 giving better agreement at higher ℎ. The increasing discrepancy at higher ℎ for larger values of is consistent with the larger relative size of the evanescent velocities, since neglecting these components will result in a larger error in these cases.
The phase of , and is shown in Fig. 15 for the case = 0.5. The case = 0 is not shown as the phase is zero for all models in this case. The phases are generally is good agreement for all three models, but the agreement is marginally closer for method 1.

Effect of inertial coefficient
In Section 4.4 the it was shown that the analytic solutions fall on near-semicircular arcs in the complex plane, with the position on the arc determined by the value of . Runs with the BEM model were conducted to verify this behaviour for the numerical solution. The analysis presented in Section 4.4 considered the behaviour of the solution for constant values of the coefficient . The coefficient is a function of ℎ, and , but with the dependence on ℎ differing between methods 1 and 2. The runs with the BEM model were conducted for ℎ = 1 and ℎ = 4, which represents solutions for long and short waves respectively. For each value of ℎ, runs were conducted for = 0.1, = 0.1 and = 0, 0.5, 1, 2, 5, 10. The results are shown in Fig. 16.
For the long wave case with ℎ = 1, the results from the BEM model and analytic solutions agree very well. There is negligible difference between the two analytic solutions, since the functions ( ℎ) and ( ℎ) have the same asymptotic behaviour for low ℎ (see Fig. 4).
Moreover, as discussed above, at low ℎ the evanescent component of the potential is smaller, meaning that the analytic solutions will be good approximations.
For the short wave case with ℎ = 4, there is a small difference between the two analytic solutions when > 0, which is due to the difference in the functions ( ℎ) and ( ℎ). When = 0 these functions are not part of the solution, hence the two analytic solutions coincide. The analytic solution for = 0 is the exact solution to the linear potential flow problem, so in theory the BEM solution should exactly coincide. In practice, there are some small differences due to numerical effects related to discretisation of the problem. However, the agreement is good overall. For the cases with > 0 the BEM solution lies somewhere between the two analytical solutions, but again the agreement is good overall.

Conclusions
Two analytical solutions for wave interaction with a vertical porous barrier have been presented. The first extends the solution of Kriebel (1992) to include the inertial effects of the porous barrier and the second simplifies the solution of Kim (1998). The two methods were shown to both be the solution of the same quadratic equation, with the difference occurring in the depth dependence for the quadratic coefficient. The analytical solutions provide a relationship between the behaviour of reflection and transmission coefficient and forces on the porous barrier with the wave conditions and porosity parameters.
It was shown that the range of values that the reflection and transmission coefficient can take is bounded by a semicircle in the complex plane of radius 1/2, centred on the real axis at = 1∕2, with the semicircle lying in the first quadrant of the plane for the reflection coefficient and in the fourth quadrant for the transmission coefficient. It was also shown that the maximum phase angle for a given value of | | or | | is bounded by the relations 0 ≤ ≤ arccos(| |) and − arccos(| |) ≤ ≤ 0, where and are the phase angles of and respectively. The analytical solution was also used to derive the asymptotic behaviour of the solution for long and short waves for the cases of constant amplitude and constant steepness. For the case of constant amplitude waves, the reflection and transmission coefficients tend to constant values for long waves. Similarly for constant steepness waves, the reflection and transmission coefficients tend to constant values for short waves. The cases of constant amplitude in the short wave limit and constant steepness in the long wave limit both correspond to non-physical cases of either infinite steepness or infinite amplitude respectively. The analytical solutions are approximate, in that they both neglect the evanescent components in the wave field. Results from an iterative BEM model demonstrate that for longer waves, this assumption is well justified, but for shorter waves the amplitude of the velocity through the porous wall from the evanescent modes can be up to 25% of the velocity from the progressive modes at the free surface. Despite neglecting the evanescent components of the wave field, both analytical solutions agree well with the BEM model in the prediction of the reflection and transmission coefficients and forces and moments on the porous wall. For the reflection coefficient, the solution based on Kim's method has slightly better agreement with the BEM model than the solution based on Kriebel's method, with Kriebel's method leading to a small underprediction of the reflection coefficient at higher frequencies. However, the solution based on Kriebel's method gives better agreement in terms of forces and moments.

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

Acknowledgement
This work was funded under EPSRC, UK grant EP/R007519/1.

Appendix. Matrix equations for BEM model
The boundary of each domain is discretised with panels on the porous wall, panels on the free surface, panels on the control surface and panels on the seabed. Panels in the up-wave domain are numbered sequentially in an anticlockwise direction starting from the seabed at the porous wall. Similarly, panels in the down-wave domain are numbered clockwise, starting at the same point.
Define matrices where and are the centre and length of panel , respectively. After gathering terms, the system of equations for the potentials in each domain can be written as where 1 = The vector w is given by The matrices , , and are sub-matrices of S, corresponding to the terms multiplying the potentials on the porous, surface, control and bottom panels respectively, so that = [ ] . The notation • denotes the Hadamard product (elementwise multiplication) and is a matrix the same size as , with rows corresponding to the porous parameter on each panel on the porous wall (updated on each iteration), with = ( ). The matrix is the same size as , with terms given by