Development of An Integrated Numerical Model for Simulating Wave Interaction with Permeable Submerged Breakwaters Using Extended Navier–Stokes Equations

An integrated two-dimensional vertical (2DV) model was developed to investigate wave interactions with permeable submerged breakwaters. The integrated model is capable of predicting the flow field in both surface water and porous media on the basis of the extended volume-averaged Reynolds-averaged Navier–Stokes equations (VARANS). The impact of porous medium was considered by the inclusion of the additional terms of drag and inertia forces into conventional Navier–Stokes equations. Finite volume method (FVM) in an arbitrary Lagrangian–Eulerian (ALE) formulation was adopted for discretization of the governing equations. Projection method was utilized to solve the unsteady incompressible extended Navier–Stokes equations. The time-dependent volume and surface porosities were calculated at each time step using the fraction of a grid open to water and the total porosity of porous medium. The numerical model was first verified against analytical solutions of small amplitude progressive Stokes wave and solitary wave propagation in the absence of a bottom-mounted barrier. Comparisons showed pleasing agreements between the numerical predictions and analytical solutions. The model was then further validated by comparing the numerical model results with the experimental measurements of wave propagation over a permeable submerged breakwater reported in the literature. Good agreements were obtained for the free surface elevations at various spatial and temporal scales, velocity fields around and inside the obstacle, as well as the velocity profiles.


Introduction
Porous structures are constructed along coasts for shoreline protection and beach erosion prevention against wave attack. Permeable submerged breakwater is one of the artificial structures that reduce the impact of waves and currents to beaches without affecting aesthetic aspects of coastal zones. Permeable submerged breakwaters have several advantages over impermeable obstacles or emerging structures. They provide lower construction costs and, unlike impermeable barriers, the wave interaction with porous structure takes place both inside and outside of the obstacle. This results in additional energy dissipation due to the friction created by porous medium. Although the emerged structures act as complete barriers, these low-crested breakwaters do not obstruct the overtopping water circulation, fish passage, or ship navigation. Therefore, understanding the wave behavior passing over permeable submerged structures is essential for many coastal engineering applications. The interactions between the waves and porous structures have been studied by many researchers. First, research in the scope of flow motion through porous medium were conducted by Sollitt and Cross [1]. They proposed the momentum equations by including additional terms of inertia and nonlinear resistance forces for the pore flow motion. They used their model to investigate wave transmission and reflection due to a permeable breakwater. Since then, a number of improved numerical models have been developed to investigate the interactions between wave and coastal constructions. These models are based on the mild-slope equations [2][3][4][5], shallow-water equations [6][7][8][9], and Boussinesq equations [10][11][12][13][14]. Because these simplified models are unable to predict some basic flow features such as nonlinearity, frequency dispersion, and wave breaking, recent studies have therefore concentrated on the modelling of flow motion using volume-averaged Reynolds-averaged Navier-Stokes (VARANS) equations. Sakakiyama and Kajima suggested extended Navier-Stokes equations, which include the impacts of porous medium by drag and inertia force components, and studied the wave transformation interacting with a permeable breakwater [15]. Van Gent utilized a model using Navier-Stokes equations to predict the wave motion both inside and outside the porous structures [16]. Liu et al. presented a model that simulates the interactions between waves and permeable structures [17]. The model calculates the flow field around the structure on the basis of the Reynolds-averaged Navier-Stokes (RANS) equations using an improved k-ε turbulence model. They utilized the spatially averaged Navier-Stokes equations for modeling pore flow motion. In their model, they assumed the turbulence within the porous medium negligible and resistance to flow through a porous medium was taken into account by considering linear and nonlinear frictional forces as suggested by Van Gent [16]. Hsu et al. improved the model by adding a k-ε turbulence model for the permeable structure, which makes the model applicable for modeling turbulent pore flow motion [18]. The newly developed model was employed for describing the wave motions around a coastal structure, which could be either a rigid obstacle, a permeable structure, or a composite breakwater. Huang et al. developed a numerical model in order to study the solitary wave interaction with a submerged structure [19]. They solved the unsteady Navier-Stokes equations for the flow outside the breakwater and another Navier-Stokes-based equation for solving the flow motion inside the structure. They added the convective inertia force and the viscous force into the suggested equations by Sollitt and Cross [1]. Then, they studied the effects of parameters such as the aspect ratio of the breakwater, porosity, and the wave height on the wave-submerged breakwater interactions. Unlike the earlier models, the numerical solution proposed by Karunarathna and Lin represents the fact that the flow resistance depends on the Reynolds number, which results in covering a wide range of flows [20]. They applied their model to investigate wave behavior passing over a permeable bed. Del Jesus et al. developed the multiphase VARANS equations that consider the spatial variation of porosity to simulate wave-structure interaction problems [21]. Wu and Hsiao studied the non-breaking solitary wave propagation over a permeable submerged structure using experimental and numerical models [22]. They employed the particle image velocimetry (PIV) method for measuring free surface elevation and velocity fields around a permeable structure instantaneously. They used VARANS equations including a nonlinear k-ε turbulence closure model as expressed by Hsu et al. [18]. Ma et al. presented a non-hydrostatic numerical model that simulates wave interactions with permeable structures [23], using the VARANS equations suggested by Del Jesus et al. [21]. To consider the temporally varying porosity at the computational cells which is due to surface displacements, they considered the porosity inside the derivatives and calculated the porosity values at each time step. Furthermore, the boundary element method (BEM) can be used for modeling wave interaction with solid obstacles, as described in [24,25]. In the present paper, a non-hydrostatic model presented by Hejazi et al. [26] was deployed, and the model was further developed to study the problems involving wave interactions with a permeable submerged breakwater. The extended Navier-Stokes equations derived by Sakakiyama and Kajima [15] were used to calculate the flow field in and outside the porous structure. Testing of the model's accuracy was performed by comparing the model results with experimental measurements for the surface displacement, velocity fields, as well as velocity profiles around the structure.

