Using a porous-media approach for CFD modelling of wave interaction with thin perforated structures

This work presents the use of a porous-media approach for computational fluid dynamics (CFD) modelling of wave interaction with thin perforated sheets and cylinders. The perforated structures are not resolved explicitly but represented by a volumetric porous zone where a volume-averaged pressure gradient in the form of a drag term is applied to the Navier–Stokes momentum equation. The horizontal force on the structures and the free-surface elevation at wave gauges around the cylinder model have been analysed for a range of porosities and regular wave conditions. The CFD results are verified against results from a linear potential-flow model and validated against experimental results. The applied pressure gradient formulation produces good agreement for all porosity values, wave frequencies and wave steepnesses investigated. It is demonstrated that an isotropic macroscopic porosity representation used for large volumetric granular material can also be used for thin perforated structures. This approach offers greater flexibility in the range of wave conditions that can be modelled compared to approaches based on linear potential-flow theory and requires a smaller computational effort compared to CFD approaches which resolve the flow through the openings. The approach can therefore be an efficient alternative for engineering problems where large-scale effects such as global forces and the overall flow-behaviour are of the main interest.


Introduction
Modelling wave interaction with structures consisting of thin porous or perforated elements is of interest in various contexts. Examples are fixed or floating breakwaters, cages for aquaculture or tuned liquid dampers with slotted baffles.
In general, numerical modelling of fluid-structure interaction with thin perforated barriers is often based either on potential-flow theory or on the Navier-Stokes (NS) equations. The first assumes the fluid to be inviscid, incompressible and irrotational. The latter, which is commonly referred to as computational fluid dynamics (CFD) accounts for the viscosity of the fluid, can be used with both compressible or incompressible fluid properties and generally offers various levels of turbulence modelling.
In the context of potential-flow modelling for porous structures, two main types of linearisation can be employed. One concerns the wave condition, the other concerns the implementation of macroscopic porosity representation.
Regarding the wave condition, there is a large volume of work on both linear and fully non-linear potential-flow models for wave interaction with impermeable structures. Potential-flow models for wave interaction with porous structures have mainly been based on linear wave theory. This assumes that the wave steepness and body motions are small. Some examples of existing work are studies on a simplified representation of a Jarlan-type breakwater, done by Fugazza and Natale (1992), work on a vertical porous barrier, done by Mei et al. (1974), Hagiwara (1984) and Bennett et al. (1992), or modelling of a series of vertical porous plates in a flume by Molin and Fourest (1992). Relatively little work has been done on higher-order or non-linear potential-flow models for porous structures. Interaction of cnoidal waves with an array of vertical concentric porous cylinders has been investigated by Weng et al. (2016) and with a concentric cylindrical structure with an arc-shaped outer cylinder by Zhai et al. (2020). Solitary waves have been used to study interaction with a concentric porous cylinder system, done by Zhong and Wang (2006), or with vertical wall porous breakwaters, done by Lynett et al. (2000).
Regarding the porosity representation, most potentialflow models apply the porous pressure-drop in a linearised way and only few models use a non-linear formulation. Examples for work that use a linear relationship between pressure-drop and flow velocity are studies on a nearly vertical porous wall, done by Chakrabarti and Sahoo (1996), or for investigations on the effects of bottom topography, done by Kaligatla and Sahoo (2017) and Kaligatla et al. (2018). A quadratic pressure-drop formulation has for instance been used by Mei et al. (1974) in combination with shallow water theory, by Bennett et al. (1992) for thin vertical barriers and by Liu and Li (2017) for an iterative boundary element model.
Potential-flow models based on linear wave theory have the advantage of short computation times but their fundamental assumptions can limit them in their application such as in the presence of steep or breaking waves. Non-linear potential-flow models can overcome the limitations in terms of wave condition but require longer computation times and are still not able to capture viscous effects such as viscosityinduced drag forces.
In contrast, CFD modelling is very versatile and can account for all non-linearities but requires longer computation times, which can further increase in the presence of porous or perforated structures. One approach to CFD modelling of fluid interaction with a perforated structure is to resolve the microstructural geometry explicitly. This approach is necessary when the detailed fluid flow behaviour needs to be analysed. It is often used for simulating smaller domains. For example, Filho et al. (2007) simulated water flow through perforated sheets as part of a pressurized water reactor. Mentzoni and Kristiansen (2019) developed a twodimensional NS solver for perforated plates in oscillating flow to study the sensitivities of geometrical porosity and motion parameters. This approach is also used to derive porosity coefficients from the model as done by Valizadeh et al. (2018), Chen et al. (2019) or Poguluri and Cho (2020) in the context of vertical perforated sheet breakwaters. It is also used when the perforations are large and the flow behaviour is highly anisotropic, for instance to study the interaction of irregular waves with a caisson breakwater with circular perforations as done by Lee et al. (2017). In these cases, a fine mesh with a large number of mesh cells is required for both a sufficient geometry resolution as well as to achieve high mesh-quality levels in order to support solver and scheme stability. This can lead to a very high computational demand and can make the microscopic approach prohibitive for larger domains and when the perforations are small and numerous.
A macroscopic approach can offer a more efficient alternative for engineering problems where large-scale force effects and the overall flow behaviour are of interest rather than details of the flow near the structure. Here, the micro-scale geometry of the perforated structure is not resolved explicitly but represented by a volumetrically averaged macro-scale model in a geometrically defined porous region. The porous structure can either be represented as a porous region with additional momentum source terms or as a surface with a pressure-jump condition. A macroscopic porous-media CFD approach addresses the limitations of both linear potentialflow models and microscopic CFD models. This approach does not require the restrictive assumptions of linear wave theory. It also simplifies mesh generation and reduces the number of cells required compared to microscopic CFD approaches.
To date-to the best of the authors' knowledge-the porous-media approach has been used for large porous structures with granular material characteristics such as rubble or block mound type break waters or dams (Higuera et al. 2014b;Lara et al. 2012;del Jesus et al. 2012;Jensen et al. 2014). In this paper we explore its use for thin perforated structures by using the volume-averaged Navier-Stokes equations as implemented in OlaFlow/IHFoam (Higuera et al. 2014a). These follow averaging procedures of Whitaker (1967), Slattery (1967) and Gray (1975) which assume the porous structure is homogeneous and with isotropic characteristics.
The closure terms that are needed as a result of the volume averaging process depend on the characteristics of the structure and its porosity coefficients are usually derived (semi-) empirically, often based on potential-flow or Euler assumptions. This area of research is still very active (Madsen 1974;Mei et al. 1974;Chwang and Chan 1998;McIver 1998;Li et al. 2006;Crowley and Porter 2012;Mentzoni et al. 2018) and the present work uses an existing concept rather than investigating an alternative formulation. A theoretical model based on the work of Molin (2011) is applied in the form of a turbulent drag term with a drag coefficient as a function of porosity n and a constant discharge coefficient δ, independent of flow conditions. This assumption is based on work by Tait et al. (2005) and Hamelin et al. (2013) and is further described in Sect. 2.3.
A two-dimensional (2D) model with a vertical perforated sheet and a three-dimensional (3D) model with a circular vertical perforated cylinder have been simulated and the horizontal force on the structures as well as the free-surface elevation at wave gauges around the cylinder have been analysed for a range of porosities and regular wave conditions. The CFD results are verified against results from a linear potential-flow model (Mackay and Johanning 2020) and validated against experimental results.
The paper is structured as follows: The governing equations including the pressure-drop model are explained in Sect. 2. A brief description of the experiments that have been conducted for validation is covered in Sect. 3. The specific setup of the CFD models is described in Sect. 4 for the 2D model and in Sect. 5 for the 3D model. Section 6 presents the CFD results and their verification and validation. A discussion of the results and the limitations of the model is given in Sect. 7, followed by conclusions in Sect. 8.

