*2.1. Flow Model*

In this paper, the flow model is based on one of the solvers in OpenFOAM®, IHFOAM. Recently, some other open-source codes have been developed for modeling wave propagation and wave–structure interaction problems [39,40] as well. To simulate the coastal, offshore, and hydraulic engineering process, this model solves the three-dimensional Volume Averaged Reynolds Averaged Navier–Stokes (VARANS) equations with regard to two incompressible phases of air and water. The fluid model adopts the finite volume discretization and the volume of fluid (VOF) method [37]. Including continuity and momentum conservation equations; the VARANS equations as the governing mathematical expressions in this model can be expressed as:

$$\frac{\partial \langle u\_i \rangle}{\partial x\_i} = 0,\tag{1}$$

*J. Mar. Sci. Eng.* **2019**, *7*, 369

$$\frac{\partial \rho\_f \langle u\_i \rangle}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} \left[ \frac{1}{n} \rho\_f \langle u\_i \rangle \langle u\_j \rangle \right] = -n \frac{\partial \langle p^\* \rangle^f}{\partial \mathbf{x}\_i} + n \rho\_f \mathbf{g}\_i + \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu\_{eff} \frac{\partial \langle u\_i \rangle}{\partial \mathbf{x}\_j} \right] - [\mathbf{C}T]\_\mathbf{\prime} \tag{2}$$

where hi is Darcy's volume averaging operator, while hi*<sup>f</sup>* is the intrinsic averaging operator; *ρf* represents the density which is computed by *ρ<sup>f</sup>* = *αρwater* + (1 − *α*)*ρair*; *α* is the indicator function defined in (5); *u<sup>i</sup>* is the velocity vector while *n* is the porosity; *p* ∗ is the pseudo-dynamic pressures; **g**, *g<sup>i</sup>* is the gravitational acceleration; *µe f f* is the efficient dynamic viscosity, defined as *µe f f* = *µ* + *ρ<sup>f</sup> νturb*, in which *µ* is the molecular dynamic viscosity and *νturb* is the turbulent kinetic viscosity, given by the chosen turbulence model. The *k* − *e* turbulence model is adopted in this study; the last term in (2) corresponds to the resistance of porous media, which is shown as:

$$\left[\left|\mathbf{C}T\right] = A\left + B\left|\left\right|\left + \mathcal{C}\frac{\partial\left}{\partial t},\tag{3}$$

compared to factors *A* and *B*, factor *C* is less significant. Moreover, 0.34 (kg/m<sup>3</sup> ) is often adopted for the value of *C* by default [41].

In present fluid model, the values of *A* and *B* are derived by Engelund [42]'s formulae, which also employed in Burcharth and Andersen [43].

$$A = E\_1 \frac{(1-n)^3}{n^2} \frac{\mu}{d\_{50}^2}, \text{ and } B = E\_2 (1 + \frac{7.5}{\text{KC}}) \frac{1-n}{n^2} \frac{\rho\_f}{d\_{50}}.\tag{4}$$

in which *d*<sup>50</sup> is the medium grain diameter of the materials; *KC* is the Keulegan Carpenter number, which is defined as *KC* = *umT*0/(*nd*50); *u<sup>m</sup>* is the maximum oscillating velocity while *T*<sup>0</sup> is the period of the oscillation. *E*<sup>1</sup> and *E*<sup>2</sup> are parameters that define the linear and non-linear friction terms. The default values of these parameters (*E*<sup>1</sup> = 50 and *E*<sup>2</sup> = 1.2) are used in this study [44].

A two-phased fluid mixture containing air and water is taken into consideration in each cell, which can be controlled by an indicator function *α*. In the present model, *α* is defined as the quantity of water per unit of volume, which varies from 0 (air) to 1 (water):

$$a = \begin{cases} 1, & \text{water} \\ 0, & \text{air} \\ 0 < a < 1, & \text{free surface} \end{cases} \tag{5}$$

As an phase function, *α* can represent any variation of fluid properties considering the mixture properties, such as density and viscosity:

$$
\Phi = \mathfrak{a}\Phi\_{water} + (1 - \mathfrak{a})\Phi\_{air\prime} \tag{6}
$$

where Φ*water* and Φ*air* stand for the properties of water and air, separately, such as density of the fluid.

The fluid movement could be traced by solving the advection equation as [44]:

$$\frac{\partial \alpha}{\partial t} + \frac{1}{n} \frac{\partial \langle u\_l \rangle \alpha}{\partial x\_l} + \frac{1}{n} \frac{\partial \langle u\_{cl} \rangle \alpha (1 - \alpha)}{\partial x\_l} = 0,\tag{7}$$

where |*uc*| = *min* [*cα*|*u*|, *max*(|*u*|)], in which the default value of *c<sup>α</sup>* is 1, however, the user can specify a greater value to enhance the compression of the interface, or zero to eliminate it.

The solving algorithm used in flow model is PIMPLE, which is a combination of PISO (Pressure Implicit with Splitting of Operators) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithms. In this study, the *κ* − *e* RAS (Reynolds-averaged simulation) turbulence model is adopted to model the turbulent viscosity *νturb* as:

$$\nu\_{turb} = \frac{\mathbb{C}\_{\mu} \kappa^2}{\mathfrak{e}} \tag{8}$$

where *C<sup>µ</sup>* is an empirical constant; *κ* is the turbulence kinetic energy while *e* is the turbulence energy dissipation rate, separately.

The IHFOAM implements the wave generation and active wave absorption in the fluid domain, which introduce several boundary conditions: (i) the inlet boundary condition allows the generation of wave according to different wave theories as well as adding different steady current flow; (ii) the outlet boundary condition applies an active wave absorption theory to prevent the re-reflection of incoming wave; (iii) the slip boundary condition (zero-gradient) is applied on the bottom of the fluid domain and the lateral boundary of the numerical wave flume; (iv) the top boundary condition is set as the atmospheric pressure. The details of IHFOAM could be found in Higuera et al. [37].