Governing Equations
The basic equations in this paper are the extended continuity and Navier-Stokes equations presented by Sakakiyama and Kajima [15], which solve the flow motion in both porous medium and fluid region, simultaneously. The extended continuity and Navier-Stokes equations in their conservative form are expressed as follows, which consider the effects of porous medium by including the additional terms of drag and inertia forces: where u and w are the flow velocity components in the x and z directions, respectively; t is the time; g is the gravitational acceleration; w g is the vertical grid velocity; v T is the dynamic viscosity; ρ is the local density; ρ r is the reference density; P * = ∆P/ρ r , where ∆P is the dynamic pressure that is obtained by reduction of hydrostatic pressure from the total instantaneous pressure (∆P = P t − P g ); γ v is the volume porosity; γ x and γ z are the surface porosity components in the x and z directions, respectively. The parameters λ v , λ x , and λ z in the above equations are defined as follows: where C M is the inertia coefficient. Furthermore, the drag force components in the above-mentioned equations are described by R x and R z as follows: where C D represents the drag coefficient.

Discretization of the Computational Domain
Finite volume method (FVM) in an arbitrary Lagrangian-Eulerian (ALE) system was applied to discretize the governing equations. The arbitrary Lagrangian-Eulerian formulation was developed to combine the advantages of the Lagrangian-based and Eulerian-based formulations, while minimizing their corresponding drawbacks as far as possible. When using the ALE formulation, the nodes on the boundaries and interfaces of the computational mesh can move along with fluid particles in normal Lagrangian fashion to precisely track the boundaries and interfaces, while the nodes of the computational mesh inside the domain can move arbitrarily in a normal Eulerian manner to optimize the shapes of volumes. Because of this freedom in nodes movement was allowed in the ALE formulation, significant distortions of the computational mesh can be reduced in comparison with a purely Lagrangian approach and more resolution can be provided than that allowed by a purely-Eulerian method. Figure 1 illustrates the grid generation system, points of action of scalar and vector quantities, as well as their corresponding control volumes. As shown in Figure 1, scalar variables such as density, pressure, and dynamic viscosity are defined at the center of an individual computational grid, and velocity components are considered at the center of the faces of each control volume. Furthermore, the bed level and free surface elevation in a water column are defined at the center of the lower and the upper side of the control volumes in the base layer and top layer, respectively.
in the ALE formulation, significant distortions of the computational mesh can be reduced in comparison with a purely Lagrangian approach and more resolution can be provided than that allowed by a purely-Eulerian method. Figure 1 illustrates the grid generation system, points of action of scalar and vector quantities, as well as their corresponding control volumes. As shown in Figure  1, scalar variables such as density, pressure, and dynamic viscosity are defined at the center of an individual computational grid, and velocity components are considered at the center of the faces of each control volume. Furthermore, the bed level and free surface elevation in a water column are defined at the center of the lower and the upper side of the control volumes in the base layer and top layer, respectively. The volume ( v) and surface porosity components ( x, z) are defined at the center and the faces of each computational grid. A schematic diagram of a partially filled control volume and the definitions of the associated parameters are demonstrated in Figure 2. By considering the temporally varying free surface, the porosity components are updated at each time step by calculating the fraction of each computational grid open to water and the total porosity of porous medium, as follows: The volume (γ v ) and surface porosity components (γ x , γ z ) are defined at the center and the faces of each computational grid. A schematic diagram of a partially filled control volume and the definitions of the associated parameters are demonstrated in Figure 2. By considering the temporally varying free surface, the porosity components are updated at each time step by calculating the fraction of each computational grid open to water and the total porosity of porous medium, as follows:  It should be noted that for a control volume full of water, the volume and surface porosities are equal to 1 (γv = γx = γz = 1), whereas for a control volume full of porous medium, the porosity components are obtained as total porosity (γv = γx = γz = n). It should be noted that for a control volume full of water, the volume and surface porosities are equal to 1 (γ v = γ x = γ z = 1), whereas for a control volume full of porous medium, the porosity components are obtained as total porosity (γ v = γ x = γ z = n).