Numerical method
For the CFD modelling, the finite-volume based open source code OpenFOAM ® (The OpenFOAM Foundation v5) is used in combination with the OlaFlow/IHFoam (Higuera et al. 2014a) toolbox. OlaFlow provides boundary conditions and solvers tuned for wave modelling and incorporates a macroscopic porosity implementation within a Reynolds-Averaged Navier-Stokes (RANS) approach. Waves2Foam (Jacobsen et al. 2012;Jensen et al. 2014) would be an alternative tool with similar capabilities.

Governing equations
In the porous-media approach, the porous or perforated structure is represented by its volume-averaged properties via a geometrically defined porous region of volumetric dimension. In this continuous and homogeneous porous region, resistance source terms are applied as closure terms to the momentum equation. The source terms represent the macroscopic effects of the porous structure on the flow.
The volume averaging procedure in OlaFlow/IHFoam (and also in Waves2Foam) follows Whitaker (1967), Slattery (1967) and Gray (1975) and the comprehensive derivation can be found in Higuera et al. (2014a) or Jensen et al. (2014). The resulting Volume-Averaged RANS (VARANS) equations for incompressible, immiscible two-phase flow formulated for the intrinsic velocity U/n (the averaged velocity of the fluid inside the voids of the porous media) are the mass conservation equation, and the momentum equation, where the last three terms on the right-hand side are the closure terms that represent the volume-averaged pressure-drop across the structure. a, C f and c are porosity coefficients, further explained in Sect. 2.3. U is the averaged volumetric flow rate per unit area, sometimes referred to as Darcy-velocity, n is the porosity, X is the position vector in Cartesian coordinates, g the vector of gravitational acceleration and p * the pseudo-dynamic pressure. For two-phases, air and water in this case, ρ represents the weighted averaged density for each finite volume cell, calculated from the water and air densities ρ w and ρ a via ρ = αρ w +(1−α)ρ a . μ eff represents the effective dynamic viscosity and is defined as μ eff = μ + μ turb , where μ is the molecular viscosity and μ turb is the turbulent eddy viscosity. The latter represents small-scale stochastic turbulent effects by means of a turbulence model. Same as with the density, the viscosity μ represents the weighted average, calculated from the water and air viscosities μ w and μ a via μ = αμ w + (1 − α)μ a . Surface tension effects are represented with σ as the surface tension coefficient and κ as the curvature of the interface. α represents the phase-fraction field that is used to capture the interface as part of the volume of fluid (VOF) method. α = 1 corresponds to a cell full of water and α = 0 to a cell full of air. Hence, cells with α-values between 0 and 1 are part of the phase interface and free water surface, respectively. The present work uses the algebraic VOF method [with the multidimensional universal limiter for explicit solution (MULES) for boundedness preservation] adapted for a porous structures and formulated for the intrinsic velocity. The advection of α follows: where U c is an artificial compression velocity as part of an interface compression term (Berberovic et al. 2009). Inside the porous zone, the pressure-drop due to the porous medium is applied as a momentum resistance source and with a limited fluid amount. Outside the porous zone, the VARANS-equations are the same as the standard RANS equations and OlaFlow's porosity solver olaFlow reduces to OpenFOAM ® 's standard solver for incompressible, immiscible two-phase flow including an algebraic VOF interface capturing method interFoam. Both solve the unsteady (VA)RANS equations with the PIMPLE algorithm. Further information on those algorithms can for instance be found in Moukalled et al. (2016).
The pressure-drop coefficients a, C f and c are input parameters and are defined beforehand. They can either be obtained by experiments or from theoretical models. The present work uses a theoretical model, described in Sect. 2.3.
It is worth mentioning that in OlaFlow a, C f and c are implemented as scalar values which implies the porous media to be of isotropic nature with a constant pressuredrop in all directions within one porous zone. Although this is not strictly appropriate for the physical situation of a thin perforated sheet, this simplification is considered as a valid approximation for small sheet thicknesses and is in accordance with other simplifications such as the assumptions made for the derivation of the theoretical pressure-drop model. Furthermore, the computational effort is smaller compared to an anisotropic implementation where the porosity coefficients are implemented in tensor form.

Flow scales and turbulence-free modelling
In the present context of wave interaction with thin perforated sheets and cylinders, two main flow scales are of interest, the flow through the sheet openings on a small scale, and the flow around the perforated structure on a large scale. The characteristic non-dimensional numbers here are the Reynolds-number (Re), the Keulegan-Carpenter-number (KC) and the Stokes-parameter β. Re is defined as where u m is the amplitude of the oscillating flow velocity, L is a characteristic length and ν is the kinematic viscosity. Re can be expressed as Re = βKC, where and T is the oscillation period (equal to the wave period). The parameters KC, Re and β are calculated separately for the small and large scale flow. The relevance of the smallscale values comes in for the theoretical pressure-drop model in Sect. 2.3. In this case, the characteristic length scale L is set to the sheet thickness d. For the large-scale wave flow around a cylinder, L is set to the cylinder diameter D. For linear waves in deep water, the amplitude of the horizontal fluid velocity at the free surface is u m = ωH /2, where H is the wave height and ω is the wave angular frequency. The wave period is T = 2π/ω, substituting this into (5) gives: Although the waves in the physical experiment and simulations are not linear and the water depth is finite, Eq. (7) can still be used as an indicator of the flow regime.
In this context, large-scale turbulent effects such as horseshoe or lee-wake vortex generation must be considered. For an impermeable cylinder under wave loads, Fredsoe and Sumer (2006) note that the KC number has dominant influence on the flow regime, with β having a smaller effect. Furthermore, Chakrabarti (1987) notes that the relative importance of drag and inertial forces depends on the ratio cylinder diameter, D, to the wavelength, λ.
For the present setup, with D = 0.5m and the wave conditions listed in Table 1, we have 46946 ≤ β ≤ 149374, 0.19 ≤ KC ≤ 1.6 (calculated with (7)) and 0.048 ≤ D/λ ≤ 0.27. For these conditions, wave forces on an impermeable cylinder are inertia-dominated, with small or negligible drag. In particular, Sumer et al. (1997) found that there was no significant lee-wake vortex generation for KC < 4 and no horse-shoe vortex generation for KC < 6. For the case of a perforated cylinder, vortex generation is expected to be lower than for an impermeable cylinder.
Corresponding to the differentiation made between largescale and small-scale flow, a distinction is made between large-scale and small-scale turbulence. The former concerns anisotropic effects such as horse-shoe or lee-wake vortex shedding around the cylinder which can be simulated directly to a large extent by a sufficiently fine temporal and spatial resolution. The latter refers to isotropic turbulence which in the RANS-approach is not simulated but modelled and represented by a turbulence model. For the specific present conditions neither large-scale nor small-scale turbulence are expected to have significant effects. Therefore, no turbulence closure model is applied and μ eff = μ in (2) (μ turb =0), respectively, for all models. This is elaborated further as follows.
The effects of local turbulence generation and related losses caused by flow separation across the porous barrier since well as increased shear stresses around the fluidstructure interface are taken into account by the theoretical pressure-drop model and the structure's representation as continuous medium. The application of a turbulence model would represent an extra contribution to the applied momentum resistance. In general (although this is not relevant for the present work), this gets even more important when the applied porosity coefficients in the model are obtained from experiments since they already include all physical turbulent effects. This was highlighted by Jensen et al. (2014) for the CFD modelling of porous dams. Furthermore, the wakes close to the structure are quickly regularized and homogenized away from the porous barrier since the openings can be assumed to be small (Molin 2011). This means that in the clear flow region (i.e. away from the structure and sea bed), turbulence can be neglected for the propagation of non-breaking waves. Also, for the present case of wave propagation over a smooth, flat bottom, no significant bottom boundary layer is expected to develop, and a free-slip condition is applied at the bottom boundary in all present models.
Since the contribution of large-scale turbulent forces are expected to be negligible and the pressure-drop model takes small-scale turbulence effects into account, no turbulence model is applied to the governing equations. In other words, the full Navier-Stokes equations are used but due to the relatively coarse spatial resolution it cannot be considered as Direct Numerical Simulation (DNS) here. The use of the full Navier-Stokes equations with a relatively coarse mesh is considered as sufficient for the present conditions and possibly also for other marine engineering problems where the turbulence levels are minimal. However, a turbulence model may be important for more general cases of wave-structure interaction such as wave breaking, cases where significant boundary layer effects exist or when large-scale turbulent vortex generation and shedding must be accounted for.

Theoretical pressure-drop model
The pressure-drop, p, across a porous barrier can be written as a function of the flow velocity (usually formulated with the Darcy-velocity U): where a, C f and c are porosity or pressure-drop coefficients. The first linear term of (8) represents viscous friction effects, typically dominant for low Re-numbers. The second term represents turbulent drag of high Re-number regimes. The third transient term accounts for the acceleration of the fluid through the voids or openings. The relative importance of the three components depend on the flow regime. In the context of granular material these terms are typically referred to as Darcy-term (Darcy 1858) for the linear term, as Forchheimerterm (Forchheimer 1901) for the nonlinear drag component and as Polubarinova-Kochina-term (Polubarinova-Kochina 1962) for the transient term, with corresponding formulations for the porosity coefficients.
In the context of gravity wave induced flow through thin porous structures, turbulent losses dominate over viscous losses and the linear term can be neglected (Sollitt and Cross 1972). Furthermore, the transient term tends to zero when the size of the openings is small relative to the wavelength (Molin 2011). Consequently, the linear term and the transient term in (8) are set to zero in this work and the pressure-drop across the perforated barrier is represented by a drag term with the drag or friction coefficient C f .
For this work, a formulation by Molin (2011) is applied for C f , defined as where δ is an empirically defined discharge coefficient, usually in the range of 0.4-0.5. For oscillatory flow the discharge coefficient δ and hence the drag coefficient C f vary with the KC-number, calculated for the small-scale flow across the structure. Tait

Wave modelling
For the wave modelling, OlaFlow provides wave generation boundary conditions and an active wave absorption treatment based on common potential-flow wave theories. The active wave absorption method (AWA) minimizes wave reflections via correcting the incident free-surface elevation by imposing the incident velocity profile in the opposite direction to the wave propagation (Higuera et al. 2013). This condition aims to cancel out the incoming waves at both the wave generation boundary as well as the pure absorption boundary. In the initial version, the correction velocity profile was implemented based on shallow water theory (h < λ/20 with h as the water depth and λ as the wavelength) where a constant velocity profile is applied. As this does not account for the variable velocity profiles for intermediate (λ/20 < h < λ/2) and deep water ranges (h > λ/2), the performance of this method decreases with increasing water depth. Recently, an improved version was released that enables the use of a more general formulation for any water depth (Higuera 2019) for the pure absorption boundary. Still, at the generation boundary a constant velocity profile is assumed for absorption. The initial Fig. 1 Photos of the experiments: a wave interaction with the perforated sheet which is mounted onto a rigid steel frame (brown color) facing towards the wavemaker, b the cylinder model in the empty wave flume with the separation wall on the right (yellow) and c under wave interaction. The direction of wave propagation is from the left to the right in the case of the sheet and from the right to the left in the case of the cylinders version is referred to as shallow-water AWA (SW-AWA) and the latter as extended-range-AWA (ER-AWA). Both the SW-AWA and ER-AWA are used in the present work.