Numerical Implementation
In this study, projection approach suggested by Chorin [27] and Temam [28] was used to solve the extended Navier-Stokes equations in two main steps. In the first fractional step, the pressure fields were excluded from the Navier-Stokes equations and the momentum equations in the absence of pressure gradients were solved to compute an intermediate velocity field (u * , w * ) that did not satisfy the continuity constraint. The momentum equations in the absence of pressure terms are as follows: where u * and w * are intermediate velocity components, and the superscript n denotes the time level, where t n = n.∆t. The terms on the left-hand side of Equations (9) and (10) are advective terms, the first and second terms on the right-hand side of the above equations are diffusive terms, and the third is the drag force component. Therefore, the first step was split into three sub-fractional steps. The first sub-fractional step consists of solving the advective terms and computing the corresponding velocity fields (u A , w A ) as follows: where the superscript A stands for advection and shows the time level after completing the advection process. It should be noted that Equations (11) and (12) were further split into other three sub-fractional steps. A locally one-dimensional method was applied for solving the above equations and a fifth-order-accurate scheme was utilized to obtain more accurate results for flow characteristics. The second sub-fractional step consists of solving the diffusive terms and calculating the associated velocities (u D , w D ) as follows: where the superscript D stands for diffusion and shows the time level after completing the diffusion process. θ D is the implicit weighting factor, where θ D = 1 denotes a fully implicit approach, whereas θ D = 0 implies a fully explicit approach. Equations (13) and (14) were split into other two sub-fractional steps, and the Crank-Nicolson method (θ D = 0.5) with a central discretization scheme was employed to solve them. The third sub-fractional step consists of solving the drag force components and computing the intermediate velocity fields (u * , w * ) as follows: In the second fractional step, the extended continuity and momentum equations in the absence of advective and diffusive terms and drag force components were solved for the computational cells below the top layer as follows: According to the Gauss divergence theorem, Equation (17) can be written as: where A i,k can be calculated as: By considering Equations (18) and (19), the parameters u n+1 i±1,k+2 , u n+1 i±1,k , u n+1 i±1,k−2 , and w n+1 i,k±1 can be expressed as follows: Substituting Equations (22)- (25) into Equation (20), the pressure Poisson equation (PPE) can be derived as follows: Using the Gauss divergence theorem for pressure derivatives on the velocity locations, the pressure Poisson equation can be written according to the pressure values in the 15 adjacent cells.

Pressure Equation at Free Surface Layer
The surface equation is calculated by integrating the extended continuity equation over the depth and considering the kinematic boundary conditions at the impermeable bed level and free surface, as follows: where ξ is defined as the surface elevation above the mean water level and z b is the impermeable bed level above the datum, which is equal to zero for a fixed bed. The mass continuity equation for the ith water column can be obtained as: where kc min and kc max are the reference numbers of the bottom and the top layers, respectively. The hydrostatic pressure assumption at the top layer may result in inaccurate predictions of the wave phase and water surface elevation. Therefore by applying the non-hydrostatic pressure distribution at the top layer, more accurate results can be obtained for the flow features. In this study, the method proposed by Yuan and Wu was utilized to calculate the pressure equation at free surface layer [29,30]. For the ith water column, the momentum equation in the vertical direction from the center of the top layer to free surface level can be written as follows: where the subscripts T and S represent the center of the top half-grid and the free surface level, respectively. The parameters w n+1 i,T and w * i,T are the vertical velocity components at the distance of 1/4∆z i from the free surface level and can be approximated by the following equations: In Equation (29), p * n+1 i,S and p * n i,S are the pressure at free surface level (kc max + 1) and are calculated as follows: By substituting Equations (30)-(33) into Equation (29), the momentum equation in the vertical direction and at the top half-grid can be written as follows: The volume porosity and density are defined at the center of each grid. Hence, by considering the assumption that (ρ/ρ r ) i,T = (ρ/ρ r ) i,kc max , γ v i,T = γ v i,kcmax , and λ v i,T = λ v i,kcmax , and by multiplying Equation (34) by 4∆t, the following equation can be obtained: According to Dean and Dalrymple, the kinematic free surface boundary condition is expressed as [31]: By substituting Equation (27) into Equation (36), the above equation at surface level can be discretized as: The parameter w n+1 i,kc max −1 can be obtained by writing the momentum equation in the vertical direction, as follows: By substituting ξ n+1 i , w n+1 i,kc max +1 , and w n+1 i,kc max −1 from Equations (28), (37), and (38), respectively, Equation (35) can be written as: By substituting horizontal velocity components from Equation (23) into Equation (39), the pressure equation at surface layer is obtained as: By applying the finite difference method (FDM) for discretization of pressure derivatives at the top layer and gauss divergence theorem for pressure derivatives below the top layer, the pressure equation at the free surface layer can be obtained according to the pressure values at ith column and two adjacent water columns.

Numerical Solution of the Pressure Poisson Equation (PPE)
By applying Equation (40) for cells that are located at the top layer of the numerical domain and pressure Poisson Equation (26) for cells below the top layer, the system of partial differential equations (PDEs) is written as follows: where ic is the reference number of the center of computational grids in the x direction. The above equation can be rewritten as follows: where p * n+1 ic is the unknown pressure vector at icth column and is expressed as: The coefficient matrix in Equation (42) forms a block tri-diagonal matrix, which has the dimension of 1/2 ic max − ic min + 1 × 1/2 ic max − ic min + 1 , and each block takes the following form: with the exception of top layer, it is a five-diagonal matrix with the dimension of 1/2 kc max − kc min + 1 × 1/2 kc max − kc min + 1 . In the presented model, the system of partial differential equations (Equation (42)) was solved by applying a direct approach of the Thomas algorithm, and pressure fields were calculated for the whole computational domain. After calculating the pressure fields, the velocity components at the new time level (u n+1 , w n+1 ) were obtained by solving Equations (18) and (19). Then, the free surface elevation was updated by using Equation (28) and the new grid system was generated to fit the new surface level and to initiate the next time step.