Physical experiments
Experiments were conducted at Dalian University of Technology, China to measure wave loads on thin perforated sheets and cylinders. The wave flume used is 60 m long and 4 m wide, has a single piston wavemaker at one side and a beach at the other end. A water depth of h = 1 m was used for all studies. Along the central tank section a vertical wall with a length of 13.20 m was positioned in order to split the tank into two separate test sections where one contained the sheet and the other one the cylinder models. The section with the sheet had a width of 1 m and the section with the cylinder a width of 3 m. Load cells were mounted at the bottom and top ends of the structures. The flat sheets occupied the full width of the channel and full height of the water column and were mounted on a rigid frame. Figure 1 shows photos of the experiments for both the sheet and cylinder. An example of the wave interaction of the porous sheet during the experiments in the wave flume is shown in Fig. 1a where one can see how the sheet is mounted onto the brown rigid frame. The configuration of the cylinder model in the empty tank is shown in Fig. 1b and during wave tests in Fig. 1c.
The perforated sheets and cylinders had circular holes arranged in a regular square grid of side s and hole radius r h , such that the porosity is defined by n = πr 2 h /s 2 . In the experiments, a number of geometrical parameters such as sheet thicknesses (3 mm ≤ d ≤ 10 mm), hole separation distances (25 mm ≤ s ≤ 100 mm), outer cylinder diameters (0.375 m ≤ D ≤ 0.750 m) and porosities (0.1 ≤ n ≤    Table 1.The input waves cover nondimensional wave numbers of 0.61 < kh < 3.34 and wave steepnesses of 0.05 ≤ k A ≤ 0.20 where k = 2π/λ is the wave number and A is the wave amplitude. All conditions listed in Table 1 were used for the two-dimensional model with the porous sheets. For the three-dimensional (2D) model with the porous cylinder only a subset with k A = 0.1 (B03-B06) was used.
The experimental data are suspected to be unreliable for the wave conditions A09, A10, A12 and A13 for both the sheet and cylinder tests for all configurations. For the sheet with porosity of n = 0.1 the wave condition A05 is suspected additionally. For these conditions there is a large deviation between experimental results and numerical results from both the CFD and potential-flow model (see Fig. 9a). It is suspected that the discrepancy is caused by a raised section of the floor in the physical tank which was required to submerge the load cells beneath the structures. The raised section had a total length of 11 m and the side facing the wave maker had a gentle slope with the height increasing over the first 5-0.17 m above the tank floor. The rear sloped section was 2 m in length and the central flat section was 5 m long. Below the structure, a water depth of h = 1 m was maintained (next to the raised floor the water was 1.17 m deep). A sketch of the raised floor can be found in Feichtner et al. (2019). Since differences in the measured wave height between the deeper section and the raised section have been observed, it is hypothesized that it generates a diffracted wave field that interacts with the reflected wave from the porous sheet and cylinders. The effect is discussed further by Mackay et al. (2019) who also found discrepancies with a potential-flow model.

CFD model development in 2D
This section presents the specific setup of a 2D model containing a vertical thin perforated sheet. This model was used to evaluate the minimum mesh requirements and to identify the most influential settings more extensively. The outcome of the model development for the 2D model is then used for the setup of the 3D model with the perforated cylinder.

Computational domain and boundary conditions
The numerical wave tank was set to a length of 26 m in order to contain a minimum of approximately 2.5 wavelengths based on the longest wavelength of λ = 10.36 m (A13). The domain height is set to 1.3 m and the water depth to h = 1.0 m in accordance with the experiment. The porous sheet was positioned at the flume centre at x = 13 m in order to achieve the same amount of time for reflected waves travelling back from both ends of the domain. Despite variable sheet thicknesses in the experiments, a constant sheet thickness of d = 10 mm was used for simplicity for all CFD cases and the pressure-drop coefficient C f has been calculated via (9) correspondingly for porosities n = 0.1, 0.2 and 0.3. It was verified that the CFD results are independent of the sheet thickness d in the CFD model for a range of 5 mm ≤ d ≤ 20 mm as long as the drag coefficient per unit length has been adapted accordingly. For all comparative models, a minimum of 16 cells per thickness (N x /d) was applied as a result of the mesh convergence study below in Sect. 4.2, see the summary in Table 3. Consequently, the cell size across the sheet has changed. The top boundary was set to an atmospheric condition which allows both air and water to flow out but only air to flow back into the domain. As stated in Sect. 2.2, a free-slip wall condition was applied at the tank bottom. This corresponds to an unrestricted tangential velocity component and a normal velocity component set to zero. Second order Stokes-waves were imposed at the generation boundary, defined by the freesurface elevation η(t) and the horizontal and vertical velocity profiles u x and u z respectively. η(t) is given by where θ = kx −ωt +ψ is the wave phase, k the wave number, ω is the angular frequency, ψ is the wave phase shift, x is the horizontal coordinate and t is the time. u x and u z are defined as and where z is the vertical coordinate of the water column. The waves were generated at the left tank boundary and propagate towards the pure absorption boundary on the right-hand side. The wave input parameters cover the whole range of regular wave conditions as shown in Table 1. The active wave absorption (AWA) treatment was activated for both the generation and the pure absorption boundary. A sketch of the model configuration including all boundary conditions and wave gauge (WG) positions is shown in Fig. 2 (the sketch includes the cylinder setup as well). The position of the WGs relative to the centre of the structure in the model correspond to the WG position in the physical tank. Further information on the applied numerical solver and scheme settings, mesh cell numbers, simulation and execution times can be found in Appendix A.
The following sensitivity studies were carried out for a shorter wave (A02) and a longer wave (A09) from the range of wave conditions with k A = 0.05 as shown in Table 1.

Spatial discretisation and mesh independence
Accurate wave propagation and a low level of wave reflections is recognised as an important starting point for the present wave-structure interaction problem. However, as a wide range of wave conditions with wavelengths in the range Mesh convergence studies were carried out in two steps for the waves A02 and A09 (see Table 1). Firstly, mesh independence of the free-surface elevation was examined for WG00 at the centre of the empty flume. This is at an x-position of 13 m where the porous sheet would be placed. Secondly, the convergence of the force on the porous sheet was examined concerning the mesh size around the sheet in terms of number of cells per thickness N x /d. The baseline mesh was generated using a uniformly structured mesh with regular hexahedral cells with a size of l x = l z = 20 mm, where l x is the cell length in horizontal x-direction and l z the cell height in vertical z-direction. In the course of the following mesh independence study the mesh was successively refined locally and separately for the free-surface region and the area around the sheet, as explained further below. Figure 3 shows a close-up of the mesh and wave interaction with the sheet.
For transient wave-structure interaction modelling the achievement of monotonic convergence is considered as a challenging requirement since an exact periodically steady behaviour is rarely observed. However, the mesh convergence metric was specified as the averaged wave amplitude A at WG00 in the empty tank as well as the averaged amplitude of the horizontal force on the sheet F. For the averaging process, the initial wave ramp-up time and first transient section of the time series were removed and the nearly periodically steady time series section was then cropped to a whole number of at least 10 wave periods. The average amplitude was then calculated with the standard deviation (STD) of the time series for both the wave elevation and force on the sheet. The averaged wave amplitude is obtained via A = √ 2 * STD(η(t)), where η(t) is the free-surface elevation and the averaged force amplitude on the sheet via is the force amplitude over time. For most cases, the CFD time series contain wave reflections, as the reflected waves have travelled back before 10 wave periods have passed WG00 at the tank centre.

Free-surface elevation in the empty tank
Starting off from the baseline mesh with the uniform cell dimensions stated above, the free-surface region was successively refined locally. The refinement region at the freesurface covers the maximum target wave height H of 0.26 m of the range of wave conditions (see Table 1). Table 2 shows a summary of the mesh independence study for the freesurface region looking at the surface elevation at WG00 in terms of cells per wave height (CPH) and cells per wavelength (CPL). The measured averaged wave amplitudes A are presented including the successive percentage change between the meshes. The values are compared to the target wave amplitudes based on Stokes 2nd-order input A input . It can be observed in Table 2 that no monotonic convergence could be achieved for the averaged wave amplitude A. The differences between the target amplitudes A input based on Stokes 2nd-order input and the measured A are 4-9% for the wave A02 and around 1% for the wave A09. As the absolute wave height of A02 is relatively small compared to the cell heights, the accuracy was deemed to be acceptable nonethe- Table 2 Mesh convergence study for the free-surface region for the wave amplitude A (mm) at WG00 in terms of cells per wave height (CPH) and cells per wavelength (CPL)-selected mesh is mesh 3 with l x = 10 mm and l z = 5 mm less since the mesh aims to cover the whole range of wave conditions. Mesh 3 with l x = 10 mm and l z = 5 mm was selected for the following studies. This corresponds to a minimum number of 6 CPH for the smallest wave height (A01) and 52 CPH for the largest wave height (C04) and 188 CPL for the shortest wave (A01) and 1036 CPL for the longest wave (A13), see Table 1. The reasons for differences between the measured and target conditions are identified as a combination of highfrequency ripples at the air-water interface, numerical scheme settings, wave reflections and the ratio of wavelength to domain length as well as plate position. These influencing factors are discussed further below and their effects are incorporated within an uncertainty estimation in the Appendix A.

Force on the porous sheet
The second step of the mesh independence study was focused on the force on the porous sheet which is the main parameter of interest. The force was calculated as the difference of the integrated pressure values at both sides of the sheet, evaluated at the cell face patches that confine both sides of the porous zone that represents the sheet. Therefore, the mesh was locally refined close to the porous sheet, with the size increasing gradually away from the sheet to achieve a smooth transition in cell aspect ratios. The cell lengths in x-and zdirection were kept equal (l x = l z ) and the refinement region around the free-surface was extended in all horizontal and vertical directions to cover possible wave-structure interaction. Table 3 shows a summary of the mesh convergence study for the waves A02 and A09 for the averaged force amplitude F in terms of absolute cell size l x = l z and number of cells per sheet thickness N x /d . Table 3 shows that no clear monotonic convergence could be achieved for the averaged force on the sheet F. However, it is clear that one cell per sheet width (N x /d=1) results in significantly different forces than finer meshes for both wave conditions. For the wave condition A02, the effect of the number of cells per sheet thickness is rather oscillatory but for the longer wave A09 the force increases nearly monotonically with a successive decrease in relative change with increasing refinement. Excluding the results of the model with mesh 1 (N x /d = 1), the mean values of the averaged force amplitude are F = 0.635 N for wave A02 and F = 5.55 N for wave A09 and the maximum deviation is 4.1% for A02 and 6.4% for A09. Sections of the force time series were inspected visually. The converging behaviour around the minimum and maximum amplitudes of the time series indicates mesh convergence towards mesh 5 with N x /d = 16.
. The accuracy from mesh 5 was deemed acceptable and a mesh using 16 cells per thickness equalling a cell size of l x = l z = 0.625 mm for the present sheet thickness of d = 10 mm was selected for all further models. Mesh 5 is shown in Fig. 3.

Temporal discretisation
As often used for CFD modelling of wave flumes, automatic time stepping was employed which adjusts the time step to give a maximum Courant-Friedrichs-Lewy (CFL)-number.
The CFL-number is calculated as where t is the time step size, V is the cell volume and φ the volume flux across the cell faces. It is interpreted as the number of cells that a scalar quantity traverses during one time step. For instance, for a maximum CFL-number of 0.5 the fluid crosses a maximum of half a mesh cell during one time step. In this case, the spatial and temporal convergence studies are directly linked since the time step decreases correspondingly with a decreasing cell size for a constant maximum CFL-number.
For the current studies with the porous sheet it was found that a maximum CFL-number of 0.3 is beneficial for a stable solver run. Obviously, this maximum CFL is strongly influenced by the fine mesh region around the sheet. This entails that the CFL-numbers in the free-surface region are automatically around 0.1 which is considered to be sufficient for accurate wave propagation. This magnitude is in agreement with similar work on propagating waves, for example work by Roenby (2017) or Larsen et al. (2019). As this CFLnumber is considered as small already, no separate temporal convergence study was conducted as a consequence.

Factors of influence
As observed above, no monotonic convergence behaviour could be achieved for both the averaged wave and force amplitudes A and F in terms of spatial discretisation. The main causes are identified as numerical scheme and solver settings, wave reflections at the boundaries and highfrequency ripples at the air-water interface in the pure wave propagation region. In particular, wave reflections at the boundaries combined with unfortunate ratios of wavelength to domain length as well as sheet position can have significant effects on the results. In this context, the wave absorption method is a crucial factor.
For each of these factors a sensitivity analysis and a corresponding uncertainty estimation was performed in the Appendix A. An estimated overall mean uncertainty of U = ±16.4% was obtained for the CFD results for the horizontal force on the porous sheet.

CFD cylinder model setup
Based on the development of the 2D model in Sect. 4 and its outcomes, a 3D model containing a vertical perforated cylinder was created. In order to minimize the computational effort for the 3D simulations, the domain size as well as the number of mesh cells have been reduced. The length of the 3D numerical flume is 16 m, the domain width is 3.0 m, the domain height is 1.3 m and the water depth is set to h = 1.0 m. A porous cylinder is positioned at the centre of the flume at an x-position of 8.0 m and y-position of 1.5 m. The cylinder was set to a constant outer diameter of D = 0.5m (and is hollow in the middle). Although the cylinder thickness was 3 mm in the experiments, the thickness in the numerical model was set to d = 5 mm for mesh generation reasons and the porosity coefficient was calculated correspondingly for porosity values of n = 0.2 and 0.3.
As in the 2D model, the top boundary was set to an atmospheric condition and the bottom boundary to a free-slip wall. The front and back boundaries were set to a free-slip wall condition in the 3D model as well. From the wave conditions in Table 1, only a subset (B03-B06) was used in the 3D cases. The ER-AWA was used at both the wave generation boundary at the left-hand side and the pure absorption boundary at the right-hand side. The reflection coefficients were below 0.5% for equivalent empty 2D models for all conditions, estimated using a three-gauge method (Mansard and Funke 1980).
The computational mesh is generated solely out of hexahedra and set up to meet the basic requirements for the specific conditions. Since it was determined in Sect. 2.2 that no large-scale turbulent vortex generation is expected, refinement regions were employed along the air-water interface and in the vicinity of the porous cylinder but not in possible horse-shoe or lee-wake regions. Following the 2D model setup, the cell dimensions were set to l x =l z = 20 mm in the clear flow region away from the structure. Along the water surface the cell size was set to l x = 10 mm and l z = 5 mm which corresponds to a minimum of 21 CPH and 335 CPL (B03). In contrast to the 2D model, the local mesh refinement around the structure was only performed in horizontal direction in order to reduce the total mesh cell number. The number of cells per sheet thickness in radial direction were set to N r /d = 16 (l r = 0.3 mm), the cell dimension in vertical direction was kept to l z = 5 mm along the free water surface and l z = 20 mm for the rest of the water column (in contrast to l x = l z = 0.625 mm for the 2D model). As for the 2D model, automatic time stepping was employed which adjusts the time step to give a maximum CFL of 0.3.
Illustrations of the simulation for the case with the wave B06 and a cylinder porosity of n = 0.2 are shown in Fig. 4. The water surface is shown for a point in time of t = 40 s for the whole domain (Fig. 4a) and a close-up of the cylinder including the mesh (Fig. 4b).

Results
The horizontal force on the structures is specified as the main parameter of interest and has been analysed in two ways for both the sheet and cylinder. Firstly, the force time series f (t) Fig. 4 3D cylinder model for the wave condition B06 and porosity n = 0.2, showing the water surface and the cylinder. The water isosurface (α = 0.5) is colored in the magnitude of velocity within the limits of 0 and 1 m/s from the CFD results and the experimental results have been compared for the sheet and cylinder. Since the lengths of the time series are different due to the different physical and numerical tank lengths, the experimental time series have been cropped to correspond with the CFD time series. Secondly, the variation of the force with wave frequency and steepness is examined, for a range of values of porosity. The force on the vertical 2D sheet has been calculated for a wide range of parameter combinations in order to assess the validity of the theoretical pressure-drop model for various porosities and wave conditions. The results of the cylinder models are assessed for a subset of conditions with the focus on investigating the validity of the present approach for 3D structures.
For the 3D model containing the cylinder, the free-surface elevation around the cylinder was compared for a number of WGs to verify the capability of the model to replicate the mean flow behaviour.

Time domain results
The time series of the horizontal force on the structures f (t) are compared for the CFD models and experiments. Differences are expected for two reasons. Firstly, the lengths of the physical and the numerical wave tank are different which produce differences in the interaction between the incident and reflected waves. All following figures contain a vertical line that indicates the point in time when the wave has travelled across the numerical wave tank, been reflected at the tank end and reached the tank centre again. Secondly, the distance from the wave generation boundary to the sheet/cylinder in the CFD model is different to that in the physical tank. Moreover, the ramping-up of the waves differs between the CFD and the experiments. This results in a difference in the initial transient part of the record. Consequently, the differences in the initial section of the free-surface elevation directly results in deviations of the force response. However, both the CFD and experimental results reach an approximately periodically steady state. Figure 5 shows a comparison between the CFD and experimental results of the time series of the horizontal force f (t) on the perforated sheet for the cases n = 0.3 with wave A02, n = 0.2 with wave A07 and n = 0.1 with wave A08. Since the physical sheet in the experiments had a width of 1 m, the experimental results were scaled accordingly.

2D sheet model
Differences between the CFD and experimental time series for the initial transient sections can be observed and are more significant for the longer waves, in particular for the case with n = 0.1 with wave A08. Also the differences in the reflection behaviour due to different tank lengths can be observed for all example cases. Reflections in the experimental results can be observed for the case with n = 0.1 with wave A08 and noticeable reflections in the CFD results are As this is not related to the porosity representation in the model, the agreement between the CFD and experimental results is considered as relatively good, in particular for the case with n = 0.2 and wave A07 after about t = 30s.
The shape of the measured force time series is well reproduced by the CFD model for all combinations of porosity n and wave condition, clearly showing the non-sinusoidal profile, due to the quadratic pressure-drop.

Cylinder model
For the 3D model, both the force and wave elevation time series are analysed comparing the CFD and experimental results. Figure 6 shows the time series of the horizontal force on the perforated cylinder. Results are presented for the cases with n = 0.2 with wave B03, n = 0.3 with wave B04, and n = 0.3 with wave B06.
The CFD and experimental force results agree very well with the CFD model being capable of reproducing the shape and amplitude of the horizontal force on the cylinder. The differences are related to the same factors affecting the 2D model, discussed above. The initial transient period differs significantly due to the differences in the wave ramp-up. The reflections are more noticeable for the experimental results, in particular for the second half of the time series between about t = 25 s and t = 50 s. In contrast, the CFD results exhibit a nearly periodically steady force response with negligible reflections from about t = 15 s onwards.
Furthermore, the time series of the surface elevation, η(t) and averaged amplitudes A are analysed for the set of WGs shown in Fig. 2. Firstly, the time series are presented for the WGs closest to the structure which are WG03 at x = 6.4 m (1.6 m before the cylinder centre) and WG04 at x = 9.21 m (1.21 m after the cylinder centre). Figure 7 shows a selection of time series and amplitude spectra from the CFD model and experiments. The amplitude spectrum is defined as where S( f ) is the variance density spectrum and f is the frequency resolution.
The agreement is very good for all cases shown. The surface elevation was also accurately reproduced at the WGs further away from the structure (not shown). The amplitude spectra indicate that the measured fundamental frequencies and higher harmonics are well captured by the CFD model. Small deviations are present with the computational results tending to underpredict the experiment. The devia- Next, the normalized averaged wave amplitudes, A/A input , where A input is the target CFD wave input, are analysed for all WG positions. The comparison between experimental and CFD results is shown in Fig. 8 for the cases with n = 0.2 with the waves B03 and B04, and n = 0.3 with wave B06. The spatial variation in the wave amplitude agrees reasonably well between the CFD model and experiments. The agreement is better for the shorter wave conditions B03 and B04. The longest wave B06 exhibits the largest deviations at WG01 which is located furthest away from the cylinder and closest to the wavemaker. The experimental results do not contain wave reflections since the averaged results were estimated from the portion of the record after the transient response and before the reflected waves had arrived back from the beach. The CFD results are estimated from the portion of the time series including reflected waves, but due to the low reflection coefficient in all the 3D CFD models (R < 0.5% using ER-AWA), the influence of reflections on the CFD results is minimal.
The good agreement between the numerical and experimental WG results indicates that the present porosity representation is able to capture the large-scale characteristics of the wave-structure interaction over time. It is shown that the present approach is valid not only in 2D but also for 3D structures.

Frequency domain results
In this section we examine the ability of the CFD model to predict the variation of the force amplitude with wave frequency and steepness. For reference, the CFD and experimental results are also compared to results from a potentialflow model (Mackay and Johanning 2020). The potentialflow model is based on linear wave theory. It uses the same quadratic pressure-drop formulation as the CFD model, but the time dependence in the velocity-squared term is linearised using Lorenz's principle of equivalent work. The non-linearity of the wave conditions increases with decreasing relative water depth kh and increasing wave steepness k A. The accuracy of the linear potential-flow model is expected to decrease as the wave non-linearity increases.
The results are presented in terms of a normalised force, f x , defined as the ratio of the measured force, F porous , to the force on the equivalent solid structure in the same wave conditions: The value of F solid is based on linear wave theory. For the flat sheet it is and for the cylinder-based on MacCamy and Fuchs (1954)it is where w is the sheet width and r is the cylinder radius. H (1) 1 is the Hankel function of the first kind of order 1 and the prime denotes differentiation with respect to the arguments. Figure 9a shows the normalised force on the sheet for a constant target normalized wave steepness of k A = 0.05 and for porosity values of n = 0.1, 0.2 and 0.3. For n = 0.1 the whole range of wave conditions was simulated, for n = 0.2 and Overall, the averaged normalized force amplitudes f x are in reasonably good agreement between the experimental, potential-flow and CFD results. Generally, the potential-flow and CFD results agree well for most frequencies. The experimental results exhibit more scatter.

2D sheet model
The results match particularly well for the larger porosity values (n = 0.2 and n = 0.3) with the theoretical pressuredrop model applied in the CFD models being capable of reproducing the variation with kh to a satisfactory extent. The deviations increase for the smaller porosity value of n = 0.1. Overall, the results indicate that the theoretical model with a constant discharge coefficient δ is capable of reproducing the change in the force with wave frequency and steepness as well as sheet porosity with sufficient accuracy.
The potential-flow and experimental results are not influenced by wave reflections. The experimental results are derived from the portion of the record after the initial transient and before the reflected waves arrive back from the beach. In contrast, the averaged CFD results are affected by wave reflections as well as high-frequency ripples at the air-water interphase. These and other factors are taken into account by an error and uncertainty analysis in Appendix A, resulting in a mean uncertainty of the CFD results of U C F D = ±16.4%.
The absolute mean uncertainty of the normalized experimental results was estimated as U f x = 2.1%, based on analysis of repeat tests where repeatability was proven (the analysis is outside the scope of this paper). The uncertainties of the experimental and CFD results are represented by error bars in Fig. 9.
Noticeably, the CFD results lead to larger f x values compared to both the potential-flow and experimental results for nearly all cases. The difference ranges between − 3.5 and + 21.9% relative to the potential-flow results and between − 10.5% and + 21.9% relative to the experimental results (after removing the unreliable data). Figure 10 shows the normalised force on the perforated cylinders for a constant wave steepness of k A = 0.10 and porosities n = 0.2 and 0.3. The error bars indicate the level of uncertainty of the CFD results that was estimated for the 2D model (U = ±16.4%).

Cylinder model
The agreement between the numerical and experimental results is very good, in particular between the potential-flow and CFD results with deviations ranging between − 1.5 and + 2.9% of the CFD results relative to the potential-flow results. The experimental results exhibit larger scatter and the differences between the CFD results and the experimental results range between + 1.9 and + 12.4% of the CFD relative to the experimental results.
The good agreement between the CFD, potential-flow and experimental results supports the idea that the present theoretical pressure-drop formulation with the use of a discharge coefficient δ independent of any wave parameters gives good results not only in 2D but also for 3D wave-structure interaction.

Discussion
The overall agreement between the CFD, potential-flow and experimental results is reasonably good, with the CFD model being able to reproduce the horizontal force on both the 2D and 3D structures as well as the influence of the structure on the surrounding wave field. The limitations of the present model and its porous-media representation are discussed below.
The porous-media approach, by its nature, is not capable of reproducing the detailed flow field and free-surface interface very close to the porous structure. A volume-averaged porous zone results in a smoother flow across the porous barrier. Small-scale effects of the flow close to the openings in the perforated barrier, such as water jets or vortex shedding, are not replicated directly in the CFD model. Instead, the effects of these small-scale features are parameterized in terms of the pressure-drop model. The flow patterns near the perforated structure in the CFD model can be observed in Fig. 3 which shows a snapshot of the 2D numerical tank as well as a close-up of the wave interaction with the sheet. The flow patterns of the CFD model can be compared to the real flow patterns as shown in a photo of the wave interaction with the sheet in the experiments in Fig. 1a. The pressuredrop model is derived under the assumption of steady flow through the openings, so it is not expected to be able to repli-cate the effect of the porous barrier on the flow around the free-surface when there is a large difference in water level on either side. However, for the wave conditions considered, the pressure-velocity approximation used here appears to capture the large scale effects of the flow (quantified in terms of the horizontal force on the porous sheet) reasonably well.
In the current models, porosity is approximated with isotropic characteristics which means that the pressure-drop is applied volumetrically with homogeneous and constant properties in all directions. More realistically, the pressuredrop could be applied proportional to the component of the flow normal to the porous boundary only. Whether an anisotropic volumetric porosity representation would avoid the observed over-prediction of the force for the 2D models, will be investigated in future work. However, an anisotropic implementation would increase the computational cost and the isotropic simplification can be considered advantageous in this context.
With a macroscopic porosity representation no actual boundary exists in the model where a boundary layer condition could be applied. Consequently, no boundary layer effects can be taken into account with default boundary conditions. For the present conditions, viscous boundary layer effects seem to be negligible, and the results are reasonably good. However, for more general cases a different approach may have to be considered.
For the evaluation of the 2D model, the same domain size and spatial discretisation were used for the whole range of wave conditions. Dependent on the wavelength, the domain length can have a significant effect due to wave reflections at the boundary and phase differences between the incident and reflected waves. In this work, the key effects are taken into account by uncertainty values and the overall results are considered as satisfactory. If higher accuracy is required, the model setup and numerical settings could be optimized separately.
The current porous-media approach gives good results for the wave interaction with simple perforated structures in two and three dimensions. It needs to be highlighted that the present models are specifically set up for the geometrical parameters and wave conditions that have been present in the experiments. For more general cases, the setup would need to be adapted dependent on the dominating physics.

Conclusions
CFD modelling using a simplified porous-media approach is presented for wave interaction with thin perforated structures. Perforated sheets and cylinders are represented by their macro-scale effects on the fluid as volume-averaged homogeneous volumetric porous zones. The volume-averaged NS-equations (as in Higuera, 2017) are used where poros-ity is taken into account with a limited amount of fluid inside the structure and the application of a resistance source term as closure term. For the present case of a thin perforated barrier, the theoretical model for the porosity resistance is applied as a quadratic drag term with a drag coefficient as a function of porosity following Molin (2011) and a constant discharge coefficient (thus independent of any wave parameters).
The current work aims to model large-scale effects and the mean flow behaviour of the wave-structure interaction. In this context, two simplifications are made in the setup. The porous media is assumed to be of isotropic nature and no turbulence closure is applied in the governing equations.
The assessment focuses on the horizontal force on the perforated structures and the free-surface elevation for a number of wave gauges near the perforated cylinder. The investigations covered a range of regular wave conditions with various wave frequencies and wave steepnesses as well as multiple sheet porosities. The CFD results have been compared against results from a linear potential-flow model as well as against experimental results, in both the time and frequency domain.
Overall, the experimental, potential-flow and CFD results show good agreement. The CFD models are capable of reproducing the horizontal force response as well as the mean flow behaviour with satisfactory accuracy. A slight overprediction of the force amplitude was observed for the 2D model of the perforated sheet, with the over-prediction being larger for lower porosities.
The good agreement indicates that the volume-averaged NS-equations without the use of a turbulence model are sufficiently accurate when no large-scale turbulent effects such as horse-shoe or lee-wake vortex generation are expected.
The results indicate that a porous-media approach with the assumption of isotropic material characteristics (as implemented in OlaFlow/IHFoam) is capable of reproducing the large-scale interactions between waves and thin perforated structures. It is shown that the applied theoretical pressuredrop model as volume-averaging closure term is capable of replicating the characteristic quadratic pressure-drop of the flow across thin perforated barriers for the range of regular wave frequencies, wave steepnesses as well as sheet porosities considered.
As with general wave-structure interaction problems in CFD, the investigations conducted here highlighted that aspects of the model setup, such as domain length, wave absorption method and spatial resolution have a significant influence on the results.
This work serves as a base for future extensions to more complex fixed and floating offshore structures under both operational and extreme wave conditions. It can provide a starting point for engineering applications where global forces and the mean flow behaviour are of the main interest.

A.2 Sensitivity study and uncertainty estimation
A sensitivity study was performed for the 2D model with the sheet. The main factors that influence the results are wave reflections, the ratio of wavelength to domain length and sheet position, solver and scheme settings, and highfrequency ripples at the air-water interface. The ripples are considered to be related to the numerical settings and are not analysed separately. The other aspects and their effect on the results are estimated below.

A.2.1 Influence of wave reflections
To assess the influence of the wave absorption method, the reflection coefficients R were compared for models set up with the SW-AWA and ER-AWA. This was performed on a subset of wave conditions (A02, A03, A04, A09) with k A = 0.05. The reflection coefficients R have been calculated with a three-gauge method (Mansard and Funke 1980) using WG02, WG00 and WG05 as shown in Fig. 2 for an empty numerical flume. The results are shown in Table 5. As expected, the performance of the SW-AWA decreased for the shorter waves and increasing non-dimensional wave numbers kh, respectively. The ER-AWA method improved the absorption performance for all wavelengths down to reflection coefficients of R = 1-2%. The improvement increases with decreasing wavelength.

A.2.2 Influence of the domain length
The size of the CFD domain is a trade-off between computational cost and the accuracy of results, with shorter domains resulting in shorter durations before reflected waves influence results. For the current studies, the baseline domain length was set to 26 m for the 2D model. To estimate the influence of the domain length, a few simulations with tank lengths of 30 m, 34 m and 38 m were performed for the wave conditions A02 and A09 including the porous sheet. The averaged force amplitudes F for various tank lengths are shown in Table 6. The deviation in relation to the base length of 26 m as well as the ratio X /λ for each domain length, X , and wave period are included in the table.
The deviations of the force amplitude F are less than 1.7% for the shorter wave (A02), but up to 14.1% for the longer wave (A09). This behaviour is not related to the reflection coefficient R which is larger for shorter waves (R A02 = 21%) and smaller for longer waves (R A09 = 6%) but to the position of the structure within the standing wave that is created in the numerical wave tank. Larsen et al. (2019) reported a relatively strong influence of the numerical setup in the course of the assessment of OpenFOAM ® 's solver interFoam for the use of simulating pure wave propagation. Among other aspects such as the maximum CFL-number, they identified the time discretisation scheme (ddt), the convection scheme (div(rhoPhi,U)), the surface normal gradient scheme (snGrad, e.g. used in the discretisation of the Laplacian) and the interphase-capturing method as critical. Although these factors were identified for wave propagation without a structure, their effect was estimated for the present work dealing with wave-structure interaction. Variations of time discretisation and convection schemes as well as the use of an alternative free-surface capturing method [isoAdvector by Roenby et al. (2016)] were used for comparisons and were incorporated in the uncertainty estimation below. Similar to Larsen et al. (2019), high-frequency ripples in the air-water interphase have also been observed in the present study.

A.2.3 Influence of solver settings and numerical schemes
It may be possible to optimise the solver and scheme settings for each wave condition. However, for consistency the settings were kept constant over all model runs. The settings used are listed in Table 4.

A.2.4 Uncertainty calculation
As a comprehensive error analysis is complex with many subcomponents and usually requires monotonic convergence, a simpler approach is adopted here. In addition to the factors discussed above (wave reflections, domain length, scheme settings), the unsteady nature of the solution is considered as an additional source of uncertainty. To assess the influence of the unsteady solution on the results, the mean force and wave amplitude were calculated over various windows of length 10 wave periods. For wave A09, which has a longer period and consequently fewer periods in the record, the uncertainty was estimated by comparing results averaged over 5 or 10 periods.
The overall uncertainty U in the simulations is estimated as where u i , (i = 1, . . . 4) is the uncertainty and v i is the mean value of each uncertainty related to each of the four factors discussed above. The individual mean values v i were calculated as where S max i and S min i are the largest and smallest values of the variable of interest in the corresponding sensitivity studies. Following the approach from Stern et al. (2001) for oscillatory convergence behaviour, the individual uncertainties were calculated as where F s = 1.1 is a safety factor.
A summary of the uncertainty analysis is shown in Table 7. The overall uncertainty was estimated to ±U = 17.4 % for wave A02 and U = 15.4 % for wave A09. The average over the two cases, U = 16.4 % has been applied in the visualisation of the CFD results in Figs. 9 and 10 .