Initial Boundary Conditions
Two types of boundary conditions, namely, Dirichlet and Neumann were used in the present model. The Dirichlet boundary condition assumed zero for vertical velocity components at impermeable bed, and the Neumann boundary condition was imposed by assuming zero for normal gradient of the horizontal velocity components at impermeable bed and the vertical velocity components at the left side of the numerical domain. Furthermore, horizontal velocity components at the upstream boundary were specified as inflow conditions at x = 0. A free exit for the flow was considered at the downstream boundary of the domain by considering the dynamic pressure equal to zero.

Results and Discussions
Before proceeding to the study of wave interaction with a permeable coastal structure, the accuracy of the model in the absence of the porous medium must be evaluated. In this regard, the numerical model is firstly validated by comparing the numerical predictions with the analytical solutions proposed for the determination of the Stokes and solitary wave propagation in a wave tank. It can be derived from Equations (7) and (8)-for a control volume full of water the porosity components are equal to 1 (γ v = γ x = γ z = 1) and the governing extended Navier-Stokes equations become identical to the conventional equations. This implementation was considered in the two following sections. Then, the model was utilized to investigate the solitary wave interaction with a permeable submerged breakwater and the model results were compared with the experimental measurements.

Progressive Stokes Waves Propagation
The capability of the newly developed model for simulating the non-linear progressive wave train in a deep water was firstly evaluated. An incident wave with the height of H = 0.5 m and wave period of T = 5 s propagated in a water depth of h = 15 m. According to the validity range of various wave theories presented by Le Méhauté [32], the second-order Stokes wave was obtained for the above features. For second-order Stokes waves, the linear dispersion Equation (45) is valid [33]; therefore, the wave number can be obtained as k = 0.16344 rad/m.
Using equations k = 2π/L and C = L/T, the wave length and wave celerity are calculated as 38.44 m and , respectively. The numerical domain with a length of L = 1500 m was discretized with uniform grids of ∆x = 1 m in the horizontal direction and 15 layers in the vertical direction. The increment of dimensionless time was chosen to be ∆t = 0.05 s and a total simulation time of t = 140 s was adopted. For second-order Stokes waves, the analytical equations for surface elevation and velocity components are obtained by [34]: where z is the elevation above the datum and σ = 2π/T is the angular frequency. Comparison between numerical result and analytical solution for free surface displacements at t = 140 s, is shown in Figure 3. As shown in Figure 3, the predicted results agree well with analytical solution for free surface evolution.  Figure 3, the predicted results agree well with analytical solution for free surface evolution.

Solitary Wave Propagation
A numerical wave tank with the length of L = 2000 m was used and a solitary wave with the amplitude of H = 1 m propagated in a constant water depth of h = 10 m. A uniform grid system with ∆x = 2 m in the x direction and 10 layers in the z direction was deployed for the discretization of the computational domain. The numerical time interval was set at ∆t = 0.1 s, and calculations were continued for t = 180 s. According to Boussinesq equations, the analytical free surface displacements and velocity components are derived as proposed by Lee et al. [33]: where c = g(h + H) is the celerity of the solitary wave which can be obtained as c = 10.388 m/s. By assuming the upstream inflow boundary at x = −300 m, the solitary wave in the computational domain will start at its highest point.           The figures show that the newly developed model provides a high accuracy level in predicting the free surface distribution and velocity components.

Solitary Wave Interaction with a Permeable Submerged Breakwater
In this section, the newly developed model was deployed to investigate the solitary wave propagation over a submerged permeable breakwater. Numerical outputs were compared with the experimental data presented by Wu and Hsiao (2013). The experiment was conducted in a wave tank with the length of L = 25 m and a permeable submerged structure with the length of a = 13 cm and height of b = 6.5 cm was situated at the bottom. The breakwater was made of glass spheres with a diameter of D = 1.5 cm, yielding a total porosity of n = 0.52. A solitary wave with a height of H = 4.77 cm propagated in a water depth of h = 10.6 cm. The PIV method was utilized for the experimental measurements of flow field around the porous structure. The origin of the coordinate system was assumed to be at the intersection of the impermeable bed and the left side of the breakwater. The time series of the surface displacement was measured by two capacity gauges located at x = −1.8 m and x = 1.8 m. In calculations, the time t = 0 s was defined as the time when the solitary wave crest reached the first gauge. Figure 7 shows the schematic description of the numerical wave tank.
The length of wave tank was discretized with uniform grids of ∆x = 0.005 m, and 20 layers were considered over the depth of calculation domain. The total simulation time was t = 25 s and the time interval was chosen to be ∆t = 0.001 s. The coefficients C D and C M that produced the best agreement with the experimental measurements were considered as the calibration parameters in the numerical simulation. A number of values for each coefficient were considered and the effect of each individual parameter was studied by varying one coefficient while keeping another one unchanged and calculating the root mean square error (RMSE) for each set. In this study, the set of C D = 3.5 and C M = 0.5 resulting in the minimum RMSE was utilized for the numerical calculations. Figure 8 illustrates the comparison between experimental and numerical time histories of surface displacement at the location of measurement gauges.

Solitary Wave Interaction with a Permeable Submerged Breakwater
In this section, the newly developed model was deployed to investigate the solitary wave propagation over a submerged permeable breakwater. Numerical outputs were compared with the experimental data presented by Wu and Hsiao (2013). The experiment was conducted in a wave tank with the length of = 25 m and a permeable submerged structure with the length of = 13 cm and height of = 6.5 cm was situated at the bottom. The breakwater was made of glass spheres with a diameter of = 1.5 cm, yielding a total porosity of = 0.52. A solitary wave with a height of = 4.77 cm propagated in a water depth of ℎ = 10.6 cm. The PIV method was utilized for the experimental measurements of flow field around the porous structure. The origin of the coordinate system was assumed to be at the intersection of the impermeable bed and the left side of the breakwater. The time series of the surface displacement was measured by two capacity gauges located at = −1.8 m and = 1.8 m. In calculations, the time = 0 s was defined as the time when the solitary wave crest reached the first gauge. Figure 7 shows the schematic description of the numerical wave tank. The length of wave tank was discretized with uniform grids of ∆ = 0.005 m, and 20 layers were considered over the depth of calculation domain. The total simulation time was = 25 s and the time interval was chosen to be ∆ = 0.001 s. The coefficients and that produced the best agreement with the experimental measurements were considered as the calibration parameters in the numerical simulation. A number of values for each coefficient were considered and the effect of each individual parameter was studied by varying one coefficient while keeping another one unchanged and calculating the root mean square error (RMSE) for each set. In this study, the set of = 3.5 and = 0.5 resulting in the minimum RMSE was utilized for the numerical calculations. Figure 8 illustrates the comparison between experimental and numerical time histories of surface displacement at the location of measurement gauges. As shown in Figure 8, the computed outputs agreed fairly well with the measured outputs for the incident, reflected, and transmitted waves. Figures 9-18 show the comparison between numerical predictions and experimental measurements for the spatial free surface distributions, velocity fields, and velocity profiles at five different time levels. According to the modelled and measured results, when the solitary wave front reached the permeable submerged breakwater, due to the hydraulic jump formation, the flow was separated from the top of the obstacle and small clockwise vortices were generated at the leading edge of the obstacle (Figure 9). When the solitary wave passed over the structure, a primary vortex was generated on the lee side of the structure (Figure 11). This vortex grew in size and penetrated into the deeper layers of water by increasing time (Figure 13). On the basis of the last two phases, the primary vortex gradually reached the free surface, lost its strength, and finally faded away due to the As shown in Figure 8, the computed outputs agreed fairly well with the measured outputs for the incident, reflected, and transmitted waves. Figures 9-18 show the comparison between numerical predictions and experimental measurements for the spatial free surface distributions, velocity fields, and velocity profiles at five different time levels. According to the modelled and measured results, when the solitary wave front reached the permeable submerged breakwater, due to the hydraulic jump formation, the flow was separated from the top of the obstacle and small clockwise vortices were generated at the leading edge of the obstacle (Figure 9). When the solitary wave passed over the structure, a primary vortex was generated on the lee side of the structure (Figure 11). This vortex grew in size and penetrated into the deeper layers of water by increasing time (Figure 13). On the basis of the last two phases, the primary vortex gradually reached the free surface, lost its strength, and finally faded away due to the diffusive effects (Figures 15 and 17). It is also seen that due to the resistance forces of the porous medium, the velocity inside the structure was substantially smaller than the velocity values on the top of the obstacle.                  Comparisons between numerical model predictions and experimental measurements by Wu and Hsiao [23] show the capability of the newly developed model for the simulation of wave interactions with permeable submerged structures.  Comparisons between numerical model predictions and experimental measurements by Wu and Hsiao [23] show the capability of the newly developed model for the simulation of wave interactions with permeable submerged structures.  Comparisons between numerical model predictions and experimental measurements by Wu and Hsiao [23] show the capability of the newly developed model for the simulation of wave interactions with permeable submerged structures.

Conclusions
Extended Navier-Stokes equations with a non-hydrostatic pressure distribution were solved numerically in the whole computational domain to investigate the wave interactions with permeable submerged structures. The impact of porous medium was considered by including additional terms for drag and inertia forces into general Navier-Stokes equations. Finite volume method (FVM) in an arbitrary Lagrangian-Eulerian (ALE) formulation was utilized for the discretization of the numerical domain, in which the number of layers was constant during the whole calculations and the new mesh was generated after updating the surface elevation at each time step. The porosity components were updated by calculating the fraction of cell open to water and the total porosity of porous medium. Projection method was used to solve the extended Navier-Stokes equations in two main steps. In the first step, the Navier-Stokes equations in the absence of pressure gradients were solved to calculate an intermediate velocity field, which did not satisfy the continuity equation. In the second step, the continuity and momentum equations in the absence of advection, diffusion, and drag force terms were solved, and the pressure Poison equations were obtained for the computational grids below the top layer. The pressure equation for the computational cells located at the free surface layer was obtained by integrating the continuity equation over the depth and considering appropriate boundary conditions for the impermeable bed and free surface. The newly developed model in the absence of porous medium was firstly validated against analytical solutions of second-order Stokes and solitary wave propagation under non-breaking conditions. Numerical results agreed well with the analytical solutions. The model was then deployed for the simulation of a solitary wave passing over a rectangular porous structure. Numerical predictions for the time histories and spatial distribution of the free surface elevations, velocity fields, and velocity profiles were compared with the reported experimental measurements using the PIV method. Results showed that when the solitary wave reached the permeable submerged breakwater, a hydraulic jump was formed and the flow was separated from the top-left corner of the breakwater. When the solitary wave passed over the obstacle, a primary vortex was created behind the breakwater, which grew in time and penetrated into the bottom of the wave tank. This vortex gradually reached the upper layers of water and decreased in strength and finally disappeared due to the diffusive impacts. Comparisons between numerical results and experimental data proved the capability of the newly developed model in simulating the integrated fluid-porous medium problems and the interaction of wave with submerged breakwaters.
Author Contributions: P.P. contributed to the methodology, investigation, validation, visualization and writing the manuscript of the paper, while H.K. contributed to the conceptualization, resources, project administration and supervision. Authorship has been limited to those who contributed substantially to the work reported. All authors have read and agreed to the published version of the manuscript.